R语言分析α、β多样性（TAXA分组）

抽平、TAXA分组

flatten&abundance_calc.R
rm(list = ls())
library(vegan)
input_otu <- "otu_corrected"

# 读取OTU表
otu_table <- read.table(paste0(input_otu,".csv"), header = TRUE, row.names = 1, sep = ",")

# 计算每个样本的总序列数
total_sequences <- colSums(otu_table)

#使用该代码进行抽平
otu_flattened = as.data.frame(t(rrarefy(t(otu_table), min(colSums(otu_table)))))
write.table (otu_flattened, , col.names = NA,file =paste0(input_otu,"_flattened.csv"),sep =",") #结果导出

# 计算相对丰度
relative_abundance <- sweep(otu_flattened, 2, total_sequences, FUN = "/")

# 定义分组函数
classify_taxa <- function(row) {
  max_abundance <- max(row)
  min_abundance <- min(row)
  if (max_abundance <= 0) {
    return("None")  # 全为0行
  } else if (max_abundance <= 0.0001) {
    return("RT")  # 稀有类群 rare taxa
  } else if (min_abundance >= 0.01) {
    return("AT")  # 丰富类群 abundant taxa，AT
  } else if (min_abundance > 0.0001 & max_abundance < 0.01) {
    return("MT")  # 中间类群 moderate taxa，MT
  } else if (min_abundance <= 0.0001 & max_abundance < 0.01) {
    return("CRT")  # 条件稀有类群 conditionally rare taxa，CRT
  } else if (min_abundance > 0.0001 & max_abundance >= 0.01) {
    return("CAT")  # 条件丰富类群 conditionally abundant taxa，CAT
  } else if (min_abundance <= 0.0001 & max_abundance >= 0.01) {
    return("CRAT")  # 条件稀有或丰富类群 conditionally rare or abundant taxa，CRAT
  } else {
    return("Unclassified")  # 未分类
  }
}

# 应用分组函数
taxa <- apply(relative_abundance, 1, classify_taxa)

# 添加taxa列
otu_flattened$taxa <- taxa
write.table(otu_flattened, paste0(input_otu,"_flattened_taxa.csv"), col.names = NA, sep = ',', quote = FALSE)

# 生成taxaCombined列
otu_flattened$taxaCombined <- ifelse(otu_flattened$taxa %in% c("AT", "CAT", "CRAT"), "AT",
                        ifelse(otu_flattened$taxa %in% c("RT", "CRT"), "RT&CRT", otu_flattened$taxa))

# 创建仅包含原始数据和taxaCombined列的新数据框
otu_flattened_combined <- otu_flattened[, !names(otu_flattened) %in% "taxa"]

# 输出文件
write.table(otu_flattened_combined, paste0(input_otu,"_flattened_taxaCombined.csv"), col.names = NA, sep = ',', quote = FALSE)

α多样性

代码说明

1. 读取数据并提取分组信息：
   • 使用 read.csv 读取 final_resample_otu_table_sediment_taxa.csv 文件，得到包含 OTU 数据和分组信息的完整数据框。
   • 提取最后一列 taxa 作为分组信息，并获取所有不同的 taxa 值。
2. 定义计算 alpha 多样性的函数：
   • 定义 alpha 函数，计算多种 alpha 多样性指数，包括 Richness, Shannon, Simpson, Pielou, Chao1, ACE, goods_coverage。
3. 遍历每个 taxa 组进行分析：
   • 对于每个 taxa 值，筛选出该组的数据，并移除 taxa 列。
   • 转置数据，使每列代表一个样本。
   • 使用 alpha 函数计算该组的 alpha 多样性指数，并添加组信息。
   • 将结果存储在一个列表中，并保存为独立的 CSV 文件。
4. 合并所有组的数据：
   • 使用 do.call(rbind, all_alpha_data) 将所有组的 alpha 多样性数据合并成一个数据框。
5. 绘制各指数的箱线图：
   • 定义需要绘制的指数列表 indices。
   • 对每个指数，使用 melt 函数将数据重塑为长格式。
   • 使用 ggplot2 创建箱线图，x 轴为组名，y 轴为指数值，使用不同颜色填充表示不同组。
   • 将所有图形存储在一个列表中，并使用 wrap_plots 函数将它们组合成一个 1 行 5 列的布局。
   • 保存组合图像为 SVG 文件。

