批次效应的识别：定量指标

前面我们聊了怎么用可视化工具“看”出批次效应。但说实话，光靠眼睛看是不够的。

你想想看，一张PCA图说两个批次分开了，那到底分开了多少？这个“分开了”能不能算个具体的数字出来？

嗯，这就是我们今天要聊的核心——定量指标。我个人习惯用两个指标来衡量批次效应的严重程度：批次差异指数和基因表达变异。

有了数字，我们才好做决定。到底是“轻微影响，不用管”，还是“重大干扰，必须校正”。

指标一：批次差异指数

批次差异指数，英文叫 Batch Divergence Index，简称 BDE。

这玩意儿怎么理解呢？我习惯把它看作是一个“距离感”的度量。

计算思路其实不复杂：

• 先算同一个批次内部，样本之间的平均相似度。
• 再算不同批次之间，样本的平均差异度。
• 用一个比率把这两个结果糅合在一起。

BDE值越大，说明批次之间的结构差异越明显。曾经我手头有个单细胞数据，跑完BDE一看，数值高达0.89，那基本不用怀疑，批次效应严重到必须校正。

具体的计算公式，我们用一个简化版来理解：

# 伪代码：BDE 简化计算思路
batch_divergence = 0
for each gene:
    between_batch_var = variance_of(gene_means_per_batch)
    within_batch_var  = mean_of(variance_within_each_batch)
    gene_score = between_batch_var / within_batch_var
    batch_divergence += gene_score

# 最终 BDE 通常需要归一化到 [0, 1] 区间
BDE_final = normalize(batch_divergence, across_all_genes)

你可能会问，BDE有没有现成的包？

有的。比如 R 包 variancePartition 里面的 fitExtractVarPartModel 函数，就能算出类似的指标。我们稍后会展示具体用法。

指标二：基因表达变异

第二个指标，我把它称为“基因表达变异”，核心是计算变异系数 CV。

为什么用CV？因为基因表达量本身有均值的差异，高表达基因的方差天然就大。直接用方差比较不公平。

CV = 标准差 / 均值，它消除了均值的干扰，纯粹反映波动程度。

具体操作步骤：

1. 对每个基因，分别计算它在每个批次内的CV。
2. 把所有基因在每个批次内的CV算出来，画个箱线图。
3. 看不同批次之间，CV的中位数或者分布是否有明显偏移。

举个例子：

# Python 伪代码：计算基因级 CV 并比较批次
import numpy as np
import pandas as pd

def compute_gene_cv(expression_matrix, batch_labels):
    # expression_matrix: 基因 x 样本
    # batch_labels: 样本的批次标签
    results = []
    for batch in np.unique(batch_labels):
        batch_data = expression_matrix[:, batch_labels == batch]
        cv_per_gene = np.std(batch_data, axis=1) / np.mean(batch_data, axis=1)
        results.append({
            'batch': batch,
            'median_cv': np.median(cv_per_gene),
            'mean_cv': np.mean(cv_per_gene)
        })
    return pd.DataFrame(results)

# 调用示例
# df_cv = compute_gene_cv(expr_data, batch_vec)
# print(df_cv)

输出的表格类似这样：

你看，如果Batch_2的中位CV明显高于其他批次，那说明这个批次要么存在技术偏差，要么样本质量有问题。我在一个癌症转录组项目中就碰到过这种情况——Batch_2的RNA提取时间比别的批次晚了两个月，结果CV值直接飙上去。

结合BDE和CV，我们就能对批次效应做个全面的定量评估。BDE看宏观结构差异，CV看微观基因波动。两者互相印证，决策起来心里才有底。