由于高通量测序的价格降的越来越低，现在很多人的研究已经从早期几个、十几个样品的研究发展到了几十、几百、甚至上千个样品，这种确实在以扩增子测序为基础的研究中越发明显。
样品数目多了很多的时候也就意味着分组数目的增加，在文章必备的beta多样性分析结果图中，很有可能会造成不同组样品之间的互相重叠，直观无法有效的看出分组聚类。
今天给大家带来一个beta多样性分析与统计学检验结合的图像，整体的效果就想下面这幅图一样。
在图中整体上使用PCoA来展示样品的beta多样性分析结果，之后分别在上方和右侧添加相须图来展示不同组样品在PC1和PC2轴的分布情况，并给出统计学检验结果，最后在右上角的位置给出PERMANOVA的检验结果。
这样大家就可以直观的看到分组是否对样品具有显著的影响，以及不同组间是否具有明显差异，样品之间的变化规律是什么样子的。
绘图过程
绘图前准备
首先载入分析和绘图所需要的R包。
代码语言：javascript复制library(vegan)library(ape)library(ggplot2)library(grid)library(dplyr)library(multcomp)library(patchwork)本文使用R语言自带的iris数据进行示例展示，请根据下方代码运行得到的data和groups文件格式自行准备输入文件。
代码语言：javascript复制data("iris")
data <- iris[,1:4]
groups <- data.frame(rownames(iris),iris[,5])
colnames(groups) <- c("V1","V2")定义一些绘图所需要用到的参数。
代码语言：javascript复制length=length(unique(as.character(groups$V1)))
times1=length%/%8res1=length%%8times2=length%/%5res2=length%ol1=rep(1:8,times1)
col=c(col1,1:res1)
pich1=rep(c(21:24),times2)
pich=c(pich1,15:(15 res2))
cbbPalette <- c("#B2182B","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5")PCoA分析
使用下方代码进行PCoA分析并建立PCoA绘图数据文件。
代码语言：javascript复制data <- vegdist(data)
pcoa<- pcoa(data, correction = "none", rn = NULL)
PC1 = pcoa$vectors[,1]
PC2 = pcoa$vectors[,2]
plotdata <- data.frame(rownames(pcoa$vectors),PC1,PC2,groups$V2)
colnames(plotdata) <-c("sample","PC1","PC2","Group")
pc1 <-floor(pcoa$values$Relative_eig[1]*100)
pc2 <-floor(pcoa$values$Relative_eig[2]*100)
plotdata$Group <- factor(plotdata$Group,levels = c("setosa","versicolor","virginica"))⚠️正常的微生物丰度数据应该是列是样本、行是物种，分析需要的数据是行是样本、列是物种，所以有些数据可能需要转置一下。
⚠️另外就是最后一行需要根据你自己的分组名自定义排列顺序，关系到分组在图中显示的先后顺序。
PC1和PC2的显著性检验
使用下方代码分别对上一步得到的PCoA结果中PC1和PC2进行组间差异检验，所以这一步一定要在PCoA分析之后运行。
代码语言：javascript复制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.1yd2$Max <- yd2$Max max(yd2$Max)*0.1fit1 <- 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,levels = c("setosa","versicolor","virginica"))⚠️这里同样在自定义一下分组的顺序，顺序与上一步保持一致。
相须图绘制
特别强调，一定要先画上方和右侧的相须图！！！
这里有一个细节，就是因为相须图是添加了差异检验字母的，就会导致相须图和PCoA散点图的坐标轴范围不一致，如果直接合并的话会导致图像扭曲，箱子无法准确对应PCoA中点的分布。
所以一定要先画相须图，然后在后面PCoA图的绘制过程中调用两个相须图的坐标轴范围，以达到4个图的完美匹配。
代码语言：javascript复制p1 <- ggplot(plotdata,aes(Group,PC1))
geom_boxplot(aes(fill = Group))
geom_text(data = test,aes(x = Group,y = yd1,label = PC1),
size = 7,color = "black",fontface = "bold")
coord_flip()
scale_fill_manual(values=cbbPalette)
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.y=element_text(colour='black',size=20,face = "bold"),
axis.text.x=element_blank(),
legend.position = "none")
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")
scale_fill_manual(values=cbbPalette)
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=20,angle = 45,
vjust = 1,hjust = 1,face = "bold"),
axis.text.y=element_blank(),
legend.position = "none")⚠️由于每个人分组名称的长短不一样，为了保持图像美观，在出图之后可能需要根据结果返回来调整一下相须图中分组名的字号大小。
PCoA结果图绘制
使用下方代码进行PCoA结果图绘制。
代码语言：javascript复制p2<-ggplot(plotdata, aes(PC1, PC2))
geom_point(aes(fill=Group),size=8,pch = 21)
scale_fill_manual(values=cbbPalette,name = "Group")
xlab(paste("PC1 ( ",pc1,"%"," )",sep=""))
ylab(paste("PC2 ( ",pc2,"%"," )",sep=""))
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=30))
geom_vline(aes(xintercept = 0),linetype="dotted")
geom_hline(aes(yintercept = 0),linetype="dotted")
theme(panel.background = element_rect(fill='white', colour='black'),
panel.grid=element_blank(),
axis.title = element_text(color='black',size=34),
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=34,vjust = 7),
axis.title.y=element_text(colour='black', size=34,vjust = -2),
axis.text=element_text(colour='black',size=28),
legend.title=element_text(size = 24,face = "bold"),
legend.text=element_text(size=20),
legend.key=element_blank(),legend.position = c(0.88,0.13),
legend.background = element_rect(colour = "black"),
legend.key.height=unit(1,"cm"))
guides(fill = guide_legend(ncol = 1))⚠️这里已经帮你自动匹配的相须图的坐标轴范围。
⚠️由于图例是放在图像内部的，所以需要根据出图结果调整图例的位置、行列数目、标题、字号大小等等，这样才能不遮挡图中的点。
⚠️由于分组名的长度可能会发生变化，需要根据出图结果自行调整横纵坐标轴标题的vjust数值，不然可能会出现重叠或者中间空很大的现象。
PERMANOVA分析
使用下方代码对数据进行PERMANVOA分析并绘制图像。
代码语言：javascript复制otu.adonis=adonis(data~V2,data = groups,distance = "bray")
p4 <- ggplot(plotdata, aes(PC1, PC2))
geom_text(aes(x = -0.5,y = 0.6,label = paste("PERMANOVA:ndf = ",otu.adonis$aov.tab$Df[1], "nR2 = ",round(otu.adonis$aov.tab$R2[1],4), "np-value = ",otu.adonis$aov.tab$`Pr(>F)`[1],sep = "")),
size = 7)
theme_bw()
xlab("") ylab("")
theme(panel.grid=element_blank(),
axis.title = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank())图像拼接
使用patchwork包将4幅图拼在一起。
代码语言：javascript复制p5 <- p1 p4 p2 p3
plot_layout(heights = c(1,4),widths = c(4,1),ncol = 2,nrow = 2)结果导出
将绘制好的图像保存为pdf和png格式。
代码语言：javascript复制pdf("PCoA.pdf",height=12,width=15)
p5
png(filename="PCoA.png",res=600,height=7000,width=9000)
p5
dev.off()
