如何在ggplot2绘制的PCA双标图中添加分组符号？

需求背景 #

拥有70个池塘的蝾螈存在-缺失（1=存在，0=缺失）数据及13种重金属浓度数据，目标是通过PCA双标图，用不同符号区分蝾螈存在/缺失的池塘，分析重金属与蝾螈分布的关联。此前使用vegan包生成了PCA双标图，但未实现分组符号标注，尝试ggplot2也未解决。

解决方案 #

方法1：基于vegan自定义绘图 #

利用vegan的rda模型提取样本得分和重金属载荷，手动绘制双标图并添加分组符号：

# 加载包
library(vegan)

# 提取重金属数据（沿用你的代码）
heavymetals <- cbind(Newts[,24:38], Newts[,40:46])
# 提取蝾螈存在/缺失数据（请根据实际列名调整，示例假设列名为newt_presence）
newt_status <- Newts$newt_presence

# 运行PCA（无约束rda等价于PCA）
heavymetals_model <- rda(heavymetals, scale = TRUE)

# 提取样本得分和重金属载荷
sample_scores <- scores(heavymetals_model, display = "sites")
metal_loadings <- scores(heavymetals_model, display = "species")

# 初始化绘图
plot(sample_scores, type = "n", main = "PCA双标图（蝾螈存在/缺失分组）")
# 添加重金属箭头
arrows(0, 0, metal_loadings[,1], metal_loadings[,2], col = "darkred", lwd = 1.5)
text(metal_loadings, rownames(metal_loadings), col = "darkred", pos = 3)
# 添加样本点，用不同形状区分蝾螈状态
points(sample_scores[newt_status == 1,], pch = 19, col = "blue", cex = 1.2) # 存在用实心圆
points(sample_scores[newt_status == 0,], pch = 24, col = "orange", cex = 1.2) # 缺失用空心三角
# 添加图例
legend("topright", legend = c("蝾螈存在", "蝾螈缺失"), 
       pch = c(19,24), col = c("blue","orange"), bty = "n")

方法2：基于ggplot2绘制灵活双标图 #

ggplot2更便于自定义样式，步骤如下：

library(ggplot2)
library(vegan)

# 提取数据（同前）
heavymetals <- cbind(Newts[,24:38], Newts[,40:46])
# 转因子便于映射，根据实际列名调整
newt_status <- factor(Newts$newt_presence, labels = c("蝾螈缺失", "蝾螈存在"))
heavymetals_model <- rda(heavymetals, scale = TRUE)

# 整理样本得分数据框
sample_df <- as.data.frame(scores(heavymetals_model, display = "sites"))
sample_df$newt_status <- newt_status

# 整理重金属载荷数据框
metal_df <- as.data.frame(scores(heavymetals_model, display = "species"))
metal_df$metal <- rownames(metal_df)

# 绘制PCA双标图
ggplot() +
  # 添加重金属箭头
  geom_segment(data = metal_df, aes(x = 0, y = 0, xend = PC1, yend = PC2),
               color = "darkred", arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = metal_df, aes(x = PC1, y = PC2, label = metal),
            color = "darkred", hjust = 0.5, vjust = -0.5) +
  # 添加样本点，形状映射到蝾螈状态
  geom_point(data = sample_df, aes(x = PC1, y = PC2, shape = newt_status, color = newt_status),
             size = 3) +
  # 自定义坐标轴标签（包含方差解释率）
  labs(title = "PCA双标图：蝾螈存在/缺失与重金属关联",
       x = paste0("PC1 (", round(heavymetals_model$CA$eig[1]/sum(heavymetals_model$CA$eig)*100, 1), "%)"),
       y = paste0("PC2 (", round(heavymetals_model$CA$eig[2]/sum(heavymetals_model$CA$eig)*100, 1), "%)")) +
  # 自定义形状和颜色
  scale_shape_manual(values = c(24, 19)) +
  scale_color_manual(values = c("orange", "blue")) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

关键说明 #

• 请根据Newts数据框的实际列名，调整蝾螈存在/缺失数据的提取代码（比如列是第23列，就写Newts[,23]）。
• 两种方法均可实现分组符号标注，ggplot2更适合后续调整样式（如颜色、形状、图例位置等）。

内容的提问来源于stack exchange，提问作者Seppe