第4章：PCA在降维中的应用：从高维数据到二维/三维投影

这一章，我们来聊聊 PCA（主成分分析）。我敢说，这是生物信息学里用得最多的降维方法之一。说白了，它就是把一堆乱麻似的高维数据，拧成几条清晰的主干。

我之前跟一个基因表达数据分析项目时，看到同事直接拿原始表达矩阵画热图，结果图上一片噪点，根本看不出样本分组。我当时说：「试试 PCA 吧。」几分钟后，二维散点图上三组样本清清楚楚分开了。这就是 PCA 的威力——把噪音滤掉，把结构露出来。

4.1 为什么要用PCA降维？

先理解一个简单问题：我们为什么非要降维？

• 可视化：人类的眼睛只能理解2D、3D。几千维数据你怎么看？只能靠投影。
• 去噪声：高维数据里很多维度是噪音，PCA能帮我们扔掉。
• 减少过拟合：维度越高，模型越容易「记住噪声」而不是真正模式。
• 加速计算：几十维比几万维算起来快得多。

但 PCA 也有个前提：它假设数据的重要结构存在于线性子空间里。如果你的数据是非线性的（比如时空折叠的那种流形），那 PCA 就有点力不从心了。那时候你可能得考虑 UMAP 或者 t-SNE。

4.2 PCA的数学直觉——三句话讲清楚

我不会在这里推导矩阵特征值分解，毕竟我们不是数学课。但核心思想我可以用三句话概括：

1. 找方向：在数据里找到一个方向，使得所有样本在这个方向上的投影方差最大。这个方向叫「第一主成分」。
2. 找下一个：再找一个与上一个方向正交的方向，让剩余方差最大。这叫「第二主成分」。
3. 重复：直到你拿到足够解释数据方差的主成分，通常选前2个或3个拿来画图。

你看，其实就是一步步地「旋转坐标轴」，让新坐标轴的方向和数据本身最「合拍」。

嗯，这里要注意：PCA 对数据的尺度非常敏感。我之前遇到过一个问题：一个基因的表达值范围在 0~10，而另一个基因因为检测方法不同，范围在 10000~20000。如果不做标准化，PCA 会被那个大尺度的基因完全带偏。所以跑 PCA 前，一定要先对数据做 Z-score 标准化（每个基因减去均值，再除以标准差）。

4.3 实操：用Python做PCA降维

别嫌啰嗦，我带你一步步写。先准备一个基因表达矩阵假数据：

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# 假装我们有100个样本，2000个基因
np.random.seed(42)
X = np.random.randn(100, 2000)
# 人为加入两组差异：前50个样本模拟处理组，后50个模拟对照组
X[:50, :100] += 1.5
X[50:, 100:200] += 1.0

# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCA降到2维
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# 画图
plt.figure(figsize=(6,5))
colors = ['red']*50 + ['blue']*50
plt.scatter(X_pca[:,0], X_pca[:,1], c=colors, alpha=0.7, edgecolors='k')
plt.xlabel('PC1 (explained variance ratio: {:.2%})'.format(pca.explained_variance_ratio_[0]))
plt.ylabel('PC2 (explained variance ratio: {:.2%})'.format(pca.explained_variance_ratio_[1]))
plt.title('PCA of gene expression data')
plt.grid(alpha=0.3)
plt.show()

运行这段代码，你应该能看到两组样本在 PC1 方向上分开了。这就是 PCA 的功劳——把两千维数据里的主要差异提取出来，浓缩在一张二维图上。

4.4 解释方差——选几个主成分才够？

你可能会问：「我选2个主成分够吗？3个呢？」这个问题靠「解释方差比」来回答。

每个主成分都能解释原始数据总方差的一部分。通常我们画一个「碎石图」：

pca_full = PCA()
pca_full.fit(X_scaled)
explained_variance = pca_full.explained_variance_ratio_

