RNA-seq 转录组分析流程

引言

RNA-seq 是应用最广的测序技术之一：它测量「细胞在某个状态下表达哪些基因、表达多少」。与基因组测序（测一次就够）不同，转录组是动态的——同一细胞在不同处理、不同时间点的表达谱差异巨大，这既是它的价值（能反映状态），也是它的挑战（必须用统计方法区分真实差异与噪声）。

RNA-seq 的分析难点集中在「统计」而非「算法」。上游的比对、定量相对标准化，但差异表达分析涉及大量统计概念：为什么必须做生物学重复？为什么不能用「表达量变化 2 倍」当唯一判据？为什么 p 值要校正？多重检验校正会不会太保守？这些问题没有工程直觉可依赖，必须理解统计原理。

工程上的陷阱也很典型：一是把技术重复当生物学重复，导致假阳性；二是混淆 CPM、TPM、FPKM 三种归一化，用错了导致结论错误；三是忽略批次效应，把批次差异当成处理效应；四是 p 值不校正，报告出成百上千个「显著」基因却大多不可复现。

本文按「实验设计 → 定量路线 → 表达矩阵 → 归一化 → 差异表达 → 统计陷阱 → 富集分析 → 剪接分析」的顺序展开。工具基于 STAR 2.7、Salmon 1.10、DESeq2 1.38。读完你应该能独立设计并完成一个 RNA-seq 分析，且理解每一步的统计依据。

目录

1. RNA-seq 与 DNA 测序的本质差异
2. 实验设计：重复、批次与深度
3. 比对路线与无比对定量
4. 表达矩阵的生成
5. 归一化：CPM、TPM 与 FPKM
6. 差异表达分析
7. 多重检验校正与统计陷阱
8. 功能富集分析
9. 可变剪接与融合基因

1. RNA-seq 与 DNA 测序的本质差异

RNA-seq 与 WGS/WES 的根本区别在于「测量对象」和「动态性」：

表达量的动态范围是 RNA-seq 的核心特征。一个细胞里，高表达基因（如肌动蛋白 ACTB）可能占所有 mRNA 的百分之几，而低表达基因（如某些转录因子）可能只有几个拷贝。同一份数据里，表达量可以跨越 5-6 个数量级。这意味着：

• 测序深度被高表达基因「吃掉」：高表达基因占据大量读段，留给低表达基因的读段很少。想检测低表达基因的差异，需要更深度的测序。
• 归一化至关重要：样本间的总读段数不同、高表达基因占比不同，直接比较原始计数是错误的。这就是 CPM/TPM 等归一化方法存在的理由。
• 零膨胀（zero inflation）：低表达基因在部分样本里计数为 0，这不是「不表达」，而是「测序没测到」。统计模型必须处理这种稀疏性。

另一个区别是 RNA 的剪接：mRNA 成熟时内含子被切除，所以 RNA 读段在基因组上可能跨越外显子边界（跨内含子比对）。这要求比对工具具备剪接感知能力，也是 RNA-seq 必须用 STAR/HISAT2 而非 BWA 的原因。

2. 实验设计：重复、批次与深度

实验设计决定了分析的统计功效（power），事后无法弥补。三个核心要素：

生物学重复（biological replicates）：独立生物学样本（不同个体、不同培养皿）。至少 3 个，推荐 4-5 个。重复的作用是估计「组内变异」——没有组内变异，就无法判断「组间差异」是否显著。技术重复（同一文库测多次）不能替代，因为它不包含生物学变异。

批次（batch）：不同日期、不同操作员、不同试剂批次的系统性差异。设计铁律是批次平衡：处理组和对照组要均匀分布在各个批次，不能「一批全是对照、一批全是处理」。如果批次与处理完全混淆（confounded），统计方法无法分离两者。

错误设计（批次与处理混淆）：
  批次1: 对照1 对照2 对照3
  批次2: 处理1 处理2 处理3
  → 无法区分「处理效应」与「批次效应」

正确设计（批次平衡）：
  批次1: 对照1 处理1 对照2
  批次2: 处理2 对照3 处理3
  → 处理与批次可分离

测序深度：RNA-seq 的深度以「读段数」计。常规差异表达每个样本 20-30M 读段（PE100）；检测低表达基因或做剪接分析需要 50-100M。深度不足会导致低表达基因计数为 0，降低检测灵敏度。

