批次效应检测方法：可视化与统计检验

做批次校正之前，你得先确认一个问题：到底有没有批次效应？

我刚开始做多批次的单细胞数据整合时，上来就一股脑地跑Harmony，结果校正完之后发现生物信号被抹掉了不少。后来才明白——检测比校正更重要。你连有没有效应都不清楚，怎么敢随便调参数？

我个人习惯是先用可视化方法瞄一眼全局，再用统计检验给个量化结论。两个层面的证据对上号了，再动手校正。这一章我就把这两类方法拆开讲讲。

可视化方法：先看个大概

可视化方法的本质是降维。把高维的基因表达数据压到2D或3D空间，然后按批次着色，看看不同批次的样本是不是混在一起。如果它们各自抱团、分得清清楚楚，那批次效应八九不离十。

常用的三个工具：PCA、t-SNE、UMAP。各有各的脾气。

1. PCA（主成分分析）

PCA是最老实的一个。它是线性降维，速度快，解释性强。适合先看一眼数据的整体结构。

2. t-SNE（t分布随机邻域嵌入）

t-SNE是非线性的，擅长把局部结构展示出来。但它有个大问题：随机性。每次运行结果都不一样，而且对困惑度（perplexity）参数敏感。

我记得有次分析一个免疫细胞数据集，明明PCA看着批次效应不大，结果t-SNE画出来，各批次的T细胞居然各自缩成一团。后来调整了困惑度，又正常了。你看，这东西不能全信。

3. UMAP（均匀流形逼近与投影）

UMAP是目前最流行的选择。它保留了更多全局结构，速度也比t-SNE快，结果相对稳定。我个人习惯用UMAP做主图。

不过UMAP也有参数坑。比如n_neighbors设太小，局部细节会被放大，批次效应容易被"画"出来。我一般设为15～30，先跑一次，再用不同参数验证。

统计检验方法：用数字说话

可视化是感性认识，统计检验才能给出"是否显著"的硬证据。

我最常用的两个指标：ANOVA的F统计量和R²（决定系数）。

1. ANOVA（方差分析）

ANOVA的核心思想：把总方差分解成"批次造成的方差"和"随机误差"。如果批次间差异远大于批次内差异，F值就大，p值就小。你可以对每个基因都跑一个单因素ANOVA，然后看有多少基因的p值小于0.05。

# 对每个基因跑单因素ANOVA（伪代码示意）
import numpy as np
from scipy import stats

genes = expression_matrix.T  # 基因×样本
batch_labels = meta['batch']

for gene in genes:
    groups = [gene[batch_labels == b] for b in unique_batches]
    f_stat, p_val = stats.f_oneway(*groups)
    # 记录f_stat, p_val

2. R²（决定系数）

R²衡量"批次这个变量能解释多少方差"。你把批次作为一个分类变量，对每个基因做线性回归，看加入批次变量后R²增加了多少。或者更直接地，用批次做响应变量？不对，应该是将基因表达量作为响应变量，批次作为自变量。

我个人习惯计算每个基因的批次相关方差占比：

# 伪代码：计算每个基因的R²
from sklearn.linear_model import LinearRegression

r2_list = []
for gene in genes:
    X = pd.get_dummies(batch_labels)  # one-hot
    y = gene
    model = LinearRegression().fit(X, y)
    r2 = model.score(X, y)
    r2_list.append(r2)
mean_r2 = np.mean(r2_list)

如果平均R²超过0.1，就说明批次效应不可忽视。我见过一些数据平均R²能到0.4，那基本就是批次为主、生物学差异为辅了。

方法对比总结