📚 课程目录第4章批次效应校正概述

批次效应校正概述：常见策略分类

做组学数据分析的人，迟早都会碰上批次效应。我自己第一次遇到它，是在处理两个月的RNA-seq数据时。结果一聚类，样本按采样日期分成了两群，和生物学变量毫无关系。当时我对着PCA图愣了十分钟，那种感觉就像你明明做的是水稻，结果长出来的全是稗草。

批次效应说白了就是：明明同一种处理、同一个条件，就因为实验日期、操作员、试剂批次不一样，数据里就硬生生多出了一套“假信号”。你想想看，这对后续差异分析、聚类、分类的影响有多大。

好在业界已经总结出了几类成熟的校正策略。我个人习惯把它们分成三大流派：基于模型的、基于位置的、基于学习的。这三类思路不一样，适用场景也不同。下面我一个个来讲。

基于模型：最常用的“统计控制法”

这类方法的核心思路是：把批次效应当作一个额外的协变量，塞进统计模型里去估计，然后从数据里剥离掉。你看，这其实就像在多元回归里加一个“批次”的虚拟变量——它本身不参与生物学变异，但你必须把它控制住。

最出名的就是ComBat，来自约翰霍普金斯的一篇经典论文。它的思路很巧妙：先假设每个批次里的每个基因都有一个批次偏移量和方差缩放因子，然后用经验贝叶斯去估算这两个量，最后把数据调回“同一把尺子”。

我在项目中遇到过好多次跨平台的芯片数据合并，用的就是ComBat。效果确实不错，尤其是跨平台、跨实验室的场景。但有一点要注意：ComBat默认假设批次间只有平移和缩放差异。如果你的数据中某个批次连分布形态都和别的批次不一样，那ComBat可能就力不从心了。

经验贝叶斯ComBat的核心思想

假设第i个基因在第j个样本中的表达值为Y_ij，拆成三部分：

真实生物学信号（Xβ）
批次偏移（γ）
批次方差异动（δ）

ComBat利用所有基因的数据来“借力”估计γ和δ，这就是“经验贝叶斯”的含义。

另一个常见的是limma包的removeBatchEffect函数。它的做法更简单：直接线性回归，把批次效应当作一个fixed effect剔除。不过它不直接修改原始方差结构，所以如果你后续还要做差异分析，我个人建议把批次变量直接放进limma的设计矩阵里，而不是先用removeBatchEffect再跑差异。

举个例子，假如你有两个批次的数据，用R语言调ComBat是这样的：

library(sva)
# 表达矩阵：行是基因，列是样本，已做过log2变换
expr <- as.matrix(read.table("expression_matrix.txt"))
# 批次编号向量（0和1代表两个批次）
batch <- c(rep(0, 20), rep(1, 20))
# 可选：生物学协变量（比如处理组/对照组）
mod <- model.matrix(~ as.factor(group))
# 执行ComBat校正
expr_corrected <- ComBat(dat = expr, batch = batch, mod = mod)

我的小建议

跑ComBat之前最好确认一下数据已做过合适的归一化（比如log2转换或quantile normalization）。ComBat对输入数据的尺度很敏感，原始counts直接丢进去效果会大打折扣。

基于位置：让数据“站在同一起跑线”

这类方法更直接——我不去建模批次效应的分布，而是直接操作数据的“位置参数”，比如均值、中位数、分位数。

最简单的例子就是z-score标准化：对每个批次内的每个基因，减去该批次的均值，除以该批次的标准差。这样所有批次的数据就“居中”到了均值0、方差1的尺度上。但这也带来了一个问题：你失去了原来的生物学量纲，并且不同基因之间的可比性被强制做了变换。

还有一种更强力的叫分位数归一化（Quantile Normalization）。它强制让所有样本的分布完全一致。嗯，这里要注意：分位数归一化原本是用于芯片数据探针水平的归一化，但它有一个很强的假设——所有样本的全局分布应该是一样的。如果你的条件处理真的改变了部分基因的高表达分布，分位数归一化会把这些差异也给抹掉。

我曾经在整合两个不同平台的蛋白表达数据时试过分位数归一化，结果发现一个平台里某个通路的信号被压平了。后来改用ComBat才把正确的信号还原出来。所以，基于位置的方法虽然快，但适用范围有限。

避坑提醒

不要对差异分析前的数据直接做z-score标准化。因为这种方法会无差别地抹平批次效应和生物学差异，导致你后续找不到真正的差异基因。

基于学习：高维数据的新解法

近五年，随着单细胞测序数据爆发，大家发现传统方法在高稀疏、高维度的数据面前有些吃力。于是基于学习的策略开始流行起来。

这里的学习，更多指的是基于子空间映射（比如Harmony）或基于深度生成模型（比如scVI）的做法。它们不假设批次效应是简单的平移缩放，而是把数据映射到一个低维共享的“潜空间”里，在这个空间里强制让不同批次的样本混在一起。

Harmony是我自己非常喜欢的工具。它的做法很优雅：先对数据做PCA降维，然后在软聚类的基础上不断迭代调整每个细胞的嵌入坐标，直到不同批次的样本在潜空间中“你中有我，我中有你”。我处理过一个9个批次的单细胞数据整合，跑了10轮Harmony，分群结果和生物学注释基本吻合，效果比ComBat好得多。

import scanpy as sc
import harmonypy as hp

# 加载单细胞数据
adata = sc.read_h5ad("scRNA_data.h5ad")
# 前期标准化和PCA
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata)
sc.pp.pca(adata)

# Harmony校正
ho = hp.run_harmony(
    adata.obsm['X_pca'],
    adata.obs,
    vars_use=['batch']
)
# 将校正后的嵌入放回adata
adata.obsm['X_pca_harmony'] = ho.Z_corr.T

除了Harmony，基于深度学习的做法（如scVI、scAlign）走的更远。它们用变分自编码器（VAE）来拟合每个批次的条件分布，然后从一个“去除了批次标签”的条件分布里重新采样表达值。效果通常很好，但缺点是训练时间长、可解释性差，而且超参数不少。

三类策略怎么选？我自己的经验法则

场景	推荐策略	原因
少量样本（<100），2~5个批次	基于模型（ComBat）	统计稳定，假阳性可控
跨平台芯片数据合并	基于模型（ComBat/BEAGLE）	专为芯片设计，效果好
单细胞大通量数据	基于学习（Harmony/scVI）	能处理稀疏性，高维度稳定
快速查看批次问题	基于位置（z-score/分位数）	速度快，适合诊断性绘图

你想想看，这三种策略本质上是一个“假设强度”由弱到强的谱。基于位置的方法几乎不做假设，但手段粗糙。基于模型的方法做了较强的统计假设，但在满足条件时精度高。基于学习的方法最灵活，但代价是计算和可解释性。

最后补充一句：没有银弹。我个人的习惯是，对一份新数据，先跑一个快速的基于位置的方法做可视化诊断，确认批次效应确实存在，再根据数据量、平台、下游任务选具体工具。千万别一上来就想“这个场景必须用某一种方法”——做数据分析，灵活比固执重要得多。