📚 课程目录第4章 PCA在R和Python中的标准实现与常用函数

4. PCA 在 R 和 Python 中的标准实现与常用函数归纳总结

你也许会问：PCA 的理论我都懂了，那在代码里到底怎么用？
说实话，我一开始也卡在这步。书上讲的都是矩阵分解，但打开 R 或 Python，函数一堆，参数一堆。
今天我就把 R 和 Python 里最常用的 PCA 实现掰开揉碎讲清楚。后面做批次效应检测时，你会反复用到这些函数。

4.1 R 语言中的 PCA 实现

R 里实现 PCA 的主流函数有两个：prcomp() 和 princomp()。
我个人习惯用 prcomp()，因为它基于 SVD（奇异值分解），数值稳定性更好。
在我处理过的那批基因表达数据里（样本 200，基因 20000），prcomp() 从没出过报错，而 princomp() 却因为协方差矩阵奇异崩过一次。

避坑指南： 我曾经发现 princomp() 默认使用协方差矩阵分解，当变量数 > 样本数时直接报错。而 prcomp() 用 SVD 就没事。所以，大多数场景请优先选 prcomp()。

4.1.1 prcomp() 核心用法

# 标准 PCA 调用
pca_result <- prcomp(data_matrix, center = TRUE, scale. = TRUE)

参数很直观：

center：是否中心化，默认 TRUE，强烈建议保持。
scale.：是否标准化为方差 1，默认 FALSE。但如果你的变量量纲不统一（比如有的浓度、有的 pH），一定设成 TRUE。我在做代谢组数据时吃过这个亏，没标准化，结果 PCA 第一主成分被浓度最大的变量拽走了。

结果对象里常用的东西：

属性	说明	维度
`pca_result$x`	主成分得分（样本在新坐标下的坐标）	n × k
`pca_result$rotation`	载荷矩阵（变量对每个主成分的贡献）	p × k
`pca_result$sdev`	每个主成分的标准差	长度 k
`summary(pca_result)`	方差解释比例一眼看	—

4.1.2 提取方差解释比例的小技巧

# 直接算
var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)
round(var_explained[1:5], 3)   # 前 5 个主成分的解释方差

# 或者用 summary
summary(pca_result)$importance[, 1:5]

记得我最开始学 PCA 时，老以为 sdev 就是方差比例，其实它是标准差！平方后再归一化才是正解。很多新人踩过这坑。

4.2 Python 中的 PCA 实现

Python 这边最标准的是 sklearn.decomposition.PCA。它底层用的也是 SVD，稳定可靠。
在我做的单细胞 RNA-seq 批次效应项目中，几万个细胞用 PCA(n_components=50) 几秒钟就算完了，速度挺满意。

4.2.1 基础用法

from sklearn.decomposition import PCA
import numpy as np

# 假设 data 是 numpy 数组，行为样本，列为变量
pca = PCA(n_components=50, whiten=False)
pca.fit(data)

# 或者一步到位
X_pca = pca.fit_transform(data)

核心参数：

n_components：保留的主成分个数。若设为 0~1 之间的小数，表示保留方差比例（例如 0.95 保留 95% 方差），自动选 k。我常用这个，省心。
whiten：是否白化（使各主成分方差为 1）。一般批次效应检测不用，留到后续分析再考虑。

4.2.2 结果属性

属性	说明
`pca.components_`	载荷矩阵，形状 (k, p)，每行对应一个主成分
`pca.explained_variance_`	每个主成分解释的方差值（绝对值）
`pca.explained_variance_ratio_`	方差解释比例（0~1），直接拿来画 scree plot
`pca.singular_values_`	奇异值，SVD 直接产物

# 快速查看前 5 个主成分的方差比例
print(pca.explained_variance_ratio_[:5])

  个人经验： 我习惯在 fit_transform 之前先用 StandardScaler 标准化数据，保证所有变量均值 0 方差 1。尤其当变量单位不同时，标准化是必须的。我曾经偷懒没做，直接跑 PCA，结果第一主成分全是表达量最高的基因，批次效应完全被淹没。

4.3 R vs Python：PCA 函数对照速查表

为了方便你切换语言，我整理了一张对照表。

功能/属性	R（prcomp）	Python（sklearn）
核心函数	`prcomp(x, center, scale.)`	`PCA(n_components)`
标准化参数	scale. = TRUE	需额外用 `StandardScaler`
得分矩阵	`$x`	`.fit_transform(data)`
载荷矩阵	`$rotation`	`.components_`
方差比例	`summary()$importance[2,]`	`.explained_variance_ratio_`
特征值（标准差平方）	`$sdev^2`	`.explained_variance_`
自动选出 k（解释 >= 95%）	手动算累加	`n_components=0.95`

4.4 一个完整的 PCA 流程（以 Python 为例）

下面这段代码是我项目里的标准模板，直接复制改个数据就能用。

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

# 1. 加载数据（行为样本，列为变量）
df = pd.read_csv("expression_data.csv", index_col=0)

# 2. 标准化
scaler = StandardScaler()
data_scaled = scaler.fit_transform(df)

# 3. PCA 保留 95% 方差
pca = PCA(n_components=0.95)
X_pca = pca.fit_transform(data_scaled)

# 4. 检查结果
print(f"保留的主成分数量: {pca.n_components_}")
print(f"累计方差比例: {sum(pca.explained_variance_ratio_):.3f}")

# 5. 转换为 DataFrame 方便画图
pca_df = pd.DataFrame(
    X_pca,
    columns=[f"PC{i+1}" for i in range(X_pca.shape[1])],
    index=df.index
)
pca_df["batch"] = df["batch"]   # 假设原始数据里有一列 batch 信息

你可以直接把 pca_df 传给 Seaborn 画散点图或热图。我在识别批次效应时，常用前两个主成分着色看是否有 batch 相关聚类。

4.5 PCA 工作流与函数选择逻辑

下面这张图总结了从原始数据到 PCA 结果的核心步骤，以及在不同环节该用什么函数。

4.6 避坑总结

标准化不可省略 —— 我见过太多人直接拿原始 counts 跑 PCA，第一 PC 永远是表达量最高的基因，毫无意义。量纲不同的变量，必须先标准化。
R 里 prcomp 的 retx=TRUE 才是得分矩阵 —— 默认就是 TRUE，但你如果手贱关了，$x 会是空的。我以前调试了半天才发现。
Python 中 n_components 设为小数时 —— 内部会先计算全部 SVD（昂贵），再截断。如果变量数很多（>10000），建议显式指定整数，比如 50，避免计算全量 SVD。
注意 prcomp() 默认将变量放在列 —— 确保你的数据矩阵是样本×变量，这和大多数搞生物信息的人习惯一致。但有些旧的 R 包会要求转置，检查清楚。

严重警告： PCA 结果用于批次效应识别时，千万不要只看 PC1 vs PC2 的散点图。批次效应往往隐藏在后面的主成分里。我建议至少看前 10 个 PC 的热图，或者用我下一段话里提到的“方差解释比例累积曲线”来辅助判断。不过那已经是下一节的内容了，这里先记着。

好了，R 和 Python 的 PCA 核心函数就这些。说白了，你只需要记住两个函数：prcomp() 和 PCA，配合标准化步骤，就足以应对 90% 的场景。下次当你拿到一批新数据，心里没底时，先跑个 PCA 热图看看——这活儿我们已经铺垫好了。