酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

最近因为需要“恰饭”的缘故,有一个星期没有给大家更新干货了,今天它终于来了!!

由于高通量测序的价格降的越来越低,现在很多人的研究已经从早期几个、十几个样品的研究发展到了几十、几百、甚至上千个样品,这种确实在以扩增子测序为基础的研究中越发明显。

样品数目多了很多的时候也就意味着分组数目的增加,在文章必备的beta多样性分析结果图中,很有可能会造成不同组样品之间的互相重叠,直观无法有效的看出分组聚类

今天给大家带来一个beta多样性分析与统计学检验结合的图像,整体的效果就想下面这幅图一样。

酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

在图中整体上使用PCoA来展示样品的beta多样性分析结果,之后分别在上方和右侧添加相须图来展示不同组样品在PC1和PC2轴的分布情况,并给出统计学检验结果,最后在右上角的位置给出PERMANOVA的检验结果。

这样大家就可以直观的看到分组是否对样品具有显著的影响,以及不同组间是否具有明显差异,样品之间的变化规律是什么样子的。

绘图过程

绘图前准备

首先载入分析和绘图所需要的R包。

library(vegan)library(ape)library(ggplot2)library(grid)library(dplyr)library(multcomp)library(patchwork)

本文使用R语言自带的iris数据进行示例展示,请根据下方代码运行得到的data和groups文件格式自行准备输入文件。

data("iris")
data <- iris[,1:4]
groups <- data.frame(rownames(iris),iris[,5])
colnames(groups) <- c("V1","V2")

定义一些绘图所需要用到的参数。

length=length(unique(as.character(groups$V1)))
times1=length%/%8res1=length%%8times2=length%/%5res2=length%%5col1=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绘图数据文件。

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分析之后运行。

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个图的完美匹配。

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完美解决打样本量多组数据不好区分的问题!!

PCoA结果图绘制

使用下方代码进行PCoA结果图绘制。

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数值,不然可能会出现重叠或者中间空很大的现象。

酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

PERMANOVA分析

使用下方代码对数据进行PERMANVOA分析并绘制图像。

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幅图拼在一起。

p5 <- p1 + p4 + p2 + p3 + 
plot_layout(heights = c(1,4),widths = c(4,1),ncol = 2,nrow = 2)

酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

结果导出

将绘制好的图像保存为pdf和png格式。

pdf("PCoA.pdf",height=12,width=15)
p5
png(filename="PCoA.png",res=600,height=7000,width=9000)
p5
dev.off()

酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!
我就知道你在
生物医学科研方法

Nat. Genet. | 基于泛癌单细胞RNA-seq,识别22种癌症类型的细胞异质性重复程序

2020-11-18 17:57:25

生物医学科研方法

不影响结论的一图多用是否不算学术造假?对于造假,导师/通讯作者该不该来背锅?

2020-11-18 18:05:27