📚 课程目录第4章批次效应检测技术

4、批次效应检测技术：PCA、t-SNE与热图

光说要校正批次效应，但你得先知道它到底存不存在，对吧？

我个人习惯把批次效应检测比作“看病先做检查”。你不能稀里糊涂地就上模型，先得看看这批次的“病灶”有多明显。今天我们就聊聊三种最常用的检查手段：PCA、t-SNE 和热图。

4.1 主成分分析（PCA）：快速概览

PCA，说白了就是把高维数据压到我们能看的二维或三维上。你想啊，一个QSAR数据集动辄上千个分子描述符，人眼根本看不透。PCA就帮我们抓住最“主”的几个方向，看看能不能把不同批次的样本自动分开。

我在早期的一个项目中，用PCA一跑，发现前3个主成分只解释了不到35%的方差。当时我还以为模型选错了，后来才发现是批次效应太强，把真实的结构-活性关系给淹没了。

我的一点经验：先跑不加标签的PCA，看看散点图的颜色是按批次分布，还是按活性值分布。如果批次颜色分得比活性值还清楚，那基本可以确定有严重的批次效应。

下面是一个用scikit-learn实现的基础PCA分析片段，我习惯先标准化再用PCA：

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# X是分子描述符矩阵，batch_labels是批次标签列表
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

plt.figure(figsize=(8,6))
for batch in set(batch_labels):
    idx = [i for i, b in enumerate(batch_labels) if b == batch]
    plt.scatter(X_pca[idx, 0], X_pca[idx, 1], label=f'Batch {batch}')
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
plt.legend()
plt.show()

4.2 t-SNE：细节放大镜

PCA是线性的，有时候批次效应藏得比较深，比如不同批次的实验系统误差只在局部特征上有差异。这时候我会换用t-SNE。t-SNE擅长保留数据的局部结构，能把那些“纠缠”在一起但实际属于不同批次的点给揪出来。

我记得有一次分析一个跨实验室的QSAR数据集，PCA看着挺正常的，所有点混在一起。但用t-SNE一画，哇，立马分出三簇，正好对应三个实验室。你说巧不巧？这就是局部差异在作祟。

但t-SNE有个坑：它是个随机算法，每一次跑结果都可能不一样。我建议你多跑几次，或者固定一下随机种子（random_state），不然容易被个别“随机抖动”带偏判断。

注意：t-SNE对参数很敏感。perplexity太大（比如超过样本数的一半），结果会过于平滑；太小又容易形成离散的小岛。一般我从30开始试，根据样本量调整。

from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X_scaled)

plt.figure(figsize=(8,6))
for batch in set(batch_labels):
    idx = [i for i, b in enumerate(batch_labels) if b == batch]
    plt.scatter(X_tsne[idx, 0], X_tsne[idx, 1], label=f'Batch {batch}', alpha=0.7)
plt.title('t-SNE Visualization of Batches')
plt.legend()
plt.show()

4.3 热图：全局“大扫描”

PCA和t-SNE都是降维后的可视化，有时候你更想知道：到底是哪些描述符或者哪些样本受批次效应影响最大？这时候就轮到热图出马了。

我常做的是“批次-描述符相关性热图”。简单来说，就是对每个描述符计算它与批次标签之间的相关性（比如用ANOVA的F值，或者直接用Kruskal-Wallis检验的p值）。然后把结果画成一个矩阵，一目了然：颜色越深，说明该描述符受批次影响越大。

避坑指南：我只曾用过一次把所有描述符不筛选就直接画热图，结果图太密，根本看不出规律。后来我学会了先按方差排个序，或者只取p值最小的top 50个描述符，图才清晰起来。

import seaborn as sns
import numpy as np
import pandas as pd

# 假设你有描述符矩阵X（n_samples, n_features）和批次标签batch_labels（0,1,2,...）
# 这里用Kruskal-Wallis检验示例
from scipy.stats import kruskal

def batch_correlation_matrix(X, batch_labels):
    p_values = []
    for i in range(X.shape[1]):
        groups = [X[batch_labels == g, i] for g in np.unique(batch_labels)]
        _, p = kruskal(*groups)
        p_values.append(-np.log10(p + 1e-10))  # 取负对数方便可视化
    return np.array(p_values)

# 取top 50个描述符
n_top = 50
top_indices = np.argsort(batch_correlation_matrix(X_scaled, np.array(batch_labels)))[-n_top:]

# 构建数据框用于热图
heatmap_data = X_scaled[:, top_indices].T
heatmap_df = pd.DataFrame(heatmap_data, 
                          columns=[f'Sample_{i}' for i in range(X_scaled.shape[0])],
                          index=[f'Desc_{i}' for i in top_indices])

sns.heatmap(heatmap_df, cmap='viridis', xticklabels=False, yticklabels=True)
plt.title('Top 50 Descriptors: Batch Effect Profile')
plt.show()

4.4 综合检测流水线

实际干活的时候，我很少只用一种方法。我的策略是这样的：

先用PCA快速扫一眼：如果前两主成分就能明显分开批次，直接确认有批次效应，不用再多犹豫。
再用t-SNE做补充：尤其当PCA结果模糊、样本点混在一起但你又怀疑隐藏的批次差异时，t-SNE能帮你还原局部结构。
最后用热图定位：把那些“重灾区”描述符揪出来，看看是不是某些化学特征在跨批次时测量不稳定。

下面这张图总结了这个流程，你可以跟着走一遍：

4.5 实战：一个完整的检测案例

我给你一个我最近做的实际例子。数据集来自两个实验室的同一份化合物库，各测了一批活性值。我一开始不放心，就用PCA看：

PCA结果：PC1 + PC2解释了约55%方差。散点图上，两个实验室的样本在PC1方向上有轻微偏移，但没有完全分开。
t-SNE结果：设置perplexity=40，两个实验室的点几乎成了两个独立的簇。这就确定了，肯定有批次效应。
热图分析：我对所有描述符做了Kruskal-Wallis检验，发现约30%的描述符在实验室间存在显著差异（p<0.01）。主要集中在分子量、LogP和一些拓扑特征上。这就解释了为什么活性值会受批次影响——跨实验室的测量系统对某些化学空间敏感。

所以你看，这套组合拳打下来，不管问题多么隐蔽，基本都能发现。

一句话总结：PCA给你全局感，t-SNE给你细节，热图给你证据链。三条腿走路，批次效应无所遁形。

4.6 常见误区与避坑

误区	正确做法
只看PCA就下结论	PCA线性假设强，批次效应若非全局线性时会被掩盖。必须结合t-SNE或热图。
t-SNE跑一次就信	t-SNE是随机初始化，单次结果不可靠。至少跑3次，或者固定随机状态观察稳定性。
热图颜色深浅靠肉眼判断	最好加上p值筛选，或者用显著性标记（如星号），避免主观解读。
样本量太小还硬上t-SNE	样本数小于50时，t-SNE的perplexity范围很窄，效果差。建议换成UMAP或直接PCA。

以上这些，就是我个人在项目中常用的批次效应检测三板斧。你下次拿到一个新数据集，不妨也按这个流程走一遍，我相信你也会有收获的。