R：PCA（第三版）

# 清除所有变量
rm(list = ls()) 
# 设置工作目录
setwd("C:\\Users\\Administrator\\Desktop\\新建文件夹\\PCA_Pathway") 

# 1. 加载所需的库
library(vegan)
library(tidyverse)
library(ggpubr)
library(patchwork)
library(ggforce)


# 2. 读取和处理数据
otu <- read.table("./1.txt", row.names = 1, sep = "\t", header = TRUE) %>% as.data.frame()
map <- read.table("./group.txt", sep = "\t", header = TRUE) 
colnames(map)[1] <- "ID" 
row.names(map) <- map$ID 
idx <- rownames(map) %in% colnames(otu) 
map1 <- map[idx,]
otu <- otu[, rownames(map1)]

# 3. 绘制PCA图
otu_centered <- scale(t(otu), scale = TRUE)
pca <- prcomp(otu_centered)
summary_pca <- summary(pca)
points <- as.data.frame(pca$x) %>% dplyr::rename(x = "PC1", y = "PC2")
points <- cbind(points, map1[match(rownames(points), map1$ID),]) 
colors <- c("B73_DAS28"="#8FC9E2","B73_DAS42"="#8FC9E2","B73_DAS56"="#8FC9E2","B73_DAS70"="#8FC9E2","Mo17_DAS28"="#ECC97F","Mo17_DAS42"="#ECC97F","Mo17_DAS56"="#ECC97F","Mo17_DAS70"="#ECC97F")
shapes <- c("B73_DAS28"=24, "B73_DAS42"=22, "B73_DAS56"=21, "B73_DAS70"=23, "Mo17_DAS28"=24, "Mo17_DAS42"=22, "Mo17_DAS56"=21, "Mo17_DAS70"=23)
levels_order <- c("B73_DAS28", "B73_DAS42", "B73_DAS56", "B73_DAS70", "Mo17_DAS28", "Mo17_DAS42", "Mo17_DAS56", "Mo17_DAS70")
points$Group <- factor(points$Group, levels = levels_order)

# 在ggplot中使用这些形状和颜色进行绘制
p1 <- ggplot(points, aes(x = x, y = y, fill = Group, shape = Group)) + 
  geom_point(alpha = .7, size = 6) + 
  scale_shape_manual(values = shapes) + 
  scale_fill_manual(values = colors) +
  labs(x = paste("PC1 (", format(summary_pca$importance[2, 1] * 100, digits = 4), "%)", sep = ""),
       y = paste("PC2 (", format(summary_pca$importance[2, 2] * 100, digits = 4), "%)", sep = "")) +
  geom_mark_ellipse(aes(fill = Group, label = Group), alpha = 0.1, color = "grey", linetype = 3) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text = element_text(color = "black", size = 12),
        axis.title = element_text(size = 16),
        legend.text = element_text(size = 14),
        legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 2),
        axis.ticks = element_line(linewidth = 2),
        legend.key.size = unit(1, "cm")) +
  coord_cartesian(xlim = c(-max(abs(points$x)) * 1.1, max(abs(points$x)) * 1.1), ylim = c(-max(abs(points$y)) * 1.1, max(abs(points$y)) * 1.1)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")

# 显示绘制的图
p1

# 保存PCA图
ggsave(filename = "PCA_plot.png", plot = p1, width = 10, height = 8, units = "in", dpi = 600)