# 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()
# 对factors中的每一对执行adonis分析
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)
pairw.res = data.frame(pairs, F.Model, R2, p.value, p.adjusted)
return(pairw.res)
}
# 3. 读取和处理数据。
setwd("C:\\Users\\Administrator\\Desktop")
otu <- read.table("./otu_table.txt", row.names = 1, sep = "\t", header = TRUE) %>% as.data.frame()
map <- read.table("./metadata3-1.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")
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 ", R2_value, " p: ", 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"="#00BFFF","Mo17"="#FF4500")
p1 <- ggplot(points, aes(x = x, y = y, fill = Group)) +
geom_point(alpha = .7, size = 5, pch = 21) +
scale_shape_manual(values = points$Group) +
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) +
geom_mark_ellipse(aes(x = x * n, y = y * n, 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")
# 存储为文本文件
write.table(as.data.frame(pair_bray_adonis), "table.txt", sep = "\t", quote = FALSE, row.names = FALSE)
tab2 <- ggtexttable(pair_bray_adonis, rows = NULL)
p2 <- tab2
p2
# 使用ggsave保存PCoA图为PNG格式
ggsave(filename = "PCoA_plot.png", plot = p1, width = 10, height = 8, units = "in", dpi = 300)
以下是对给定的R代码的逐行解释:
# 1. 导入所需的库。
这部分导入了一系列R包,用于数据处理、统计分析和绘图。
library(vegan) # 用于生态多样性分析和多元统计分析。
library(tidyverse) # 包括一系列数据处理和可视化的包。
library(ggalt) # 提供了一些ggplot2中没有的图层、坐标系统等。
library(car) # 为线性模型提供了一些高级函数和数据集。
library(ggforce) # ggplot2的扩展,提供额外的几何对象和统计变换。
library(ggpubr) # 用于绘制出版质量的ggplot2图。
library(patchwork) # 允许将多个ggplot图组合成一个图。
# 2. 定义所需的函数。
这部分定义了一个函数 pairwise.adonis1
,用于进行成对的adonis分析(一种多元回归方法),并返回统计结果。
# 3. 读取和处理数据。
这部分的代码首先设置工作目录,然后读取两个数据文件:otu_table.txt
(OTU表)和metadata3-1.txt
(样本元数据)。接着,代码将OTU表和元数据根据共有的样本名称进行对齐。
# 4. 进行adonis分析并计算统计值。
首先,代码计算了基于Bray-Curtis距离的样本间距离矩阵。然后,使用adonis分析测试了一个组(从元数据中提取)是否可以解释距离矩阵中的变化。最后,代码计算了adonis分析的R²值和p值,并为后面的图表标题生成了一个字符串。
# 5. 绘制PCoA图。
这部分的代码首先进行了主坐标分析(PCoA),然后使用ggplot2绘制了PCoA图,显示了每个样本在第一和第二主坐标上的位置,并根据其组标记颜色。此外,还为每个组绘制了一个置信椭圆。
# 6. 输出pairwise adonis结果。
这部分的代码首先使用之前定义的
pairwise.adonis1
函数计算了成对的adonis结果,然后将这些结果保存为一个名为table.txt
的文本文件。接着,它使用ggtexttable函数将统计结果转化为一个ggplot对象,并显示它。最后,它使用ggsave函数将PCoA图保存为一个PNG文件。