alpha_taxa.r
rm(list = ls())

## 加载必要的包
library(vegan)
library(ggplot2)
library(reshape2)
library(patchwork)

filename <- "otu_corrected_flattened"
Legend_column <- "Group" #分组依据是csv中的哪一列

## 读入分组文件
groups <- read.table(paste0(filename, "_group.csv"), header = TRUE, row.names = 1, sep = ",")
groups <- groups[, Legend_column, drop = FALSE]
## 读入数据并进行转置
otu <- read.table(paste0(filename, "_taxaCombined.csv"), sep = ",", header = TRUE, row.names = 1)

## 提取最后一列的taxa值
taxa_values <- unique(otu$taxa)
# 去除 "None" 值
taxa_values <- setdiff(taxa_values, "None")

## 定义 alpha 多样性计算函数
alpha <- function(x, groups, base = exp(1)) {
  est <- estimateR(x)
  Richness <- est[1, ]
  Chao1 <- est[2, ]
  ACE <- est[4, ]
  Shannon <- diversity(x, index = "shannon", base = base)
  Simpson <- diversity(x, index = "simpson")
  Pielou <- Shannon / log(Richness, base)
  goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
  
  result <- data.frame(SITE = rownames(x), GROUP = groups, Richness, Shannon, Simpson, Pielou, Chao1, ACE, goods_coverage)
  colnames(result) <- c("SITE", "GROUP", "RICHNESS", "SHANNON", "SIMPSON", "PIELOU", "CHAO1", "ACE", "GOODS_COVERAGE")
  result
}

## 创建一个列表用于存储所有图形
plot_list <- list()

## 针对每个taxa值进行处理
for (taxa_value in taxa_values) {
  subset_data <- otu[otu$taxa == taxa_value, -ncol(otu)]  # 去除最后一列 taxa
  subset_data <- t(subset_data)
  
  ## 使用已读入的分组文件中的组信息
  subset_groups <- groups[rownames(subset_data), , drop = FALSE]
  alpha_all <- alpha(subset_data, subset_groups[[1]], base = 2)
  
  ## 保存为 CSV 文件
  #write.csv(alpha_all, file = paste0("alpha_all_", taxa_value, ".csv"), row.names = FALSE, quote = FALSE)
  
  ## 绘制各指数箱线图
  indices <- c("SHANNON", "RICHNESS", "PIELOU", "SIMPSON", "CHAO1")
  
  # 创建一个空列表来存储图形对象
  plots <- list()
  
  # 使用循环创建每个指数的箱线图
  for (i in seq_along(indices)) {
    index <- indices[i]
    data_melted <- melt(alpha_all, id.vars = c("GROUP", "SITE"), measure.vars = index)
    colnames(data_melted)[which(colnames(data_melted) == "value")] <- "value"
    p <- ggplot(data_melted, aes(x = GROUP, y = value, fill = GROUP)) +
      geom_boxplot() +
      labs(x = paste(index, "Index"), y = "") +
      theme_bw() +
      theme(axis.text.x = element_blank())
    if (i == length(indices)) {
      p <- p + labs(fill = paste0("Group-",taxa_value)) # 仅在最后一个图表上添加图例
    } else {
      p <- p + guides(fill = FALSE) # 隐藏其余图表的图例
    }
    plots[[index]] <- p
  }
  
  # 使用patchwork将所有图表组合成一个1行5列的布局
  p_combined <- wrap_plots(plots, ncol = 5)
  
  # 将合并后的图表添加到列表中
  plot_list[[taxa_value]] <- p_combined
  
  # 输出处理完成信息
  cat("Processed taxa:", taxa_value, "\n")
}

# 将所有taxa_value的图合并为一个图
final_plot <- wrap_plots(plot_list, ncol = 1)

# 保存合并后的图像
ggsave(paste0('output\\diversity\\alpha_',filename,'.png'), final_plot, width = 6 * 0.6 * 5, height = 4 * length(plot_list))

# 输出处理完成信息
cat("End of processing.\n")

