Transcriptomic Analysis Based on Bulk RNA-seq Data

Summary

The present protocol establishes a complete pipeline for analyzing the process of bulk RNA-seq from raw data to functional enrichment analysis.

Abstract

Nonalcoholic fatty liver (NAFL) is usually considered a benign condition; however, once it progresses to non-alcoholic steatohepatitis (NASH), patients face a significantly elevated risk of developing end-stage liver disease. Many studies are attempting to elucidate the molecular mechanism underlying the transition from NAFL to NASH. High-throughput sequencing technologies (such as bulk RNA-seq) have provided researchers with a deeper understanding by examining the transcriptome, revealing the expression of molecules, activation of signaling pathways, and other factors associated with disease progression. There is a wealth of open-source data available for researchers to analyze in order to identify potential targets for disease treatment. However, related research is limited by the lack of an efficient and reliable process for upstream analysis of the transcriptome. Here, a highly reproducible and user-friendly upstream analysis and subsequent related differential gene analysis pipeline is provided to achieve standardized processing and deep parsing of private or public data. The pipeline is divided into four steps: (1) quality control of data; (2) gene mapping; (3) differential gene analysis; and (4) functional analysis. This process aims to uncover the molecular mechanisms of disease transformation and assist researchers in screening potential drug targets and therapeutic approaches through the analysis of Bulk RNA-seq data.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is the most prevalent chronic liver disease globally, affecting more than a quarter of the population. Its incidence has increased dramatically in recent decades^1,2,3. The growing disease burden, especially its more advanced form, non-alcoholic steatohepatitis (NASH), poses a major global health challenge and a heavy economic burden⁴. The first stage of NAFLD is non-alcoholic fatty liver (NAFL), which is accompanied by inflammation and fibrosis that can progress to NASH. The latter significantly increases the risk of progression to end-stage liver disease, including cirrhosis and hepatocellular carcinoma (HCC)^5,6,7. HCC incidence and mortality are associated with an increase in NASH^8,9, and it is expected that NAFLD/NASH will become the leading indication for liver transplantation by 2030¹⁰. However, the clinical progression of NAFLD is highly heterogeneous¹¹, which severely hampers the development of relevant drugs¹², making it particularly important to precisely explore the molecular mechanisms involved.

Bulk RNA-seq-based acquisition of cellular compositional information can significantly elucidate the pathogenesis of various diseases. In recent decades, numerous bulk RNA-seq studies have been conducted in model organisms and humans to elucidate gene expression differences in NASH progression^13,14,15, to identify new therapeutic targets for intervention. Based on bulk RNA-seq analysis, Xiong et al. found that nonparenchymal cells (NPCs) in the liver are involved in processes such as extracellular matrix formation and cell adhesion, which contribute to the progression of NASH¹⁶. Li et al. demonstrated that hepatic Wilms' tumor 1-associating protein (WTAP) in hepatocytes regulates ectopic lipid accumulation and inflammation, thereby promoting NASH formation¹⁷. Although bulk RNA-seq analysis is a powerful tool for elucidating the mechanisms of NASH, its results are highly sensitive to the quality of upstream data. The heterogeneity of upstream experimental operations and analysis processes can seriously impair the reliability of data, thereby masking true biological information and interfering with the accuracy of subsequent analyses. Therefore, it is important to establish a set of standardized upstream analysis procedures.

Compared with single-cell RNA sequencing (scRNA-seq), bulk RNA-seq offers several distinct advantages in both experimental design and practical applications. While scRNA-seq enables the identification of cellular heterogeneity at the single-cell level and allows precise analysis of cell type-specific transcriptional features, it is associated with high cost, complex data processing requirements, and limited sensitivity for detecting low-abundance transcripts¹⁸. In contrast, bulk RNA-seq provides higher sequencing depth, lower cost, and greater sample throughput, making it particularly suitable for population-level differential gene expression analyses and the exploration of molecular mechanisms¹⁹. Therefore, when guided by standardized analytical workflows, bulk RNA-seq remains an efficient, cost-effective, and robust approach for investigating the molecular basis of complex diseases.