设计完成后，统计功效可以预估（power analysis）。工具如 RNASeqPower 或 pwr 包能计算「给定重复数，能检测到多大倍数的差异」。3 个重复通常能可靠检测 2 倍以上的差异，检测 1.5 倍差异需要更多重复。

3. 比对路线与无比对定量

RNA-seq 的定量有两条技术路线：

比对路线（alignment-based）：先把读段比对到基因组（剪接感知），再统计每个基因/外显子的读段数。

# STAR 比对（剪接感知，需要先建索引）
STAR --runMode genomeGenerate --genomeDir star_idx \
  --genomeFastaFiles ref.fa --sjdbGTFfile annot.gtf \
  --runThreadN 16

STAR --runThreadN 16 --genomeDir star_idx \
  --readFilesIn R1.fq.gz R2.fq.gz --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate --outFileNamePrefix sample_

# 计数（featureCounts）
featureCounts -T 8 -p -a annot.gtf -o counts.txt sample_Aligned.sortedByCoord.out.bam

无比对路线（alignment-free）：直接把读段比对到转录本序列，用概率模型估计每个转录本的丰度。

# Salmon 定量（快，10 倍于比对）
salmon index -t transcripts.fa -i salmon_idx -k 31
salmon quant -i salmon_idx -l A \
  -1 R1.fq.gz -2 R2.fq.gz -p 8 \
  --validateMappings -o sample_quant

两条路线的对比：

选择逻辑：只做差异表达 → Salmon（快且准）；需要剪接分析、变异检测、可视化 → STAR + featureCounts。两者在基因水平定量的相关性通常 > 0.98，所以常规差异表达用哪个都行。

4. 表达矩阵的生成

无论哪条路线，最终都汇总成「表达矩阵」：行是基因，列是样本，值是计数。

gene_id          ctrl_1  ctrl_2  ctrl_3  treat_1  treat_2  treat_3
ENSG00000000003  1234    1100    1300    2450     2380     2500
ENSG00000000005  0       2       1       0        1        0
ENSG00000000007  5678    5600    5800    1200     1150     1250

矩阵生成的几个要点：

• 用原始计数（raw counts）进入统计：DESeq2/edgeR 需要整数计数，不要用已归一化的值。归一化在统计模型内部完成。
• 基因 ID 一致：确保所有样本用同一套注释版本，否则基因对应不上。
• 过滤低表达基因：计数全为 0 或极低的基因会干扰统计（增加多重检验负担）。通常过滤「至少在一半样本中计数 ≥ 10」的基因。

import pandas as pd
# 读入表达矩阵并过滤
counts = pd.read_csv('counts.txt', sep='\t', index_col=0, skiprows=1)
counts = counts.iloc[:, 5:]           # 去掉 featureCounts 的注释列
keep = (counts >= 10).sum(axis=1) >= (counts.shape[1] / 2)
filtered = counts[keep]
print(f"过滤后保留 {filtered.shape[0]} / {counts.shape[0]} 个基因")

用 pandas 处理表达矩阵是标配技能，参见 pandas 数据分析 。矩阵的行列对齐、缺失值处理、分组聚合都依赖它。

5. 归一化：CPM、TPM 与 FPKM

归一化（normalization）解决「样本间可比性」问题。三种常用指标，用途完全不同：

CPM 只校正测序深度：CPM = 计数 / 总计数 × 1e6。它不校正基因长度，所以不能比较不同基因的表达（长基因天然有更多读段）。

FPKM/RPKM 同时校正深度与基因长度：FPKM = 计数 / (基因长度kb × 总计数百万)。但它有个致命缺陷——样本间不可比。因为每个样本的总数归一化后，各基因的 FPKM 之和不为常数，导致不同样本的 FPKM 不能直接比较。这就是它被 TPM 取代的原因。

TPM 修正了 FPKM 的缺陷：先除以基因长度得到「每碱基读段率」，再归一化到百万。这样每个样本的 TPM 之和恒为 1e6，样本间可比。

TPM 计算：
  1. rate_i = counts_i / length_i
  2. TPM_i = rate_i / sum(rate) × 1e6
  → 所有基因的 TPM 之和 = 1e6（跨样本可比）