其中，final_resample_otu_table_sediment_taxa.csv的格式如：

OTUID,D109_1,D109_2,D109_3,D109_4,D109_5,D109_6,D109_7,D109_8,D111_1,D111_2,D111_3,D111_4,D111_5,D111_6,D111_7,D112_1,D112_2,D112_3,D112_4,D112_5,D112_6,D112_7,D112_8,D113_1,D113_2,D113_3,D113_4,D113_5,D113_6,D113_7,D113_8,M1_1,M1_2,M1_3,M2_1,M2_2,M2_3,D148_1,D148_2,D148_3,D148_4,D148_5,D148_6,D148_7,D148_8,D149_1,D149_2,D149_3,D149_4,D149_5,D149_6,D149_7,D149_8,D150_1,D150_2,D150_3,D150_4,D150_5,D150_6,D150_7,D150_8,D151_1,D151_2,D151_3,D151_4,D151_5,D151_6,D151_7,D151_8,D152_1,D152_2,D152_3,D152_4,D152_5,D152_6,D152_7,D152_8,S01_1,S01_2,S01_3,S01_4,S01_5,S01_6,S01_7,taxa
OTU_1,1,1,0,1599,4341,4409,2869,3226,0,0,0,1891,11086,2327,3318,0,0,0,2275,5208,1449,5727,2546,0,0,1,8505,10716,10983,8403,10830,5,1710,8734,1,0,6153,713,1036,465,1129,404,304,5541,598,831,384,564,1419,563,340,677,758,8153,3084,2262,11689,3769,1656,2713,5313,1575,3188,4494,5807,615,6821,11074,6589,2693,2243,2141,1246,2662,7719,10479,2439,473,505,680,272,1592,162,2212,CRAT

β多样性

基于OTU表 (bray、jaccard)

nmds_bray&jaccard_taxa.R
library(vegan)
library(ape)
library(ggplot2)
library(grid)
library(ggalt)
library(dplyr)
library(multcomp)
library(patchwork)

rm(list = ls())

method <-"bray" #bray jaccard #https://rdocumentation.org/packages/vegan/versions/2.6-4/topics/vegdist 

filename<-"otu_corrected_flattened"
output_url <- paste0("beta_", filename, "_NMDS_", method)

groups = read.table(paste0(filename,"_group.csv"), header=T, row.names=1, sep=",")
Legend_column <- "Group" #分组依据是csv中的哪一列
data = read.table(paste0(filename,"_taxaCombined.csv"), header=T, row.names=1, sep=",")

taxa_values <- unique(data$taxa)
# 去除 "None" 值
taxa_values <- setdiff(taxa_values, "None")

# 创建一个列表用于存储所有图形
plot_list <- list()
adonis_label_list <- list()  # 创建一个列表用于存储所有标签


