第四章 基于统计检验的检测方法：ANOVA、Kruskal-Wallis、ComBat

各位做药物批次效应检测，最怕什么？

最怕你不知道有批次效应，或者知道了却拿不准它的严重程度。

这一章咱们就来聊聊怎么用统计检验给批次效应做个"体检"——ANOVA、Kruskal-Wallis和ComBat。这三招我在项目里反复用过，各有各的脾气。

4.1 为什么先用统计检验打头阵？

说实话，很多人一上来就上深度学习模型，结果数据本身就有结构性偏差，模型学到的全是批次差异。我个人习惯是：可视化之后，立刻上统计检验，先定性再定量。

统计检验能帮你回答三个核心问题：

  1. 有没有差异？——批次间的均值或分布是否显著不同
  2. 差异多大？——效应量的大小
  3. 哪些特征受影响？——具体是哪些基因/化合物被批次"污染"了

你想想看，如果没有这一步就直接上校正，很容易过度校正或校正不足。我在一个RNA-seq项目里就吃过这个亏，当时跳过ANOVA直接跑ComBat，结果把真实生物学差异也给抹掉了，后来重跑才补救回来。

4.2 ANOVA：参数检验的"老大哥"

ANOVA（方差分析）是检验多个批次均值是否相等的经典方法。它假设数据服从正态分布且方差齐性——说白了就是各批次的"身材"差不多，只是"站位"可能不同。

适用场景：连续型数据、批次数量≥2、正态分布假设基本成立

ANOVA的核心逻辑是这样的：

  • 总变异 = 批次间变异 + 批次内变异
  • 如果批次间变异远大于批次内变异，F值就大，p值就小
  • p < 0.05 说明至少有一个批次和其他批次不一样

我一般先用ANOVA做个全局扫描。代码很简单：

import pandas as pd
import numpy as np
from scipy import stats

# data: 行为样本，列为特征，batch: 批次标签
def batch_effect_anova(data, batch_labels):
    results = []
    for col in data.columns:
        groups = [data[col][batch_labels == b] for b in np.unique(batch_labels)]
        f_stat, p_val = stats.f_oneway(*groups)
        results.append({'feature': col, 'F_stat': f_stat, 'p_value': p_val})
    return pd.DataFrame(results).sort_values('p_value')

# 举个实际例子
df_anova = batch_effect_anova(expression_data, batch)
print(df_anova.head(10))

跑完之后，我会重点看p值最小的那些特征——它们被批次效应影响最严重。

注意：ANOVA对异常值非常敏感。如果你的数据里有离谱的离群点，先做截尾处理或改用非参数方法。我曾经在一个代谢组学数据上直接跑ANOVA，结果被两个异常样本带偏了，差点误判。

4.3 Kruskal-Wallis：非参数检验的"稳健派"

遇到数据不满足正态分布，或者样本量小、方差不齐，怎么办？

这个时候Kruskal-Wallis检验就派上用场了。它不关心数据的具体分布，只关心数据的秩（排名）。说白了，就是把所有数据排个队，看各批次的排名分布是否差不多。

适用场景：连续型或有序数据、分布未知或非正态、样本量差异大

Kruskal-Wallis的H统计量计算：

  • 把所有数据混合排序，赋秩次
  • 计算每个批次的平均秩次
  • 如果各批次的平均秩次差异大，H值就大，p值就小
from scipy.stats import kruskal

def batch_effect_kruskal(data, batch_labels):
    results = []
    for col in data.columns:
        groups = [data[col][batch_labels == b] for b in np.unique(batch_labels)]
        h_stat, p_val = kruskal(*groups)
        results.append({'feature': col, 'H_stat': h_stat, 'p_value': p_val})
    return pd.DataFrame(results).sort_values('p_value')

# 小样本场景下我更喜欢用这个
df_kw = batch_effect_kruskal(expression_data, batch)
print(df_kw[df_kw['p_value'] < 0.05].shape[0], "个特征存在显著批次效应")

Kruskal-Wallis的优点是稳健，缺点是灵敏度比ANOVA稍低——如果数据真的满足正态分布，ANOVA的检验力更强。我个人习惯是：样本量大且正态时用ANOVA，否则用Kruskal-Wallis。

个人技巧：两个方法都跑一遍，看看结果的一致性。如果ANOVA和Kruskal-Wallis同时指向同一个特征，那这个特征的批次效应基本上是"实锤"了。

4.4 ComBat：经验贝叶斯的"矫正大师"

好，前面两个方法都是在"检测"问题。那检测出来之后呢？

ComBat就是目前最流行的批次效应校正方法之一。它由Johnson和Li在2007年提出，用经验贝叶斯估计每批次的均值和方差偏移，然后把这些偏移"掰正"。

