RNA-seq差异表达分析：DESeq2、edgeR与limma-voom全解析

拿到featureCounts输出的count矩阵那一刻RNA-seq分析的下半场就开始了。前面我们跑完了比对、定量手里这张基因×样本的表达矩阵只是一堆数字真正要回答的问题是处理组和对照组之间哪些基因的表达量显著变了、变了多少、方向如何。这套差异表达分析业内基本绕不开三个R包——DESeq2、edgeR、limma。这三个包我在过去几年里的每个转录组项目几乎都会用到今天就完整拆一遍它们各自的统计原理是什么、该怎么用、输入数据要准备成什么样以及跑完结果不一致时怎么处理。1. 三大工具的分工DESeq2、edgeR、limma各管哪一段很多刚接触RNA-seq的同学会问一个很实在的问题这三个包都是做差异表达的到底选哪个我的答案是与其纠结选哪个不如先搞清楚它们各自的出身和定位。DESeq2和edgeR是专门为RNA-seq的count数据量身定做的。RNA-seq的原始定量结果是整数计数比如某个基因在样本A里比对到100条reads、在样本B里比对到380条reads。这类数据的特点是离散且方差大直接用普通的t检验或线性模型会出大问题。DESeq2和edgeR都基于负二项分布建模核心思路是用离散分布去拟合计数数据的波动。它们的差别在于内部参数估计的方式DESeq2用中位数比值法估计标准化因子用经验贝叶斯收缩离散度edgeR则用TMM标准化 经验贝叶斯估计离散度两个包思路相近但细节不同。limma的出身要古老一些它最早是给基因表达芯片microarray设计的一套线性模型工具核心方法是经验贝叶斯调节的t检验。芯片数据是连续信号limma在里面表现极好但RNA-seq计数数据不符合正态假设。后来limma团队提出的voom方法改变了局面先把count数据转换成log2-CPM每百万reads的log2计数再用lowess拟合均值-方差关系为每个基因赋予一个精度权重之后直接用线性模型作答。这等于给RNA-seq数据穿上了一件芯片数据的外衣把计数问题转化成大家都很熟悉的线性模型问题。从应用场景上区分如果你的实验设计相对简单就是两组比较DESeq2和edgeR都是稳妥选择如果你要从头建模、有多个分组、要加入协变量比如批次、性别、年龄或者想统一使用一套线性模型的灵活性limma-voom会更顺手。DESeq2当然也支持多因素设计但limma在矩阵设计上的自由度确实更高。实际项目里我从来不做非此即彼的选择题而是把三个工具都跑一遍。原因很简单三个工具算法不同、假设不同、细节处理不同同时显著的结果更容易在审稿人那里站得住脚。后面我会专门说结果怎么交汇。2. 算法本质负二项分布、经验贝叶斯与均方差模型到底在算什么理解一个工具最忌讳的就是只记代码不记原理。代码会过时但统计思路是通用的。我把三个工具的核心算法拆开讲一遍尽量说人话。先回答一个基础问题为什么RNA-seq的count数据不能直接用泊松分布描述泊松分布的经典假设是均值等于方差。对测序技术本身来说技术重复下同一个基因的reads计数确实接近泊松分布但生物学重复引入了样本间个体差异导致count数据的方差远大于均值这就是所谓的过离散overdispersion。如果硬用泊松模型会把生物学波动误判成真实差异产生大量假阳性。负二项分布多加了一个离散度参数允许均值小于方差正好能同时覆盖技术波动和生物学波动。DESeq2的流程是这样的第一步用每个样本所有基因的几何均值构造一个参考样本然后计算每个基因的表达量相对参考的比例取中位数得到每个样本的size factor。这个size factor就是你熟悉的上桌时每人分到的菜量不一样——有的样本测序深度高、有的低size factor就是用来把大家拉到同一饭量的。第二步估计每个基因的离散度方式是用样本均值去拟合离散度的经验贝叶斯先验分布再把每个基因的离散度向全局趋势收缩。这一步对低表达基因尤其重要因为低表达基因的样本间波动往往不可靠收缩可以让估计更稳。第三步拟合负二项广义线性模型用Wald检验比较两组差异输出log2 fold change和经过Benjamini-Hochberg校正的p值。此外DESeq2还会做一步independent filtering自动滤掉那些几乎不表达、不参与检验的基因这也是它经常报告的差异基因数量看起来比edgeR略少的原因之一。edgeR的思路和DESeq2是同源异流。它同样用负二项分布建模差别集中在标准化和离散度处理上。edgeR默认的标准化方法是TMMtrimmed mean of M-values思路是选一个参考样本然后对每个样本取其与参考样本的M值log2倍数和A值平均表达的加权均值经过截尾后估计出一个缩放因子。相比DESeq2的size factorTMM对高表达差异基因更不敏感两组共同高表达的基因不会把整个样本的比例带偏。离散度方面edgeR先估计一个共同的离散度再通过经验贝叶斯方法把每个基因的离散度向共同值收缩得到基因特异的离散度。检验环节有两套一是经典的exactTest适用于没有复杂协变量的两组比较二是更推荐的glmQLFTest适用于多组设计或多因素模型它先做负二项广义线性模型拟合再对每个基因做似然比检验或准似然F检验。limma-voom是完全不同的思路。voom先把count矩阵转成log2-CPM然后对表达量排序用lowess曲线拟合均值-方差关系发现低表达基因的方差大、高表达基因的方差小。voom会为每个基因、每个样本算出一个precision weight相当于在告诉线性模型这个数据点的可靠性有多高方差大的给低权重方差小的给高权重。之后用加权最小二乘拟合线性模型再走limma标志性的经验贝叶斯收缩——把所有基因的方差估计向一个公共值收缩得到moderated t统计量。这一步其实是limma的精髓它利用了全基因组数万个基因的信息来增强单个基因的方差估计哪怕某个基因只有两个重复也能获得一个相对稳定的方差估计而不是直接用那两个重复算出来的极不稳定的样本方差。下面这个表可以帮你快速记住三者差异工具数据假设标准化方法核心检验最大优势DESeq2负二项分布size factor中位数比值法Wald检验 / LRT小样本下离散度估计稳健自动过滤低表达基因edgeR负二项分布TMMexactTest / 准似然F检验低表达基因处理精细经典老牌limma-voom均值-方差加权线性模型TMM voom权重经验贝叶斯moderated t灵活度高支持复杂协变量计算速度快看到这里你应该明白这三个工具不是在谁算得更准上有本质差别而是在如何估计数据中的噪声上走了三条路。你手里那批批次效应明显、样本量不大不小的真实数据最终谁表现好真的只有跑过才知道。3. 同一份count数据三个工具各自要求的输入长什么样我见过太多人在这一步卡住代码复制下来了数据读进去报错原因全是输入格式没对齐。三个工具的共同起点是一张count矩阵但各自对列名、行名、样本信息表的格式都有要求我逐个说并且附上我踩过坑之后总结出的检查清单。你需要准备的核心文件有两份第一份是count矩阵通常来自featureCounts、HTSeq或STAR生成的基因定量结果。要求是行名为基因名Ensembl ID或SYMBOL均可列名为样本名单元格为整数计数。所有样本的基因必须完全一致、顺序可以不同反正读取时按行名对齐。featureCounts输出的文本文件里会有几列注释信息Geneid、Chr、Start、End等真正要读入的只是后几列样本计数。第二份是样本信息表colData描述每个样本所属的分组和你想纳入模型的协变量。它是data.frame行名必须与count矩阵的列名完全一一对应。这个完全一一对应是血泪教训的重灾区count矩阵列名叫Ctrl_1样本信息表行名叫ctrl1看起来没啥问题但R会认为它们是不同的样本直接报错或者更糟——静默地错位。我现在每次跑分析前固定做一步all(colnames(countData) rownames(colData))的检查返回TRUE才继续。group这一列最好直接设置成factor并且要主动指定水平顺序。R默认的factor水平是按字母序排列的假设你的处理组叫Treat、对照组叫Ctrl字母序会把Ctrl当成第一个水平做出来的log2FC方向就会反过来。别觉得这是小事我见过有人把处理组和对照组跑反了上游富集分析全做在错误方向上的案例。对DESeq2来说它还需要知道设计公式design。常见的设计是~ condition意思是表达量只受分组影响。如果你有批次信息设计公式要写成~ batch condition。注意DESeq2的设计公式不包含截距项这一点和limma正好相反——写代码的时候别把习惯带乱。对edgeR和limma-voom来说输入方式反而更灵活。edgeR先读入DGEList对象group参数可以是分组的因子如果要做多组比较需要自己构建设计矩阵model.matrix(~ condition)。limma则完全基于设计矩阵工作group变量要手动转成model matrix。第三个容易被忽略的是基因注释。三个工具本身都不需要注释文件但如果你计划把Ensembl ID换成基因名做后续分析建议在DEG分析前就准备好ID映射表。我习惯用biomaRt包做注释或者直接用clusterProfiler自带的bitr函数一步到位转成SYMBOL和ENTREZID。最后提醒一点实际数据分析中我通常还会把基因长度和GC含量存下来备用不是DEG必需但如果后续要做可视化或某些特定分析比如用TPM值画热图有这些信息会省很多事。4. 实战代码一套数据跑完DESeq2、edgeR、limma-voom全流程理论讲完直接上实战。我模拟一个最常见的实验设计6个样本3个正常对照Ctrl_1、Ctrl_2、Ctrl_33个药物处理Treat_1、Treat_2、Treat_3。假设你已经有了一个count矩阵文件counts.txt行是基因列是样本整数计数。4.1 读入数据与QC检查# 读取count矩阵 countData - read.table(counts.txt, header TRUE, row.names 1, check.names FALSE) # 读取样本信息表 colData - read.table(colData.txt, header TRUE, row.names 1, check.names FALSE) # 强制检查行名是否完全一致 all(colnames(countData) rownames(colData)) # 把分组转为factor并明确指定水平 colData$condition - factor(colData$condition, levels c(Ctrl, Treat))在跑差异分析之前强烈建议先做一次样本层面QC。我最常用的方式是PCA看看处理后样本是否和对照样本在PC1或PC2方向明显分开顺便看有没有离群样本。如果对照组里有一个样本跑到处理组那边去了别急着跑DEG先回去查这个样本的建库或测序情况。4.2 DESeq2完整流程library(DESeq2) # 建立DESeqDataSet对象 dds - DESeqDataSetFromMatrix(countData countData, colData colData, design ~ condition) # 预过滤保留至少在3个样本中count 10 的基因 keep - rowSums(counts(dds) 10) 3 dds - dds[keep, ] # 运行差异分析 dds - DESeq(dds) # 提取结果指定比较方向和阈值 res - results(dds, contrast c(condition, Treat, Ctrl), alpha 0.05) summary(res) # 转成data.frame方便输出 res_df - as.data.frame(res) res_df$gene - rownames(res_df)三个细节需要说明contrast参数的顺序非常重要c(condition, Treat, Ctrl)表示计算的是Treat相对于Ctrl的log2FC正值代表在Treat组上调。alpha参数是独立的FDR阈值不仅仅是提取结果的阈值results()会依据它重新进行独立过滤。如果你不只是想看某一对比较而是想看处理组内多个时间点可以用results()配合lfcThreshold或者用DESeq里的LRT检验做多水平测试。DESeq2我一直觉得最省心的是它对低表达基因做了自动处理——DESeq()内部会做独立过滤输出结果里NA的基因通常就是被滤掉的低表达基因。它另一个好处是可以直接做log2FC的收缩lfcShrink对排名基因列表、后续做排序GSEA特别有用。我自己的规则是报告用的log2FC用原始值做基因排序或可视化时用shrink值。4.3 edgeR完整流程edgeR有两个层次的分析路线经典路线exactTest适合简单两组比较更通用的是广义线性模型路线glmQLFitglmQLFTest支持多因素设计。我推荐直接用后者不仅灵活而且对复杂设计不容易出错。library(edgeR) # 建立DGEList对象 dge - DGEList(counts countData, group colData$condition) # 低表达过滤edgeR自带filterByExpr非常省事 keep - filterByExpr(dge) dge - dge[keep, , keep.lib.sizes FALSE] # TMM标准化 dge - calcNormFactors(dge, method TMM) # 构建设计矩阵注意这里可以加入批次等协变量 design - model.matrix(~ colData$condition) # 估计离散度 dge - estimateDisp(dge, design) # 拟合广义线性模型并进行准似然F检验 fit - glmQLFit(dge, design) qlf - glmQLFTest(fit, contrast c(0, 1)) # 提取结果 topTags(qlf, n 20)这里有个容易踩的坑glmQLFTest里的contrast参数顺序和design矩阵的列对应。我的design是model.matrix(~ condition)默认生成两列第一列是截距所有样本的基线第二列编码Treat与Ctrl的差异所以contrast c(0, 1)代表提取第二列的系数。如果你加入了其他协变量设计矩阵的列会变多contrast也要随之调整。我通常会在跑之前用colnames(design)确认一下列名再决定contrast怎么写这比凭记忆瞎写稳妥得多。keep.lib.sizes FALSE这一步的作用是在过滤低表达基因后重新计算library size。如果不加这步过滤掉的基因仍然会占据library size的比例标准化的SCF会偏高影响后续的CPM计算虽然差异基因结果影响不大但会有细节偏差。4.4 limma-voom完整流程library(limma) # 复用edgeR的DGEList对象 # 如果没有可以用 DGEList(counts countData, group colData$condition) # 设计矩阵 design - model.matrix(~ colData$condition) # voom转换加权重 v - voom(dge, design, plot TRUE) # 线性模型拟合 fit - lmFit(v, design) fit - eBayes(fit) # 提取结果 topTable(fit, coef 2, number 20, sort.by P)voom的plot TRUE会画一张均方差关系图我建议每次跑都瞄一眼图中的曲线应该在低表达端有明显抬升说明低表达基因方差大voom正在给它们分配较小的权重。如果曲线平直的像一条直线说明你的数据可能已经预处理过比如已经做了某种标准化要回头查一查。limma的优势在这种场景还不明显等你的实验有批次效应时就体现出来了。把设计矩阵改成下面这样即可design - model.matrix(~ batch condition)batch列变量可以是因子limma天然支持。DESeq2也有类似机制设计公式里加batch但limma在线性模型框架下处理协变量更自然对连续型协变量比如RNA Integrity Number也支持得很好。5. 结果比较三个工具跑出不同的基因列表怎么调和跑完三个工具最常见的困惑来了三个包报告的显著差异基因数量不同基因列表也在重叠之外各有差异到底以哪个为准先讲为什么会有差异。以我的经验数量差别主要来自三点一是低表达基因的处理策略不同DESeq2的独立过滤会把一部分本来就不表达的基因排除在检验之外edgeR的filterByExpr则更偏向保留那些虽然低表达但信噪比高的基因二是对离散度的估计方式不同间接影响了p值的大小三是标准化步骤的细节不同TMM和size factor计算出的倍数变化本身就有细微差别。举个例子我最近处理的一个植物胁迫转录组项目3个对照 3个处理在FDR0.05、log2FC绝对值1的阈值下DESeq2筛出2143个edgeR筛出2381个limma-voom筛出2529个三者的交集是1729个基因。这个格局非常典型limma-voom通常报出的基因数最多DESeq2最少edgeR居中。原因不难理解limma-voom的经验贝叶斯收缩利用了全基因组信息让许多本来方差估计不稳的基因也能得到较小的p值因此检出力更强DESeq2的收缩更保守假阳性控制更严格。怎么调和我的工作流是这样的如果只是要一个可靠的DEG列表用于后续功能富集我会优先取三者交集再把DESeq2的结果单独导出一份用于下游GSEA排序。理由很实际交集是三个工具共同支持的铁板钉钉的显著基因拿去富集很少有争议而GSEA这类排序型分析对log2FC的绝对值更敏感DESeq2自带收缩后的log2FC排序效果好天然适合。这三个工具分别跑出的log2FC其实高度一致假如你发现某个基因在DESeq2里显著上调、在edgeR里却是显著的负值这基本说明出bug了不是统计差异而是对比方向设反了。遇到这种情况回去查各自的contrast定义特别是factor水平顺序。下面这份表格是我习惯用来汇总三个工具结果对比的格式你也可以照这个思路做工具显著上调基因数显著下调基因数交集内上调交集内下调DESeq211201023896833edgeR12511130896833limma-voom138611438968336. 真实项目里的坑低表达过滤、批次效应和对比方向这一节我想集中说说实战里最常翻车的几个问题每个都是我亲眼见过或亲手埋过的雷。低表达基因过滤这事看起来不起眼影响却很大。如果不做任何过滤几万个基因里有大量在所有样本中计数都是0或极低的基因它们会带来两个麻烦一是离散度估计会受到极端值干扰尤其是那些偶尔冒出一个read的基因方差被低估或高估直接影响p值二是多重检验校正时上万个基因的检验次数会增加FDR的压力把真正显著的基因挤下去。当然像DESeq2有内置的独立过滤可以帮忙兜底但我仍然建议在三个工具之前做一次统一过滤——至少保证三个工具处理的是同一批基因比较才有意义。常见做法是保留至少在最小样本组中有一定计数的基因filterByExpr就是一个不错的通用标准。批次效应是个大坑。RNA-seq实验往往分多批建库测序批次之间会有系统性偏差如果不控制差批次可能被检测为差异。处理批次效应的方式有几层最理想的是实验设计时做随机化把各个分组均匀分布在批次里数据层面可以在模型里加协变量这是limma和DESeq2都支持的还可以用专门的工具比如ComBat-seq、RUVSeq、sva包做显式批次校正。我的经验是除非你打算做复杂的数据校正否则优先在模型里加协变量这样最不容易引入额外偏差。ComBat-seq这类方法虽然好用但会把数据过度拉平反而压低真实差异需谨慎。再说对比方向。前面反复强调factor水平顺序再补充一个相关的坑results()默认采用的比较方向是design公式里第二个水平的factor相对于第一个水平。如果你用的是DESeq2且levels写成c(Treat, Ctrl)那默认结果就是Ctrl相对于Treat反向的。所有工具都提供显式指定对比的参数我建议永远别用默认方向每次明确对比顺序宁可多打几个字。处理组和对照组跑反了后面所有下游分析和生物学解释全废这种事故一旦发生重跑成本极高。多分组比较是另一个容易混乱的场景。假设你有Ctrl、Low、High三个剂量组想找Low和High各自相对于Ctrl的差异。边缘情况还好但如果想找Low组独有、High组没有的差异基因许多人会直接把三者同时放进设计矩阵再跑结果发现用的对比系数根本说不清。常见的做法是两两比较Ctrl vs Low单独跑一次Ctrl vs High单独跑一次再把两个DEG列表做Venn或UpSet图找交集和差集。这么做代码多跑两遍但逻辑清楚得多下游解释也方便。还有一个我反复跟人说的问题生物学重复数量。三个工具在只有2个重复时都能计算但结果极不稳定。差异表达分析的本质是估计组内方差2个样本算方差能算但可信度很低经验贝叶斯收缩能拉一把却救不回设计上的先天不足。最少3个生物学重复是我能做建议的下限4-6个更舒服。如果你的重复样本其实来自同一只动物的不同测序文库那是技术重复统计上等于1个生物学重复用三个工具里的任何一个做组间比较都会得到一堆假阳性。7. 选型与组合我的推荐工作流文章最后给出一套我经过多个项目验证的工作流你可以直接拿去用。第一步比对定量拿到count矩阵后先做PCA和样本聚类用pheatmap画样本相关热图检查离群样本和批次结构。这一步不是DEG分析的一部分但决定后面所有分析的走向。如果PCA图里处理组和对照组完全混在一起要么是样本分组问题要么是处理效应太弱先别急着跑差异分析。第二步用统一的过滤标准整理基因列表保证三个工具分析的基因集一致。然后并行跑DESeq2、edgeR和limma-voom各自输出完整结果。第三步对比三个DEG列表。我一般先看三个列表的韦恩重叠率如果重叠率低于70%通常说明数据里有没处理干净的问题优先查批次效应和离群样本。重叠率正常时把交集基因作为可靠的DEG集合用于Venn、热图和功能富集。第四步根据你的实验设计选择主报工具。如果是标准的两组比较、3-5个重复我建议主报DESeq2它的独立过滤和log2FC收缩能帮你省掉很多解释上的麻烦审稿人也最眼熟。如果有明确的批次信息或其他协变量可以主报limma-voom毕竟线性模型处理协变量最顺手。edgeR更适合那些基因表达非常稀疏、大量基因在多个样本里计数极低的特殊场景它对低表达基因的过滤更精细。第五步差异表达只是起点。拿到DEG列表后接着做GO/KEGG富集clusterProfiler或者用fgsea做GSEA。GSEA对明确阈值的依赖小能利用所有基因的表达变化趋势我最近两个项目都靠它补出了单纯看DEG列表看不到的生物学信号。最后再分享一个我个人的习惯整套分析过程一定要留好脚本和sessionInfo。做生信分析的时间跨度长今天跑的结果可能三个月后要重跑或补分析R包版本一变结果就可能对不上了。我把每次分析的R版本、Bioconductor版本、三个工具版本都记录在案回头复盘时省了无数事。做RNA-seq差异分析工具是固定的数据是流动的真正拉开差距的永远是分析者对原理的理解和对数据的敏感度。希望这篇把DESeq2、edgeR、limma的底层逻辑和实践细节讲清楚之后你在自己的项目里能少踩几个坑。