基于大量RNA测序数据的转录组分析

摘要

本方案建立了一个完整的分析流程，用于从原始数据到功能富集分析的批量RNA测序（bulk RNA-seq）研究。

摘要

非酒精性脂肪肝（NAFL）通常被认为是一种良性疾病；然而，一旦进展为非酒精性脂肪性肝炎（NASH），患者发展为终末期肝病的风险将显著增加。目前许多研究正致力于阐明NAFL向NASH转变的分子机制。高通量测序技术（如批量RNA测序，bulk RNA-seq）通过分析转录组，揭示了与疾病进展相关的分子表达、信号通路激活及其他因素，使研究人员获得了更深入的理解。目前已有大量开源数据可供研究人员分析，以识别潜在的疾病治疗靶点。然而，相关研究受限于缺乏高效且可靠的转录组上游分析流程。本文提供了一种高度可重复且用户友好的上游分析流程及后续差异基因分析流程，以实现对私有或公共数据的标准化处理与深度解析。该流程分为四个步骤：（1）数据质量控制；（2）基因比对；（3）差异基因分析；（4）功能分析。该流程旨在揭示疾病转变的分子机制，并通过批量RNA-seq数据分析，帮助研究人员筛选潜在的药物靶点和治疗策略。

引言

非酒精性脂肪性肝病（NAFLD）是全球最常见的慢性肝病，影响超过四分之一的人口。近几十年来，其发病率显著上升^1,2,3。该疾病的负担日益加重，尤其是其更严重的形态——非酒精性脂肪性肝炎（NASH），已成为重大的全球健康挑战，并带来沉重的经济负担⁴。NAFLD 的第一阶段为非酒精性脂肪肝（NAFL），常伴随炎症和纤维化，可进一步发展为 NASH。后者显著增加了进展为终末期肝病（包括肝硬化和肝细胞癌（HCC））的风险^5,6,7。HCC 的发病率和死亡率与 NASH 的增加密切相关^8,9，预计到 2030 年，NAFLD/NASH 将成为肝移植的首要指征¹⁰。然而，NAFLD 的临床进展具有高度异质性¹¹，严重阻碍了相关药物的研发¹²，因此，精确探究其涉及的分子机制显得尤为重要。

基于批量RNA测序（bulk RNA-seq）获取细胞组成信息可显著阐明多种疾病的发病机制。近年来，已有大量针对模式生物和人类的批量RNA-seq研究，用于揭示非酒精性脂肪性肝炎（NASH）进展过程中的基因表达差异^13,14,15，并识别可用于干预的新治疗靶点。根据批量RNA-seq分析结果，Xiong等人发现肝脏中的非实质细胞（NPCs）参与了细胞外基质形成和细胞黏附等过程，从而促进NASH的进展¹⁶。Li等人证明，肝细胞中肝源性Wilms'肿瘤1相关蛋白（WTAP）可调控异位脂质积聚和炎症反应，进而促进NASH的形成¹⁷。尽管批量RNA-seq分析是解析NASH机制的有力工具，但其结果对上游数据质量高度敏感。上游实验操作和分析流程的异质性可能严重损害数据的可靠性，从而掩盖真实的生物学信息，并干扰后续分析的准确性。因此，建立一套标准化的上游分析流程至关重要。

与单细胞RNA测序（scRNA-seq）相比，批量RNA测序（bulk RNA-seq）在实验设计和实际应用中具有若干显著优势。尽管scRNA-seq能够在单细胞水平上识别细胞异质性，并实现对细胞类型特异性转录特征的精确分析，但其成本高昂、数据处理复杂，且在检测低丰度转录本方面灵敏度有限¹⁸。相比之下，bulk RNA-seq具有更高的测序深度、更低的成本以及更高的样本通量，因此特别适用于群体水平的差异基因表达分析和分子机制的探索¹⁹。因此，在标准化分析流程的指导下，bulk RNA-seq仍然是研究复杂疾病分子基础的一种高效、经济且稳健的方法。

本方案专为源自人类组织且RNA完整性较高（RIN ≥ 7.0）并具有足够输入量RNA（每样本≥ 500 ng）的批量RNA测序（bulk RNA-seq）数据集设计。为确保比对和定量步骤的可靠执行，建议使用配备至少10核CPU、32 GB内存以及不少于200 GB可用磁盘空间的本地工作站。在满足上述要求的基础上，本方案提供了一套高效且易于使用的分析流程，包含详细的操作说明和标准化的参数配置，以满足研究人员分析大规模转录组数据的需求。

方案

