R : PCOA

library(vegan)
library(tidyverse)
pairwise.adonis1 <-function(x,factors,p.adjust.m)
{
x = x %>% as.matrix()
co = as.matrix(combn(unique(factors),2))
pairs = c()
F.Model =c()
R2 = c()
p.value = c()
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.adjusted =p.adjust(p.value,method=p.adjust.m)
pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)
return(pairw.res)
}
setwd("D:\\Desktop")

#---导入数据
#otu =read.csv("./otu.csv",row.names = 1) %>% t() %>% as.data.frame()
otu =read.table("./OTU_table.txt",row.names = 1,sep = "\t",header=T) %>% as.data.frame()
head(otu)

#map = read.csv("./map.csv")
map = read.table("./metadata3-1.txt", sep = "\t",header=T)
head(map)
colnames(map)[1] = "ID"
row.names(map) = map$ID

idx = rownames(map) %in% colnames(otu)
map1 = map[idx,]
map1
otu = otu[,rownames(map1)]
head(otu)

bray_curtis =vegan:: vegdist(t(otu),method = "bray", na.rm=TRUE)

pcoa <- otu %>%
t() %>%
vegan:: vegdist(method = "bray", na.rm=TRUE) %>%
cmdscale(k=2, eig=T) ###数据矩阵的多维尺度

points <- pcoa$points %>% as.data.frame() %>%
# as.tibble() %>%
dplyr::rename(.,x = "V1",y = "V2" )
head(points)

eig = pcoa$eig
eig
points = cbind(points, map1[match(rownames(points), map1$ID), ])
ado = adonis(bray_curtis~ map1$Group,permutations = 999,method="bray")
a = round(as.data.frame(ado$aov.tab[5])[1,1],3)
R2 <- paste("adonis:R ",a, sep = "")
b = as.data.frame(ado$aov.tab[6])[1,1]
p_v = paste("p: ",b, sep = "")
title = paste(R2," ",p_v, sep = "")
title
# anosim.result<-anosim(bray_curtis, map1$Group,permutations = 999)
# summary(anosim.result)

pair_bray_adonis = pairwise.adonis1(bray_curtis,map1$Group, p.adjust.m= "bonferroni")
head(pair_bray_adonis)

library("ggalt")
library(car)
library(ggforce)

n = 0.85
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=c("#E41A1C","#E6A429","#44AC41","#377EB8","#B055B0"))+ ##底色设置
scale_color_manual(values=c("#BD3C29","#0172B6","#E18727","#44AC41","#B055B0"))+
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())+
theme(axis.text = element_text(color = "black",size =9))

p1

library(ggpubr)
library(patchwork)
tab2 <- ggtexttable(pair_bray_adonis, rows = NULL)

p2 = p1 | tab2
p2

##############################################################################################################

这段R语言代码主要用于分析生态数据的差异和相似性。以下是代码的主要步骤和功能：

1. 导入所需的R包（vegan、tidyverse等）。
2. 定义了一个名为pairwise.adonis1的函数，用于执行多个配对的adonis分析，并计算F模型值、R方和p值。函数中使用了adonis函数来执行adonis分析。
3. 设置工作目录，读取并导入数据文件（OTU表和元数据表）。
4. 对OTU表进行初步处理，保留与元数据匹配的样本。
5. 计算样本间的Bray-Curtis距离矩阵。
6. 对Bray-Curtis距离矩阵进行PCoA分析，提取前两个主坐标轴的值。
7. 将PCoA结果与元数据合并，用于后续的可视化和统计分析。
8. 使用adonis函数执行整体的adonis分析，并计算R方和p值。
9. 进行配对的adonis分析，计算配对间的F模型值、R方、p值，并进行多重检验校正。
10. 使用ggplot2包进行PCoA的可视化，根据元数据中的分组信息对样本进行着色，并绘制群落组成的椭圆。
11. 使用ggtexttable函数将配对的adonis结果以表格形式展示。
12. 使用patchwork包将PCoA图和配对的adonis结果表格合并为一个图形。

###############################################################################################################

library(vegan)
library(tidyverse)

这两行代码用于加载所需的R包（vegan和tidyverse）。

这段代码定义了一个名为pairwise.adonis1的函数，用于执行多个配对的adonis分析。函数中的循环用于逐对比较不同组合之间的差异，并计算F模型值、R方和p值。