PCA主成分分析实战和可视化 | 附R代码和测试数据

摘要原文证据摘要文章给出从表达矩阵读入、可选 MAD 过滤、prcomp() PCA，到碎石图、样本图、分组置信椭圆、变量贡献图和 biplot 的一套完整示例；同时讨论中心化与 scale 的选择。

页面核验信息#

标题： PCA主成分分析实战和可视化 | 附R代码和测试数据
公众号/作者： 生信宝典 / 陈同
发布时间： 2018-12-17 08:20
原文链接： https://mp.weixin.qq.com/s/4R14xJkQVPtaufaoXOcPIw
浏览器核验结果： 页面包含“数据初始化”“数据读入”“数据预处理（可选）”“PCA分析”“PCA结果展示”和“PCA分析注意事项”等章节。

原文流程摘录#

数据预处理与 PCA#

原文表达矩阵约定为“行是基因、列是样本”，先去除全零行和零方差行，再转置为“行是样本、列是变量”：

data <- read.table(exprData, header = TRUE, row.names = NULL, sep = "\t")
rownames_data <- make.names(data[, 1], unique = TRUE)
data <- data[, -1, drop = FALSE]
rownames(data) <- rownames_data
data <- data[rowSums(data) > 0, ]
data <- data[apply(data, 1, var) != 0, ]

mads <- apply(data, 1, mad)
data <- data[rev(order(mads)), ]

data_t <- t(data)
pca <- prcomp(data_t[, 1:variableL], scale = TRUE)

原文还提示：PCA 前通常需要中心化；是否 scale 要结合变量单位、量级和噪声判断。对于 RNA-seq，建议先使用合适的库深度校正和 log 变换/方差稳定变换，再进行 PCA。

结果展示#

原文使用 factoextra 生成以下结果：

fviz_eig(pca, addlabels = TRUE)
fviz_pca_ind(pca, repel = TRUE)
fviz_pca_ind(
  pca,
  col.ind = data_t$conditions,
  mean.point = FALSE,
  addEllipses = TRUE,
  legend.title = "Groups",
  ellipse.type = "confidence",
  ellipse.level = 0.95
)
fviz_pca_var(pca, select.var = list(cos2 = 0.99), repel = TRUE,
             col.var = "cos2", geom.var = c("arrow", "text"))
fviz_pca_biplot(
  pca,
  fill.ind = data_t$conditions,
  addEllipses = TRUE,
  ellipse.type = "confidence",
  ellipse.level = 0.95,
  mean.point = FALSE,
  col.var = "contrib",
  select.var = list(contrib = 10),
  ggtheme = theme_minimal()
)

关键解释#

碎石图用于查看每个主成分的方差贡献率。
样本图的点表示样本，点间距离表示在选定 PCA 空间中的相对距离。
分组椭圆是组内分布的可视化辅助，不等于组间显著性检验。
cos2 表示变量在所选坐标轴上的表示质量；contrib 表示变量对主成分的贡献。
scale = TRUE 会让不同量级变量的权重趋于可比，但也可能放大噪声；参数必须写入分析记录。

复用时需要修正的地方#

原文示例含有公众号排版导致的换行和变量拼接，不能直接当作可运行脚本；应重新整理输入、样本信息合并和输出路径。
variableL 必须在转置后明确计算，不能依赖未定义变量。
分组列应来自样本信息表，并按表达矩阵样本名重新排序，不能假设两个文件顺序天然一致。
正式分析建议把椭圆阈值写成显式参数；本仓库脚本默认 ellipse.level = 0.95，不把椭圆当作统计检验。

与综合方法页的关系#

本页保留原文方法证据；参数化的完整 Bulk RNA-seq PCA 流程见远程仓库 04_Bulk转录组/4.2_PCA/02_prcomp_factoextra_Bulk_PCA完整流程.R。