企业官网建设流程全解析

作为一个常年泡在转录组数据里的生信人，我几乎每天都要跟差异表达分析打交道。过去几年，不管是用芯片数据还是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) # 设计矩阵 design <- model.matrix(~ 0 + group) colnames(design) <- levels(group) # voom转换 v <- voom(dge, design, plot = TRUE)

voom输出的v$E就是log2-CPM表达矩阵，可以直接交给后面的lmFit。

注意：这里用了~ 0 + group而不是~ group。两者的区别在于，~ 0 + group会生成每个分组一列的design矩阵（没有截距列），这种形式在做多组比较时更直观，不容易搞混contrast的定义。如果R²新手习惯用带截距的形式，则在构建contrast矩阵时容易出错，建议统一用“无截距”的设计。

2.2 构建合理的实验设计矩阵

设计矩阵是limma分析的骨架。构建时我习惯用model.matrix函数，以因子的形式传入分组信息。

group <- factor(c("Control", "Control", "Control", "Low", "Low", "Low", "High", "High", "High", "Recovery", "Recovery", "Recovery")) design <- model.matrix(~ 0 + group) colnames(design) <- levels(group)

这时的design矩阵长这样（简略示意）：

每一行是一个样本，每一列是一个分组。矩阵的含义非常明确：某个样本属于哪个组，就在对应列取1，其余列为0。

如果实验设计中还有批次信息，可以将其加入模型，作为协变量校正：

batch <- c("B1", "B1", "B2", "B2", ...) design <- model.matrix(~ 0 + group + batch)

但如果样本量不大，加入太多协变量需要谨慎，因为每一个协变量都要消耗自由度，而经验贝叶斯虽然能缓解方差估计不稳定的问题，但自由度太少仍然会严重影响结果可靠性。

2.3 用makeContrasts构造多组比较的组合

设计矩阵准备好之后，下一步就是定义我们关心的比较。这里强烈推荐用makeContrasts，因为它是专门干这个的，而且代码可读性高。

contr.matrix <- makeContrasts( LowvsControl = Low - Control, HighvsControl = High - Control, RecoveryvsControl = Recovery - Control, HighvsLow = High - Low, levels = colnames(design) )

这里levels = colnames(design)表示contrast矩阵的列名基于design矩阵的列来定义。makeContrasts内部会解析Low - Control这样的字符串，生成对应的数值向量。

需要注意的一点是，contrast矩阵中的每个contrast，本质上是对design矩阵中某一组系数的线性组合。由于design矩阵是“每列代表一组”的形式，所以Low - Control就是“Low组系数减去Control组系数”，这正好对应两组间的log2 fold change。

后面拟合模型和做检验的代码就是标准的limma流程：

fit <- lmFit(v, design) fit <- contrasts.fit(fit, contr.matrix) fit <- eBayes(fit) results <- decideTests(fit) summary(results)

decideTests输出的是一个矩阵，行是基因，列是各个contrast，值表示该基因在对应比较中是否显著以及上下调方向（1为上调，-1为下调，0为不显著）。默认的判定标准是调整后P值小于0.05且logFC绝对值大于1，可以通过adjust.method和p.value参数调整。

2.4 整体F检验与两两比较的区别

多组比较中，很多人容易忽略一个问题：如果只做两两比较，那么“四组之间是否存在整体差异”这个问题并没有被直接回答。比如四个组之间差异都不显著，但在某些contrast组合下，可能单个两两比较都达不到显著阈值，而整体检验却可能显著。

limma提供了两种方式来看整体差异。一种是在没有做contrasts.fit之前直接对fit对象做eBayes，然后看每个基因的F统计量。这里的F检验对应的是“这个基因在所有分组之间的表达均值是否存在显著差异”。另一种是使用decideTests的global模式，将多个contrast的检验结果综合考虑，判断该基因是否在任意一个contrast中显著。

fit_all <- eBayes(fit) topTable(fit_all, number = 10, sort.by = "F")

这个F检验在多组比较中非常有用，它可以帮助我们快速筛选出那些“在任何一个组间比较中可能有差异”的基因，作为后续两两比较的重点关注对象。很多时候，我会先用F检验做一个初筛，再针对有整体差异的基因集去细看每个contrast的方向和大小。

但也要注意：F检验显著不等于每个两两比较都显著，它只说明至少有一组和其他组不同。具体是哪两组不同，还是要看contrast的结果。

