limma多组差异表达分析：从设计矩阵到contrast矩阵的完整指南

作为一个常年泡在转录组数据里的生信人，我几乎每天都要跟差异表达分析打交道。过去几年，不管是用芯片数据还是RNA-seq数据，只要涉及多分组比较（比如对照组、处理组、时间点系列），我第一反应就是用limma包。虽然现在也经常用DESeq2或edgeR处理count矩阵，但如果数据是芯片表达矩阵、或者已经标准化好的log2表达量，limma始终是稳定性和速度方面的首选。尤其是当你的实验设计不是简单的“处理vs对照”两组比较，而是涉及到三组甚至更多分组时，limma凭借其线性模型框架和强大的contrast矩阵设计，处理起来非常顺手。

这篇文章我就把用limma做多组差异表达分析的完整思路、实操代码和踩坑记录全部梳理一遍。内容不局限于“运行一下limma就行”，而是会讲清楚设计矩阵怎么构建、contrast矩阵怎么写、多组比较的结果怎么解读、以及如何避免那些常见的“看起来跑通了但实际上是错的”问题。适合正在处理多分组转录组数据、想要用R做差异分析的研究生和科研人员参考。

1. 多组差异表达的整体设计与思路拆解

1.1 多组比较和两组比较的本质差异

很多人一开始学limma都是从两组比较入门的，比如“肿瘤组 vs 正常组”，代码逻辑很简单：样本分成两组，直接跑 lmFit 和 eBayes ，然后 topTable 拉出显著基因列表。这个流程很好理解，本质上就是为每个基因拟合一个简单的线性模型，然后看处理组相对对照组的log2 fold change是否显著。

但多组比较就不一样了。假设你有四个分组：正常对照（Control）、低剂量处理（Low）、高剂量处理（High）、恢复期（Recovery）。这时候面临的已经不是“一个系数”的问题，而是需要回答以下几类问题：

• 四个组之间是否存在任何显著的表达差异？
• 哪个处理组相对于对照组表达谱发生了显著变化？
• 高剂量和低剂量处理之间的差异是否显著？
• 恢复期样本是否回归到了接近对照的状态？

这些问题如果拆开来做两两t检验，每一次只取两组数据跑一遍limma，虽然也能得到结果，但会带来两个很麻烦的后果：一是多次比较导致假阳性累积；二是每次只用两组数据，相当于放弃了其他分组的样本信息，统计效力下降。

limma解决这个问题的思路是：把所有样本一次性纳入一个线性模型，用一个“整体模型”去描述所有分组，然后再通过contrast矩阵在这个模型基础上做任意组间比较。这样做的好处是，估计基因表达方差时用到的是所有样本的信息，而不是只限于当前比较的两个组。

1.2 为什么limma适合处理多组比较

limma的核心优势在于它的经验贝叶斯（empirical Bayes）方法：它会借用所有基因的表达波动信息，去调整单个基因的方差估计，避免那些表达量低、波动大的基因因为偶然性被判定为显著。这个思路在多组比较中尤其重要，因为分组越多，每个组内的样本量可能越少，方差估计越不稳定，经验贝叶斯正好能起到“借力”的作用。

另外，limma的线性模型框架非常灵活，支持任意的实验设计。它不仅能够处理单因素多水平设计（也就是我们常说的多组比较），还能处理双因素设计（比如“处理×时间”的交互）、协变量校正（比如批次效应、性别、年龄）以及配对的blocking设计。这意味着你在分析前不需要把数据拆成小块分别跑，而是可以构建一个完整的模型，一次性回答多个生物学问题。

注意：limma最初是为芯片数据设计的，但后来通过 voom 函数扩展到了RNA-seq count数据。如果你手头是RNA-seq的count矩阵，建议用 voom 把counts转换为log2-CPM并估计均值-方差关系，然后再走limma的流程。下面会给出具体操作。

1.3 方案选型的考量：contrast矩阵是核心

多组比较最简单的落地方式是利用 makeContrasts 函数构造contrast矩阵。它的作用是在拟合好的线性模型上定义我们关心的比较。

还是举四组数据的例子：Control、Low、High、Recovery。你可以定义这些比较：

• Low vs Control：低剂量处理效应
• High vs Control：高剂量处理效应
• Recovery vs Control：恢复情况
• High vs Low：剂量依赖性效应

这些contrast本质上就是设计矩阵各列系数的线性组合。比如设计矩阵中，Control是截距列，Low、High、Recovery分别是各组的指示变量，那么“High vs Control”对应的contrast向量就是 High - Control ，在 makeContrasts 里写成 HighvsControl = High - Control 。

这样做的好处是，所有比较共享同一个拟合模型，方差估计一致，结果之间可比性强。而且 decideTests 函数可以一次性给出所有比较的显著基因上下调情况，方便做后续的Venn图、热图等可视化。

2. 核心细节解析与实操要点

2.1 输入数据的格式和预处理

在用limma之前，第一步是把表达矩阵准备好。这个步骤看似简单，但坑最多。

对于芯片数据，表达矩阵通常已经经过了标准化（如RMA、MASS），一般不需要额外处理，直接读入R即可。但需要检查数据中是否有缺失值、是否有重复的探针或基因名、是否有明显的批次效应。

对于RNA-seq数据，原始数据最好是基因水平的count矩阵，行是基因，列是样本。分析前建议过滤掉在绝大多数样本中表达量都很低的基因，具体标准可以根据数据情况调整，比较常用的做法是保留至少在某个比例的样本中CPM大于1的基因。这一步很重要，因为低表达基因的counts波动很大，不仅会拖慢计算速度，还会在后续多重检验校正中引入大量无意义的检验，降低统计效力。

过滤后，用 voom 做转换：

library(limma)
library(edgeR)

# 假设count_matrix是基因×样本的count数据，group是分组因子
dge <- DGEList(counts = count_matrix)
dge <- filterByExpr(dge, group = group)
dge <- calcNormFactors(dge)

# 设计矩阵
desig