R语言批次效应建模与自动化评估 - 第4节

4.1 批次效应移除的线性模型思想

先聊一个基础问题：批次效应到底怎么“移掉”？你直接减去平均数？那肯定不行。我早年刚接触这个功能时，还真试过按批次算均值再扣掉——结果后续分析全乱了。

实际上，limma 的 removeBatchEffect 背后是个标准的线性模型。简单说就是：

Y = Xβ + 批次效应 + 误差

其中 Y 是基因表达量矩阵，X 是设计矩阵（包含你关心的生物学变量），批次效应 是批次协变量对应的系数的乘积。函数做的事是：先拟合这个模型，估计出批次项的系数，然后从原始表达量中减去批次项贡献的值，只留下 Xβ + 误差。

我习惯把这种操作理解为“在控制生物学差异的前提下，把批次带来的偏移拉平”。说白了，它不是简单地“去噪”，而是用线性回归把可解释的批次分量揪出来并减去。

从设计矩阵到批次因子

那设计矩阵长什么样？假设你有 6 个样本，分两个批次，每组都有处理 vs 对照：

样本  批次  条件
S1     1    处理
S2     1    对照
S3     1    处理
S4     2    对照
S5     2    处理
S6     2    对照

对应的设计矩阵（不含截距）可能是：

        条件处理  批次2
[1,]        1      0
[2,]        0      0
[3,]        1      0
[4,]        0      1
[5,]        1      1
[6,]        0      1

拟合模型后，批次2的系数就是批次2相对批次1的表达偏移。这个偏移被估计出来后，removeBatchEffect 就会把它从所有样本中减掉。

你想想看，如果批次和条件完全混在一起（比如所有处理都在批次1，对照都在批次2），那就没办法区分了——这叫完全混杂。我遇到过这种“完美”设计，数据拿到手直接傻眼。

4.2 removeBatchEffect() 函数详解

函数在 limma 包里，签名长这样：

removeBatchEffect(x, batch, batch2 = NULL, covariates = NULL, design = NULL, ...)

我列个表格，把几个关键参数说清楚：

参数	用途	注意事项
x	表达量矩阵（行基因，列样本）	通常需 log2 转换后的数据，如 logCPM
batch	批次分组因子（向量）	必须与样本列一一对应
batch2	第二个批次因子（可选）	比如日期与操作员两个混杂因素
covariates	其他需要保留的连续协变量	比如年龄、GC含量等，不会被移除
design	完整的设计矩阵（包含生物学变量）	如果提供了，函数会用这个设计矩阵拟合，而不是自动构造

基本用法示例

# 假设我们已经有了 logCPM 表达矩阵 expr，批次信息 batch
# 以及一个分组因子 group（处理 vs 对照）
design <- model.matrix(~ group)
expr_corrected <- removeBatchEffect(expr, batch = batch, design = design)

这里有个坑：design 参数一定要包含你关心的生物学变量。我见过有人只传 batch 不传 design，结果函数内部会自动假设一个 intercept-only 模型——这样批次效应估计时就没有考虑生物学差异，相当于把生物学信号也擦掉了一部分。

与差异分析的关系

这个问题值得单拎出来说。我经常在论坛上看到有人问：“我用了 removeBatchEffect 做了批次矫正，然后跑 limma 差异分析，为什么结果不对劲？”

原因很简单：差异分析中的批次效应应该作为协变量留在线性模型里，而不是先一步去掉。

你把批次预移除后，原本属于批次项的自由度被“吞掉”了，残差自由度变大，检验会偏宽松——说白了就是假阳性。反过来，如果批次效应不明显，你又提前移除了，也可能把微弱的真实信号一起丢掉。

正确的做法是：

• 用 removeBatchEffect 的结果做探索性分析（PCA、聚类）
• 差异分析时，直接在 lmFit 的设计矩阵中包含批次变量

# 差异分析的正确姿势（无需预矫正）
design_full <- model.matrix(~ batch + group)
fit <- lmFit(expr, design_full)
fit <- eBayes(fit)
results <- topTable(fit, coef = "group")

我刚开始做分析时，为了可视化美观，先把批次去掉，然后顺手就把矫正后的矩阵拿去做差异分析。结果跟朋友的标准流程对比，假阳性基因多了整整一倍。从那以后，我总会在差异分析脚本开头加一句注释：# 此处不要用 pre-corrected data。

4.3 实战案例：真实数据演示

用手边的小例子走一遍。我们用 GSE5859 数据的一个子集（13 男 7 女，两个批次）。

library(limma)
library(GEOquery)

# 加载数据（此处省略完整下载，仅示意流程）
gset <- getGEO("GSE5859", GSEMatrix=TRUE)[[1]]
expr <- log2(exprs(gset))
batch <- as.factor(phenoData(gset)$"batch:ch1")
sex <- as.factor(phenoData(gset)$"Sex:ch1")

# 用于可视化的移除批次
design <- model.matrix(~ sex)
expr_adj <- removeBatchEffect(expr, batch = batch, design = design)

# PCA 对比
par(mfrow=c(1,2))
plotPCA(expr,    col=batch, main="原始数据")
plotPCA(expr_adj, col=batch, main="批次矫正后")

我在项目里用这个流程检查过多次矫正效果。如果 PCA 图上批次聚集现象明显减少，说明矫正有效。

4.4 避坑指南

我这些年踩过的坑，整理成列表，你能少走弯路：

• 忘记传 design：等于只做了批次平均校正，生物学差异被当作批次抹掉。
• 批次与条件完美混杂：这种设计本身无解，removeBatchEffect 也救不了。只有重新设计实验。
• 对非 log2 数据使用：函数默认数据是正态分布的，如果你用的是原始 count，最好先用 cpm 或 voom 转换。
• 批次因子顺序敏感：我曾经因为 batch 和 batch2 顺序写反，导致第二批次被当作主要批次，结果矫正后批次效应更明显。
• 直接用于差异分析：前面已经说过，这是最常见的错误。最好在脚本里画个图提醒自己。