微生信生物新年放大招：一条代码完成堆叠柱状图-冲击图的操作-终结版

00. 写在前面

眼下2019年的余额不足一天了，2020年最终还是要到来，首先祝愿大家元旦快乐，其次新的一年祝大家心想事成。说实在的最近老是想出去玩，哈哈哈，文章还没写完，不可以！不可以！

这只是大招的前奏，我们扩增子一体化分析exe系统已经包括这条功能了：敬请期待：微生信生物-扩增子V0.1。

  • 堆叠柱状图用于表征高维数据，例如微生物群落，代谢物组成是一种常见的图表类型，所以使用的次数也非常多，理所应当被封装成函数，一键绘制。

  • 冲击图其实添加了堆叠柱状图之间的连线，这对于柱子之间的横向比较会更加容易一些，同堆叠柱状图互为替代。

  • 堆叠柱状图在R语言中实现似乎不难，但是有许多的细节需要把握：整条柱子的子柱块排序问题，是否在柱子上添加标签，过滤低丰度柱子等在我们这条代码中都会解决。甚至前一段时间大家需要在堆叠柱状图上添加误差线，也都在这条代码里的。

01. 一条函数和参数解释


result = barplot(otu = "./otutab.txt",tax = "./taxonomy.txt",map = "./metadata.tsv",j = "Phylum",k= 0.01,rep = 6,axis_ord = "KO-OE-WT",label = FALSE,sd = TRUE)

  • otu = "./otutab.txt" ：OTU表格

  • tax = "./taxonomy.txt" ：物种注释文件

  • map = "./metadata.tsv"：样本分组文件

  • j = "Phylum" ：设置按照门为分类水平合并OTU，当然可以选择不同的分类水平；j = "Class"；j = "Order"；j = "Family"；j = "Genus"；甚至是OTU，如果你不嫌图例多的话。

  • k= 0.01：按照丰度过滤掉第丰度的物种，减少图例信息。

  • rep = 6 ：测定每个处理重复数量。

  • axis_ord = "KO-OE-WT" ：柱子在x轴额排布顺序。

  • label = FALSE：是否在柱子上添加标签。

  • sd = TRUE :是否对堆叠柱状图添加标准

02. 处理本地输出图片外我们还可以通过result提取图片

我设置函数返回图形结果的目的是方便大家根据自己需求修改颜色等图形属性。

#提取堆叠柱状图
p = result[[1]]
p
#提取冲击图
p = result[[2]]
p

3. 附：主函数

注意一下几个包，大家都看看装上没有，务必都安装后在运行，方可成功。

