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