limma 在批次校正中的应用：线性模型方法
说实话，我觉得 limma 是生信批次校正领域的“瑞士军刀”。它原本是为基因表达差异分析设计的，但它的线性模型框架实在太灵活了。我今天想跟你聊聊，怎么用这把刀去处理批次效应。
limma 的核心思想很简单：把批次效应看作线性模型里的一个协变量。你想想看，如果你的数据里包含了多个批次，理论上每个批次都会给表达量加一个（或乘一个）固定偏移。那我们只要估计出这个偏移，把它减去不就完事了？
limma 是怎么工作的？
limma 的全称是 Linear Models for Microarray Data。它最早是为芯片数据设计的，但现在广泛应用于 RNA-seq、蛋白质组学等各种组学数据。它的核心管线分三步：
- 构建一个全模型：把处理效应和批次效应都塞进设计矩阵里。
- 估算系数：拟合模型，求出批次效应的系数。
- 剥离批次效应：从原始表达量减去批次效应部分，保留残差和处理效应。
这里最妙的地方在于——它用的不是普通的线性模型。limma 通过经验贝叶斯方法（eBayes）来收缩估计的方差。什么意思呢？就是它会“借用”所有基因的信息来稳住每个基因的方差估计。当你样本量很少的时候，这个特性特别关键。
我早期做的一个项目里，前两个批次只有各 3 个生物学重复。如果用普通的线性模型去分离批次效应，方差简直飘得像风筝。但 limma 的 eBayes 一出手，方差估计一下就稳住了。这是 limma 的招牌本领。
怎么在 R 里实现？
代码非常简单。假设你有一个表达矩阵 expr（行为基因，列为样本），以及一个数据框 pdata 包含批次信息和处理分组：
# 步骤1：构建设计矩阵
design <- model.matrix(~ 处理分组 + 批次, data = pdata)
# 步骤2：拟合并用经验贝叶斯
fit <- lmFit(expr, design)
fit <- eBayes(fit)
# 步骤3：提取仅包含批次效应的系数矩阵
# 假设批次因子有2个水平，批次系数列名为"批次2"
batch_coef <- fit$coefficients[, "批次2", drop = FALSE]
# 如果有多个批次列，用 cbind 合并
# 步骤4：校正后的表达量 = 原始值 - 批次效应贡献
corrected_expr <- expr - (design[, "批次2", drop = FALSE] %*% t(batch_coef))
# 或者使用 limma 自带的 removeBatchEffect()
corrected_expr <- removeBatchEffect(expr, batch = pdata$批次)我的个人习惯：我喜欢显式构建设计矩阵去做校正，而不是直接调用 removeBatchEffect()。为什么呢？因为手动构建设计矩阵，我可以明确看到哪些系数被当作批次效应剥离了，避免误伤处理效应。
适用场景：什么时候用 limma？
| 场景 | 推荐程度 | 说明 |
|---|---|---|
| 有明确的处理分组与批次分组 | ★★★★★ | limma 能同时估计处理效应和批次效应，是标准做法 |
| 批次效应恒定，但与处理效应有重叠 | ★★★★☆ | 线性模型可以区分，只要设计矩阵非退化 |
| 样本量少（每批次3~5个样本） | ★★★★☆ | 经验贝叶斯对样本量小的批次特别友好 |
| 事后批量整合多个独立数据集 | ★★★☆☆ | 如果各批次处理标签不一致，需要先做标签对齐 |
| 批次效应非线性或样本间有依赖关系 | ★★☆☆☆ | 线性模型假设批次效应是加性常数，非线性时不够用 |
说白了，limma 最适合的场景是你的实验设计比较干净——每个批次里都包含所有处理组，样本来自同一个总体。它对于“批次效应是加性常数”这个假设非常依赖。如果批次效应是乘性的或者在少数样本中剧烈突变，那线性模型就慢慢失灵了。
核心局限：你必须知道的坑
局限一：对非线性批次效应无能为力
limma 假设批次效应是以加性常数存在的。如果你的数据来自不同测序平台（比如有的用 Illumina，有的用 Nanopore），这两种平台带来的偏差远不止一个常数那么简单。我曾经在整合 RNA-seq 与 microarray 数据时硬上 limma，结果校正后的数据出现了一大堆伪差异表达基因。回头看，这种非线性效应应该用 ComBat 或 Harmony 才能处理。
我曾经踩过的一个坑：在一个项目里，前两个批次跑在旧版试剂盒上，第三个批次换成了新版试剂盒。我直接套用了 limma 的线性模型去校正，结果发现新版试剂盒带来的表达量整体上移了约 20%，但某些基因却出现了下移。其实这是因为新试剂盒对不同基因有差异化的扩增偏好。线性模型无法捕捉这种基因特异性的非线性偏差。
局限二：不能处理非平衡设计
假设批次1 里有处理组 A、B、C，但批次2 里只有处理组 A 和 B，缺少 C。这种情况下，你的设计矩阵会变得严重共线性，批次效应和处理效应无法分离。limma 会报错，或者估计出极度不稳定的系数。
我的建议：碰到非平衡设计时，别硬用 limma。可以考虑用 sva 里的 fsva，或者干脆把缺少处理组的样本去掉，重新设计分析。当然，这种事后调整需要提前在论文的方法部分写清楚。
局限三：假设样本独立
如果同一个病人有不同时间点的样本，或者样本之间来自同一个家系，它们之间的相关性被 limma 忽略了。不过这个问题可以通过 duplicateCorrelation 函数来解决，limma 生态里是有配套工具的，只是很多人不知道。
# 处理相关样本（同一病人的重复测量）
corfit <- duplicateCorrelation(expr, design, block = pdata$病人ID)
fit <- lmFit(expr, design, block = pdata$病人ID, correlation = corfit$consensus.correlation)
fit <- eBayes(fit)与 ComBat 的对比
很多同学会问我：limma 和 ComBat 到底选哪个？我通常这么回答：
| 维度 | limma（线性模型） | ComBat（经验贝叶斯） |
|---|---|---|
| 模型假设 | 批次效应为加性常数 | 批次效应可加性+可乘性（调整均值和方差） |
| 对处理效应的处理 | 明确保留处理效应 | 可指定保留处理效应（通过mod参数） |
| 对样本量的要求 | 每批次3~5个样本即可 | 每批次建议5个以上，否则经验贝叶斯收缩不稳定 |
| 处理非线性偏差 | 差 | 中等（只能调整均值和方差，不能处理更复杂的非线性） |
| 速度 | 快，数秒可处理万基因 | 较慢，对大数据集可能需要几十秒 |
我个人工作流里的选择逻辑：
- 如果批次数量少（2~3个批次），且每个批次内处理分组完整 → 优先用 limma，速度快，模型透明。
- 如果批次数量多（5个以上），且批次间差异不仅在于均值，还在于方差 → 用 ComBat。
- 如果数据来自不同平台或不同实验类型 → 两个都不太好，考虑数据融合前的标准化预处理。
总结一下关键点
- limma 通过线性模型加经验贝叶斯，稳健分离批次效应和处理效应。
- 代码极其简单，可以用
removeBatchEffect()一行搞定，但我更推荐手动构建设计矩阵。 - 适用范围是有明确分组的平衡设计，批次效应近似加性常数。
- 局限很明显：不能处理非线性偏差、非平衡设计、样本间相关性（除非用
duplicateCorrelation）。
嗯，我建议你拿到数据后先做一个简单的 PCA 看看批次有没有明显的分离。如果 PCA 上不同批次的样本点只是整体平移（位置不同），那 limma 基本上够用了。如果它们分散在不同的方向或形状不同，那你可能就需要考虑 ComBat 或者其他非参数方法了。