# 针对每个 taxa 值进行处理
for (taxa_value in taxa_values) {
  # 筛选出对应 taxa 值的行
  subset_data <- data[data$taxa == taxa_value, -ncol(data)]  # 去除最后一列 taxa
  subset_data = t(subset_data) # 每行为一个样品
  dist_matrix <- vegdist(subset_data, method = method)# 计算距离矩阵
  
  # 进行 NMDS 分析
  nmds <- metaMDS(dist_matrix, distance = method, k = 2)
  NMDS1 <- nmds$points[, 1]
  NMDS2 <- nmds$points[, 2]
  
  # 创建数据框包含 NMDS 坐标和分组信息
  plotdata <- data.frame(sample = rownames(nmds$points), NMDS1, NMDS2, groups)
  colnames(plotdata)[1:3] <- c("sample", "NMDS1", "NMDS2") # 重命名1~3列
  names(plotdata)[names(plotdata) == paste0(Legend_column)] <- 'Group' # 重命Group列
  plotdata$Group <- factor(plotdata$Group)
  
  p1 <- ggplot(plotdata, aes(Group, NMDS1)) +
    geom_boxplot(aes(fill = Group)) +
    coord_flip() +
    theme_bw() +
    theme(axis.ticks.length = unit(.4, "lines"), 
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_text(size = 15),
          axis.text.x = element_blank(),
          legend.position = "none")
  
  p3 <- ggplot(plotdata, aes(Group, NMDS2)) +
    geom_boxplot(aes(fill = Group)) +
    theme_bw() +
    theme(axis.ticks.length = unit(0.4, "lines"), 
          axis.ticks = element_line(color = 'black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x = element_text(colour = 'black', size = 15, angle = 45,
                                     vjust = 1, hjust = 1),
          axis.text.y = element_blank(),
          legend.position = "none")
  
  p2 <- ggplot(plotdata, aes(NMDS1, NMDS2)) +
    geom_point(aes(fill = Group, color = Group), size = 5, pch = 21, stroke = 2, alpha = .8) +
    labs(fill = Legend_column) + labs(color = Legend_column) +
    stat_ellipse(level = .8, aes(colour = Group)) + # 画置信圈, alpha是透明度
    xlab("NMDS 1") + 
    ylab("NMDS 2") +
    theme_classic() +
    theme(axis.ticks.length = unit(0.4, "lines"), axis.ticks = element_line(color = 'black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x = element_text(colour = 'black', size = 18, vjust = 3),
          axis.title.y = element_text(colour = 'black', size = 18, vjust = -1),
          axis.text = element_text(colour = 'black', size = 14),
          legend.title = element_text(size = 14),
          legend.text = element_text(size = 12),
          legend.key = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1),
          legend.background = element_rect(colour = "black", fill = alpha("white", 0.4)))
  
  otu.adonis <- adonis2(subset_data ~ groups[[paste0(Legend_column)]], data = groups, distance = method)
  
  # 创建标签字符串
  adonis_label <- paste0(method, "-", taxa_value, "\ndf = ", otu.adonis$Df[1],  
                         "\nR² = ", round(otu.adonis$R2[1], 4),  
                         "\nP-value = ", otu.adonis$`Pr(>F)`[1])
  label_grob <- textGrob(adonis_label, x = 0.5, y = 0.5, gp = gpar(fontsize = 12, fontfamily = "serif"))
  
  p4 <- wrap_elements(grid::grobTree(label_grob))+
    theme(panel.border = element_rect(colour = "black", fill = NA, size = 1))  # 添加边框
  
  p5 <- p1 + p4 + p2 + p3 + 
    plot_layout(heights = c(1, 4), widths = c(4, 1), ncol = 2, nrow = 2)
  
  # 将生成的图添加到列表中
  plot_list[[taxa_value]] <- p5
  #ggsave(paste0("output\\",output_url,"_",taxa_value,".svg"), p5, width = 12, height = 12/sqrt(2))
  # 输出处理完成信息
  cat("Processed taxa:", taxa_value, "\n")
}

# 将所有图形合并为一个图
final_plot <- wrap_plots(plot_list,ncol = length(taxa_values))

# 保存最终图像
ggsave(paste0("output\\diversity\\",output_url, ".png"), final_plot, width = 12*length(taxa_values), height = 12/sqrt(2),limitsize = FALSE)

# 输出处理完成信息
cat("End of processing.\n")

PCOA_bray&jaccard_taxa.R
library(vegan)
library(ape)
library(ggplot2)
library(grid)
library(ggalt)
library(dplyr)
library(multcomp)
library(patchwork)

rm(list = ls())

method <-"bray" #bray jaccard #https://rdocumentation.org/packages/vegan/versions/2.6-4/topics/vegdist 

filename<-"final_resample_otu_table_sediment_taxa"
output_url <- paste0("beta_", filename, "_PCOA_", method)

groups = read.table(paste0(filename,"_group.csv"), header=T, row.names=1, sep=",")
Legend_column <- "Group" #分组依据是csv中的哪一列
data_all = read.table(paste0(filename,".csv"), header=T, row.names=1, sep=",")

taxa_values <- unique(data_all$taxa)
# 去除 "None" 值
taxa_values <- setdiff(taxa_values, "None")

# 创建一个列表用于存储所有图形
plot_list <- list()
adonis_label_list <- list()  # 创建一个列表用于存储所有标签