library(phyloseq)
library(tidyverse)
library(vegan)
library(reshape2)
library("plyr")
library(ggalluvial)
library(ggplot2)
barplot = function(otu = "./otutab.txt",tax = "./taxonomy.txt",map = "./metadata.tsv",j = "Phylum",k= 0.01,rep = 6,axis_ord = "KO-OE-WT",label = TRUE){
library(phyloseq)
library(tidyverse)
library(vegan)
library(reshape2)
library("plyr")
library(ggalluvial)
library(ggplot2)
axis_order = strsplit(basename(axis_order), "-")[[1]]
#导入otu表格
otu = read.delim(otu,row.names = 1)
head(otu)
otu = as.matrix(otu)
str(otu)
#导入注释文件
tax = read.delim(tax,row.names = 1)
head(tax)
tax = as.matrix(tax)
# taxa_names(tax)

#导入分组文件
map = read.delim(map,row.names = 1)
head(map)

# #导入进化树
# tree = read.tree("./otus.tree")
# tree

ps <- phyloseq(otu_table(otu, taxa_are_rows=TRUE),
sample_data(map) ,
tax_table(tax)
# phy_tree(tree)

)
ps

ps1 = ps
i = ps1

colnames(tax_table(ps1))
##这里我们过滤一定阈值的otu，会出现最后堆叠柱状图总体丰度高于100%的情况，这是合理的
###########绘制不同分类等级的柱状图
Taxonomies <- i %>%
tax_glom(taxrank = j) %>% # agglomerate at Class level Class
transform_sample_counts(function(x) {x/sum(x)} )%>%# Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance >= k) %>% # Filter out low abundance taxa
arrange(Phylum)

# head(Taxonomies)
# dim(Taxonomies)
colbar <- dim(unique(select(Taxonomies, one_of(j))))[1]

colors = colorRampPalette(c( "#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(colbar)

# mi = colorRampPalette(c( "#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD",
# "#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
# "#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(colbar)
#
#

Taxonomies$Abundance = Taxonomies$Abundance * 100
Taxonomies$Abundance = Taxonomies$Abundance/rep
# head(Taxonomies)

#按照分组求均值
colnames(Taxonomies) <- gsub(j,"aa",colnames(Taxonomies))
by_cyl <- group_by(Taxonomies, Group,aa)
zhnagxu2 = dplyr :: summarise(by_cyl, sum(Abundance))
#colnames(zhnagxu2) = c("group", j,"Abundance")
# head(zhnagxu2)

##确定因子，这里通过求和按照从小到大的顺序得到因子
##长变宽

# head(Taxonomies)

Taxonomies2 = dcast(Taxonomies,aa ~ Sample,value.var = "Abundance")
head(Taxonomies2)
Taxonomies2[is.na(Taxonomies2)] <- 0
aa = Taxonomies2
# head(aa)

n = ncol(aa)
#增加一行，为整列的均值，计算每一列的均值，2就是表示列
aa[n+1]=apply(aa[,c(2:ncol(aa))],1,sum)
colnames(aa)[n+1] <- c("allsum")
# str(aa)
bb<- arrange(aa, allsum)
# head(bb)
bb = bb[c(1,ncol(bb))]
cc<- arrange(bb, desc(allsum))
# head(cc)
##使用这个属的因子对下面数据进行排序

head(zhnagxu2)
colnames(zhnagxu2) <- c("group","aa","Abundance")
zhnagxu2$aa = factor(zhnagxu2$aa,order = T,levels = cc$aa)
zhnagxu3 = plyr::arrange(zhnagxu2,desc(aa))
# head(zhnagxu3)
##制作标签坐标，标签位于顶端
# Taxonomies_x = ddply(zhnagxu3,"group", transform, label_y = cumsum(Abundance))
# head(Taxonomies_x )
#标签位于中部
Taxonomies_x = ddply(zhnagxu3,"group", transform, label_y = cumsum(Abundance) - 0.5*Abundance)
head(Taxonomies_x,20 )
Taxonomies_x$label = Taxonomies_x$aa
#使用循环将堆叠柱状图柱子比较窄的别写标签，仅仅宽柱子写上标签
for(i in 1:nrow(Taxonomies_x)){
if(Taxonomies_x[i,3] > 3){
Taxonomies_x[i,5] = Taxonomies_x[i,5]
}else{
Taxonomies_x[i,5] = NA
}
}

##普通柱状图
p4 <- ggplot(Taxonomies_x , aes(x = group, y = Abundance, fill = aa, order = aa)) +
geom_bar(stat = "identity",width = 0.5,color = "black") +
scale_fill_manual(values = colors) +
theme(axis.title.x = element_blank()) +
theme(legend.text=element_text(size=6)) +
scale_y_continuous(name = "Abundance (%)")+
scale_x_discrete(limits = axis_order)

if (label == TRUE) {
p4 = p4 +
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")
}

# print(p4)

# install.packages("ggalluvial")
p4 = p4+theme_bw()+
scale_y_continuous(expand = c(0,0))+

theme(

panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
text=element_text(face = "bold"),
plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14,),
axis.text.y = element_text(colour = "black",size = 14),

legend.text = element_text(size = 15)
#legend.position = "none"#是否删除图例

)
p4

FileName1 <- paste("./a2_",j,k,"_bar",".pdf", sep = "")

ggsave(FileName1, p4, width = 12, height =8 )

##柱状图冲击图
#stratum定义堆叠柱状图柱子内容，以weight定义柱子长度，alluvium定义连线
head(Taxonomies_x )
cs = Taxonomies_x $aa
# head(cs)
# as.factor(Taxonomies_x $Genus)
# cs = as.character(Taxonomies_x $Genus)
# cs1 = as.factor(cs)
cs1 = cs
#提取真正的因子的数量
lengthfactor = length(levels(cs1))
#提取每个因子对应的数量
cs3 = summary (as.factor(cs1))
cs4 = as.data.frame(cs3)
cs4$id = row.names(cs4)
#对因子进行排序
df_arrange<- arrange(cs4, id)
#对Taxonomies_x 对应的列进行排序
Taxonomies_x1<- arrange(Taxonomies_x , aa)
head(Taxonomies_x1)
#构建flow的映射列Taxonomies_x
Taxonomies_x1$ID = factor(rep(c(1:lengthfactor), cs4$cs3))

#colour = "black",size = 2,,aes(color = "black",size = 0.8)

p3 = ggplot(Taxonomies_x1,
aes(x = group, stratum = aa, alluvium = ID,
weight = Abundance,
fill = aa, label = aa)) +
geom_flow(stat = "alluvium", lode.guidance = "rightleft",
color = "black",size = 0.2,width = 0.3,alpha = .2) +
geom_bar(width = 0.45)+
geom_stratum(width = 0.45,size = 0.2) +
#geom_text(stat = "stratum", size = 3,family="Times New Roman",fontface = "bold.italic") +
#theme(legend.position = "none") +
scale_fill_manual(values = colors)+
#ggtitle("fow_plot")+
scale_x_discrete(limits = axis_order)+
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")+
labs(x="group",
y="Relative abundancce (%)",
title="")
# p3
if (label == TRUE) {
p3 = p3 +
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")
}

p3 =p3+theme_bw()+
scale_y_continuous(expand = c(0,0))+
#geom_hline(aes(yintercept=0), colour="black", linetype=2) +
#geom_vline(aes(xintercept=0), colour="black", linetype="dashed") +
#scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+
theme(

panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
text=element_text(face = "bold"),
plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14),
axis.text.y = element_text(colour = "black",size = 14),

legend.text = element_text(size = 15,face = "bold.italic")
#legend.position = "none"#是否删除图例

)
p3

FileName2 <- paste("./a2_",j,k,"_bar_flow",".pdf", sep = "")

ggsave(FileName2, p3, width = 12, height =8)

return(list(p4,p3))

}

