limma多组差异表达分析全流程：设计矩阵、对比矩阵与voom实战

这篇内容打算好好写一下。多组差异表达分析是组学数据里绕不开的活尤其当你手上有“对照组 多个处理组”这种设计时很多入门教程还停留在两组比较的思维里一上来就跑循环 t 检验。limma 这套线性模型加经验贝叶斯的框架才是处理多组问题最顺手的那把刀。这里把完整流程、设计矩阵的原理、对比矩阵的写法、RNA-seq 数据要过的 voom 这一关以及我实际踩过的坑都梳理一遍。1. 为什么多组比较绕不开 limma从一次两两 t 检验翻车说起先讲个真实场景。我早年间帮一个课题组处理一批微阵列数据实验设计是“野生型对照 三个突变体”每组三个生物学重复一共 12 张芯片。当时课题组的分析流程还是一个一个脚本拼起来的里面最核心的差异表达部分用的是 for 循环每两个组之间跑一次 t 检验p 值小于 0.05 就算显著。翻车出在哪呢三个突变体都要跟野生型对照比那就是三组 t 检验每组大概两万多个基因最后全部基因算下来 p 值小于 0.05 的有七千多个课题组拿着这个名单去做富集分析出来的通路五花八门完全没法解释。问题不是实验做得差而是多重检验这个坑完全没考虑进去——两万次独立检验里按 0.05 的阈值算光随机误差本身就能产生一千个左右的假阳性。三组 t 检验一叠加假阳性数量就更多了。1.1 两两 t 检验的三大硬伤第一多重检验问题。这是最致命的一点。每次 t 检验都是独立决策基因数量越大假阳性越失控。你可以用 Bonferroni 或者 BH 校正但如果整个分析流程建立在“循环 t 检验 手动校正”的基础上代码又长又容易漏。第二方差估计不稳定。每组只有三四个重复的时候单个基因自身方差的估计噪音非常大。有的基因本身表达波动小有的基因波动大t 检验对每基因独立处理相当于完全无视其他基因的信息波动大的基因很容易被误判为显著。第三无法处理多组之间的整体比较。t 检验天然只能处理两组。如果你的实验设计里有三个处理组你想回答“这三组的表达谱整体上有没有差异”t 检验就干不了这件事你必须用方差分析ANOVA那种思想先做一个全局 F 检验再去看哪些组之间有差异。1.2 limma 解决多组比较的核心逻辑先整体建模再局部对比limma 的全称是 Linear Models for Microarray Data但现在已经完全兼容 RNA-seq 的计数数据。它的做法概括起来就三步把每个基因的表达量当作因变量把所有实验组的信息编码成设计矩阵里的一列列指示变量为每个基因拟合一个线性模型。这一步不是只比较两个组而是把所有组的数据一次性放进模型里拟合出一个整体的表达模式。然后如果你想回答“处理组 A 对比对照组哪些基因变化”就用一个对比向量把相应的回归系数相减得到一个对比。如果你想回答“三个处理组整体上有没有差异”就做一个 F 检验同时检验多个系数是否全为零。最后经验贝叶斯方法会“借力”——它把所有基因的方差估计往一个全局的基准上收缩方差大的基因拉回来一点方差小的基因放松一点。这个收缩过程把单个基因方差不稳定的问题大幅削弱即使每组只有三个重复也能得到比较稳健的统计推断。一句话总结limma 不是一组一组去比而是把所有数据一次建模再通过对比矩阵切出你关心的那些比较。这才是多组比较的正确打开方式。2. 准备工作表达矩阵、分组因子和设计矩阵很多教材讲 limma 上来就直接给代码跳过了最关键的“为什么设计矩阵长这样”的部分。实际上多组分析里最容易出错的并不是 limma 本身而是设计矩阵里栏位跟你实验分组之间的对应关系。这里我把从数据到步入模型的完整思路捋一遍。2.1 表达矩阵和分组因子的组织方式假设你有四组样本对照ctrl、处理 AtreatA、处理 BtreatB、处理 CtreatC每组三个重复共 12 个样本。那你的表达矩阵expr_mat至少应该是这样行是基因或探针列是样本列名与样本一一对应表达量如果是微阵列数据一般是 log2 处理后、已经过归一化的值如果是 RNA-seq 计数数据后续要通过 voom 处理。分组因子group是 limma 建模的核心。它的长度必须等于样本数值必须是因子类型水平顺序就是你希望的设计矩阵列的顺序。这一点很多人不重视结果model.matrix出来之后列名完全不是自己想要的顺序后面看对比矩阵时一脸懵。group - factor( c(ctrl, ctrl, ctrl, treatA, treatA, treatA, treatB, treatB, treatB, treatC, treatC, treatC), levels c(ctrl, treatA, treatB, treatC) )levels的顺序要刻意设定别让 R 帮你按字母排。这个顺序会直接映射到设计矩阵的列也映射到你后面写对比矩阵时引用的列名。2.2 设计矩阵的两种编码方式为什么 -0 group 更常用设计矩阵用model.matrix生成典型写法有两种# 写法一带截距项treatment contrast design1 - model.matrix(~ group) # 写法二无截距项cell means model design2 - model.matrix(~ 0 group)写法一默认把第一个水平这里就是ctrl设成基线估计出来的系数是“处理组相对对照组的变化量”——第一列是截距代表对照组均值后面每一列是某个处理组均值减对照组均值。写法二把每个组的均值直接作为一列共 4 列。初次接触的人可能会问为什么 limma 教程里几乎都用写法二原因特别直白写法二让列名和组名一一对应起来后面写对比矩阵时你能用ctrl、treatA、treatB这样直观的名字去做加减法而不是记着“第 2 列是 treatA 对 ctrl 的变化第 3 列是 treatB 对 ctrl 的变化……”。尤其在组多的时候写法二极大降低坐标错位带来的风险。提示两种写法最终得到的差异表达结果在数学上完全等价区别只是系数含义和写对比矩阵的直觉性。我强烈建议多组场景一律用~ 0 group。2.3 列名就是对比矩阵里的变量名这一步必须确认生成设计矩阵后立刻做一件事检查列名。design - model.matrix(~ 0 group) colnames(design) # [1] groupctrl grouptreatA grouptreatB grouptreatC注意列名是groupctrl而不是ctrl带着group前缀。后面写对比矩阵时要么带上这个前缀要么直接把列名改干净colnames(design) - levels(group) # [1] ctrl treatA treatB treatC这一步非常值得做。对比矩阵里引用变量名时用干净的组名比用默认前缀要省心很多也方便看输出结果。3. 对比矩阵多组比较的核心“说明书”设计矩阵只是把实验结构表达了出来limma 还不知道你要做哪些比较。对比矩阵就是告诉模型“我想看哪几个对比每个对比里谁减谁”。3.1 用 makeContrasts 创建多组对比假设你的核心科学问题是“三个处理组分别相对对照组有哪些差异表达基因”那对比矩阵这样建contr_matrix - makeContrasts( treatA_vs_ctrl treatA - ctrl, treatB_vs_ctrl treatB - ctrl, treatC_vs_ctrl treatC - ctrl, levels colnames(design) )如果你还想知道“处理组之间是否有整体差异”这种全局问题limma 也支持直接用设计矩阵的系数做 F 检验这个后面专门讲。对比矩阵的本质是系数的线性组合。treatA_vs_ctrl treatA - ctrl相当于把设计矩阵里treatA这列的系数减去ctrl列的系数得到的值就是两组均值差。3.2 不只是两两比较多组场景下的 F 检验我见过不少人在多组比较里犯一个思维惯性错误认为差异表达分析就等于“两两比较”。实际上如果你的问题是“三个处理组到底有没有任何一组与对照组不同”可以先做一个整体检验这就是 F 检验。在 limma 里实现很简单。先不急着写 makeContrasts。直接用lmFit拟合完整模型然后利用eBayes对系数做检验再用topTable的coef参数指定多个系数同时检验fit - lmFit(expr_mat, design) fit - eBayes(fit) topTable(fit, coef c(treatA, treatB, treatC), number 10)这里coef传了三个系数进去limma 会自动做一个 F 检验检验“这三个系数是否同时为零”。如果某个基因的 F 检验显著说明它在至少一个处理组里相对对照组有差异。看到显著基因后再走对比矩阵做两两比较定位到底哪个处理组贡献了差异。先全局后局部这是多组比较里非常稳妥的分析策略。3.3 对比矩阵的系数与回归系数的关系再解释一个坑设计矩阵是 cell means 模型时每个系数是组均值。对比矩阵里treatA_vs_ctrl treatA - ctrl算出来的新系数就是两组均值的实际差。而contrasts.fit做的核心工作就是把这组线性变换应用到模型的拟合值和方差上得到一个新的拟合对象。有个容易踩的坑如果你用了带截距的设计矩阵~ group那么设计矩阵里的grouptreatA这一列已经表示“treatA 对 ctrl 的差值”了你写对比矩阵时如果还写treatA - ctrl就要用到截距列的名字也就是“对比矩阵里的系数名”跟“设计矩阵列名”必须对得上。如果用的是带截距项的设计矩阵列名变成(Intercept), grouptreatA, grouptreatB, grouptreatC对比矩阵就得按这些名字来。两种方式的坑点不一样我个人认为~ 0 group在设计上最容易自查。4. 核心分析流程逐步拆解从 lmFit 到 topTable前面背景讲完这部分给一套可以直接复用的完整代码。假设你已经准备好了expr_mat行是基因列是样本和group因子微阵列数据不需要 voom。4.1 第一步lmFit 拟合线性模型library(limma) design - model.matrix(~ 0 group) colnames(design) - levels(group) contr_matrix - makeContrasts( treatA_vs_ctrl treatA - ctrl, treatB_vs_ctrl treatB - ctrl, treatC_vs_ctrl treatC - ctrl, levels colnames(design) ) fit - lmFit(expr_mat, design)lmFit做的事情很直白对每个基因把表达量向量当作因变量12 个样本对应 12 行用最小二乘法拟合出设计矩阵对应的系数。注意这里拟合的是“线性模型”不是那些弯弯绕绕的机器学习模型。它假设每个基因在处理组间的关系可以通过组均值的组合来表达这个假设在绝大多数标准实验设计下是合理的。lmFit对每个基因是完全独立的因此这一步通常很快哪怕几万个基因也就几秒的事。4.2 第二步contrasts.fit 应用对比fit2 - contrasts.fit(fit, contr_matrix)这一步把线性模型的系数从“每个组的均值”转换到“组与组之间的差值”并且同步更新方差-协方差矩阵。这里要特别提醒contrasts.fit不会改变拟合的数学本质它只是做一个线性变换。如果你之前的设计矩阵列是组均值那这里生成的新拟合对象的系数就是对比矩阵定义的那些对比值。后续eBayes和topTable都在这个新对象上操作。4.3 第三步eBayes 经验贝叶斯方差收缩fit3 - eBayes(fit2)这是 limma 区别于普通线性建模的灵魂步骤。对于每个基因limma 先估计一个原始方差然后利用所有基因的方差分布情况把这个估计向全局均值方向收缩。直观理解就是某个基因的原始方差如果异常小可能是偶然就把它往上微调一点如果异常大就往下压一点。这样得到的方差估计比单基因独立估计稳健得多尤其在重复数少每组 3 个的时候效果极其明显。收缩程度由先验自由度控制limma 内部会自动估计。eBayes之后每个基因对应一组修正后的 t 统计量、p 值以及对整个实验的 F 统计量。4.4 第四步topTable 提取结果results - topTable( fit3, coef treatA_vs_ctrl, number Inf, sort.by none, adjust.method BH )coef指定你要看哪个对比number Inf表示返回全部基因不要只返回默认的前 10 个adjust.method BH是 Benjamini-Hochberg 校正也是差异表达分析里的默认主流选择。很多人在这一步漏掉一个细节topTable返回的adj.P.Val是针对你指定的coef这一列对比的校正 p 值。如果你想做多个对比比如三个处理组分别对比对照组通常需要把每个对比单独跑一次topTable或者用decideTests做一个联合决策。4.5 decideTests多组比较里被低估的函数当你关心三组对比时一个现实问题浮出水面如何综合三个对比的结果做基因筛选是三个对比各自 BH 校正后各筛各的还是三个合在一起筛limma 提供了一个漂亮的工具decideTeststest_results - decideTests(fit3, method global, adjust.method BH) summary(test_results)method global意味着 limma 会同时考虑所有对比结果用一个全局错误率控制标准来标记每个基因在每个对比中是否显著。输出是一个矩阵行是基因列是对比1 代表上调、-1 代表下调、0 代表不显著。用summary(test_results)可以看到每个对比里上调/下调基因的数量。这对多组比较特别有用因为你可以一眼看出哪个处理组的影响最大。5. RNA-seq 数据先用 voom 把计数转成适合线性模型的权重limma 早期是为微阵列设计的但 RNA-seq 时代它照样活跃在生信前线。关键在于一个函数voom。5.1 为什么计数数据不能直接进 lmFitRNA-seq 的本质输出是 read count 整数。这个数据有一个通性期望值越高方差也越大——低表达基因的计数之间的方差小高表达基因的计数之间的方差大而且这种均值-方差关系不是线性的。直接把 count 填进lmFit相当于假设所有基因的方差同质这严重违背实际。有人可能会想那我先把 count 转成 log2 CPMcounts per million不就行了转完之后数值范围倒是正常了但低表达基因的 log2 CPM 离散度跟高表达基因仍然不一样只是程度改善了并没有根本解决方差随均值变化的问题。5.2 voom 的工作原理与标准用法voom解决这个问题的思路很聪明它先对 count 做 log2 CPM 的归一化再拟合一个“均值-方差关系”的经验曲线然后根据每个基因、每个样本的表达水平给每个数据点算出一个权重inverse variance weight。权重信息随线性模型一起使用等于给高方差点降权给低方差低表达的点提权。用法非常直接library(edgeR) counts - read.delim(your_count_table.tsv, row.names 1, check.names FALSE) group - factor(...) design - model.matrix(~ 0 group) colnames(design) - levels(group) v - voom(counts, design, plot TRUE) fit - lmFit(v, design) fit2 - contrasts.fit(fit, contr_matrix) fit3 - eBayes(fit2) results - topTable(fit3, coef treatA_vs_ctrl, number Inf)注意这里counts的格式行是基因列是样本。进入voom之前请确保数据已经过滤掉低表达基因用filterByExpr或手动设定 CPM 阈值都能完成存在可选的归一化需求时用calcNormFactors取得 TMM 归一化因子再传给voom生物学重复尽量保持平衡设计voom 对不平衡样本量理论上能处理但统计功效会明显下降。plot TRUE这个参数值得用上。它会画一张均值-方差关系图。如果看到一个平滑下降的曲线说明 voom 的权重方案运作正常如果看到完全平坦的散点说明数据可能已经接近微阵列特性或者过滤基因太狠了。5.3 voom 的多组处理注意点多组场景下voom的设计矩阵是全局设计矩阵不要按组别分开展开。这意味着设计矩阵的行数是全部样本数列数是全部组数。voom会一次估计整个数据的均值-方差关系然后用这个关系为每个样本点赋权。如果你手工把每组切片单独跑 voom会破坏全局方差结构的估计后面统计推断会系统性地偏温和或偏保守。提示RNA-seq 多组分析里另一个常见流派是edgeR的 QLF 检验或DESeq2。对多组全局检验DESeq2 的 LRT似然比检验也提供了不错的方案。如果你更习惯 limma 的灵活性、对比矩阵的直观性和速度那么 voom limma 是一个成熟度非常高的选择两种方法得出的主流结论通常高度一致。6. 结果可视化与差异基因筛选策略统计推断做完下一步就是面对几千个“显著基因”怎么筛、怎么展示的问题。这部分我给出我在多组项目里常用的筛选策略和可视化组合。6.1 多组比较结果的__筛选__标准不能只用一个固定阈值多组比较时建议先整体后局部。先看全局 F 检验显著的基因数量。如果 F 检验显著基因很少说明这个实验的整体效应很小后面就算两两比较筛出几百个基因大概率也是边际效应。再看各对比的显著基因数量。用decideTests的 summary 快速掌握分布。通常我会记录每个对比里上调/下调基因数再统计“只在一个处理组差异”和“在所有处理组都差异”的基因数。对于最终的差异基因列表我的常用标准是abs(logFC) 1且adj.P.Val 0.05。如果你的实验比较的是强扰动处理比如敲除某个核心转录因子阈值放松一点也没关系如果处理效应微弱这可能过于严格。实操中往往要先跑一遍看看结果分布再定阈值而不是从第一秒就锁死参数。6.2 热图多组样本表达模式的最佳展示方式多组数据最适合用热图看全局模式。用显著基因的表达量画热图行按基因聚类列按样本分组排列。下面是最小可用的示例library(pheatmap) sig_genes - rownames(results_subset) expr_hm - expr_mat[sig_genes, ] annotation_col - data.frame(group group) rownames(annotation_col) - colnames(expr_hm) pheatmap( expr_hm, annotation_col annotation_col, scale row, show_rownames FALSE, clustering_method ward.D2 )scale row是必要的因为不同基因的表达量基数差异非常大无论按原始表达量还是 log2 CPM 直接画高表达基因都会把热图的颜色范围压扁。按行标准化后每个基因在样本间的相对变化才能体现出来热图才有阅读价值。画完之后仔细看列的分组聚集情况。如果对照组和处理组在热图上完全分不开那即使统计上跑出一堆显著基因也要对实验质量和批次效应打个问号。6.3 火山图和对比之间交集分析单个对比的火山图是最直观的差异表达展示手段plot(results$logFC, -log10(results$adj.P.Val), pch 20, col ifelse(results$adj.P.Val 0.05 abs(results$logFC) 1, red, grey80), xlab log2 fold change, ylab -log10 adjusted p-value)在多组项目里我还特别建议画一个 contrast 之间的 UpSet 图或维恩图。用一个布尔矩阵标记每个基因是否在某个对比中显著然后看交集和并集。这能回答一些很有价值的生物学问题某个处理组特异影响的基因有哪些所有处理组共同响应的核心基因集是什么用UpSetR包就能非常方便地完成这件事。library(UpSetR) sig_mat - data.frame( treatA rownames(res_A) %in% sig_A_genes, treatB rownames(res_B) %in% sig_B_genes, treatC rownames(res_C) %in% sig_C_genes ) upset(sig_mat, nsets 3)这类图在组会汇报和论文里都相当加分——评审人一眼就能看出你的分组策略和组间关系是经过逻辑设计的。7. 实操总结多组差异分析的完整流程与我的建议最后把这套完整流程串起来给你一个可以直接套用的总览。7.1 完整流程速览一次标准的多组差异表达分析我认为可以拆成这几个环节数据预处理完成质控、过滤低表达基因、TMM 归一化RNA-seq构建设计矩阵model.matrix(~ 0 group)列名设为组名构建对比矩阵makeContrasts定义你要回答的科学问题拟合模型lmFit-contrasts.fit-eBayes全局检验查看 F 检验显著基因数量局部比较逐个 contrast 跑topTable联合决策用decideTests或交集分析综合多个对比结果可视化热图、火山图、UpSet 图下游分析根据差异基因列表做富集分析。7.2 我在实操中踩过的三个坑第一个坑是设计矩阵的列名问题。早期我用model.matrix(~group)后直接写自定义对比结果列名是grouptreatA这种东西写对比时引用变量名报错而且报错提示还不太直观。现在我一律用~ 0 group再手动改列名整个流程清爽很多。第二个坑是不能直接对 count 数据跑lmFit。哪怕只是试跑、看流程通不通也别这么干。我见过有人用log2(count 1)糊弄过去虽然跑得通但结果的排序和显著性跟 voom 后的结果差别很大后续的生物学结论完全可能被带偏。第三个坑跟筛选阈值有关。我早期喜欢一上来就adj.P.Val 0.05 abs(logFC) 2一股脑筛结果在微弱效应的实验里筛出来的基因少得可怜富集分析根本跑不动。后来学乖了先看全局检验和火山图分布再根据数据形态决定阈值。阈值是工具不是目的你要保证最终拿到的基因列表能够回答你的生物学问题。7.3 一个额外提醒生物学重复的数量很影响经验贝叶斯的效果经验贝叶斯的核心优势是“信息借力”但组内重复数也不能少到离谱。每组两个重复虽然技术上能跑但方差估计太粗即使有收缩可靠性也有限。我的个人建议是至少每组三个重复四个更好。在预算有限时把重复从两组加到三组带来的统计功效提升往往比增加测序深度更明显这一点在 do 组学实验设计时常被忽略。多组差异表达分析本身并不复杂把设计矩阵、对比矩阵和 limma 的经验贝叶斯原理理解到位你对分析流程的掌控感会完全不一样。之后哪怕遇到更复杂的设计比如含协变量、含批次效应、时间序列实验也能基于同一套建模逻辑信手拈来。