重要提醒：差异表达分析不要用 TPM/FPKM。DESeq2/edgeR 的统计模型基于「原始计数」的离散分布（负二项分布），它们内部有自己的归一化（如 DESeq2 的 median-of-ratios）。喂给它们 TPM 值会破坏统计假设。TPM 适合可视化、聚类、跨样本展示，不适合差异检验。

6. 差异表达分析

差异表达（Differential Expression, DE）是 RNA-seq 的核心分析。主流工具是 DESeq2 和 edgeR，都基于负二项分布模型：

library(DESeq2)
# 构建 DESeqDataSet
dds <- DESeqDataSetFromMatrix(
  countData = counts, colData = coldata,
  design = ~ batch + condition)      # 把批次纳入模型
dds <- dds[rowSums(counts(dds)) >= 10, ]

# 运行分析
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treat", "ctrl"),
               alpha = 0.05, lfcThreshold = 0)  # 显著性水平
res <- res[order(res$padj), ]

# 导出
write.csv(as.data.frame(res), "DE_results.csv")

关键输出字段：

判读标准：通常用 padj < 0.05 且 |log2FoldChange| > 1（即变化 ≥ 2 倍）。但**「2 倍」是约定而非真理**——它取决于数据的噪声水平和生物学意义。噪声大的数据可能要求 3 倍，而精心设计的实验能可靠检测 1.5 倍。

为什么用负二项分布？ 因为 RNA-seq 计数的方差大于均值（过离散，overdispersion）。泊松分布假设方差等于均值，会低估变异、产生过多假阳性。负二项分布引入了额外的离散参数，更符合实际。这也是为什么不能用简单的 t 检验或卡方检验。

收缩（shrinkage）：DESeq2 的 lfcShrink 会把低表达基因的 log2FC 向 0 收缩（因为它们估计不稳定），减少假阳性。这在排序和可视化时特别有用。

7. 多重检验校正与统计陷阱

RNA-seq 要同时检验上万个基因，多重检验问题无法回避。如果对每个基因用 p < 0.05，检验 20000 个基因就会产生约 1000 个「假阳性」（纯随机）。校正方法：

BH 法是 RNA-seq 的标准，padj 就是 BH 校正后的 p 值。它控制「被判定为显著的基因中，假阳性比例不超过 5%」——这个目标比 Bonferroni 的「完全不犯错」更实用。

常见统计陷阱：

1. 不做重复或把技术重复当生物学重复：无法估计组内变异，任何差异都「显著」；
2. 用 TPM 做差异检验：破坏负二项模型假设，p 值不可信；
3. 忽略批次效应：把批次差异当处理效应，需要把批次写进 design 公式；
4. 用原始 p 值而非 padj：报告大量假阳性；一律用 padj；
5. 只按 log2FC 排序：忽略统计显著性，低表达基因的极端 log2FC 往往是噪声；
6. p 值做「反向筛选」：用 p 值筛掉「不显著」的基因再分析，是数据窥探（p-hacking）。

批次校正 的两种时机：设计阶段（批次平衡，最有效）和统计阶段（把批次作为协变量，design = ~ batch + condition）。事后校正工具（如 ComBat、limma removeBatchEffect）只在无法设计平衡时使用，且有风险——它们可能过度校正掉真实的生物学信号。

8. 功能富集分析

差异表达给出「哪些基因变了」，富集分析回答「这些基因共同参与了什么功能」。三类方法：

ORA（Over-Representation Analysis）：给定差异基因列表，看哪些功能条目（GO term、KEGG pathway）里差异基因比例异常高。

library(clusterProfiler)
# GO 富集
ego <- enrichGO(gene = deg_genes, OrgDb = org.Hs.eg.db,
                keyType = "ENSEMBL", ont = "BP",
                pAdjustMethod = "BH", pvalueCutoff = 0.05)
# KEGG 富集
ekegg <- enrichKEGG(gene = deg_genes, organism = "hsa",
                    pvalueCutoff = 0.05)

GSEA（Gene Set Enrichment Analysis）：不需要「切阈值选差异基因」，而是用全部基因的排序（按 log2FC 或统计量）来检验功能集是否集中在排序的顶部或底部。它更敏感，能发现「整体轻微上调」的功能。