# 针对每个 taxa 值进行处理
for (taxa_value in taxa_values) {
  # 筛选出对应 taxa 值的行
  data <- data_all[data_all$taxa == taxa_value, -ncol(data_all)]  # 去除最后一列 taxa
  data = t(data) # 每行为一个样品
  
  length=nrow(groups)
  times1=length
  res1=length
  times2=length
  res2=length
  col1=rep(1:8,times1)
  col=c(col1,1:res1)
  pich1=rep(c(21:24),times2)
  pich=c(pich1,15:(15+res2))
  
  data <- vegdist(data, method = method)
  pcoa<- pcoa(data, correction = "none", rn = NULL)
  PC1 = pcoa$vectors[,1]
  PC2 = pcoa$vectors[,2]
  plotdata <- data.frame(rownames(pcoa$vectors),PC1,PC2,groups)
  colnames(plotdata)[1:3] <-c("sample","PC1","PC2") #重命名1~3列
  names(plotdata)[names(plotdata) == paste0(Legend_column)] <- 'Group' #重命Group列
  pc1 <-floor(pcoa$values$Relative_eig[1]*100)
  pc2 <-floor(pcoa$values$Relative_eig[2]*100)
  plotdata$Group <- factor(plotdata$Group)
  
  yf <- plotdata
  yd1 <- yf %>% group_by(Group) %>% summarise(Max = max(PC1))
  yd2 <- yf %>% group_by(Group) %>% summarise(Max = max(PC2))
  yd1$Max <- yd1$Max + max(yd1$Max)*0.1
  yd2$Max <- yd2$Max + max(yd2$Max)*0.1
  fit1 <- aov(PC1~Group,data = plotdata)
  tuk1<-glht(fit1,linfct=mcp(Group="Tukey"))
  res1 <- cld(tuk1,alpah=0.05)
  fit2 <- aov(PC2~Group,data = plotdata)
  tuk2<-glht(fit2,linfct=mcp(Group="Tukey"))
  res2 <- cld(tuk2,alpah=0.05)
  test <- data.frame(PC1 = res1$mcletters$Letters,PC2 = res2$mcletters$Letters,
                     yd1 = yd1$Max,yd2 = yd2$Max,Group = yd1$Group)
  test$Group <- factor(test$Group)
  
  p1 <- ggplot(plotdata,aes(Group,PC1)) +
    geom_boxplot(aes(fill = Group)) +
    geom_text(data = test,aes(x = Group,y = yd1,label = PC1),size = 7) +
    coord_flip() +
    theme_bw()+
    theme(axis.ticks.length = unit(.4,"lines"), 
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_text(size=15),
          axis.text.x=element_blank(),
          legend.position = "none")
  yrange <- ggplot_build(p1)$layout$panel_scales_y[[1]]$range$range
  p1
  
  p3 <- ggplot(plotdata,aes(Group,PC2)) +
    geom_boxplot(aes(fill = Group)) +
    geom_text(data = test,aes(x = Group,y = yd2,label = PC2),
              size = 7,color = "black",fontface = "bold") +
    theme_bw()+
    theme(axis.ticks.length = unit(0.4,"lines"), 
          axis.ticks = element_line(color='black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.x=element_text(colour='black',size=15,angle = 45,
                                   vjust = 1,hjust = 1),
          axis.text.y=element_blank(),
          legend.position = "none")
  yrange <- ggplot_build(p3)$layout$panel_scales_y[[1]]$range$range
  p3 <- p3
  p3
  
  p2<-ggplot(plotdata, aes(PC1, PC2)) +
    geom_encircle(aes(fill=Group,colour=Group),alpha=.1, show.legend=F,size=2,expand=0.05)+ #同组阴影底色
    #stat_ellipse(level=.6,aes(colour=Group))+ #画置信圈,alpha是透明度
    geom_point(aes(fill=Group,colour=Group),size=5,pch =21, stroke=2, alpha=.8)+ #点
    #geom_text(aes(label = Mark), size = 3)+
    xlab(paste0("PC1 ( ",pc1,"%"," )")) + 
    ylab(paste0("PC2 ( ",pc2,"%"," )"))+
    xlim(ggplot_build(p1)$layout$panel_scales_y[[1]]$range$range) +
    ylim(ggplot_build(p3)$layout$panel_scales_y[[1]]$range$range) +
    theme(text=element_text(size=18))+
    geom_vline(aes(xintercept = 0),linetype="dotted")+
    geom_hline(aes(yintercept = 0),linetype="dotted")+
    theme_classic()+
    theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x=element_text(colour='black', size=18,vjust=3),
          axis.title.y=element_text(colour='black', size=18,vjust=-1),
          axis.text=element_text(colour='black',size=14),
          legend.title=element_text(size = 14),
          legend.text=element_text(size=12),
          legend.key=element_blank(),legend.position = c(1, 1),legend.justification = c(1,1),
          legend.background = element_rect(colour = "black",fill=alpha("white",0.4)),
    )
  p2
  
  otu.adonis=adonis2(data~groups[[paste0(Legend_column)]],data = groups,distance = method)
  adonis_label <- paste0(method, "-", taxa_value, "\ndf = ", otu.adonis$Df[1],  
                         "\nR² = ", round(otu.adonis$R2[1], 4),  
                         "\nP-value = ", otu.adonis$`Pr(>F)`[1])
  label_grob <- textGrob(adonis_label, x = 0.5, y = 0.5, gp = gpar(fontsize = 12, fontfamily = "serif"))
  p4 <- wrap_elements(grid::grobTree(label_grob))+
    theme(panel.border = element_rect(colour = "black", fill = NA, size = 1))  # 添加边框
  
  p5 <- p1 + p4 + p2 + p3 + 
    plot_layout(heights = c(1,4),widths = c(4,1),ncol = 2,nrow = 2)
  p5
  plot_list[[taxa_value]] <- p5
  
  # 输出处理完成信息
  cat("Processed taxa:", taxa_value, "\n")
}

