生信批次效应统计检测与多算法校正 - 第1节
一、什么是批次效应？为什么我总提它？
做组学数据分析的朋友，或多或少都遇到过这种场景：明明是同一种样本，只是分了两批送测，结果一画PCA图，两批样本各自抱团，离得老远。这可不是什么好现象，这就是批次效应，一个让我头大了好多年的问题。
说白了吧，批次效应就是实验过程中引入的系统性技术偏差。它不是生物学差异，却常常和生物学差异搅在一起，干扰我们的判断。
1.1 它到底长什么样？
我习惯从三个维度去观察它。这张图概括了批次效应的全貌：
比方说，你上午测一批样本，下午换了个操作员接着测。两批数据放在一起一比较，差异可能根本不是生物学上的，而是操作员习惯造成的。这就是批次效应。
1.2 三类常见的批次效应
| 类型 | 典型来源 | 我遇到过的例子 |
|---|---|---|
| 实验批次 | 不同的实验日期、试剂批次、仪器、操作员 | 两批RNA-seq建库用了不同批次的接头 |
| 技术平台 | 芯片类型、测序平台、质谱型号 | 一批样本用Illumina，另一批用MGI |
| 时间效应 | 几个月甚至几年的数据集整合 | 公共数据库中不同年份的芯片数据 |
1.3 为什么我们必须要重视它？
好多新手朋友问过我：好像有点偏差而已，问题到底有多大？我来给你算笔账。
- 假阳性飙升： 我见过一个项目，不校正批次直接做差异分析，有60%的显著性基因其实是批次标签在起作用，和生物学本身没关系
- 假阴性误判： 更隐蔽的是，真正的差异信号可能被批次效应掩盖掉，让你漏掉关键基因
- 模型崩塌： 你要是用这批数据训练分类器，放到新批次上，准确率直接跌到50%以下，跟抛硬币一样
二、识别批次效应——先确认敌人长什么样
知道了概念还不够，得学会怎么识别出它来。我个人习惯从上手最快的方法开始。
2.1 第一步：可视化的直觉判断
最常用的办法就是PCA图。画一个PCA，按批次上色。如果不同批次在PC1或PC2上明显分开，那基本上可以断定存在批次效应。
# 一个简单的PCA检查，用Python的scikit-learn
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
# X 是样本x特征矩阵，batch_labels 是批次标签
pca = PCA(n_components=2)
pca_result = pca.fit_transform(X)
plt.figure(figsize=(8,6))
for b in set(batch_labels):
idx = batch_labels == b
plt.scatter(pca_result[idx,0], pca_result[idx,1], label=f'Batch {b}', alpha=0.6)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend()
plt.show()同理，也可以用热图观察样本聚类树，看是否先按批次分叉。
2.2 第二步：量化评价——用统计指标说话
光看图还不够，万一感觉不准呢。我习惯用两种方式量化：
- ANOVA / Kruskal-Wallis： 对每个特征（基因），看它在不同批次间的差异是否显著。如果太多基因的p值小于0.05，说明批次效应普遍存在
- 批次效应比例： 计算PCA前几主成分中有多少变异可以归因于批次标签——用线性模型 R² 来衡量
# 计算每个基因的批次效应p值，示意
from scipy.stats import kruskal
import numpy as np
pvals = []
for gene in range(X.shape[1]):
groups = [X[batch_labels == b, gene] for b in set(batch_labels)]
_, p = kruskal(*groups)
pvals.append(p)
# 看有多少基因的批次p < 0.01
prop_significant = np.mean([p < 0.01 for p in pvals])
print(f'受批次显著影响的基因比例: {prop_significant:.2%}')三、校正的大方向——不是所有方法都万能
先别急着上代码。我记得当年第一次处理批次效应时，直接找了个ComBat函数就跑，结果把生物学差异也给抹掉了，气得我直跺脚。所以搞清楚策略方向很重要。
3.1 按设计阶段分
| 策略 | 适用场景 | 关键点 |
|---|---|---|
| 预防性设计 | 新实验规划阶段 | 随机化样本、交叉分组、统一操作 |
| 统计校正 | 已有数据整合 | 选择合适的模型，避免过校正 |
预防永远比事后校正省心。但现实是，我们拿到的很多是公共数据，或实验已经做完了，那就只能走统计校正这条路。
3.2 按数据特性分
不同的数据类型，要选不同的刀。你想想看，RNA-seq的count数据分布，和芯片的连续信号能一样吗？
- 芯片/连续信号数据： 我常用 ComBat (基于经验贝叶斯)、Limma 中的批次校正
- RNA-seq 计数数据： 推荐 RUVseq (基于无监督因子)、SVA、或者用 limma-voom + 批次作为随机效应
- 甲基化/ChIP-seq： 常通过分位数归一化配合线性模型处理
四、一点小建议——记在最后
我见过太多人直接拿着数据就开始跑校正，连基本的PCA检查都没做。嗯，这不比做菜，不能上来就把所有调料倒进去。
我个人习惯是每拿到一个多批次数据集，先走一遍这个流程：
- 画PCA热图，看有没有明显分层
- 量化批次效应比例
- 判断严重程度——不严重的也许不需要校正，很多分析对轻微批次效应是鲁棒的
- 选择合适的校正方法，保留原始数据备份
- 校正后对比验证：生物差异是否还在？批次差异是否减小？
这门课接下来的章节，我会一步步拆解每一步的统计学原理和代码实现。从PCA的细节讲到ComBat的贝叶斯推导，再到多平台整合的进阶技巧。每一步，我都会把我在项目里摔过的跤、绕过的路讲给你听。