3. 实操过程与核心环节实现

3.1 从表达矩阵到差异基因的完整R代码

下面我给出一个可以直接改着用的完整脚本，假设你已经有了一个count_matrix，列是样本，行是基因，并且有一个sample_info数据框包含样本的分组信息。

# ===================== 加载包 ===================== library(limma) library(edgeR) # ===================== 读入数据 ===================== # count_matrix: 行为基因，列为样本 # sample_info: 必须包含样本ID和分组信息两列 count_matrix <- read.csv("count_matrix.csv", row.names = 1, check.names = FALSE) sample_info <- read.csv("sample_info.csv") # 确认列名和样本信息顺序一致 stopifnot(all(colnames(count_matrix) == sample_info$sample_id)) # ===================== 构建分组因子 ===================== group <- factor(sample_info$group, levels = c("Control", "Low", "High", "Recovery")) # ===================== 过滤低表达基因 ===================== dge <- DGEList(counts = count_matrix) keep <- filterByExpr(dge, group = group) dge <- dge[keep, , keep.lib.sizes = FALSE] dge <- calcNormFactors(dge) # ===================== 设计矩阵 ===================== design <- model.matrix(~ 0 + group) colnames(design) <- levels(group) # ===================== voom转换 ===================== v <- voom(dge, design, plot = TRUE) # ===================== 构建contrast矩阵 ===================== contr.matrix <- makeContrasts( LowvsControl = Low - Control, HighvsControl = High - Control, RecoveryvsControl = Recovery - Control, HighvsLow = High - Low, levels = colnames(design) ) # ===================== 线性拟合与检验 ===================== fit <- lmFit(v, design) fit <- contrasts.fit(fit, contr.matrix) fit <- eBayes(fit) # ===================== 结果输出 ===================== results <- decideTests(fit) summary(results) # 输出每个contrast的top基因表 for (contrast_name in colnames(contr.matrix)) { top <- topTable(fit, coef = contrast_name, number = Inf, sort.by = "P") write.csv(top, paste0("topTable_", contrast_name, ".csv")) }

这段代码跑完之后，你会得到每个contrast对应的差异基因表CSV，里面包含了logFC、AveExpr、t统计量、P值、调整后P值（adj.P.Val）和B统计量等列。

3.2 结果表的解读：logFC、P值和B统计量的使用

差异基因表里最常看的几列是logFC、AveExpr、t、P.Value和adj.P.Val。其中：

• logFC：log2倍变化。正值表示该基因在比较的第一组（如Low）相对于第二组（如Control）表达上调，负值表示下调。实际生物学意义中，logFC绝对值大于1通常意味着表达量变化超过2倍。
• AveExpr：该基因在所有样本中的平均表达量。这一列可以用来判断差异是否出现在低表达基因中，低表达基因的差异往往可靠性较差。
• t： moderated t统计量。limma用经验贝叶斯调整了基因特异的方差，因此这里的t统计量比普通t检验更稳健。
• P.Value和adj.P.Val：P值和多重检验校正后的P值。差异基因筛选时应当使用adj.P.Val，而不是P.Value。

筛选显著差异基因时，我常用的标准是：

sig_genes <- topTable(fit, coef = "HighvsControl", number = Inf) sig_genes <- sig_genes[abs(sig_genes$logFC) > 1 & sig_genes$adj.P.Val < 0.05, ]

如果想一次性从多个contrast中提取共同的显著基因，可以利用decideTests返回的矩阵，配合VennDiagram包画韦恩图，看看不同比较间显著基因的重叠情况。

3.3 使用Venn图和多维标度图辅助结果展示

多组比较的分析不能只停留在输出表格上，可视化是理解和解释结果的关键。

多维标度图（MDS plot）：相当于PCA的另一种表现形式，用来检查样本之间的总体相似度。通常在跑正式差异分析之前我就会先看一下MDS图，确认样本是否按照分组自然聚类。

plotMDS(v, col = as.numeric(group), labels = group)

如果MDS显示样本没有按照分组聚类，那么后面的差异分析结果可能并不可靠，需要回到数据质量本身去排查。

Venn图：在多组比较中，Venn图适合看不同contrast之间显著基因的重叠和差异。

