R语言PCA可视化进阶:用ggplot2给不同分组加上置信椭圆(附完整代码)
在生物信息学和多元统计分析中,主成分分析(PCA)是最常用的降维技术之一。但仅仅展示散点图往往难以直观呈现不同组别间的分布差异与重叠程度。本文将深入探讨如何利用ggplot2的stat_ellipse()函数,为PCA结果添加具有统计意义的置信椭圆,让数据可视化兼具科学严谨性与美学表现力。
1. PCA基础与数据准备
1.1 理解PCA结果的数据结构
典型的PCA分析会生成两个核心结果:特征向量(旋转矩阵)和主成分得分(样本在新坐标系的投影)。以经典的iris数据集为例:
# 加载数据并计算PCA data(iris) pca_result <- prcomp(iris[,1:4], scale. = TRUE) # 查看PCA结果结构 str(pca_result)输出结果包含:
$rotation: 各原始变量对主成分的贡献(载荷)$x: 样本在主成分空间的坐标$sdev: 各主成分的标准差
1.2 数据框的构建技巧
为便于ggplot2绘图,我们需要将PCA得分与原始分组信息合并:
library(tibble) pca_df <- as_tibble(pca_result$x[,1:2]) %>% mutate(Species = iris$Species, PC1 = PC1 * -1, # 可选:调整坐标轴方向 PC2 = PC2 * -1)提示:主成分的符号是任意的,有时需要根据生物学意义手动调整方向
2. 置信椭圆的统计原理
2.1 椭圆背后的数学基础
置信椭圆是基于多元正态分布的假设,表示给定置信水平下数据点的可能分布范围。关键参数包括:
| 参数 | 说明 | 典型值 |
|---|---|---|
| level | 置信水平 | 0.95 |
| type | 分布假设 | "norm"或"t" |
| segments | 椭圆平滑度 | 51 |
2.2 分布类型的选择
stat_ellipse()支持两种分布假设:
- norm:多元正态分布(默认)
- t:多元t分布(更保守,适合小样本)
# 比较两种分布假设 ggplot(pca_df, aes(PC1, PC2, color = Species)) + stat_ellipse(type = "norm", linetype = 2) + stat_ellipse(type = "t", linetype = 1) + geom_point()3. ggplot2高级可视化技巧
3.1 基础椭圆绘制
最简单的置信椭圆只需一行代码:
library(ggplot2) ggplot(pca_df, aes(PC1, PC2, color = Species)) + stat_ellipse() + geom_point()3.2 样式深度定制
通过调整geom参数,可以实现多种视觉效果:
# 填充式椭圆 ggplot(pca_df, aes(PC1, PC2, color = Species)) + stat_ellipse(aes(fill = Species), geom = "polygon", alpha = 0.2, show.legend = FALSE) + geom_point(size = 3) + scale_fill_brewer(palette = "Set1") + theme_minimal()关键参数组合:
geom = "path":仅显示轮廓geom = "polygon":填充区域alpha:透明度控制level:调整置信区间大小
3.3 多图层叠加策略
为增强可视化效果,可以组合多种几何对象:
ggplot(pca_df, aes(PC1, PC2)) + stat_ellipse(aes(fill = Species), alpha = 0.3) + geom_point(aes(color = Species), size = 2.5) + geom_text(aes(label = rownames(iris)), check_overlap = TRUE, size = 3, vjust = -1) + labs(title = "PCA of Iris Dataset with 95% Confidence Ellipses", x = paste0("PC1 (", round(summary(pca_result)$importance[2,1]*100, 1), "%)"), y = paste0("PC2 (", round(summary(pca_result)$importance[2,2]*100, 1), "%)")) + theme_bw(base_size = 12)4. 实战案例:微生物组数据分析
4.1 处理真实科研数据
以微生物组学常用的OTU表为例:
# 模拟微生物组数据 set.seed(123) microbiome_data <- matrix(rpois(100*50, 5), nrow=100) rownames(microbiome_data) <- paste0("Sample", 1:100) colnames(microbiome_data) <- paste0("OTU", 1:50) groups <- rep(c("Control", "Treatment1", "Treatment2"), c(30, 40, 30)) # Hellinger转换 microbiome_hel <- sqrt(microbiome_data / rowSums(microbiome_data)) # 执行PCA pca_micro <- prcomp(microbiome_hel)4.2 复杂分组可视化
当存在多个分组因素时,可以使用分面(facet)展示:
# 添加实验批次信息 batch <- rep(c("Batch1", "Batch2"), each=50) pca_micro_df <- data.frame( PC1 = pca_micro$x[,1], PC2 = pca_micro$x[,2], Group = groups, Batch = batch ) ggplot(pca_micro_df, aes(PC1, PC2, color = Group)) + stat_ellipse(aes(linetype = Batch), level = 0.9) + geom_point(size = 2) + scale_linetype_manual(values = c(2, 3)) + labs(title = "Microbiome PCA with Group and Batch Effects")4.3 发表级图形优化
最终调整图形细节以满足期刊要求:
final_plot <- ggplot(pca_df, aes(PC1, PC2, color = Species)) + stat_ellipse(aes(fill = Species), geom = "polygon", alpha = 0.15, size = 0.8, level = 0.95) + geom_point(size = 3, alpha = 0.8) + scale_color_manual(values = c("#E69F00", "#56B4E9", "#009E73")) + scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73")) + labs(x = "Principal Component 1 (92.5%)", y = "Principal Component 2 (5.3%)") + theme_classic(base_size = 14) + theme(legend.position = "right", panel.grid.major = element_line(color = "grey90"), axis.line = element_line(color = "black"), plot.title = element_text(hjust = 0.5)) # 保存高分辨率图片 ggsave("PCA_with_ellipses.tiff", plot = final_plot, device = "tiff", dpi = 600, width = 8, height = 6)