📚 课程目录第3章 PCA用于批次效应检测

第三节：PCA用于批次效应检测

聊完了理论基础和原理，咱们来点实际的。数据拿到手，最怕什么？最怕你辛辛苦苦分析出来的差异，其实是批次差异。

我个人习惯，在整合数据之前，先跑个PCA看看。这一步就像体检前的“血常规”，能快速帮你判断数据有没有明显的“炎症反应”。说白了，PCA能把高维数据压缩到二维或三维，让我们用肉眼直接看到样本的“亲疏远近”。

核心概念：PCA通过最大方差方向投影，将主成分表达矩阵可视化。如果同批次样本聚在一起、不同批次样本分散开，那就是典型的批次效应。

3.1 为什么选择PCA做初筛

我不是说PCA是万能的。但相比t-SNE或UMAP，PCA有个巨大的优点：它不依赖随机初始化，结果完全可重现。

你想想看，你做一个分析，结果每次跑出来都不一样，怎么跟领导汇报？PCA就没这个问题，给定同样数据，每次投影结果一模一样。

另外PCA的计算速度极快。我手头有个二十万细胞的数据集，PCA算下来也就几十秒。换成UMAP, 嗯，至少得等个几分钟。

特性	PCA	UMAP / t-SNE
计算速度	极快（秒级）	较慢（分钟级）
结果确定性	完全确定	每次运行可能不同
批次检测能力	强（全局结构）	强（局部结构）
适用场景	初筛、全局趋势	精细聚类可视化

个人经验：在正式批次校正前，先用PCA画一张图。如果批次已经混在一起，那校正的必要性就不大；如果批次完全分离，那就要想清楚用什么方法校正了。

3.2 实际操作：跑一个PCA可视化

咱们用Python写个简单例子。假设你手头已经有一个合并后的表达矩阵，行是基因，列是样本。

import scanpy as sc
import matplotlib.pyplot as plt

# 读取数据
adata = sc.read_h5ad('merged_data.h5ad')

# 基础预处理：高变基因筛选 + 标准化 + PCA
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata.raw = adata
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=30)

# 用批次信息着色
sc.pl.pca(adata, color='batch', 
          title='PCA colored by batch',
          legend_loc='right margin')

跑完之后，你会看到一张类似散点图的画面。每个点代表一个样本或一个细胞，颜色代表不同批次。

3.3 怎么解读PCA图？

这里有个经验之谈。我刚开始做的时候，看见散点有重叠就松一口气，看见分离就紧张。但实际上没那么简单。

我建议你从三个维度去判断：

PC1 vs PC2 是否分开：这是主成分轴，解释方差最多。如果不同批次在PC1方向就分开了，说明批次效应对全局变异贡献极大。
批次内部是否聚集：同批次样本挤在一起，形成小团块，说明批次内一致性很好，但很可能引入了批次特异信号。
生物学分组是否跨批次混合：这个最关键。比如你有一个处理组和对照组，各分布在两个批次中。如果处理组样本跨批次聚在一起，说明生物学信号大于批次效应，问题不大；如果处理组样本按批次分开成两坨，那麻烦就大了。

说白了，你要把PCA图当成一个“侦察报告”——看看到底是批次在驱动分组，还是生物学变量在驱动分组。

⚠️ 常见误判：我曾遇到一个数据集，PCA图上批次看起来混在一起，我就松懈了。结果下游差异分析发现，大量假阳性基因其实都是批次相关基因。后来仔细查才发现，那批数据PCR深度不一致，在PC3/PC4上才暴露出来。所以只看PC1/PC2是不够的，至少要看前五个主成分。

3.4 批次分离的判定标准

咱们不能光凭感觉，得有个量化参考。我个人总结了一套简单的判定体系：

观测现象	严重程度	处理建议
批次完全重叠，无明显分群	✅ 没问题	可直接合并分析
批次略有偏移，但生物学分组更清晰	⚠️ 需关注	可用线性校正，如Combat
批次明显分离，盖过生物学分组	❌ 严重	必须校正，建议用深度学习方法
批次完美分离成两团，无任何重叠	🚫 不可用	检查实验设计是否出了问题

还有一种更严谨的方法：用ANOVA检验每个主成分上传款与批次标签的相关性。简单点，计算批次标签在不同主成分上的解释方差比例（R²）。

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

# 假设pca_pcs是你的PCA结果前10个主成分，batch_labels是批次标签
r2_scores = []
for i in range(10):
    model = LinearRegression()
    model.fit(pd.get_dummies(batch_labels), pca_pcs[:, i])
    r2_scores.append(model.score(pd.get_dummies(batch_labels), pca_pcs[:, i]))

# 打印每个PC上批次解释的方差比例
for i, r2 in enumerate(r2_scores):
    print(f"PC{i+1}: R² = {r2:.3f}")

R²接近1，说明该主成分完全由批次驱动。如果前几个主成分的R²都大于0.8，那就要高度重视了。

嗯，PCA检测批次效应这招，说白了就是“让数据自己说话”。不做这一步直接往下走，很容易掉坑里。

画完PCA图，如果看到批次分离严重，不要慌。后面的章节会覆盖常见的校正方法。但至少，你现在知道第一步该做什么了。

我踩过的坑：有一回分析一个scRNA-seq数据集，PCA图上批次看起来混得还行，我就直接跑差异分析了。结果发现所有差异表达基因都是核糖体蛋白，一看是不同批次的细胞裂解效率不一样导致的。后来重新跑PCA时我选择了更多高变基因，问题才暴露。所以记得，高变基因的选择会影响PCA的结果，不妨多试几组参数。