This protocol is designed specifically for bulk RNA-seq datasets derived from human tissues with high RNA integrity (RIN ≥ 7.0) and sufficient input RNA (≥ 500 ng per sample). To ensure reliable execution of alignment and quantification steps, a local workstation equipped with at least a 10-core CPU, 32 GB of RAM, and a minimum of 200 GB of free disk space is recommended. Building on these requirements, the protocol provides an efficient and user-friendly analytical workflow, including detailed operational instructions and standardized parameter configurations, to meet the needs of researchers analyzing large-scale transcriptomic data.

Protocol

For demonstration purposes, the publicly available dataset PRJNA1023502 generated by Lan Bai et al. was used to illustrate each step of both upstream and downstream analyses²⁰. As this dataset originates from the open-access NCBI SRA database, no additional permissions or ethical approvals are required. See the Table of Materials to verify all required software and R-package versions. The publicly available dataset PRJNA1023502 comprises 6 not-NASH, 6 NAFL, and 6 NASH liver RNA-seq samples. In this protocol, the dataset was used to demonstrate all steps of the bulk RNA-seq workflow, including data retrieval from the SRA database, quality control (fastp), alignment (HISAT2), quantification (featureCounts), and downstream differential expression and functional enrichment analyses.

1. SRA toolkit installation

1. Visit the SRA Toolkit official website and download version 3.2.1.

2. Public data download

1. Obtain the SRA number.
   1. Access it in the article supplement²⁰, the Data Availability section, or by searching for keywords in the NCBI SRA database.
   2. Type prefetch <SRA number> in the terminal to download it.

3. Generation of gene count matrix

