ComBat 方法概述
好,咱们进入正题。ComBat 这个名字你可能听过,全称是 Combined Batches,说白了就是一个专门用来消除批次效应的统计方法。我在处理多中心 RNA-seq 数据时,几乎每次都会用上它。你想想看,不同批次的数据放在一起,如果直接分析,结果很可能被批次差异带偏了。ComBat 就是来救场的。
ComBat 的发展历史
ComBat 最早是 2007 年由 Johnson 等人提出来的,当时主要用于基因表达微阵列数据。那会儿大家发现,同一个样本在不同批次里跑,表达量能差出好几倍。传统方法无非是线性回归调整均值,但效果有限。ComBat 的创新在于引入了经验贝叶斯框架,能利用全体基因的信息来稳定估计批次参数。
我记得第一次看到那篇论文时,心想:这不就是把贝叶斯拿来做数据清洗吗?后来自己在项目里试了一下,确实比简单减法靠谱多了。后来几年,ComBat 被扩展到了甲基化、蛋白质组学、代谢组学等领域。2018 年左右,又有了支持协变量版本的 ComBat 和能处理连续变量的 LongComBat。所以说,它从微阵列一路进化到了多组学批次校正的标配工具。
| 阶段 | 时间 | 主要进展 |
|---|---|---|
| 提出 | 2007 | Johnson 等提出 ComBat,针对微阵列基因表达 |
| 扩展 | 2010–2015 | 应用扩展到甲基化、RNA-seq 等 |
| 增强 | 2018+ | 加入协变量支持、LongComBat 处理连续批次 |
| 多组学 | 2020+ | 在蛋白质组、代谢组中广泛验证 |
ComBat 的基本原理
ComBat 的核心假设是:每个基因的批次效应由两部分组成——位置效应(均值偏移)和尺度效应(方差变化)。而且不同的基因会共享部分信息,这正好可以用经验贝叶斯来估计。
具体来说,它假设每个基因在某个批次中的残差服从一个正态分布,而这个正态分布的均值和方差又是从一个共同的先验分布里抽出来的。通过最大似然和贝叶斯更新,我们能得到每个批次每个基因的位置与尺度参数的稳健估计,然后反过来把数据校准到参考批次。
个人经验:我在做 TCGA 和 GTEx 的整合分析时,遇到过批次效应特别强的案例——同一个基因在 TCGA 里的表达量整体比 GTEx 高一倍。用 ComBat 校正后,聚类图瞬间就干净了。不过要注意,如果某个基因在所有批次里都没有信息,它可能会过度收缩到全局均值,这时检查一下 residual 的方差很重要。
ComBat 的输入与输出
- 输入:基因表达矩阵(基因×样本),批次标签向量,可选协变量矩阵(如年龄、性别)
- 输出:校正后的表达矩阵,批次参数估计值
- 关键参数:prior.plots(是否画先验密度)、mean.only(仅校正均值)、par.prior(是否使用参数先验)
下面我画了一张流程图,把整个流程串起来。
关键数学直觉:ComBat 相当于把每个基因的批次效应想象成从一个“共享经验”分布中抽取出来的。如果某个基因在批次中的差异特别大,但它与其他基因的规律冲突严重,贝叶斯就会把它向整体均值拉回一点。这种“借用信息”的做法,正是它比传统方法稳定的原因。
一个简短的使用示例
在 Python 里,我们通常用 pycombat 包。我当年第一次用的时候,装这个包还费了点劲,现在直接 pip install combat 就能搞定。下面是一段核心调用:
import combat
# dat: 基因表达矩阵,行为基因,列为样本
# batch: 数组,每个样本对应的批次编号
# mod: 可选,协变量矩阵
corrected = combat.pycombat(dat, batch, mod=mod)
我曾经踩过的坑:有一次我忘记对表达数据进行 log2 变换就直接跑 ComBat,结果校正后出现了大量负值。后来翻论文才知道,ComBat 默认假设数据是近似正态的——微阵列数据通常 log2 处理过,但 RNA-seq 的 TPM 或 Count 往往高度右偏。所以我建议你先做一步 log2 或者 vst 变换,再跑 ComBat。否则效果会打折扣。
ComBat 的局限性
- 它假设批次效应是加性的,如果批次和生物学变量有交互作用,校正可能不彻底
- 当每个批次的样本量特别少(比如小于 3)时,参数估计不稳定
- 对极端离群值比较敏感,建议先做一步质量过滤
总体来说,ComBat 是批次校正领域最常用的方法之一,尤其在基因组学领域。它的成功在于用贝叶斯思想解决了“信息不足”的问题——你想想看,如果每个基因的批次效应都要单独估计,很容易过拟合。但借用全体基因的信息,就能得到更稳健的结果。这就是 ComBat 的精髓。