plt.figure(figsize=(6,4))
plt.plot(range(1, len(explained_variance)+1), explained_variance, 'o-', color='#3498db')
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Scree Plot')
plt.grid(alpha=0.3)
plt.show()

一般情况下，前2-3个主成分能解释总方差的30%-60%就算不错了。剩下的方差通常属于噪声或者单个样本的个性。

比如上表，前三个主成分解释了45%的方差，够用了。我个人习惯是：如果前两个主成分累计超过20%，就先用二维图看看；如果效果不好，再考虑三维或者用UMAP补救。

4.5 三维PCA投影——用mpl_toolkits画出来

有时候二维图分不开，三维图就能分开了。我演示一下三维PCA的绘制方法：

from mpl_toolkits.mplot3d import Axes3D

pca3 = PCA(n_components=3)
X_pca3 = pca3.fit_transform(X_scaled)

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_pca3[:,0], X_pca3[:,1], X_pca3[:,2], 
           c=['red']*50 + ['blue']*50, alpha=0.7, edgecolors='k')
ax.set_xlabel('PC1 ({:.1%})'.format(pca3.explained_variance_ratio_[0]))
ax.set_ylabel('PC2 ({:.1%})'.format(pca3.explained_variance_ratio_[1]))
ax.set_zlabel('PC3 ({:.1%})'.format(pca3.explained_variance_ratio_[2]))
plt.title('3D PCA projection')
plt.show()

三维图看着炫，但静态图像里旋转角度比较死板。我建议你在 jupyter notebook 里交互式旋转，直观感受数据在三维空间里的形状。

4.6 PCA与批次效应——第一章就提过的检测方法

批次效应是单细胞、RNA-seq 分析里的老问题。PCA 的第一主成分往往就是批次。我教你一个快速检测方法：把样本的批次信息作为颜色画在 PCA 图上，如果 PC1 分开的就是不同批次，那就说明批次效应很强，需要做矫正。

不过，PCA 只能「看到」线性差异。有些复杂的批次效应（比如某些基因只在特定批次出现非线性模式），PCA 可能看不出来。这时候你可以混合使用 UMAP 或 Harmony 等工具。

4.7 PCA的局限性——别以为它万能

• 线性假设：数据内在结构如果是环形、弧形，PCA 可能投影得乱七八糟。
• 对异常值敏感：一个离群点就能把第一主成分方向拉歪。
• 解释性困难：PC1 是几千个基因的线性组合，你很难说「这个主成分对应哪个生物学通路」。
• 无法处理缺失值：PCA 要求输入矩阵完整，有缺失值就得先插补。

我遇到过最坑的一件事：一个同事跑PCA发现PC1完美区分了男性和女性样本，他还以为发现了性别特异表达模式，后来发现是因为样本的性别比例和批次的采集时间完全重合——实际上是批次效应。所以 PCA 给出的「差异方向」不一定是生物学原因。

4.8 本章小结——一张图说清楚PCA降维的流程

下面这个SVG图总结了PCA从高维数据到二维/三维投影的完整流程，你可以对照着看。

你可以把这张图当作一个速查路线图。照着这个流程，从标准化开始，到画图为止，一步步走，基本不会出错。

4.9 我的个人实践建议

• 先用PCA做快速探索：我每次拿到新数据，第一件事就是跑PCA。如果前两个主成分就能解释20%以上方差，而且能看到明显分组，那我就对数据质量有底了。
• 别只盯着PC1和PC2：PC3有时候也会很有信息。我遇到一个数据集，PC1和PC2都是批次，PC3才是处理效应。所以我会至少看前三张图。
• 记录解释方差比：把那个比值写在图上。否则别人看图不知道你抓住了多少信息。
• 与UMAP配合：PCA给全局线性结构，UMAP给局部非线性结构。两者互补，我通常都会画。

好，这一章就到这里。你手上有数据的话，不妨马上打开 jupyter 跑一下 PCA 试试。代码我已经贴在上面了，直接复制就能用。