woodcorpse的个人博客分享 http://blog.sciencenet.cn/u/woodcorpse

博文

物种丰度排序堆积柱形图及处理间各物种差异分析

已有 3164 次阅读 2021-10-10 14:24 |个人分类:R|系统分类:科研笔记

物种丰度排序堆积柱形图及处理间各物种丰度非参数检验多组比较的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个物种

image

###进行处理间各物种非参数检验的多组比较
#数据整理与转换
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=compare
    sub_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$p
    names(dif) <- h$G
    dif
    difL <- multcompLetters(dif)
    K.labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)
    K.labels$compare = rownames(K.labels)
    K.labels$type <- i

    mean_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

image

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

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

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

猜你喜欢

写在后面

为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外5000+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
image

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

点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA



https://m.sciencenet.cn/blog-3334560-1307437.html

上一篇:RDA_环境因子_群落结构_统计检验_可视化
下一篇:南方科技大学唐圆圆组招聘环境相关领域科研人才(年薪33~50万)

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 21:12

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部