📚 课程目录第2章批次效应解析

2. 批次效应解析：什么是批次效应？常见来源及可视化检测方法

做单细胞数据分析，第一个绕不开的坑就是批次效应。我个人习惯把它叫做“数据里的隐形幽灵”——你看不见它，但它能把你的聚类结果搞得一团糟。

什么叫批次效应？简单说，就是同样的生物学样本，因为实验操作的时间、试剂批号、操作人员不同，导致测序数据出现系统性的非生物学差异。说白了，就是技术噪声被错误地当成了生物学信号。

我最早在项目中遇到这个问题，是处理三个不同时间点采集的PBMC数据。当时直接合并后做UMAP，发现两个样本完全分开了——我还以为是发现了新的细胞亚群。后来一查，原来是操作者换了试剂盒批次。嗯，从那以后我养成了一个习惯：任何多批次数据，先做批次效应检测，再谈下游分析。

2.1 批次效应的常见来源

批次效应的来源，我归纳成三类。你想想看，实验流程里每一个环节都可能引入偏差。

来源类别	具体因素	影响程度
试剂/耗材	试剂盒批号、酶活性差异、磁珠批次、水纯度变动	高——往往能解释30%以上的技术变异
操作者	手动移液误差、细胞计数差异、设置参数习惯不同	中——在跨实验室数据中尤其明显
设备/时间	测序仪型号差异、运行周期、样本存储时间、季节温度	高——时间跨度越大，批次效应越严重

我曾遇到过一个案例：同一批细胞，上午测和下午测，PC1上的分布竟然完全错开了。为什么会这样？后来排查发现，上午的样本在室温下多放了两个小时，细胞活性下降导致基因表达谱偏移。所以，时间本身就是一种隐性的批次变量。

2.2 可视化检测方法

怎么判断你的数据里有没有批次效应？我个人建议先做三步检测，就像去医院体检一样，不能只看一个指标。

方法一：UMAP/t-SNE着色

最直观的方式，就是把样本合并后降维，按批次标签着色。如果不同批次的细胞在低维空间里形成了明显分离的簇——哪怕这些簇的细胞类型相同——那就是批次效应在作祟。

关键判断准则：如果同一细胞类型（比如CD14+单核细胞）在UMAP上堆成两堆，而这两堆恰好对应两个批次，那么批次效应比你想象的严重。

代码示例如下：

import scanpy as sc
import matplotlib.pyplot as plt

adata = sc.read("your_data.h5ad")
# 假设 adata.obs 有 'batch' 和 'cell_type' 列

# 按批次着色
sc.pp.neighbors(adata, n_neighbors=15)
sc.tl.umap(adata)
sc.pl.umap(adata, color=['batch', 'cell_type'],
           ncols=2, wspace=0.4)

小技巧：把批次和细胞类型同时展示，如果‘batch’的分离模式跟‘cell_type’的分离模式完全不同，基本可以断定批次效应主导了低维结构。

方法二：基因表达相关性热图

按批次对样本做基因表达谱的相关性矩阵。理想情况下，同一生物学重复的样本应该高度相关（r > 0.95）。如果不同批次的样本相关系数明显低于批次内部，说明有系统性偏移。

import seaborn as sns
import numpy as np

# 按批次计算平均表达谱
batch_means = []
for batch in adata.obs['batch'].unique():
    sub = adata[adata.obs['batch'] == batch]
    batch_means.append(np.mean(sub.X.toarray(), axis=0))
batch_means = np.array(batch_means)

# 相关性矩阵
corr = np.corrcoef(batch_means)
sns.heatmap(corr, annot=True, xticklabels=adata.obs['batch'].unique(),
            yticklabels=adata.obs['batch'].unique(),
            cmap='coolwarm', vmin=0, vmax=1)

注意：如果批次内部相关系数本身就低（比如 < 0.9），那可能不仅仅是批次效应，还涉及样本质量或细胞数量不足的问题。这时候先别急着调校正参数，先回头检查上游数据预处理。

方法三：PCA变异性分析（我偏爱的“硬核”方法）

用主成分分析，看前几个PC的载荷和批次标签的相关性。如果一个PC的高度变异完全由批次解释（比如PC1的方差中80%以上来自批次标签），那这个PC就是“批次PC”。

from sklearn.decomposition import PCA
import pandas as pd

pca = PCA(n_components=10)
pcs = pca.fit_transform(adata.X.toarray())

