3. 常见的差异表达分析工具 (DESeq2, edgeR, limma)
聊到差异表达分析，圈里最常用的就是这三个工具：DESeq2、edgeR 和 limma。说实话，我刚入行那会儿也纠结过该选哪个。后来做过的项目多了，慢慢摸清了它们的脾气。下面我结合自己的踩坑经历，把这三个家伙给你们讲透。
一句话总结：bulk RNA-seq 差异表达的主流工具就这三板斧。DESeq2和edgeR走的是负二项分布路线，limma原本是为芯片设计的，后来也能用 voom 处理RNA-seq数据。
3.1 DESeq2 —— 稳健派的选择
DESeq2 是我个人最常推荐给新手的工具。为什么？因为它对异常值和低表达基因的处理相当保守，不容易瞎报差异。它的核心思路是：先估计每个基因的离散度，然后通过一个「收缩」步骤把极端值拉回来——说白了就是防止假阳性。
# R 代码示例 - DESeq2 基本流程
library(DESeq2)
# 假设你已经有 counts 矩阵和 colData
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = coldata,
design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treat", "control"))
# res 里有 baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
我的一点体会：DESeq2 的 counts 输入必须是整数，千万别提前做 FPKM 或者 TPM 归一化！它内部有自己的归一化步骤（median of ratios）。我曾经手贱做过一次，结果 padj 全是 NaN —— 血的教训。
3.2 edgeR —— 灵活派的利器
edgeR 和 DESeq2 是亲兄弟，都用负二项分布。但 edgeR 更「灵活」——它允许你自定义归一化方法，而且计算速度比 DESeq2 快不少。我做过一个 100 个样本的项目，DESeq2 跑了快半小时，edgeR 十分钟出头就出结果了。
# R 代码示例 - edgeR 基本流程
library(edgeR)
y <- DGEList(counts = counts, group = group)
y <- calcNormFactors(y) # TMM 归一化，可换成其他方法
y <- estimateDisp(y)
fit <- glmQLFit(y)
qlf <- glmQLFTest(fit, contrast = c(-1, 1))
topTags(qlf)
| 特性 | DESeq2 | edgeR | limma+voom |
|---|---|---|---|
| 统计模型 | 负二项 + 收缩 | 负二项 + 经验贝叶斯 | 线性模型 + 权重 |
| 归一化方式 | median of ratios | TMM (默认) | voom 内部隐含 |
| 对小样本的稳健性 | 很强 | 较强 | 中等 |
| 计算速度 | 较慢 | 快 | 快 |
| 支持多因素设计 | 是 | 是 | 是 |
3.3 limma + voom —— 老牌劲旅的跨界
limma 本来是做芯片差异表达的，属于线性建模的老祖宗。后来有人搞了个 voom 变换，可以把 RNA-seq 计数转换成类似芯片的信号，这样 limma 就能上场了。实际效果呢？我对比过几次，limma+voom 的结果和前面两个工具的重合度很高，而且它处理复杂实验设计（比如 time course）的能力特别强。
注意：limma+voom 要求 counts 矩阵里不能有太多零。如果有大量低表达基因，建议先用 filterByExpr 筛掉。否则 voom 的权重会出问题，假阳性飙升。我遇到过一回，筛完基因 padj 从 0.9 掉到 0.05 以下……
# R 代码示例 - limma+voom 基本流程
library(limma)
library(edgeR) # 用于 DGEList
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
v <- voom(dge, design = model.matrix(~ group))
fit <- lmFit(v, design)
fit <- eBayes(fit)
topTable(fit, coef = "grouptreat")
3.4 怎么选？我的个人经验
其实这三个工具你只要熟练一个就够了，结果不会差太多。但我建议你这样选：
- 样本量少（≤3个重复） —— 优先 DESeq2，它那个收缩步骤能救命。
- 样本量中等（4~6个） —— edgeR 更灵活，而且快。
- 实验设计很复杂 —— 比如有批次、协变量、多个分组，limma+voom 的线性模型最顺手。
- 你是芯片转 RNA-seq 的老人 —— 直接 limma+voom，别折腾了。
嗯，这里还要插一句：千万不要同时跑三个工具然后取交集。我见过有人这么干，觉得取交集更保险。但其实每个工具的假阳性机制不同，取交集只会漏掉大量真正的差异基因。用一个工具把参数调好，比三个工具拼凑靠谱得多。
避坑指南：我曾经在分析中忽略了批次效应，结果 DESeq2 和 edgeR 都给出了几十个显著基因，limma+voom 却一个都没有。后来发现是批次和条件完全混杂，导致 limma 的线性模型直接把批次效应吸收后，条件效应就没了。所以实验设计一定要先画清楚再跑代码。
3.5 共同输出与解读
不管用哪个工具，最后拿到的都是类似的一张表：
| gene | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
|---|---|---|---|---|---|---|
| GeneA | 1245.6 | 2.34 | 0.45 | 5.20 | 1.2e-07 | 3.5e-05 |
| GeneB | 78.2 | -1.87 | 0.62 | 2.99 | 0.0028 | 0.045 |
重点看 padj（校正后的 p 值）和 log2FoldChange。通常 padj < 0.05 且 |log2FC| > 1（就是倍数变化 2 倍以上）作为筛选标准。不过别死抠阈值，我以前做的一个项目里只筛到 30 个基因，后来放宽到 padj < 0.1 才找到关键的调控因子。数据说话，别教条。