R：PCoA（第四版）

rm (list = ls ()) #清除所有变量
# 1. 导入所需的库。
library(vegan)
library(tidyverse)
library(ggalt)
library(car)
library(ggforce)
library(ggpubr)
library(patchwork)

# 2. 定义所需的函数。
pairwise.adonis1 <- function(x, factors, p.adjust.m) { 
  x = as.matrix(x) 
  co = as.matrix(combn(unique(factors), 2))
  pairs <- F.Model <- R2 <- p.value <- c()
  for (elem in 1:ncol(co)) { 
    ad = adonis(x[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem])), 
                  factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))] ~ 
                  factors[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))], permutations = 999) 
    pairs <- c(pairs, paste(co[1, elem], 'vs', co[2, elem])) 
    F.Model <- c(F.Model, ad$aov.tab[1, 4]) 
    R2 <- c(R2, ad$aov.tab[1, 5]) 
    p.value <- c(p.value, ad$aov.tab[1, 6]) 
  }
  
  # p值调整
  p.adjusted = p.adjust(p.value, method = p.adjust.m) #使用了p.adjust.m参数指定的校正方法对p.value向量中的p值进行了校正，然后将校正后的p值赋值给了p.adjusted变量。
  pairw.res = data.frame(pairs, F.Model, R2, p.value, p.adjusted) #创建了一个数据框pairw.res，包含了每一对因子的Adonis分析结果。
  #这个数据框有五列：pairs（因子对）、F.Model（F值）、R2（R²值）、p.value（原始p值）和p.adjusted（校正后的p值）。
  return(pairw.res)#返回Adonis分析的结果
}

# 3. 读取和处理数据。
setwd("C:\\Users\\Administrator\\Desktop\\新建文件夹\\PCoA") 
otu <- read.table("./otu_table_R.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)]

# 4. 进行adonis分析并计算统计值。
bray_curtis <- vegan::vegdist(t(otu), method = "bray", na.rm = TRUE)
ado <- adonis(bray_curtis ~ map1$Group, permutations = 999, method = "bray") #进行了Adonis分析。
R2_value <- round(as.data.frame(ado$aov.tab[5])[1, 1], 3)
p_v_value <- as.data.frame(ado$aov.tab[6])[1, 1]
title <- paste("Adonis：R^2 = ", R2_value, "  P_Value = ", p_v_value, sep = "")

# 5. 绘制PCoA图。
pcoa <- cmdscale(bray_curtis, k = 2, eig = TRUE) 
points <- as.data.frame(pcoa$points) %>% dplyr::rename(x = "V1", y = "V2") 
eig <- pcoa$eig
points <- cbind(points, map1[match(rownames(points), map1$ID),])
n <- 0.85
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)) + #创建了一个ggplot对象，设置了数据框points作为数据，设置了x和y作为x轴和y轴的变量，设置了Group作为填充色和形状的变量。
  geom_point(alpha = .7, size = 5) + #添加了散点图层。alpha = .7设置了点的透明度，size = 5设置了点的大小。
  scale_shape_manual(values = shapes) +  # 设置了形状和填充色的映射关系。使用了之前定义的shapes和colors变量。
  scale_fill_manual(values = colors) +
  labs(x = paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits = 4), "%)", sep = ""),
       y = paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits = 4), "%)", sep = ""), title = title) +
  #设置了x轴、y轴和标题的标签。使用了paste函数将多个字符串连接在一起，使用了format函数将数字格式化为指定的小数位数。
  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.x = element_text(size = 14), # 自定义横坐标标签的大小
        axis.title.y = element_text(size = 14), # 自定义纵坐标标签的大小
        plot.title = element_text(size = 16, face = "bold"), # 自定义标题文字的大小和粗细
        legend.key.size = unit(1, "cm"), # 调整图例大小
        legend.title = element_blank(), # 去掉图例标题
        legend.text = element_text(size = 14), # 自定义图例标签文字大小
        panel.border = element_rect(colour = "black", fill = NA, size = 2), # 自定义边框的粗细
        axis.ticks = element_line(size = 2)) # 自定义坐标轴刻度的大小
p1

# 6. 输出pairwise adonis结果。
pair_bray_adonis <- pairwise.adonis1(bray_curtis, map1$Group, p.adjust.m = "bonferroni") #进行了成对的Adonis分析

# 存储为文本文件
write.table(as.data.frame(pair_bray_adonis), "table.txt", sep = "\t", quote = FALSE, row.names = FALSE) #将Adonis分析的结果保存为了一个文本文件。

tab2 <- ggtexttable(pair_bray_adonis, rows = NULL)
p2 <- tab2
#创建了一个表格，并将其赋值给了p2变量。ggtexttable函数是ggpubr包中的一个函数，用于在ggplot图中添加表格。
p2

# 使用ggsave保存PCoA图为PNG格式
ggsave(filename = "PCoA_plot.png", plot = p1, width = 11.5, height = 10, units = "in", dpi = 600)