# 计算每个PC与批次标签的ANOVA F统计量
import statsmodels.api as sm
from statsmodels.formula.api import ols

df = pd.DataFrame(pcs, columns=[f'PC{i+1}' for i in range(10)])
df['batch'] = adata.obs['batch'].values

for i in range(5):
    model = ols(f'PC{i+1} ~ C(batch)', data=df).fit()
    anova = sm.stats.anova_lm(model, typ=2)
    print(f'PC{i+1}: F值={anova["F"]["C(batch)"]:.2f}, p值={anova["PR(>F)"]["C(batch)"]:.2e}')

如果PC1的F值大于10且p值小于0.001，那基本锁定了——批次效应已经占据了主成分的大部分方差。这时候你就得考虑使用我后面会讲到的校正方法了。

2.3 知识体系总览

下面这张SVG图帮你梳理本节的逻辑链条：从批次效应的来源，到检测方法，再到下一阶段的方向。

2.4 避坑指南 & 经验

不要只看UMAP。它参数敏感（特别受n_neighbors影响），有时候批次效应被强行吞掉了反而让你漏判。我建议配合PCA方差比一起看。
注意样本数量不平衡。如果批次A有5000个细胞，批次B只有200个，UMAP可能会把批次B压成小尾巴，让你误以为是稀有亚群。这种情况被称为“批次偏差放大现象”。
记得保存原始数据。我开始做项目时习惯直接在原始数据上做校正，后来发现不对——你应该先检测，再校正，而且保存一份未校正的备份。我曾经因为直接覆盖了原始数据，导致后续想验证校正效果时已经没有基线了。从那以后，我再也不敢省略备份步骤。

我的经验总结：批次效应检测不是一次性动作。你应该在数据预处理过程中至少做两次检测：第一次在QC后、合并前，第二次在合并后、校正前。两次对比才能判断哪些变异来自细胞本身，哪些来自技术噪声。

2.5 一个小工具：快速诊断函数

下面我整理了一个迷你诊断函数，适合拿到新数据后立刻跑一遍：

def batch_diagnosis(adata, batch_key='batch'):
    """快速诊断批次效应（输出三个可视化）"""
    import scanpy as sc
    import matplotlib.pyplot as plt

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # 1. UMAP按批次着色
    sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
    sc.tl.umap(adata)
    sc.pl.umap(adata, color=batch_key, ax=axes[0], show=False)
    axes[0].set_title('UMAP: 按批次')

    # 2. 相关性热图 - 批次平均谱
    import numpy as np
    import seaborn as sns
    batch_list = adata.obs[batch_key].unique()
    means = []
    for b in batch_list:
        sub = adata[adata.obs[batch_key] == b]
        means.append(np.mean(sub.X.toarray(), axis=0))
    corr = np.corrcoef(np.array(means))
    sns.heatmap(corr, annot=True, xticklabels=batch_list, yticklabels=batch_list,
                ax=axes[1], cmap='coolwarm', vmin=0, vmax=1)
    axes[1].set_title('批次间相关性')

    # 3. PCA方差比（前5PC）
    pca = PCA(n_components=10)
    pcs = pca.fit_transform(adata.X.toarray())
    from sklearn.feature_selection import f_classif
    f_stat = []
    for i in range(5):
        f, pval = f_classif(pcs[:, i].reshape(-1, 1),
                            adata.obs[batch_key].values)
        f_stat.append(f[0])
    axes[2].bar(range(1, 6), f_stat, color='coral')
    axes[2].set_xlabel('PC')
    axes[2].set_ylabel('F统计量')
    axes[2].set_title('PC vs 批次 (F值)')
    axes[2].axhline(y=10, color='gray', linestyle='--', label='阈值')
    axes[2].legend()

    plt.tight_layout()
    return fig

把这段代码保存成 batch_diag.py，拿到新数据直接调用。我个人觉得，这个函数帮你补全了肉眼判断的盲区——尤其是相关性热图和F统计量，能让你量化的看到批次效应的大小。

注意：这个函数假设你的数据已经完成了归一化和高变基因筛选。如果还没做，结果可能会被其他噪声误导。

好了，关于批次效应的定义、来源和检测方法，我们就先聊到这。记得回到你的项目里，先给数据做个“体检”。下一阶段，我们就可以开始讨论怎么把这些讨厌的技术偏差校正掉——但我不会在这里剧透。你先把手头的批次情况摸清楚，这是最重要的第一步。