物种丰度排序堆积柱形图及处理间各物种丰度非参数检验多组比较的R图形可视化

再美的可视化图形若缺少了统计检验就失去了灵魂而变得华而不实

测试数据及代码链接：https://pan.baidu.com/s/1LJzCNqJe9OcfUmhH9jfvsw
提取码：9ifa

因公众号文章不可修改，如以上链接失效，或想获取代码的更新版，请在“宏基因组”公众号后台回复本文关键字“stackplot”获取最新下载地址。

一般做16s扩增子测序，都少不了物种组成分析，但在各文献中同时将处理间各物种丰度进行统计检验R图形可视化操作的寥寥无几，若是文章中将物种组成图与处理间各物种丰度差异统计检验图同时展示将会是一道亮丽的风景。秉持着这种思考，在参考学习众多资料的情况下敲下了如下代码：

当然用的数据依然是刘永鑫老师的经典范例数据，方便需要的人共同学习进步，提出宝贵建议！

rm(list=ls())
library(tidyverse)#数据整理与数据转换包，用了一些更好用更易懂的函数
library(ggprism)
library(vegan)
otu <- read.delim('./otutab.txt',row.names = 1)
head(otu, n = 3)
#otu <-decostand(otu,'total',2)#按列标准化otu
#colSums(otu)#若想后面做成相对丰度的差异比较，可把这两行释放出来即可
tax <- read.delim('./taxonomy.txt',row.names = 1)
head(tax, n = 3)
metadata<- read.delim('./metadata.tsv')
head(metadata, n = 3)
dat <- merge(x=otu,y=tax,by='row.names')
head(dat, n = 3)
dat =dplyr::rename(dat,OTUID = Row.names)
head(dat, n = 3)
##按Phylum水平分组汇总(根据自己需求更改要展示的物种水平)
aa<-aggregate(dat[,2:ncol(otu)],by=list(dat$Phylum),FUN=sum)
head(aa)

########################三种排序方法，任选其一
#1
# aa<- mutate(aa,v=apply(aa[,c(2:ncol(aa))],1,sum))
# cc<- arrange(aa,desc(v))
# head(cc)
# cc<-select(cc,-v)
# head(cc)
# row.names(cc)=cc$Phylum
# head(cc)
# cc<-select(cc,-Phylum)
# head(cc)
#2
# row.names(aa)=aa$Phylum
# head(aa)
# aa<-select(aa,-Phylum)
# head(aa)
# cc<-aa[order(rowSums(aa),decreasing=T),]
#3
row.names(aa)=aa$Group.1
head(aa)
aa<-dplyr::select(aa,-Group.1)
head(aa, n = 3)
#根据行求和结果对数据排序
order<-sort(rowSums(aa[,2:ncol(aa)]),index.return=TRUE,decreasing=T)
#根据列求和结果对表格排序
cc<-aa[order$ix,]
head(cc, n = 3)
##只展示排名前10的物种，之后的算作Others(根据需求改数字)
dd<-rbind(colSums(cc[11:as.numeric(length(rownames(cc))),]),cc[10:1,])
head(dd, n = 3)
rownames(dd)[1]<-"Others"
head(dd, n = 3)
##再与metadata合并
bb<-merge(t(dd),dplyr::select(metadata,SampleID,Group),by.x = "row.names",by.y ="SampleID")
head(bb, n = 3)
##宽数据变长数据
kk<-tidyr::gather(bb,Phylum,Abundance,-c(Group,Row.names))
#将未注释到的Unassigned也改为Others,你也可以不改，有以下两种方式
kk$Phylum<-ifelse(kk$Phylum=='Unassigned','Others',kk$Phylum)#1
#kk[kk$Phylum=='Unassigned','Phylum']='Others'               #2
##根据Group,Phylum分组运算
hh <- kk %>%group_by(Group,Phylum) %>%dplyr :: summarise(Abundance=sum(Abundance))
head(hh, n = 3)

##更改因子向量的levels
hh$Phylum = factor(hh$Phylum,order = T,levels = row.names(dd))

