FAPROTAX预测16s扩增子群落功能信息以及R语言可视化

FAPROTAX是一款适用于对环境样本(如海洋、湖泊等)的生物地球化学循环过程(特别是碳、氢、氮、磷、硫等元素循环)进行功能注释预测的软件。我们知道微生物与土壤养分的循环利用密切相关,许多养分被植物利用过程中离不开微生物的作用,比如,氮降解途径中的亚硫酸盐呼吸。

FAPROTAX原理

首先作者先根据文献资料手动构建了联系物种分类与功能注释的FAPROTAX数据库;后又编写了联系OTU分类表与FAPROTAX数据库的python脚本;最后,只要将基于16S的OTU分类表通过python脚本就可以输出微生物群落功能注释预测结果。
另外,PICRUSt输入的物种丰度表目前只能用Greengene数据库进行物种注释,在这一点上FAPROTAX更灵活,因为它识别的是菌的属、种的名称,只要你注释到的菌己被报导有这方面的功能,无论注释数据库是Greengene、Silva和RDP,OTU是有参考还是de novo的结果都可以识别。
有关原理部分和使用部分可参考宏基因组—如何用FAPROTAX预测微生物群落功能。网址:https://mp.weixin.qq.com/s/J8EwJD_PTDhqRaD7kXlK1A
但是出现一个问题上面这篇文章给出的FAPROTAX官网网址(https://www.zoology.ubc.ca/louca/FAPROTAX/lib/php/index.php?section=Instructions)似乎失效了,无法下载该软件。

2020.3.5

FAPROTAX官网的最新网址:https://pages.uoregon.edu/slouca/LoucaLab/archive/FAPROTAX/lib/php/index.php
可能需要翻墙才能访问,没法翻墙的小伙伴也别急,关注公众号:Amoy数据分析,回复FAPROTAX即可下载最新版本的软件。
FAPROTAX官网

可以看到关于FAPROTAX的介绍,原理,使用代码,如何下载选项,点击下载选项可以看最新版本为1.2.1,点击下载。

软件的安装和使用

将下载的压缩文件上传至服务器或虚拟机。这款软件所需要的环境为Python2.7版本。

#解压缩该文件
unzip FAPROTAX_1.2.1.zip
cd FAPROTAX_1.2.1 | ls 
#使用conda软件创建名为envpython27的python=2.7环境,-n为创建环境的名称。
 conda create -n envpython27 python=2.7 -y
#激活python=2.7环境
 conda activate envpython27
#测试软件能否正确运行,如果可以运行会出现参数信息,报错说明不可用,再根据报错的信息解决。
 python /home/ubuntu/FAPROTAX_1.2.1/collapse_table.py
#使用:-i otu分类信息表:-o 输出文件名称:-g 自带的数据库名称,在FAPROTAX_1.2.1/目录下:-r 中间过程文件,可以看到那些菌属没有被数据库对比上,其它具体使用请参考官网的教程。
  python /home/ubuntu/FAPROTAX_1.2.1/collapse_table.py -i qiime1_otu.txt -o qii1_otu_tax_faprtax.txt -g /home/ubuntu/FAPROTAX_1.2.1/FAPROTAX.txt -r report.txt -v -d 'taxonomy'

输入的otu分类表文件格式:案列中的qiime1_otu.txt文件


image.png

输出文件qii2_otu_tax_faprtax.txt


image.png

使用R语言进行可视化

以下操作在R语言中进行,一共需要三个文件:上一步的输出文件qii2_otu_tax_faprtax.txt;分组文件design.txt;需要展示的功能列表文件IND.list。后两个文件自己创建即可。

#加载包和主题,我们做图的最终目的是要发表的,统一的设置一个主题可以方便我们作图,图形更加美观
library(ggplot2)
main_theme = theme(panel.background=element_blank(), panel.grid=element_blank(),
                   axis.line.x=element_line(size=.5, colour="black"), axis.line.y=element_line(size=.5, colour="black"),
                   axis.ticks=element_line(color="black"), axis.text=element_text(color="black", size=),
                   legend.position="right", legend.background=element_blank(), legend.key=element_blank(), legend.text= 
 element_text(size=),
                   text=element_text(family="sans", size=))

alpha = read.table("qii2_otu_tax_faprtax.txt", header=T, row.names=1, sep="\t", comment.char="") 
alpha =alpha[,-1]#数据第一列是NA,所以要去除
alpha = t(alpha)
alpha = alpha/100
alpha =as.data.frame(alpha)
design = read.table("design.txt", header=T, row.names=1, sep="\t")#SampleID为样品编号与qii2_otu_tax_faprtax.txt第一行一致:subspecies为分组信息,对照组和处理组:solitype代表两个不同地点。
design$group=design$groupID
sub_design = subset(design, group %in% c("MSXC","MSXT","XWZC","XWZT"))
sub_design$group  = factor(sub_design$group, levels=c("MSXC","MSXT","XWZC","XWZT"))
idx = rownames(sub_design) %in% rownames(alpha)
sub_design=sub_design[idx,]
# detach("package:ggpubr", unload=TRUE)
# library(dplyr)

list = read.table("IND.list", header=F,  sep="\t")#我们需要挑出来展示的功能,我关注的点是元素循环这块,好像还有病原菌相关的,可根据自己的需要挑选
sub_alpha=alpha[rownames(sub_design), as.vector(list$V1)]
sub_alpha$sampleID=rownames(sub_alpha) 
melted = melt(sub_alpha, id.vars = "sampleID")
melted_all = merge(melted, sub_design[,c("subspecies","solitype")], by.x="sampleID", by.y = "row.names", all.x = T ) 
x = as.data.frame(melted_all %>% group_by(variable) %>%
                    summarise(mean = mean(value)))
x = x[order(x$mean, decreasing = T),]#排序
x = head(x, 10)#选取排名前10的进行展示
melted_all = melted_all[melted_all$variable %in% x$variable,]
melted_all$variable = factor(melted_all$variable, levels = as.vector(x$variable))
melted_all$soiltype = factor(melted_all$solitype, levels=c("M", "X"))
p = ggplot(melted_all, aes(x=variable, y=value, color=subspecies)) +
  geom_boxplot(position = "dodge", alpha=1, outlier.size=0.3, size=0.5, width=0.7, fill="transparent") +
  #geom_jitter( position=position_jitter(0.17), size=0.3, alpha=0.7)+
  main_theme +
  facet_grid(soiltype ~ .)+#以soiltype 图型分页
  theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))
p
ggsave(paste("faprotax_IND_boxplot_top9..", ".pdf", sep=""), p, width = 8, height = 5)

design.txt文件格式。


image.png

结果数据未发表,就使用刘永鑫大神他们做的图Fig3

除了使用柱状图展示,还可使用热图等其它方式展示。

致谢

公众号宏基因组———如何用FAPROTAX预测微生物群落功能
R语言代码来自文章———Zhang J , Liu Y X , Zhang N , et al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice[J]. Nature Biotechnology, 2019.中Fig3中的代码部分。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容