为演示目的，本研究使用Lan Bai等人生成的公开可用数据集PRJNA1023502，以展示上游和下游分析的每一步流程²⁰。由于该数据集来源于开放获取的NCBI SRA数据库，因此无需额外的许可或伦理审批。请参见材料表以确认所有必需的软件及R包版本。公开可用数据集PRJNA1023502包含6个非NASH、6个NAFL和6个NASH肝组织RNA-seq样本。在本实验方案中，该数据集用于演示批量RNA-seq工作流程的全部步骤，包括从SRA数据库获取数据、质量控制（fastp）、比对（HISAT2）、定量（featureCounts），以及下游的差异表达分析和功能富集分析。

1. SRA 工具包安装

1. 访问 SRA Toolkit 官方网站并下载 3.2.1 版本。

2. 公共数据下载

1. 获取SRA编号。
   1. 在文章的补充材料²⁰、数据可用性部分，或通过在NCBI SRA数据库中搜索关键词来获取。
   2. 在终端中输入 prefetch <SRA编号> 以下载数据。

3. 基因计数矩阵的生成

1. 按如下说明定义路径和设置：
   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. 从官方 GENCODE 数据库（https://www.gencodegenes.org/human/）下载人类参考基因组（GRCh38 主装配）及相应的基因注释文件（gencode.v44，综合基因注释集）。
      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. 重命名 SRA 文件。
   1. 确保所有 SRA 文件均使用 .sra 扩展名，以便下游工具能够准确识别和处理。
   2. 执行以下命令：
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. 执行以下命令，将 SRA 文件转换为 FASTQ 格式：
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR -e 20; done
4. 基于参考基因组构建 HISAT2 索引。
   hisat2-build $REFERENCE $INDEX
5. 使用可视化进度监控处理 FASTQ 文件。
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. 对于双端测序数据，执行以下命令：
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. 使用 fastp 进行质量控制和过滤，并将结果输出至 $OUT_FASTP 目录。对于单端测序数据，执行以下命令：
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      对于双端测序数据，执行以下命令：
      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. 对每个样本，检查 HTML 报告以进行可视化质量评估，并查阅 JSON 报告获取结构化指标，包括读段质量、GC 分布、读段冗余度和接头污染情况。确保所有样本满足以下质量阈值：Q30 ≥ 85%，接头污染率 < 5%，且无异常 GC 模式。将所有输出文件存储在 $OUT_FASTP 目录中。
   4. 使用 HISAT2 进行序列比对，生成 SAM 文件并输出至 $OUT_HISAT2 目录。对于单端测序数据，执行以下命令：
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      对于双端测序数据，执行以下命令：
      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. 使用基于 GRCh38 主装配 FASTA 文件构建的索引，将经过质量过滤的读段比对至人类参考基因组。提供适用于单端和双端测序数据的示例命令。
   6. 将 SAM 文件转换为 BAM 文件，排序并建立索引。
      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. 为每个样本生成 SAM 文件，并使用 SAMtools 将其转换为已排序且建立索引的 BAM 文件。对于高质量的人类转录组样本（RNA 完整性数，RIN ≥ 8.0），确保总体比对率超过 85%；对于典型的批量 RNA-seq 样本（RIN ≥ 7.0），比对率 ≥ 70% 可视为可接受。
   8. 使用 featureCounts 进行基因定量。
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. 检查 featureCounts 为每个样本生成的制表符分隔输出文件（*.counts.txt）及其摘要报告（*.counts.txt.summary）。确保读段分配率达到人类批量 RNA-seq 的典型阈值 ≥70%；显著偏低的分配率可能提示链方向不匹配、注释问题或比对质量差。本非链特异性 RNA-seq 数据集使用参数 -s 0。对于链特异性文库，在命令中将 -s 0 替换为 -s 1 或 -s 2。
6. 按如下说明生成基因计数矩阵。
   1. 使用第一个样本的基因 ID 和计数初始化计数矩阵。
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. 迭代地将每个样本的计数添加到单一矩阵中。
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. 向计数矩阵添加包含样本 ID 的表头行。
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. 从 GTF 文件中提取基因长度（每个基因的外显子长度总和）。
      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. 原始计数矩阵处理与基因注释

1. R 环境初始化及加载相关 R 软件包。
   1. 使用 install.packages() 函数安装 tidyverse、ggplot2、ggrepel、RColorBrewer、ggridges、FactoMineR。其中，tidyverse 软件包用于数据处理与绘图；ggplot2 软件包用于数据可视化；ggrepel 软件包用于生成不重叠的文本标签；RColorBrewer 软件包用于提供配色方案；ggridges 软件包用于绘制山脊图；FactoMineR 软件包用于主成分分析（PCA）和多变量分析。
   2. 使用 BiocManager::install() 函数安装 biomaRt、DESeq2、clusterProfiler。其中，biomaRt 软件包用于从 Ensembl 获取基因注释信息；DESeq2 软件包用于差异表达分析；clusterProfiler 软件包用于功能富集分析。
   3. 读取步骤 2 中生成的原始计数矩阵，使用 biomaRt 将 Ensembl ID 映射为 HGNC 基因名称以用于后续分析，移除重复的 Geneid 列，并生成清洗后的计数矩阵。以原始计数矩阵（counts_matrix.csv）作为输入，其中行为 Ensembl 基因 ID，列为样本。将清洗后的计数矩阵保存为输出文件（clean_counts_SRA.csv），其中行为 HGNC 基因符号，列为样本。
      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. 基因表达定量