yanse <-c("#999999","#F781BF","#A65628","#FFFF33","#FF7F00","#984EA3","#4DAF4A","#377EB8","#74D944","#E41A1C","#DA5724","#CE50CA", "#D3D93E","#C0717C","#CBD588","#D7C1B1","#5F7FC7","#673770", "#3F4921","#CD9BCD","#38333E","#689030","#AD6F3B")#要确保颜色够用，否则会报错
##按物种丰度排序好的堆积柱形图
p1 <- ggplot(hh,aes(x = Group,y = Abundance,fill = Phylum)) + geom_bar(position="fill",stat = "identity",width = 0.5) +scale_fill_manual(values = yanse) +labs(x='Group',y='Abundance(%)')+scale_x_discrete(limits = c("KO","OE","WT"))+guides(fill=guide_legend(reverse = TRUE))+ggprism::theme_prism()+scale_y_continuous(expand = c(0,0))
p1#由于把Unassigned也算成了Others，所以只显示9个物种

###进行处理间各物种非参数检验的多组比较
#数据整理与转换
head(bb,n = 3)
cc =dplyr::select(bb,Row.names,Group,everything(),-c(Others,Unassigned))
head(cc,n = 3)
dat <- gather(cc,Phylum,v,-(Row.names:Group))
head(dat,n = 3)

(参见之前推送，一共写了五种多重比较的方法，总有一款适合你)
##非参数检验的多组比较函数
PMCMR_compare1 <- function(data,group,compare,value){library(multcompView)library(multcomp)library(PMCMRplus)library(PMCMR)a <- data.frame(stringsAsFactors = F)type <- unique(data[,group])for (i in type){g1=comparesub_dat <- data[data[,group]==i,]names(sub_dat)[names(sub_dat)==compare] <- 'g1'names(sub_dat)[names(sub_dat)==value] <- 'value'sub_dat$g1 <- factor(sub_dat$g1)options(warn = -1)k <- PMCMRplus::kwAllPairsNemenyiTest(value ~ g1,data=sub_dat)n <- as.data.frame(k$p.value)h <- n %>%mutate(compare=rownames(n)) %>%gather(group,p,-compare,na.rm = TRUE) %>%unite(compare,group,col="G",sep="-")dif <- h$pnames(dif) <- h$GdifdifL <- multcompLetters(dif)K.labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)K.labels$compare = rownames(K.labels)K.labels$type <- imean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1")names(mean_sd) <- c('compare','std','mean')a <- rbind(a,merge(mean_sd,K.labels,by='compare'))}names(a) <- c(compare,'std','mean','Letters',group)return(a)
}

##################################用函数运行输入的数据
df2 <- PMCMR_compare1(dat,'Phylum','Group','v')
df2########字母标正着标(a>b>c)(之后跑自己的数据可能出现相反的顺序，不影响结果，可在AI或PDF编辑器中调)

p2 = ggplot(dat,aes(Group,v))+geom_boxplot(aes(color=Group))+geom_text(data=df2,aes(x=Group,y=mean+2*std,label=Letters))+geom_jitter(aes(fill=Group),position = position_jitter(0.2),shape=21,size=1,color="black")+facet_wrap(~Phylum,scales = "free_y")+ labs(x='Group',y='ASVs')+ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p2

Output figure width and height
Letter纸图片尺寸为单栏89 mm，双栏183 mm，页面最宽为247 mm，推荐比例16：10，即半版89 mm x 56 mm; 全版183 mm x 114 mm

##################保存图片ggsave("./p1.pdf", p1, width = 350, height = 200, units = "mm")ggsave("./p2.pdf", p2, width = 350, height = 200, units = "mm")

参考资料

浅析R语言单因素方差分析中的多重比较

扩增子分析 | 一键式抽平和物种组成分析（1)

《R语言数据可视化之美专业图表绘制指南》

《ggplot2:数据分析与图形艺术第二版》

《R数据科学》

作者：旭日阳光、flyfly

责编：马腾飞 南京农业大学

审核：刘永鑫 中科院遗传发育所

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发NatureCell专刊肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命生命大跃进  细胞暗战 人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读