首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >酷不酷炫!想不想学!带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!

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

作者头像
DataCharm
发布2021-02-22 12:04:07
1.2K0
发布2021-02-22 12:04:07
举报

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

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

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

在图中整体上使用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结果图绘制。

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分析并绘制图像。

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)

结果导出

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

pdf("PCoA.pdf",height=12,width=15)
p5
png(filename="PCoA.png",res=600,height=7000,width=9000)
p5
dev.off()
本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-11-14,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 DataCharm 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 绘图过程
    • 绘图前准备
      • PCoA分析
        • PC1和PC2的显著性检验
          • 相须图绘制
            • PCoA结果图绘制
              • PERMANOVA分析
                • 图像拼接
                  • 结果导出
                  相关产品与服务
                  图数据库 KonisGraph
                  图数据库 KonisGraph(TencentDB for KonisGraph)是一种云端图数据库服务,基于腾讯在海量图数据上的实践经验,提供一站式海量图数据存储、管理、实时查询、计算、可视化分析能力;KonisGraph 支持属性图模型和 TinkerPop Gremlin 查询语言,能够帮助用户快速完成对图数据的建模、查询和可视化分析。
                  领券
                  问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档