# 将所有图形合并为一个图
final_plot <- wrap_plots(plot_list,ncol = length(taxa_values))

# 保存最终图像
ggsave(paste0("output\\diversity\\",output_url, ".png"), final_plot, width = 12*length(taxa_values), height = 12/sqrt(2))

# 输出处理完成信息
cat("End of processing.\n")

基于OTU表和进化树 (加权、不加权unifrac)

nmds_unifrac_taxa.R
library(vegan)
library(ape)
library(ggplot2)
library(grid)
library(ggalt)
library(dplyr)
library(multcomp)
library(patchwork)
library(ape)

# 设置工作目录
setwd("E:/ShareCache/论文/深度驱动了丰富类群和稀有类群")
rm(list = ls())

method <- "unifrac"

filename<-"otu_corrected_flattened"

output_url <- paste0("beta_", filename, "_NMDS_", method)

tree <- read.tree(paste0(filename,"_tree.nwk"))
Group = read.table(paste0(filename,"_group.csv"), header=T, row.names=1, sep=",")
Legend_column <- "Group" #分组依据是csv中的哪一列
data = read.table(paste0(filename,"_taxaCombined.csv"), header=T, row.names=1, sep=",")

taxa_values <- unique(data$taxa)
# 去除 "None" 值
taxa_values <- setdiff(taxa_values, "None")

# 创建一个列表用于存储所有图形
plot_list <- list()
adonis_label_list <- list()  # 创建一个列表用于存储所有标签

