plot34 · 代码 3 · 图 0 · 笔记 1 · 代码库专题
本目录整理了三篇 PCA 美化图教程的可复用版本，覆盖：
4.2 Bulk RNA-seq PCA 美化图（自动生成占位图）
图源：（自动生成占位图）
从专题 README 抽取；原文顺序保留，未跑过的部分照实标注。
本目录整理了三篇 PCA 美化图教程的可复用版本，覆盖：
prcomp() + ggplot2/factoextra 二维样本 PCA、分组椭圆和主题美化；ggbiplot 的分组椭圆、相关环和变量向量；pca3d 三维 PCA 的方法参考（当前目录不强制输出三维图）。top_n_mad 默认 2000，只是探索性筛选；应根据样本量和研究问题调整，并保存 selected_genes_MAD.tsv。center 控制中心化，scale_data 控制定标。RNA-seq 变换值是否 scale 取决于数据量级和分析目的，不能机械照搬教程。0.95 覆盖概率并在样本数不足时关闭椭圆。| 文件 | 用途 |
|---|---|
01_factoextra_PCA二维美化.R |
轻量二维图；输出样本 scores、方差解释率和 PNG/PDF |
02_prcomp_factoextra_Bulk_PCA完整流程.R |
完整 QC/解释流程；输出碎石图、样本图、变量贡献图和 biplot |
03_ggbiplot_Bulk_PCA美化.R |
ggbiplot 分组椭圆 biplot；可选变量向量 |
expression_transformed.tsv：
sample_id 对应；sample_metadata.tsv：
sample_id 列；condition 列，可在每个脚本顶部修改 group_col；示例表头：
sample_id condition batch
sample_01 control batch_1
sample_02 treated batch_1示例数据、.RData、.RDS、.h5ad 和生成结果不放入本仓库；请在本机 data/ 和 result/ 中准备。
top_n_mad。01 快速确认样本关系和文件格式。02 生成碎石图、样本 PCA、变量贡献和 biplot。ggbiplot 版本后，再运行 03 做风格化 biplot。常见依赖：ggplot2、factoextra；可选 ggbiplot。第一篇原文的三维部分还需要 pca3d。ggbiplot 通常需要从 GitHub 安装，安装方式和版本应以当前包状态为准。
本目录只保留可复用脚本和说明，不上传教程原始网页、PDF、压缩包或大数据对象。
3 个源码文件，按目录顺序；本机绝对路径已替换为 ~。 代码可复制，但需要按自己的数据改路径与参数。
# Bulk RNA-seq PCA 二维美化：prcomp + ggplot2/factoextra
# 输入：已变换的表达矩阵（行=基因，列=样本）和样本信息表
# 运行前修改 config 区域；不要把原始 count 或大数据对象提交到代码库。
rm(list = ls())
## ---------------------------- config ----------------------------
expr_file <- "data/expression_transformed.tsv"
metadata_file <- "data/sample_metadata.tsv"
output_dir <- "result/pca_factoextra"
group_col <- "condition"
sample_id_col <- "sample_id"
top_n_mad <- 2000L
center <- TRUE
scale_data <- FALSE
ellipse_level <- 0.95
## ---------------------------------------------------------------
required_packages <- c("ggplot2", "factoextra")
missing_packages <- required_packages[!vapply(required_packages,
requireNamespace,
logical(1), quietly = TRUE)]
if (length(missing_packages) > 0L) {
stop("请先安装 R 包：", paste(missing_packages, collapse = ", "))
}
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
expr <- read.delim(expr_file, check.names = FALSE, row.names = 1,
stringsAsFactors = FALSE)
meta <- read.delim(metadata_file, check.names = FALSE,
stringsAsFactors = FALSE)
if (!sample_id_col %in% colnames(meta)) {
stop("样本信息中缺少 sample_id_col：", sample_id_col)
}
if (!group_col %in% colnames(meta)) {
stop("样本信息中缺少 group_col：", group_col)
}
if (anyDuplicated(colnames(expr))) {
stop("表达矩阵样本名重复，请先修正")
}
if (anyDuplicated(meta[[sample_id_col]])) {
stop("样本信息 sample_id 重复，请先修正")
}
common_samples <- intersect(colnames(expr), meta[[sample_id_col]])
if (length(common_samples) < 3L) {
stop("表达矩阵与样本信息的共同样本少于 3 个")
}
expr <- expr[, common_samples, drop = FALSE]
meta <- meta[match(common_samples, meta[[sample_id_col]]), , drop = FALSE]
rownames(meta) <- meta[[sample_id_col]]
expr <- as.matrix(expr)
storage.mode(expr) <- "numeric"
if (anyNA(expr) || any(!is.finite(expr))) {
stop("表达矩阵包含 NA/Inf，请先处理")
}
expr <- expr[rowSums(expr) > 0, , drop = FALSE]
gene_mad <- apply(expr, 1, mad, na.rm = TRUE)
gene_mad <- gene_mad[is.finite(gene_mad) & gene_mad > 0]
if (length(gene_mad) < 2L) {
stop("有效且有变异的基因少于 2 个，无法进行 PCA")
}
top_n_mad <- min(top_n_mad, length(gene_mad))
selected_genes <- names(sort(gene_mad, decreasing = TRUE))[seq_len(top_n_mad)]
expr_selected <- expr[selected_genes, , drop = FALSE]
pca <- prcomp(t(expr_selected), center = center, scale. = scale_data)
variance <- pca$sdev^2 / sum(pca$sdev^2)
scores <- as.data.frame(pca$x[, seq_len(min(2L, ncol(pca$x))), drop = FALSE])
if (ncol(scores) < 2L) {
stop("PCA 只产生一个主成分，无法绘制二维图")
}
scores[[sample_id_col]] <- rownames(scores)
scores$group <- factor(meta[rownames(scores), group_col])
write.table(data.frame(component = paste0("PC", seq_along(variance)),
variance_explained = variance),
file.path(output_dir, "pca_variance.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
write.table(scores, file.path(output_dir, "pca_scores.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
write.table(data.frame(gene = selected_genes, MAD = gene_mad[selected_genes]),
file.path(output_dir, "selected_genes_MAD.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
plot_data <- scores
group_counts <- table(plot_data$group)
add_ellipse <- all(group_counts >= 3L) && nlevels(droplevels(plot_data$group)) >= 2L
p <- ggplot2::ggplot(plot_data, ggplot2::aes(PC1, PC2, color = group, fill = group)) +
ggplot2::geom_hline(yintercept = 0, color = "grey85", linewidth = 0.3) +
ggplot2::geom_vline(xintercept = 0, color = "grey85", linewidth = 0.3) +
ggplot2::geom_point(size = 3, alpha = 0.9) +
ggplot2::labs(
x = sprintf("PC1 (%.2f%%)", 100 * variance[1]),
y = sprintf("PC2 (%.2f%%)", 100 * variance[2]),
color = group_col,
fill = group_col
) +
ggplot2::theme_classic(base_size = 12) +
ggplot2::theme(legend.position = "right")
if (add_ellipse) {
p <- p + ggplot2::stat_ellipse(geom = "polygon", alpha = 0.12,
level = ellipse_level, linewidth = 0.4,
show.legend = FALSE)
}
ggplot2::ggsave(file.path(output_dir, "PCA_PC1_PC2.pdf"), p,
width = 6, height = 5, units = "in")
ggplot2::ggsave(file.path(output_dir, "PCA_PC1_PC2.png"), p,
width = 6, height = 5, units = "in", dpi = 300)
# Bulk RNA-seq PCA 完整流程：输入检查 -> MAD 筛选 -> PCA -> QC/解释图
# 表达矩阵必须是已经完成合适变换的数据（例如 logCPM 或 DESeq2 VST）。
rm(list = ls())
## ---------------------------- config ----------------------------
expr_file <- "data/expression_transformed.tsv"
metadata_file <- "data/sample_metadata.tsv"
output_dir <- "result/pca_full"
group_col <- "condition"
sample_id_col <- "sample_id"
top_n_mad <- 2000L
center <- TRUE
scale_data <- FALSE
top_n_variables <- 10L
## ---------------------------------------------------------------
required_packages <- c("ggplot2", "factoextra")
missing_packages <- required_packages[!vapply(required_packages,
requireNamespace,
logical(1), quietly = TRUE)]
if (length(missing_packages) > 0L) {
stop("请先安装 R 包：", paste(missing_packages, collapse = ", "))
}
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
expr <- read.delim(expr_file, check.names = FALSE, row.names = 1,
stringsAsFactors = FALSE)
meta <- read.delim(metadata_file, check.names = FALSE,
stringsAsFactors = FALSE)
stopifnot(sample_id_col %in% colnames(meta), group_col %in% colnames(meta))
stopifnot(!anyDuplicated(colnames(expr)), !anyDuplicated(meta[[sample_id_col]]))
common_samples <- intersect(colnames(expr), meta[[sample_id_col]])
if (length(common_samples) < 3L) stop("共同样本少于 3 个")
expr <- expr[, common_samples, drop = FALSE]
meta <- meta[match(common_samples, meta[[sample_id_col]]), , drop = FALSE]
rownames(meta) <- meta[[sample_id_col]]
expr <- as.matrix(expr)
storage.mode(expr) <- "numeric"
if (anyNA(expr) || any(!is.finite(expr))) stop("表达矩阵包含 NA/Inf")
expr <- expr[rowSums(expr) > 0, , drop = FALSE]
gene_mad <- apply(expr, 1, mad, na.rm = TRUE)
gene_mad <- gene_mad[is.finite(gene_mad) & gene_mad > 0]
if (length(gene_mad) < 2L) stop("有效且有变异的基因少于 2 个")
top_n_mad <- min(top_n_mad, length(gene_mad))
selected_genes <- names(sort(gene_mad, decreasing = TRUE))[seq_len(top_n_mad)]
pca <- prcomp(t(expr[selected_genes, , drop = FALSE]),
center = center, scale. = scale_data)
variance <- pca$sdev^2 / sum(pca$sdev^2)
scores <- as.data.frame(pca$x)
scores[[sample_id_col]] <- rownames(scores)
scores$group <- factor(meta[rownames(scores), group_col])
write.table(data.frame(component = paste0("PC", seq_along(variance)),
variance_explained = variance),
file.path(output_dir, "pca_variance.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
write.table(scores, file.path(output_dir, "pca_scores.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
write.table(data.frame(gene = selected_genes, MAD = gene_mad[selected_genes]),
file.path(output_dir, "selected_genes_MAD.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)
## 1. 碎石图
eig_plot <- factoextra::fviz_eig(pca, addlabels = TRUE, barfill = "#3C5488")
ggplot2::ggsave(file.path(output_dir, "01_scree_plot.pdf"), eig_plot,
width = 6, height = 4.5, units = "in")
## 2. 样本图：分组仅作探索性展示
sample_plot <- factoextra::fviz_pca_ind(
pca,
geom.ind = "point",
col.ind = scores$group,
mean.point = FALSE,
addEllipses = FALSE,
repel = TRUE,
legend.title = group_col,
title = ""
) +
ggplot2::labs(
x = sprintf("PC1 (%.2f%%)", 100 * variance[1]),
y = sprintf("PC2 (%.2f%%)", 100 * variance[2])
) +
ggplot2::theme_classic(base_size = 12)
ggplot2::ggsave(file.path(output_dir, "02_sample_PCA.pdf"), sample_plot,
width = 6, height = 5, units = "in")
## 3. 变量贡献和 cos2
variable_plot <- factoextra::fviz_pca_var(
pca,
col.var = "contrib",
gradient.cols = c("#2166AC", "#F7F7F7", "#B2182B"),
select.var = list(contrib = top_n_variables),
repel = TRUE
) + ggplot2::theme_classic(base_size = 12)
ggplot2::ggsave(file.path(output_dir, "03_top_variable_contribution.pdf"),
variable_plot, width = 6.5, height = 5.5, units = "in")
## 4. biplot：变量较多时只保留贡献最高的一小部分
biplot <- factoextra::fviz_pca_biplot(
pca,
geom.ind = "point",
col.ind = scores$group,
mean.point = FALSE,
addEllipses = FALSE,
col.var = "contrib",
gradient.cols = c("#2166AC", "#F7F7F7", "#B2182B"),
select.var = list(contrib = top_n_variables),
repel = TRUE,
ggtheme = ggplot2::theme_classic(base_size = 12)
)
ggplot2::ggsave(file.path(output_dir, "04_biplot.pdf"), biplot,
width = 7, height = 6, units = "in")
# Bulk RNA-seq PCA biplot：prcomp + ggbiplot
# ggbiplot 目前常需从 GitHub 安装；正式使用前请记录 R 包版本。
rm(list = ls())
## ---------------------------- config ----------------------------
expr_file <- "data/expression_transformed.tsv"
metadata_file <- "data/sample_metadata.tsv"
output_dir <- "result/pca_ggbiplot"
group_col <- "condition"
sample_id_col <- "sample_id"
top_n_mad <- 2000L
center <- TRUE
scale_data <- TRUE
ellipse_prob <- 0.95
show_variable_axes <- FALSE
## ---------------------------------------------------------------
required_packages <- c("ggplot2", "ggbiplot")
missing_packages <- required_packages[!vapply(required_packages,
requireNamespace,
logical(1), quietly = TRUE)]
if (length(missing_packages) > 0L) {
stop("请先安装 R 包：", paste(missing_packages, collapse = ", "))
}
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
expr <- read.delim(expr_file, check.names = FALSE, row.names = 1,
stringsAsFactors = FALSE)
meta <- read.delim(metadata_file, check.names = FALSE,
stringsAsFactors = FALSE)
stopifnot(sample_id_col %in% colnames(meta), group_col %in% colnames(meta))
stopifnot(!anyDuplicated(colnames(expr)), !anyDuplicated(meta[[sample_id_col]]))
common_samples <- intersect(colnames(expr), meta[[sample_id_col]])
if (length(common_samples) < 3L) stop("共同样本少于 3 个")
expr <- expr[, common_samples, drop = FALSE]
meta <- meta[match(common_samples, meta[[sample_id_col]]), , drop = FALSE]
rownames(meta) <- meta[[sample_id_col]]
expr <- as.matrix(expr)
storage.mode(expr) <- "numeric"
if (anyNA(expr) || any(!is.finite(expr))) stop("表达矩阵包含 NA/Inf")
expr <- expr[rowSums(expr) > 0, , drop = FALSE]
gene_mad <- apply(expr, 1, mad, na.rm = TRUE)
gene_mad <- gene_mad[is.finite(gene_mad) & gene_mad > 0]
if (length(gene_mad) < 2L) stop("有效且有变异的基因少于 2 个")
top_n_mad <- min(top_n_mad, length(gene_mad))
selected_genes <- names(sort(gene_mad, decreasing = TRUE))[seq_len(top_n_mad)]
pca <- stats::prcomp(t(expr[selected_genes, , drop = FALSE]),
center = center, scale. = scale_data)
variance <- pca$sdev^2 / sum(pca$sdev^2)
groups <- factor(meta[rownames(pca$x), group_col])
has_ellipse <- nlevels(droplevels(groups)) > 1L && all(table(groups) >= 3L)
plot <- ggbiplot::ggbiplot(
pca,
obs.scale = 1,
var.scale = 1,
groups = groups,
ellipse = has_ellipse,
ellipse.prob = ellipse_prob,
circle = FALSE,
var.axes = show_variable_axes,
varname.abbrev = TRUE
) +
ggplot2::labs(
x = sprintf("PC1 (%.2f%%)", 100 * variance[1]),
y = sprintf("PC2 (%.2f%%)", 100 * variance[2]),
color = group_col
) +
ggplot2::theme_classic(base_size = 12) +
ggplot2::theme(legend.position = "top")
ggplot2::ggsave(file.path(output_dir, "PCA_ggbiplot.pdf"), plot,
width = 7, height = 5.5, units = "in")
ggplot2::ggsave(file.path(output_dir, "PCA_ggbiplot.png"), plot,
width = 7, height = 5.5, units = "in", dpi = 300)
write.table(data.frame(component = paste0("PC", seq_along(variance)),
variance_explained = variance),
file.path(output_dir, "pca_variance.tsv"), sep = "\t",
quote = FALSE, row.names = FALSE)