3. 批量校正方法概述(SVA, ComBat, limma 等)
好,咱们这一章直接聊工具。做批次校正这么多年,我用过不下十种方法,但真正能在实际项目中站稳脚跟的,其实就那么几个。我自己的习惯是,先别急着调参,先把每种方法的脾气摸清楚——它们适合什么场景、有什么硬伤、什么时候千万别用。
你想想看,批量校正赛道里,主流的流派大概有三支:SVA(Surrogate Variable Analysis)、ComBat(Combatting Batch Effects)、以及limma 的线性模型族。此外还有 RUV(Remove Unwanted Variation)、Harmony 等等,但最基础、最常用、最不容易翻车的,就是前三个。我们逐一说。
3.1 SVA:用代理变量抓隐形批次
SVA 是我个人非常喜欢的一类方法。它的核心思想很简单——既然批次效应往往和实验时间、操作人员、试剂批次这些隐藏因素挂钩,那我们就先不问批次标签,直接从表达矩阵里估算出“代理变量”(surrogate variables),再把它们当作协变量回归掉。
我记得在 2020 年做过一个多中心转录组整合项目,每个中心送样时都没记录具体批次信息。我当时想,完了,没法用 ComBat 了,因为 ComBat 需要显式的批次列。后来换成了 SVA,用 sva 包的 num.sv() 先估一下代理变量个数,再用 sva() 提取,最后加到线性模型里。效果出奇地好,批次聚类在 t-SNE 上看几乎消失了。
但注意,SVA 有一个大坑——它可能会把真正的生物学差异也当成批次去掉。我曾经遇到过在估计代理变量个数时,如果设置 alpha = 0.05 太宽松,SVA 会抽出 10 多个代理变量,结果把处理组之间的差异都“校正”掉了。所以我的经验是:先用已知批次标签做验证,再决定代理变量个数。
批次标签不完整或完全未知;样本量足够(每组至少 10 个以上);生物学差异较大时效果更好。
我曾经在一个只有 6 个样本的小实验中用了 SVA,结果代理变量把几乎所有的组间差异都拿走了。小样本慎用 SVA,优先考虑更稳健的方法。
3.2 ComBat:最常用的批量校正“标杆”
ComBat 算是批次校正界的“万金油”。它的底层假设是:不同批次间的表达分布可以通过位置(均值)和尺度(方差)参数来对齐,并且这些参数可以用经验贝叶斯(Empirical Bayes)方法来稳健估计。说白了,ComBat 先估算每个基因在每个批次的偏移和缩放,再把这些偏移缩放到同一个基准。
我用 ComBat 处理过至少上百个数据集,几乎每次都能稳定地去除批次效应——只要你有明确的批次列。最舒服的是 sva 包里的 ComBat() 函数,一行代码搞定参数估计和校正。
library(sva)
corrected_matrix <- ComBat(dat = expression_matrix,
batch = batch_vector,
mod = model.matrix(~ group),
par.prior = TRUE,
prior.plots = FALSE)
要注意 mod 参数。如果你不做生物学分组回归,ComBat 会默认把组间差异也当成批次。所以我在项目中一定会在 mod 里传入实验设计矩阵,告诉 ComBat:“这些生物学变量必须保留”。
par.prior = TRUE:默认参数先验,样本量大时推荐;样本量小(每批次<5)时建议设为FALSE用非参数先验。prior.plots = FALSE:关闭先验图,节省时间。- 先运行一次
ComBat后,再用 PCA 检查批次是否混合。
3.3 limma:用线性模型“阻挡”批次
严格来说,limma 不是专门的批次校正方法,而是一个差异表达分析工具。但因为我们可以在 limma 的线性模型里直接把批次作为协变量放进去,所以它也能起到“校正”效果。做法就是把批次列转为因子,加到 model.matrix 中,然后对每个基因做 lmFit + eBayes,最后得到的残差或调整后的表达值就可以看作“校正后”的数据。
我个人的看法是:如果你只是想跑差异分析而不是做数据整合,直接用 limma 加批次协变量就足够了,不要先校正再跑差异——两步法会增加检验自由度损失。但如果你需要输出校正后的表达矩阵供下游聚类、分类用,那 limma 这一步就不太方便,因为它不会直接给你一个“清理过的矩阵”。
一个典型的 limma 批次校正写法是这样的:
design <- model.matrix(~ batch + group)
fit <- lmFit(exprs, design)
# 你可以用 getBias() 或 residuals() 提取去除批次效应后的表达
这里有个容易犯的错误:有人会把 batch 和 group 交互项也放进去,导致过拟合。除非你有充分的生物学理由认为批次和分组有交互作用,否则不要加。
3.4 其他方法简介
除了上面三个,还有几个常用的,我快速提一下:
- RUV (Remove Unwanted Variation):和 SVA 类似,也是找隐藏变量,但需要你提供“负对照基因”(稳定表达的基因)或“负对照样本”。适合单细胞或空间转录组。
- Harmony:专为单细胞整合设计,用迭代聚类的方法对齐细胞簇,不适合批量转录组。
- limma::removeBatchEffect():一个简便函数,能直接输出校正后的矩阵,但会牺牲一些统计自由度。
3.5 核心方法对比
我把最常用的三种方法整理成了一张对照表,方便你快速选择:
| 方法 | 需要批次标签? | 是否需要生物学变量? | 输出校正矩阵? | 对样本量的要求 | 核心弱点 |
|---|---|---|---|---|---|
| SVA | 否 | 建议传入 | 否(需从模型中获取) | 组内≥10 | 可能过度校正 |
| ComBat | 是 | 强烈建议 | 是(直接返回矩阵) | 每批次≥3 | 批次=0时无法使用 |
| limma | 是 | 必须 | 否(需手动提取残差) | 无严格限制 | 不适合大批次整合 |
3.6 方法选择决策流程
下面的 SVG 流程图是我在项目中常用的决策路线。你可以直接对照自己的数据情况来选:
3.7 总结:没有银弹,但有经验
本章说到这儿,你应该已经感受到:没有一种方法能通吃所有场景。我自己的建议是:
- 如果你手上有明确的批次列并且只想快速得到可用的校正矩阵→ComBat。
- 如果你不知道批次来源,但样本量足够大,想挖掘隐藏的批次因子→SVA 或 RUV。
- 如果你目标就是做差异分析,不需要输出校正矩阵→直接用 limma 在模型里控制批次。
- 千万不要在明知小样本的情况下用 SVA,也不要过度依赖后调参——我觉得多做几种方法的对比 PCA 图,比你读十篇论文都有用。
我通常会在拿到数据后,先画一个未校正的 PCA 图,标上批次颜色。如果批次有明显分离,我会先试 ComBat;如果批次混杂但看不到明显聚类,我会试 SVA。两种结果跑完后,再结合生物学意义的解释性来选择保留哪一个。记住,最终目的不是把批次全部消灭,而是让批次不影响后续的生物学发现。