本教程相关代码已经上传至 https://github.com/iMetaScience/iMetaPlot/tree/main/221129PCoA 如果你使用本代码，请引用： Yangquanwei Zhong. 2022. Differential microbial assembly processes and co-occurrence networks in the soil-root continuum along an environmental gradient. iMeta 1: e18. https://onlinelibrary.wiley.com/doi/full/10.1002/imt2.18

代码编写及注释：农心生信工作室

# 加载包
library(phyloseq)
library(vegan)
library(ape)
library(ggplot2)

set.seed(42) # 固定随机种子
#随机生成1000个otu之间的进化树
tree <- rtree(1000, rooted = T, tip.label = paste("s",1:1000, sep=""))
#随机生成OTU表 10个样品各1000个otu
otu <- matrix(sample(0:10, 10000, replace = T), nrow = 10)

#给每个otu命名，和进化树节点名称保持一致
colnames(otu) <- tree$tip.label
df <- phyloseq(otu_table(otu, taxa_are_rows = F), phy_tree(tree)) # 生成phyloseq对象
#计算非加权的uniFrac距离矩阵，衡量10个样品之间的群落谱系差异
dist.mat <- phyloseq::UniFrac(df, weighted=F)
# 用非加权的uniFrac距离矩阵做pcoa分析
df.pcoa <- pcoa(dist.mat)

# 提取pcoa结果表格，用于绘图
df.plot<-data.frame(df.pcoa$vectors)

#获得坐标轴上显示的百分比
x_label<-round(df.pcoa$values$Rel_corr_eig[1]*100, 2)
y_label<-round(df.pcoa$values$Rel_corr_eig[2]*100, 2)
#对数据进行人为分组
df.plot$group <- ifelse(df.plot$Axis.1<0,"Agropyron cristatum","Reaumuria Linn.")

#用ggplot2画图，然后添加一个表示分组的椭圆
p <- ggplot(data=df.plot,aes(x=Axis.1,y=Axis.2,color=group,shape=group))+geom_point(size=5)+theme_bw()+theme(panel.grid = element_blank())+geom_vline(xintercept = 0,lty="dashed")+geom_hline(yintercept = 0,lty="dashed")+labs(x=paste0("PCoA1",x_label,"%"),y=paste0("PCoA2",y_label,"%"))+stat_ellipse(data=df.plot,alpha=0.3)
ggsave('pcoa_plot.png', p)