07 多工具对比：DESeq2 / edgeR / limma-voom
bulk RNA-seq 差异分析有三个"标配"工具 —— DESeq2、edgeR、limma-voom。它们建模思路不同，但在绝大多数真实数据上结果高度一致。本章做一次 head-to-head 对比：同一份 airway 数据、同一条 design、同一个对比，跑三种工具，看结果是不是真的如传说那样一致。
三种工具的建模差异
| 工具 | 模型 | 特点 |
|---|---|---|
| DESeq2 | 负二项 GLM + 共享离散度估计 + Wald / LRT | 默认的归一化（median of ratios）稳；LFC shrinkage 是一大卖点 |
| edgeR | 负二项 GLM + 经验贝叶斯估计离散度 | TMM 归一化，glmQLFTest 的 F 检验统计学上最稳 |
| limma-voom | 先 voom 把 counts 变成 log-CPM + 权重，再走 limma（线性模型 + 经验贝叶斯） | 速度最快，适合大规模比较；对样本量大的数据尤其稳 |
几个经验：
- 样本量 < 10 对 10，差异不大，三者谁都能用
- 样本量 > 50，limma-voom 的速度优势很明显
- 要做复杂的 contrast（比如多因子设计），limma 的接口最灵活
- 要做单细胞数据的 pseudobulk DE，目前最主流的搭配是 edgeR
三条流水线怎么写
DESeq2、edgeR、limma-voom 的 API 有些差别，但核心都是"声明 design + 拟合模型 + 提取 contrast"：
# DESeq2
dds <- DESeqDataSetFromMatrix(cnt, colData, design = ~ cell + dex)
dds <- DESeq(dds)
res_deseq <- results(dds, contrast = c("dex", "trt", "untrt"))
# edgeR
dge <- DGEList(cnt) |> calcNormFactors(method = "TMM")
design <- model.matrix(~ cell + dex, data = coldata)
dge <- estimateDisp(dge, design)
fit_edg <- glmQLFit(dge, design)
qlf <- glmQLFTest(fit_edg, coef = "dextrt")
tt_edg <- topTags(qlf, n = Inf)$table
# limma-voom
v <- voom(dge, design, plot = FALSE)
fit_v <- lmFit(v, design) |> eBayes()
tt_lim <- topTable(fit_v, coef = "dextrt", number = Inf)
这三段 都用同一份 cnt counts 矩阵和 coldata 样本表。关键在让设计矩阵一致：都是 ~ cell + dex，都检验 dex 的 trt vs untrt 效应。
真实示例：airway 上的三工具对比
配套脚本 bulk07_tools_sci.R 在同一份 airway 数据上跑 DESeq2 / edgeR / limma-voom，做完后对比 6 个维度：
Rscript scripts/bulk07_tools_sci.R
图 1：每个工具的 p 值分布
三个工具各自的原始 p 值分布。健康的直方图都有"0 附近 peak + 剩余均匀"的形状。如果某个工具的直方图出现 U 型或者整 体左偏，说明它在这份数据上模型拟合不好。airway 是个"好数据"，三个工具都给出了教科书级的 p 值分布。
图 2：不同 padj 阈值下的 DE 基因数
同一个对比，三种工具在 padj < 0.001 / 0.01 / 0.05 / 0.1 四个阈值上的显著基因数。数值差距通常在 10% 以内 —— 这就是"三个工具大差不差"的具体含义。
值得注意的是 edgeR 的 QLF 检验（glmQLFTest）在统计学上是比 DESeq2 的 Wald 更严格的检验，所以 edgeR 的 DE 数经常比 DESeq2 略少。这不是 edgeR 少找到，而是它更严谨。
图 3：DE 基因集的重叠
三个工具的 DE 基因集（padj < 0.05）的重叠分布。最大的那一条"DESeq2 + edgeR + limma-voom"说明绝大多数 DE 基因被三个工具同时发现。只被某一个工具独有发现的基因数相对很少，通常是边界基因（接近阈值）。
实用经验：三个工具 都说是 DE 的基因，大概率是真的信号；只被一个工具说是 DE 的，要多看一眼是不是 QC 问题或边界效应。
图 4：log2FC 的散点对比
三个工具估计的 log2 fold change 两两散点，红色虚线是 y = x 的对角线。点紧贴对角线 = 三个工具对效应大小的估计几乎一致。
DESeq2 vs edgeR 的相关系数通常在 0.99 左右；DESeq2 vs limma-voom 略低一点（limma-voom 的 LFC 在低表达区会稍微更平），但整体仍 > 0.98。
图 5：p 值的 Spearman 排序相关
三个工具两两之间，对所有基因按 p 值排序的 Spearman 相关。三个相关系数都 > 0.95。意思是：如果你把"基因按显著性从高到低排名"当作输出（比如用来做 GSEA），三个工具给出的排名几乎完全一样。
这个指标比"DE 基因集重叠"更稳 —— 因为后者会被阈值人为影响，而排序相关和阈值无关。
图 6：工具独有基因的效应大小
最后一个角度：三种工具独有的那些基因有什么共同点？把每类（A only、B only、C only、三个都有）的 |log2FC| 画成小提琴图。
- Shared by all（三个都有）：效应大，中位 |LFC| 可能在 1.5 以上 —— 这些是强信号基因
- 某一个工具独有：效应小，大部分 |LFC| 在阈值附近 —— 边界基因
这张图最重要的解读：工具间的分歧主要发生在"勉强显著"的基因上。真正的强信号，三个工具都抓得到。
什么时候选哪个
实用建议（可以背下来）：
| 场景 | 推荐 |
|---|---|
| 样本量小（< 6 对 6）、新项目 | DESeq2，默认最稳 |
| 样本量大（> 30 对 30） | limma-voom，速度最快 |
| 要做 ANOVA-like 检验（LRT / 多因子交互） | DESeq2 的 LRT 或 limma 的 contrasts.fit |
| 单细胞 pseudobulk | edgeR，目前最成熟 |
| 要做严格的 FDR 控制（需要 FPR < 5%） | edgeR QLF，最保守 |
| 要做 LFC shrinkage 做下游 | DESeq2 配 apeglm |
都跑一遍是一个好习惯。在一份新数据上，先用三种工具各跑一遍，如果结果相似，报告其中一种；如果结果差很多，说明数据或设计有问题，要回头看 QC。
除了这三个
其他常见选择：
- NOISeq：非参数方法，对非常态分布友好
- sleuth：搭配 kallisto/Salmon 的 bootstrap 结果使用，能同时建模定量不确定性
- BASiCS：有 spike-in 的数据可用
- DREAM（variancePartition）：处理复杂混合效应模型（样本重复测量、配对设计）
刚上手时没必要学这些；先把 DESeq2 / edgeR / limma-voom 中一个摸熟。
下载资源
下一步
参考资源
离线资料下载
手册 HTML / PDF 已在后台预生成，点击后直接下载网站静态资源。