library(VennDiagram) venn.diagram( x = list( Low = rownames(results)[results[, "LowvsControl"] != 0], High = rownames(results)[results[, "HighvsControl"] != 0], Recovery = rownames(results)[results[, "RecoveryvsControl"] != 0] ), filename = "venn.png" )

当然，如果contrast数量超过三个，Venn图就不太直观了，这时候可以考虑用UpSetR包绘制UpSet图。

3.4 多组比较后的基因表达趋势可视化

多组比较相比两组比较，一个明显的优势是能看表达趋势。比如某个基因在Control、Low、High三组中呈现剂量依赖性的上升或下降，这是两组比较看不到的信息。

对于关注基因，我习惯用plotProfile函数或ggplot直接画每个组的表达均值和标准差：

library(ggplot2) # 提取某个基因在所有样本中的表达值 gene_name <- "ENSG00000123456" expr <- v$E[gene_name, ] plot_df <- data.frame( expression = expr, group = group ) ggplot(plot_df, aes(x = group, y = expression, fill = group)) + geom_boxplot() + geom_jitter(width = 0.2) + theme_minimal() + labs(title = gene_name, y = "log2 CPM")

这种图在文章里非常常见，而且能直观展示多组比较的生物学意义，尤其是剂量梯度实验和时间序列实验。

4. 常见问题与排查技巧实录

4.1 设计矩阵出现奇异（singular）问题

这是多组比较新手最容易遇到的问题之一。报错信息通常长这样：

Coefficients not estimable: groupLow

出现这个问题的原因，最可能是样本量和分组信息不匹配。比如某个组只有1个样本，或者design矩阵的列之间存在完全线性相关的关系。如果使用了带截距的design（~ group），列数会比无截距形式少一列，某些组别的系数会被当作基线，导致后续makeContrasts中的写法对不上。

解决方法是：优先使用无截距的设计矩阵~ 0 + group；检查每组样本量是否合理（至少3个生物学重复）；如果只有两个组，就没必要用多组比较的框架，直接两组比较更简单。

4.2 voom和limma在RNA-seq中如何配合使用

很多人会问：limma不是芯片数据的包吗，怎么用在RNA-seq上？实际上limma的voom函数就是为了让limma能够处理RNA-seq的count数据而设计的。voom的核心步骤是：先将count矩阵转换为log2-CPM，然后拟合均值-方差关系，最后为每个观测值计算一个精度权重。这个权重反映了该观测值在均值-方差关系中的可靠性，表达量越低，权重越小。

在运行voom之前，务必先做calcNormFactors。这一步是用TMM方法校正样本间的文库组成差异，不是简单的文库大小缩放，对后续差异分析的准确性有很大影响。

有一个细节需要注意：使用voom时，plot = TRUE会生成一张均值-方差关系图。如果图上趋势线显示方差随均值变化剧烈，说明数据适合用voom；如果趋势平缓，也可以考虑直接用log-CPM加lmFit的流程，但大多数情况下voom更稳妥。

4.3 多重检验校正后的显著基因数量过少

跑完decideTests后，summary结果显示显著基因数目为0或者极少，这种情况经常出现，尤其是样本量少、组内变异大的时候。

处理思路有以下几个：

• 检查数据质量：先看MDS图，如果同组样本没有聚在一起，差异信号会被组内噪声掩盖。
• 检查对照组的选择：多组比较中，基准组（baseline）的选择会直接影响contrast的生物学解释，但不应该影响显著基因的数量。如果数量差异极大，需要检查contrast定义是否正确。
• 适当调整筛选标准：有时adj.P.Val < 0.05过于严格，可以适当放宽到adj.P.Val < 0.1，或者只保留P.Value < 0.01加logFC筛选的组合。但这样做带来的假阳性风险需要自己在文中说明。
• 检查是否使用了错误的P值列：一定是用adj.P.Val，不要在差异基因筛选中只盯着P.Value。

4.4 多组比较中如何选择合适的参考组

在makeContrasts中，每一个比较都需要指定一组为参考。参考组的选择通常由生物学问题决定，比如临床研究中常用正常组织或安慰剂组作为参考，剂量实验中常用0剂量组作为参考。

这里有一个人为容易踩的坑：contrast矩阵写反方向。例如LowvsControl = Control - Low，这样写并不会报错，但结果的logFC方向就是反的。实际使用中，看到显著基因的logFC方向和预期不一致时，第一反应不是怀疑生物学，而是回去检查contrast的定义。