library(fgsea)
ranks <- sort(setNames(res$stat, rownames(res)), decreasing = TRUE)
fgsea_res <- fgsea(pathways = pathways, stats = ranks, minSize = 15)

GSVA（Gene Set Variation Analysis）：把「功能集」的活性量化为每个样本的一个分数，然后可以比较样本间功能活性的差异。适合分析功能层面的连续变化。

富集分析的常见陷阱：背景基因集选择错误。ORA 需要指定「背景」（universe），如果背景用了全部基因而非「实际检测到的基因」，会高估富集显著性。此外，GO 条目的层级结构（父子关系）会导致「父条目和子条目同时显著」，需要做冗余过滤。

9. 可变剪接与融合基因

RNA-seq 不只测表达量，还能分析剪接与融合：

可变剪接（alternative splicing）：同一个基因产生多个转录本。分析工具：

# rMATS：比较两组间的剪接差异
rmats.py --b1 treat_bams.txt --b2 ctrl_bams.txt \
  --gtf annot.gtf -t paired --readLength 150 \
  --nthread 8 --od rmats_out

# 或用 Salmon + 转录本水平定量（更简单）
salmon quant -i salmon_idx -l A -1 R1.fq.gz -2 R2.fq.gz -o quant
# 再用 tximport 汇总到基因/转录本

剪接分析的关键指标是 PSI（Percent Spliced In）：某个剪接事件中，包含特定外显子的转录本占比。PSI 的变化反映剪接调控。

融合基因（gene fusion）：两个基因的序列拼接在一起，常见于肿瘤（如 BCR-ABL、EML4-ALK）。检测工具：

# STAR-Fusion：基于 STAR 的嵌合比对
STAR-Fusion --genome_lib_dir ctat_lib \
  --left_fq R1.fq.gz --right_fq R2.fq.gz \
  --output_dir fusion_out

融合检测依赖比对工具报告「嵌合读段」（一条读段两端比对到不同基因）。这需要比对时启用嵌合检测（STAR 的 --chimSegmentMin），且对读长和深度有要求。

剪接与融合分析的数据要求比常规差异表达高：需要更深的测序（50-100M 读段）和更长的读长（PE150），因为跨剪接位点的读段数量有限。

权衡取舍

常见坑清单

1. 无生物学重复：无法估计组内变异，所有基因都「显著」；至少 3 个生物学重复。
2. 用 TPM 做差异检验：破坏负二项模型；DE 必须用原始 counts。
3. 混淆 CPM/TPM/FPKM：跨基因比较用错指标；记住 TPM 样本间可比、FPKM 不可比。
4. 批次与处理混淆：无法分离批次效应；设计阶段做批次平衡。
5. 用原始 p 值：报告大量假阳性；一律用 padj（BH 校正）。
6. 忽略低表达基因过滤：增加多重检验负担、降低功效；过滤低计数基因。
7. 只按 log2FC 排序：忽略显著性；结合 padj 与 FC 双阈值。
8. 富集背景集错误：用全部基因而非检测到的基因，高估显著性；指定正确 universe。
9. RNA-seq 用 BWA 比对：不支持剪接，大量未比对；用 STAR/HISAT2。
10. 忽略链特异性：链特异性文库用错参数导致反义链污染；确认建库类型设 -s 参数。

小结

RNA-seq 的分析价值与统计严谨性成正比。工具会跑通，但结论是否可信取决于实验设计（重复、批次平衡）、归一化正确性（counts 用于统计、TPM 用于展示）与统计方法（负二项模型、FDR 校正）。这些不是「细节」，而是「结论成立的前提」。

工程上，RNA-seq 最该建立的认知是「统计先行」：先想清楚要检测多大差异、需要多少重复、如何平衡批次，再动手做实验。事后用统计方法补救设计缺陷，效果有限且风险高。另一个要点是「区分展示与检验」：TPM、聚类、热图是展示工具，差异检验必须回到原始计数与合适的模型。

下一步可以看 单细胞测序数据分析 了解细胞水平的转录组分析，或 变异检测与 VCF 处理 了解 DNA 层面的分析。如果你的差异表达结果「显著基因太多」，先检查是否用了正确的归一化和 padj，再考虑是否过滤了低表达基因。