ComBat的核心假设：各批次的特征分布属于同一族，只是位置和尺度不同

ComBat的校正过程可以概括为：

  1. 估计每批次的均值偏移和方差偏移
  2. 用经验贝叶斯"收缩"估计值，防止过度校正
  3. 从原始数据中减去批次偏移，保留生物学变异
# 使用pycombat库
from pycombat import pycombat

# expression_data: 基因表达矩阵 (基因为行，样本为列)
# batch: 批次标签，长度与样本数一致
adjusted_data = pycombat(expression_data, batch)

# 校正前后对比
print("校正前批次间均值差异:", 
      expression_data.groupby(batch).mean().std().mean())
print("校正后批次间均值差异:", 
      adjusted_data.groupby(batch).mean().std().mean())

实战建议：ComBat对样本量较小的批次会做"收缩"，这恰恰是它的优势。如果你有一个批次只有2-3个样本，ComBat会谨慎地调整，而不是硬掰。

ComBat有很多变体，比如ComBat-Seq处理RNA-seq的count数据，ComBat-GAM处理非线性时间序列数据。我一般先用标准ComBat，如果效果不理想再换变体。

避坑指南：ComBat不能处理缺失值。如果你数据里有NaN，记得先做填充或删除。我曾经因为忘了这步，跑出来的结果是全NaN，排查了半天才发现。

4.5 三种方法的对比与选型

这三种方法我都反复用过，很多时候是组合使用的。下面这张表是我自己总结的选型指南：

维度 ANOVA Kruskal-Wallis ComBat
目的 检测差异 检测差异 校正差异
数据类型 连续型，正态 连续型/有序，非正态 连续型/离散型
批次数量 ≥2 ≥2 ≥2
样本量要求 每批≥10 每批≥5 每批≥2
输出结果 F值, p值 H值, p值 校正后数据
是否校正 否 否 是
稳健性 低（对异常值敏感） 中（基于秩次） 高（经验贝叶斯收缩）

说说我的工作流：先用ANOVA或Kruskal-Wallis做全特征扫描，筛出受批次效应影响严重的特征；然后用ComBat做整体校正；最后再用同样的检验方法验证校正效果。这样就形成闭环了。

统计检验与校正流程 统计检验与校正工作流 原始多批次数据 数据是否正态 / 方差齐？ 是 ANOVA 检测 否 Kruskal-Wallis 检测 筛选显著批次效应特征 ComBat 经验贝叶斯校正 校正后去批次数据 先检测，后校正

4.6 组合使用的实战案例

我来说一个具体的例子。之前做一个多中心药物代谢组学研究，3个实验室，每个实验室测了50个样本。我先做了主成分分析，发现PC1和PC2基本是按实验室分的——典型的批次效应。

第一步，用Kruskal-Wallis检测每个代谢物：

kw_results = batch_effect_kruskal(metabolite_data, lab_batch)
hits = kw_results[kw_results['p_value'] < 0.01]
print(f"检测到{hits.shape[0]}个代谢物有显著批次效应（p<0.01）")

第二步，用ComBat校正：

from pycombat import pycombat
corrected_data = pycombat(metabolite_data.T, lab_batch).T

第三步，用Kruskal-Wallis再次验证：

kw_after = batch_effect_kruskal(corrected_data, lab_batch)
hits_after = kw_after[kw_after['p_value'] < 0.01]
print(f"校正后剩{hits_after.shape[0]}个代谢物仍有显著批次效应")

结果：从187个显著代谢物降到了22个，校正效果很理想。那22个我进一步分析，发现其中有几个是真正的生物学差异。

经验之谈：校正前后对比p值的分布变化，比单看p值更有意义。如果校正后p值分布整体右移（变大了），说明ComBat确实在起作用。

4.7 选型建议与注意事项

总结几条我踩过坑之后得出的建议：

  • 小样本（每批<5）：首选Kruskal-Wallis做检测，ComBat做校正时注意收缩参数
  • 大样本（每批>30）：ANOVA + ComBat组合，灵敏度高
  • 高维数据（特征数>样本数）：直接上ComBat，检测步骤可以简化
  • 非连续数据（count、有序）：Kruskal-Wallis + ComBat-Seq

重要提醒：统计检验只能告诉你"有没有"，不能告诉你"该不该校正"。有些批次差异可能本身就包含了生物学意义——比如不同中心收集的病人亚型不同。你想想看，盲目校正反而会抹掉真实发现。

这一章的内容就到这里。ANOVA、Kruskal-Wallis、ComBat这三板斧，是我每次处理批次效应的标配工具。检测和校正缺一不可，组合使用才能把数据洗干净。

← 上一章 📖 返回目录 下一章 →