rm (list = ls ()) #清楚所有变量
# 1. 导入所需的库。
library(vegan) #提供了进行社区生态学分析的函数，包括多元分析、物种多样性分析等。
library(tidyverse) #一组用于数据科学的R包，包括ggplot2、dplyr、tidyr、readr、purrr和tibble等。
library(ggalt) #增强了ggplot2的功能，提供了一些额外的坐标系统、几何对象、统计变换等。
library(car) #提供了许多有用的函数，用于线性回归分析。
library(ggforce) #是ggplot2的扩展，提供了一些新的几何对象、位置调整方法、坐标系统等。
library(ggpubr) #提供了一些函数，使得在ggplot2中创建出版质量的图形变得更加容易。
library(patchwork) #允许将多个ggplot对象组合成一个图形。
# 2. 定义所需的函数。
#adonis函数是vegan包中的一个函数，用于进行多元方差分析（Permutational Multivariate Analysis of Variance，PERMANOVA）。
#这种分析可以用来测试多个群组之间的差异是否显著。
pairwise.adonis1 <- function(x, factors, p.adjust.m) { #定义了一个名为pairwise.adonis1的函数，该函数接受三个参数：x，factors和p.adjust.m
x = as.matrix(x) #将输入的数据x转换为矩阵
#创建了一个矩阵co，生成了一个包含所有可能的两两组合的矩阵
co = as.matrix(combn(unique(factors), 2))
# 初始化了四个空的向量：pairs，F.Model，R2和p.value。
# pairs：这个向量将存储每一对因子的名称
# F.Model：这个向量将存储每一对因子的Adonis分析的F值。F值是用来衡量模型的拟合优度的一个统计量。
# R2：这个向量将存储每一对因子的Adonis分析的R²值。R²值也是一个衡量模型拟合优度的统计量，表示模型可以解释的数据变异的比例。
# p.value：这个向量将存储每一对因子的Adonis分析的p值。p值是用来判断结果是否显著的一个统计量。
pairs <- F.Model <- R2 <- p.value <- c()
# 对factors中的每一对执行adonis分析
for (elem in 1:ncol(co)) { #开始了一个循环，循环的次数等于co矩阵的列数，循环的目的是对factors中的每一对因子进行Adonis分析
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[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))], permutations = 999) #这是你设置的置换次数，用于计算p值。置换次数越多，计算的p值越精确，但计算时间也会越长。
pairs <- c(pairs, paste(co[1, elem], 'vs', co[2, elem])) #将当前的一对因子的名称添加到了pairs向量中
F.Model <- c(F.Model, ad$aov.tab[1, 4]) #将当前的Adonis分析的F值添加到了F.Model向量中。ad$aov.tab[1, 4]是Adonis分析结果中的F值。
R2 <- c(R2, ad$aov.tab[1, 5]) #将当前的Adonis分析的R²值添加到了R2向量中。ad$aov.tab[1, 5]是Adonis分析结果中的R²值。
p.value <- c(p.value, ad$aov.tab[1, 6]) #将当前的Adonis分析的p值添加到了p.value向量中。ad$aov.tab[1, 6]是Adonis分析结果中的p值。
}
# 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.txt", row.names = 1, sep = "\t", header = TRUE) %>% as.data.frame()
#读取了"otu_table.txt"文件，并将其转换为数据框。read.table函数的参数指定了文件的路径、行名的位置、字段的分隔符以及是否有表头
map <- read.table("./metadata3-1.txt", sep = "\t", header = TRUE) #读取分组
colnames(map)[1] <- "ID" #这两行代码将map的第一列的列名改为"ID"
row.names(map) <- map$ID #将map的行名设置为"ID"列的值
idx <- rownames(map) %in% colnames(otu) #找出了map的行名中哪些也在otu的列名中，然后将结果保存在idx变量中
#使用idx变量对map和otu进行了子集选择，只保留了map的行名和otu的列名相同的部分。
map1 <- map[idx,]
otu <- otu[, rownames(map1)]
# 4. 进行adonis分析并计算统计值。
bray_curtis <- vegan::vegdist(t(otu), method = "bray", na.rm = TRUE)
#计算了otu数据的Bray-Curtis距离。Bray-Curtis距离是一种用于衡量样本之间差异的指标，常用于生态学中。
#vegdist函数是vegan包中的一个函数，用于计算距离矩阵。使用了"bray"方法，并且设置了na.rm = TRUE来忽略缺失值。
ado <- adonis(bray_curtis ~ map1$Group, permutations = 999, method = "bray") #进行了Adonis分析。
#使用了bray_curtis距离矩阵作为响应变量，map1$Group作为解释变量。
#设置了permutations = 999来指定置换次数，method = "bray"来指定距离计算方法。
R2_value <- round(as.data.frame(ado$aov.tab[5])[1, 1], 3)
#计算了Adonis分析的R²值。R²值是一个衡量模型拟合优度的统计量，表示模型可以解释的数据变异的比例。将R²值四舍五入到了小数点后三位。
p_v_value <- as.data.frame(ado$aov.tab[6])[1, 1]
#计算了Adonis分析的p值。p值是用来判断结果是否显著的一个统计量。
title <- paste("adonis:R ", R2_value, " p: ", p_v_value, sep = "")
#生成了一个包含R²值和p值的标题。
# 5. 绘制PCoA图。
pcoa <- cmdscale(bray_curtis, k = 2, eig = TRUE) #使用cmdscale函数进行了主坐标分析（Principal Coordinates Analysis，PCoA）
#PCoA是一种多元分析方法，可以用来探索和可视化高维数据的结构。
#使用了bray_curtis距离矩阵作为输入，设置了k = 2来指定保留的主坐标数，设置了eig = TRUE来返回特征值。
points <- as.data.frame(pcoa$points) %>% dplyr::rename(x = "V1", y = "V2") #将PCoA的结果转换为数据框，并重命名了列名.
#pcoa$points是PCoA结果中的坐标，dplyr::rename(x = "V1", y = "V2")将列名"V1"和"V2"改为"x"和"y"。
eig <- pcoa$eig
#将PCoA分析结果中的特征值赋值给了eig变量。在主坐标分析（PCoA）中，特征值可以反映每个主坐标的重要性，即它可以解释的数据变异的比例。
points <- cbind(points, map1[match(rownames(points), map1$ID),]) #将map1数据框中的元数据添加到了PCoA的结果中。
#match(rownames(points), map1$ID)找出了points的行名在map1$ID中的位置，然后使用这个结果对map1进行了子集选择。
n <- 0.85
colors <- c("B73_Week4"="#00BFFF","B73_Week6"="#00BFFF","B73_Week8"="#00BFFF","B73_Week10"="#00BFFF","Mo17_Week4"="#FF4500","Mo17_Week6"="#FF4500","Mo17_Week8"="#FF4500","Mo17_Week10"="#FF4500")
# 定义了颜色和形状的映射关系，用于后续的可视化。
shapes <- c("B73_Week4"=24, "B73_Week6"=22, "B73_Week8"=21, "B73_Week10"=23,
"Mo17_Week4"=24, "Mo17_Week6"=22, "Mo17_Week8"=21, "Mo17_Week10"=23)
# 在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 = 9))
#设置了图的主题。你使用了黑白主题，并且移除了网格线，设置了坐标轴文本的颜色和大小。
# 显示绘制的图
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 = 12, height = 10, units = "in", dpi = 600)