美文网首页扩增子分析
FAPROTAX预测16s扩增子群落功能信息以及R语言可视化

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

作者: 小三的后一位小四 | 来源:发表于2020-03-23 15:08 被阅读0次

    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中的代码部分。

    相关文章

      网友评论

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

        本文链接:https://www.haomeiwen.com/subject/indvrhtx.html