建议在每次写contrast时都加一行注释，写明“第一组相对于第二组是上调还是下调”，避免后续分析的时候把自己搞混。

4.5 结果重复性问题：批量效应和时间因素的校正

多组比较中，如果样本不是在同一个批次完成测序或芯片杂交，就极有可能引入批次效应。批次效应有时候很隐蔽，甚至会让无关的基因呈现显著的组间差异。

处理方式是在design矩阵中加入批次列，让线性模型把批次效应纳入考虑：

design <- model.matrix(~ 0 + group + batch)

但需要注意的是，加入太多协变量会让模型变复杂，样本量不大的情况下容易过度拟合。另外，加入协变量后，contrast矩阵的构建方式不变，因为在makeContrasts中我们只针对group相关的列做线性组合。

4.6 如何将limma结果与其他差异分析工具进行交叉验证

在实际项目中，我经常被问到“limma和DESeq2的结果怎么不一样”。这个问题非常正常，因为不同工具使用的统计模型不同：

• limma+voom：先转换为log2-CPM，再拟合加权线性模型，适合组间方差相似的情况。
• DESeq2：直接对counts建模，使用负二项分布，能更合理地处理低表达基因的离散度。
• edgeR：也是负二项分布模型，与DESeq2思路类似，但具体离散度估计方法不同。

在项目实操中，如果时间允许，我会用两种方法分别跑一下，然后取交集中的显著基因作为候选列表。这样做的结果更稳健，审稿人也更认可。如果交集很小，那就要回头检查数据质量、分组定义和参数选择，而不是急着选一个“看起来结果更好”的工具。

5. 实操经验总结：我从多组分析中积累的几个习惯

这些体验是我自己跑过的项目里，逐步积累出来的，不写进论文，但很影响结果的可靠性。

第一，永远先画MDS图。不管数据看起来多干净，先基于表达矩阵画一个MDS图，看看样本如何聚类。这一步能发现样本标签是否放错、是否存在离群样本、是否有明显的批次效应。如果样本没有按照预期分组聚类，我不会继续往下分析，而是先解决数据质量问题。

第二，多组比较的contrast矩阵不要写得太多。虽然makeContrasts支持定义任意多个contrast，但每多一个contrast，后续就要多输出一张表、多解读一批结果。合理做法是先根据生物学问题确定2到4个核心比较，比如“每个处理组vs对照组”，如果确实需要探讨剂量效应或交互作用，再额外添加。

第三，使用topTable时留意number = Inf。默认情况下topTable只返回前10行，这在快速查看结果时够用，但如果你要导出完整的差异基因表用于后续分析，一定要用number = Inf。这个是每次写脚本都要检查的点。

第四，decideTests的默认method是“separate”，即对每个contrast分别进行多重检验校正。如果希望从“多个contrast整体”的角度控制错误发现率，可以设置method = "global"。两种方法的结果会有差异，具体选择取决于你更关心单个比较的准确性还是所有比较整体的准确性。

第五，多组比较的显著基因筛选标准要提前定好。不要根据结果的好坏事后调整阈值，这种做法在统计上是不可接受的。我一般在分析前就会在脚本里写好adj.P.Val < 0.05 & abs(logFC) > 1，后面不做改动。

再多说一点，关于文章里展示limma结果时，除了差异基因表，一定要放一张火山图或MA图，这是审稿人很习惯看到的图表。用ggplot2自己画并不难，把logFC和adj.P.Val映射到x轴和y轴，再按阈值上色就行。

library(ggplot2) library(ggrepel) top <- topTable(fit, coef = "HighvsControl", number = Inf) top$sig <- "Not Significant" top$sig[abs(top$logFC) > 1 & top$adj.P.Val < 0.05] <- "Significant" ggplot(top, aes(x = logFC, y = -log10(adj.P.Val), color = sig)) + geom_point(size = 0.8) + scale_color_manual(values = c("grey60", "red")) + theme_minimal() + labs(x = "log2 Fold Change", y = "-log10 Adjusted P-value")

用limma做多组差异表达分析，掌握代码只是第一步，更关键的是理解线性模型的设计思路和contrast的含义。只要把design矩阵和contrast矩阵搞清楚了，多组比较的框架就可以灵活扩展到各种复杂的实验设计。希望这篇内容能帮你少走一些弯路，尤其是那些“代码没报错但结果逻辑不对”的坑，往往最难排查，也最影响结论。