注意：详细脚本请参见补充文件 1。

1. 执行以下命令以计算每百万映射读段中每千碱基转录本的片段数（FPKM）和每百万转录本数（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. 样本聚类与差异可视化

1. 通过主成分分析（PCA）考察样本聚类情况，以评估数据质量。详细脚本请参见补充文件 1。
2. 进行PCA可视化时，请使用以下代码。
   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. 差异表达分析与结果可视化

注意：详细脚本请参见补充文件 1。

1. 执行以下命令以构建 DESeq2 数据集并设置阈值。建模前仅保留所有样本中总读数 > 1 的基因。使用以下显著性阈值筛选差异表达基因（DEGs）：校正后 p 值 < 0.05 且 |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. 使用火山图可视化差异表达基因。
   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. 进行功能富集分析与可视化

注意：详细脚本请参见补充文件1。

1. 可视化 KEGG 富集结果。
   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. 可视化 GO 富集结果。根据实验要求，通过修改弧分辨率（n=500）、颜色映射和分面布局来自定义可视化效果。
   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. 可视化基因集富集分析（GSEA）结果。调整通路数量、密度山脊缩放比例和图层透明度，以适应不同的效应量分布。
   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))

结果

批量RNA测序的上游分析流程如图1A所示。该流程在Linux平台上依次执行以下关键步骤：首先，使用fastp对原始测序数据进行严格的质控，以去除低质量读段和接头序列；随后，HISAT2将高质量读段比对至参考基因组，并由Samtools转换和排序比对文件；最后，FeatureCounts进行基因水平的定量分析，生成基因表达矩阵，为下游分析提供高质量输入数据。所得表达矩阵的后续处理与统计分析在R环境中进行，相关流程及所需软件包如图1B所示。分析所用数据来自一项已发表的研究，包括6个非NASH样本、6个NAFL样本和6个NASH样本²⁰（图1C）。非NASH对照样本来自不符合肝移植标准且无NAFLD或NASH的个体。

需要注意的是，由于 Bai 等人公开提供的数据未包含明确的批次信息（例如测序批次或文库制备日期），也未提供显著差异表达基因的可下载列表，本研究无法直...

讨论

大规模RNA测序数据分析是一项跨学科任务，涉及基因组学、生物信息学、统计学和计算机科学的整合。一个完整的分析流程包括多个上游和下游步骤，如原始数据预处理、质量控制、序列比对、基因水平定量、数据标准化、差异表达分析以及生物学解释。在这些步骤中，将原始测序读段准确转化为高质量的基因表达矩阵尤为关键，因为在上游处理过程中引入的任何错误都可能传递至所有下游的生物学结论中。因此，建立透明且标准化的上游分析流程对于提高转录组学研究的可重复性至关重要。

本方案提供了一种简化的、完全基于脚本的工作流程，整合了多种广泛使用的工具，如 fastp（用于读段修剪和质量控制）、HISAT2（用于可变剪接感知的比对）以及 featureCounts（用于基因水平的定量）。这些工具已在成熟的 RNA-seq 分析框架中得到广泛应用——包括源自 Tuxedo 的分析流程和普遍采用的实验方案流程——并经过了对准确性与效率的严格验证^21,22。在这些基础方法的基础上，该工作流程通过引入明确的文件处理...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
biomaRt	Bioconductor	2.64.0	来自 Ensembl 的基因注释
clusterProfiler	Bioconductor	4.16.0	功能富集分析
DESeq2	Bioconductor	1.48.1	差异表达分析
FactoMineR	AgroParisTech	2.11.0	主成分分析和多变量分析
fastp	OpenGene	1.0.1	FASTQ 数据的质量控制与过滤
FeatureCounts	Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research	2.0.0	对映射到每个基因的读段进行计数，用于基因表达定量
ggplot2	Posit	3.5.2	数据可视化
ggrepel	Kamil Slowikowski	0.9.6	避免重叠的文本标签
ggridges	Claus O. Wilke	0.5.6	绘制山脊图
HISAT2	Johns Hopkins University	2.2.1	将过滤后的高质量读段比对至参考基因组
R	R Core Team	4.5.0	用于数据计算、分析和可视化的环境
RColorBrewer	Erich Neuwirth	1.1.3	绘图用配色方案
samtools	Large Scale Genomics work stream	1.22.0	转换和处理 SAM 文件以实现高效检索与访问
SRA Toolkit	National Center for Biotechnology Information	3.2.1	从 NCBI SRA 数据库获取并预处理原始测序数据