赞 (0)

相关推荐

  • NMDS分析

    " No one knows everything, and you don't have to."   --free傻孩子 "R数据分析"专题·第15篇   ...

  • ggplot2绘图学习 径向柱形图

    径向柱形图也被称为圆形柱形图或星图.这种图表使用同心圆网格来绘制条形图 每个圆圈表示一个数值刻度,而径向分隔线(从中心延伸出来的线)则用作区分不同类别或间隔(如果是直方图).刻度上较低的数值通常由中心 ...

  • 使用R语言的20条建议-微生信生物博主五年经验总结

    写在前面 如果说有什么理念或者习惯支撑在这几年的R语言学习中的话,我认为是这几条,如果大家将这几条能够理解大半,相信最起码会节省时间,提高效率. 注:这些建议不一定都会很好用,大家挑选适合自己的融会贯 ...

  • 微生信生物-扩增子结题报告（南京美他生物科技有限公司）

    微生物组学分析报告(南京美他生物科技有限公司) result_and_plot/Base_diversity_16s(ITS) 这部分存储微生物组部分多样性分析结果: 这部分结果使用OTU或者其他分类 ...

  • 微生信生物&根际互作生物学实验室年终总结

    阅读数量 这一年常读用户数量变化 微生信生物的前行和荆棘 2020年不平凡,有很多事情没有做,有很多事情却做了,有很多人走了,有很多人却来了,新冠疫情就像是将有限的资源进行了破碎又缝合,我认为这是新鲜 ...

  • 微生信生物历史推文集合 （持续更新）

    微生信生物历史推文集合 技巧经验思考资源 Rstudio切换挂载R版本及本地安装一些包 pubmed凉了,我们这里依然很美 ggplot版钢铁侠 当科研遇见python 学习R语言&生物信息不 ...

  • 无代码福音-微生信生物又要持续发力origin绘非典型柱状图

    上一期结束的时候留了个小问题: 一.前情回顾 首先,数据还是要分组的,因为如果放一列就是一组,最后还得一个一个改(Ctrl+鼠左双击),很麻烦. 那数据的B/C/D列一起作图会是怎样 为什么会这样,上 ...

  • 微生信生物---年中纪--2020

    2020年中纪I 抱歉占用大家一整个版面写下这个纪,毕竟什么都有个开始,有个结束,有的东西结束了,有的东西今天要开始,在此记录,铭记于心,方得始终. 2020年不平凡,经历了新冠肺炎后我们都很珍惜生命 ...

  • 0代码教程来了-来自-微生信生物-的零水平Origin制图

    写在前面 使用代码出图,R语言是最为广泛的,并且漂亮的出图和连带的分析让我们确实是受益良多.但是许多小伙伴,相信有不少人,都是没有足够的时间学习代码,因为大量的科研问题足够让我们头疼.因此,大家来看看 ...

  • 微生信生物科研爬虫项目等你来

    写在前面 微生信生物主编,最近大量任务缠身,我(五谷杂粮) 已经辅助运营了一段时间了,我们的好朋友抱起大块块毕业了.有一个月的空余时间,所以为他写了这篇推送. 前言 面对目前科研任务的多样化,越来越多 ...

  • 微生信生物学习进行时

    不知不觉,微生信生物已经陪伴我们三年了,似乎R语言和扩增子这两个主题占据了大部分的篇幅,这让这是最早涉及的内容呢?随着技术发展和研究内容改变,单一组学无法满足日益严峻的科研形式,所以我们再后来准备了代 ...