R：NMDS（第三版）

# 清除所有变量
rm(list = ls())

# 设置工作目录
setwd("C:\\Users\\Administrator\\Desktop\\新建文件夹\\NMDS")

# 导入所需的库
library(vegan)
library(ggplot2)

# 读取数据
mydata <- read.delim("otu_table_R.txt", header = TRUE, sep = "\t")
mygroups <- read.delim("group.txt", header = TRUE, sep = "\t")

# 数据处理
# 将丰度数据的第一列（分类信息）转换为行名称
rownames(mydata) <- mydata$clade_name
mydata <- mydata[, -1] # 删除分类信息列

# 计算Bray-Curtis距离
distance <- vegdist(t(mydata), method = 'bray')

# NMDS排序分析
nmds <- metaMDS(distance, k = 2)

# 获取应力值
stress <- nmds$stress

# 将NMDS结果转化为数据框
df <- as.data.frame(nmds$points)
colnames(df) <- c("MDS1", "MDS2") # 重命名列以便合并

# 将分组数据的SampleID列转换为行名称并与NMDS结果合并
rownames(mygroups) <- mygroups$SampleID
df <- merge(df, mygroups, by = "row.names")
df$SampleID <- df$Row.names
df <- df[, -1] # 删除辅助列

# 设置随机种子
set.seed(123)

# 基于Bray-Curtis距离进行PERMANOVA分析
adonis <- adonis2(as.dist(distance) ~ Group, data = df, permutations = 999)

# 应力值、Adonis R2与显著性
stress_text <- paste("Stress =", round(stress, 4))
adonis_text <- paste("Adonis: R^2 =", round(adonis$R2[1], 2), "P-value =", adonis$`Pr(>F)`[1])
title <- paste(stress_text, adonis_text)

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") #定义顺序
df$Group <- factor(df$Group, levels = levels_order)
p1 <- ggplot(df, aes(x = MDS1, y = MDS2, 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("NMDS 1", sep = ""),
       y = paste("NMDS 2", 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

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