第3章 经典方法入门——ComBat（sva包）

说到批次效应校正，ComBat 可以说是最经典、最常用的方法了。我最早接触它是在处理一批芯片数据时，几个批次间的表达值分布明显不同，用 ComBat 跑一遍后，数据就“齐整”多了。说实话，它背后的原理不算复杂，但用起来确实顺手。

原理介绍

ComBat 的全称是“ComBat for batch effects”，它假设每个基因的表达值服从正态分布，并且批次效应可以分解为两部分：整体偏移（additive）和离散度变化（multiplicative）。说白了，就是校正每个基因在每个批次里的均值和方差，让它们对齐到总体水平。

具体做法分三步：

  1. 对每个基因，用经验贝叶斯方法估计批次效应的位置参数（均值偏移）和尺度参数（方差变化）。
  2. 把这些参数从原始表达值中“减掉”并“等比缩放”。
  3. 得到校正后的表达矩阵，不同批次的数据就能放到一起分析了。

为什么要用经验贝叶斯？因为批次数量少，直接估计每个基因的批次效应会不靠谱。引入先验分布可以借用“全基因组信息”来稳定估计，这在我在一个只有3个批次的项目里尤其有用——效果比纯线性回归好不少。

ComBat批次效应校正流程 ComBat 批次效应校正流程 原始表达矩阵（多个批次） 逐基因：估计批次位置参数 γ 和尺度参数 δ 经验贝叶斯：整合先验，稳定参数估计 校正：Y_adj = (Y - γ) / δ 校正后表达矩阵

安装与加载

ComBat 在 sva 包里。安装很简单，一条命令搞定：

install.packages("sva")

如果安装遇到网络问题，我建议你试试用镜像源，比如中科大的：

install.packages("sva", repos = "https://mirrors.ustc.edu.cn/CRAN/")

加载的时候要注意，sva 包比较大，首次加载会慢一点。别急，喝口水的功夫就加载好了：

library(sva)
提示： 如果你用的是 R 4.0 以上版本，sva 包在 Bioconductor 里也能装。不过我个人更倾向于用 CRAN 版本，稳定且更新及时。

ComBat 函数使用

核心函数就是 ComBat()。它的最基本用法是传入表达矩阵和批次向量。我直接给你看代码：

# 假设 dat 是基因表达矩阵（基因行，样本列）
# batch 是与样本一一对应的批次标签向量
library(sva)

adjusted <- ComBat(dat = dat, batch = batch)

输出就是一个校正后的矩阵，行列数和原来一样。我曾在一个包含5个批次、约2万基因的数据集上跑过，大概30秒就出结果了，速度令人满意。

如果你有协变量（比如性别、年龄），想在校正时把这些因素考虑进来，可以加 mod 参数：

# 假设 covar 是协变量数据框或矩阵（样本行，协变量列）
adjusted_with_cov <- ComBat(dat = dat, batch = batch, mod = covar)
关键点： mod 里的协变量必须是“感兴趣的变量”——你想保留它们对表达值的影响，不把它们当作批次效应抹掉。比如疾病状态，你希望校正后仍能看到病例 vs 对照的差异，那就把疾病分组放进 mod 里。

参数详解

ComBat 函数里有几个参数，我整理成了一张表，方便你快速查阅：

参数名类型默认值作用
dat数值矩阵无基因表达矩阵，行是基因，列是样本
batch因子/数值/字符向量无每个样本所属批次标签，长度等于样本数
mod矩阵/数据框NULL需要保留的协变量矩阵（不含截距）
par.prior逻辑值TRUE是否使用参数先验（默认是，速度快；FALSE 则用非参数先验，更稳健但慢）
prior.plots逻辑值FALSE是否在屏幕上绘制先验分布诊断图
mean.only逻辑值FALSE是否只校正均值，不校正方差（设为 TRUE 时速度更快）
ref.batch数值/字符NULL指定一个参考批次，其他批次都向它靠拢（适合有明确基准批次的情况）

我挑几个重点参数说说我的实际经验：

  • par.prior：默认 TRUE，用的是参数法，速度很快。如果批次数量很少（比如只有2个批次），我会改成 FALSE 试试，结果有时更鲁棒，但计算时间会增加。我建议你先用默认，效果不好再换。
  • mean.only：如果批次间差异主要是均值偏移，而不是方差变化，设成 TRUE 能大幅提升速度。我处理过一个数量级差异很大的批次，用 mean.only 后校正效果完全可以接受。
  • ref.batch：以前我调试过几个 RNA-seq 项目，质检环节发现某个批次是“黄金标准”，其他批次都是后来补做的，这时候设 ref.batch 就很合理了。
注意： ref.batch 和 mod 不要同时乱用，很容易把协变量信息也拉到参考批次里，导致偏差。我曾经踩过这个坑，花了半天才排查出来——先检查 ref.batch 参数文档，再回头看的代码。

实战小贴士

最后分享几个我在项目中总结的实用技巧：

  1. 运行前检查批次标签长度：batch 向量长度必须等于 ncol(dat)，错一个程序就报错。我习惯先用 table(batch) 检查一下样本分布。
  2. 缺失值处理：ComBat 不接受缺失值（NA）。如果有，要么填充（如用行均值），要么删掉对应基因或样本。
  3. 保存中间结果：ComBat 返回的是完整矩阵，占用内存不小。如果后续要反复使用，建议用 saveRDS() 存起来，不用每次重跑。
  4. 可视化验证：校正前后画个 PCA 图或箱线图，看看批次分离是否减少。这一步我一直坚持做，毕竟“用了方法不代表就能掩盖问题”。

下次当你面对一堆批次混乱的表达矩阵，不妨先试试 ComBat。它虽然简单，但在绝大多数场景下表现都不错。

← 上一章 📖 返回目录 下一章 →