1. Define paths and settings as described below:
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. Download the human reference genome (GRCh38 primary assembly) and the corresponding gene annotation file (gencode.v44, comprehensive gene annotation set) from the official GENCODE database (https://www.gencodegenes.org/human/).
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. Rename SRA files.
   1. Ensure all SRA files use the .sra file extension to facilitate accurate recognition and processing by downstream tools.
   2. Execute the following commands：
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Execute the following command to convert SRA files into FASTQ format:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Build HISAT2 index from the reference genome.
   hisat2-build $REFERENCE $INDEX
5. Process FASTQ files with visual progress monitoring.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. For paired-end sequencing data, execute the following command
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Perform quality control and filtering using fastp and output the results to the $OUT_FASTP directory. For single-end sequencing data, execute the following command:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      For paired-end sequencing data, execute the following command:
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. For each sample, inspect the HTML report for visual quality control and consult the JSON report for structured metrics, including read quality, GC distribution, read redundancy, and adapter contamination. Ensure all samples meet the following quality thresholds: Q30 ≥ 85%, adapter contamination rate < 5%, and no abnormal GC patterns. Store all output files in the $OUT_FASTP directory.
   4. Perform sequence alignment using HISAT2, generate SAM files, and output them to the $OUT_HISAT2 directory. For single-end sequencing data, execute the following command:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      For paired-end sequencing data, execute the following command:
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. Perform HISAT2 alignment of the quality-filtered reads to the human reference genome, using an index built from the GRCh38 primary assembly FASTA file. Provide example commands for both single-end and paired-end sequencing data.
   6. Convert SAM to BAM, sort, and index.
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. Generate a SAM file for each sample and convert it to a sorted and indexed BAM file using SAMtools. For high-quality human transcriptomic samples (RNA integrity number, RIN ≥ 8.0), ensure that the overall alignment rate exceeds 85%; for typical bulk RNA-seq samples (RIN ≥ 7.0), consider alignment rates ≥ 70% acceptable.
   8. Perform gene quantification using featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Check the tab-delimited output file (*.counts.txt) and its summary report (*.counts.txt.summary) generated by featureCounts for each sample. Ensure the read assignment rate meets the typical threshold of ≥70% for human bulk RNA-seq; a markedly lower rate may indicate mismatched strand orientation, annotation issues, or poor alignment quality. Use the parameter -s 0 for this non-strand-specific RNA-seq dataset. For strand-specific libraries, substitute -s 0 with -s 1 or -s 2 in the command.
6. Generate the gene counts matrix as described below.
   1. Initialize the counts matrix with gene IDs and counts from the first sample.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Iteratively add counts from each sample into a single matrix.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Add a header line with sample IDs to the counts matrix.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Extract gene lengths from GTF file (sum of exon lengths per gene).
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. Raw count matrix processing and gene annotation

1. R environment initialisation and loading of the relevant R package.
   1. Use install.packages() function to install tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR. The tidyverse package is employed for data manipulation and plotting; the ggplot2 package is employed for visualization; the ggrepel package is employed for non-overlapping text labels; the RColorBrewer package is employed for color palettes; the ggridges package is employed for ridge plots; the FactoMineR package is employed for PCA and multivariate analysis.
   2. Use BiocManager::install() function to install biomaRt, DESeq2, clusterProfiler. The biomaRt package is employed for gene annotation from Ensembl; the DESeq2 package is employed for differential expression analysis; the clusterProfiler package is used for functional enrichment analysis.
   3. Read the original count matrix generated in step 2, use biomaRt to map Ensembl ID to HGNC gene name for subsequent analysis, remove duplicate Geneid column, and generate the cleaned count matrix. Use the original counts matrix (counts_matrix.csv) as input, with Ensembl gene IDs as rows and samples as columns. Save the cleaned counts matrix as output (clean_counts_SRA.csv), with HGNC gene symbols as rows and samples as columns.
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. Gene expression quantification

NOTE: Refer to the Supplementary File 1 for the detailed script.

1. Execute the following command to calculate fragments per kilobase of transcript per million mapped reads (FPKM) and transcripts per million (TPM).
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. Sample clustering and difference visualization

1. Assess data quality by examining sample clustering through PCA. Refer to Supplementary File 1 for the detailed script.
2. For PCA visualization, use the code below.
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. Differential expression analysis and visualization of results

NOTE: Refer to the Supplementary File 1 for the detailed script.

1. Execute the following commands to construct the DESeq2 dataset and set the threshold. Only genes with total counts > 1 across all samples were retained prior to modeling. Identify differentially expressed genes (DEGs) using the significance thresholds of adjusted p-value < 0.05 and |log2FC| > 0.5.
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. Visualize differentially expressed genes with a volcano plot.
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. Perform functional enrichment analysis and visualization

NOTE: Refer to the Supplementary File 1 for the detailed script.

1. Visualize KEGG enrichment results.
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. Visualize GO enrichment results.Customize the visualization by modifying arc resolution (n=500), color mappings, and facet arrangement according to experimental specifications.
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. Visualize gene set enrichment analysis (GSEA) results. Adjust the number of pathways, density ridge scaling, and layer transparency to accommodate different effect size distributions.
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

The upstream analysis workflow for bulk RNA-seq is illustrated in Figure 1A. This workflow sequentially executes the following key steps on a Linux platform: first, rigorous quality control of raw sequencing data is performed using fastp to remove low-quality reads and adapter sequences; subsequently, HISAT2 aligns high-quality reads to the reference genome, with Samtools converting and sorting the alignment files; finally, FeatureCounts performs gene-level quantification to generate a gene ...

Discussion

Bulk RNA-seq data analysis is characterized as an interdisciplinary task that integrates genomics, bioinformatics, statistics, and computer science. A complete analytical workflow encompasses multiple upstream and downstream steps, including raw data preprocessing, quality control, sequence alignment, gene-level quantification, data normalization, differential expression analysis, and biological interpretation. Among these steps, accurately converting raw sequencing reads into a high-quality gene expression matrix is par...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Bioconductor	2.64.0	Gene annotation from Ensembl
clusterProfiler	Bioconductor	4.16.0	Functional enrichment analysis
DESeq2	Bioconductor	1.48.1	Differential expression analysis
FactoMineR	AgroParisTech	2.11.0	PCA and multivariate analysis
fastp	OpenGene	1.0.1	Quality control and filtering of FASTQ data
FeatureCounts	Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research	2.0.0	Count the number of reads mapped to each gene for gene expression quantification
ggplot2	Posit	3.5.2	Data visualization
ggrepel	Kamil Slowikowski	0.9.6	Non-overlapping text labels
ggridges	Claus O. Wilke	0.5.6	Create ridgeline plots
HISAT2	Johns Hopkins University	2.2.1	Align the filtered high-quality reads to the reference genome
R	R Core Team	4.5.0	An environment for data computation, analysis, and visualization
RColorBrewer	Erich Neuwirth	1.1.3	Color palettes for plotting
samtools	Large Scale Genomics work stream	1.22.0	Convert and process SAM files for efficient retrieval and access
SRA Toolkit	National Center for Biotechnology Information	3.2.1	Obtain and preprocess raw sequencing data from the NCBI SRA database