# 针对每个 taxa 值进行处理
for (taxa_value in taxa_values) {
  # 筛选出对应 taxa 值的行
  subset_data <- data[data$taxa == taxa_value, -ncol(data)]  # 去除最后一列 taxa
  
  # 构建phyloseq对象
  ps <- phyloseq(otu_table(as.matrix(subset_data), taxa_are_rows = TRUE),
                 sample_data(Group),
                 phy_tree(tree))
  
  
  # 从phyloseq对象中提取元数据
  meta <- as.data.frame(sample_data(ps))
  
  # 计算加权UniFrac距离矩阵
  dist_matrix <- as.matrix(distance(ps, method = method))
  
  # 计算NMDS
  nmds <- metaMDS(dist_matrix, distance = method, k = 2)
  # 将NMDS结果转换为数据框架
  nmds_df <- as.data.frame(nmds$points)
  meta <- as.data.frame(sample_data(ps))
  otu.adonis <- adonis2(dist_matrix ~ Group, data=Group)
  
  NMDS1 = nmds$points[,1]
  NMDS2 = nmds$points[,2]
  plotdata <- data.frame(sample=rownames(nmds$points), NMDS1, NMDS2, Group=meta$Group)
  plotdata$Group <- factor(plotdata$Group) # 确保分组变量是因子类型
  
  p1 <- ggplot(plotdata,aes(Group,NMDS1)) +
    geom_boxplot(aes(fill = Group)) +
    coord_flip() +
    theme_bw()+
    theme(axis.ticks.length = unit(.4,"lines"), 
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_text(size=15),
          axis.text.x=element_blank(),
          legend.position = "none")
  yrange <- ggplot_build(p1)$layout$panel_scales_y[[1]]$range$range
  p1
  
  
  p3 <- ggplot(plotdata,aes(Group,NMDS2)) +
    geom_boxplot(aes(fill = Group)) +
    theme_bw()+
    theme(axis.ticks.length = unit(0.4,"lines"), 
          axis.ticks = element_line(color='black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.x=element_text(colour='black',size=15,angle = 45,
                                   vjust = 1,hjust = 1),
          axis.text.y=element_blank(),
          legend.position = "none")
  yrange <- ggplot_build(p3)$layout$panel_scales_y[[1]]$range$range
  p3
  
  
  # 使用ggplot2包进行绘图
  p2<-ggplot(nmds_df, aes(MDS1, MDS2, color = meta$Group)) +
    labs(color = "Group",fill="Group")+
    #geom_encircle(aes(fill=sample_data(ps)$Group,colour=sample_data(ps)$Group),alpha=.1, show.legend=F,size=2,expand=0.05)+ #同组阴影底色
    stat_ellipse(level=.8,aes(colour=sample_data(ps)$Group))+ #画置信圈,alpha是透明度
    geom_point(aes(fill=sample_data(ps)$Group,colour=sample_data(ps)$Group),size=5,pch =21, stroke=2, alpha=.8)+ #点
    #geom_text(aes(label = Mark), size = 3)+
    xlab("NMDS 1") + 
    ylab("NMDS 2")+
    xlim(ggplot_build(p1)$layout$panel_scales_y[[1]]$range$range) +
    ylim(ggplot_build(p3)$layout$panel_scales_y[[1]]$range$range) +
    theme(text=element_text(size=18))+
    geom_vline(aes(xintercept = 0),linetype="dotted")+
    geom_hline(aes(yintercept = 0),linetype="dotted")+
    theme_classic()+
    theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'),
          axis.line = element_line(colour = "black"), 
          axis.title.x=element_text(colour='black', size=18,vjust=3),
          axis.title.y=element_text(colour='black', size=18,vjust=-1),
          axis.text=element_text(colour='black',size=14),
          legend.title=element_text(size = 14),
          legend.text=element_text(size=12),
          legend.key=element_blank(),legend.position = c(1, 1),legend.justification = c(1,1),
          legend.background = element_rect(colour = "black",fill=alpha("white",0.4)),
    )
  p2
  
  # 创建标签字符串
  adonis_label <- paste0(method, "-", taxa_value, "\ndf = ", otu.adonis$Df[1],  
                         "\nR² = ", round(otu.adonis$R2[1], 4),  
                         "\nP-value = ", otu.adonis$`Pr(>F)`[1])
  label_grob <- textGrob(adonis_label, x = 0.5, y = 0.5, gp = gpar(fontsize = 12, fontfamily = "serif"))
  
  p4 <- wrap_elements(grid::grobTree(label_grob))+
    theme(panel.border = element_rect(colour = "black", fill = NA, size = 1))  # 添加边框
  
  p5 <- p1 + p4 + p2 + p3 + 
    plot_layout(heights = c(1,4),widths = c(4,1),ncol = 2,nrow = 2)
  
  plot_list[[taxa_value]] <- p5
  # 输出处理完成信息
  cat("Processed taxa:", taxa_value, "\n")
}

# 将所有图形合并为一个图
final_plot <- wrap_plots(plot_list,ncol = length(taxa_values))

# 保存最终图像
ggsave(paste0("output\\diversity\\",output_url, ".png"), final_plot, width = 12*length(taxa_values), height = 12/sqrt(2),limitsize = FALSE)

# 输出处理完成信息
cat("End of processing.\n")