查看原文
其他

功能预测之Tax4Fun

2018-01-04 郑伟 宏基因组

之前在公众号的文章《根据16S预测微生物群落功能最全攻略》阅读人数近3000人,有需求的用户还是非常多的。其中提到了4个软件,之前已经介绍了其中非常有特点的三种,分别为:

Tax4Fun简介

Tax4Fun 2015年发表于Bioinformatics,己被引用68次。

K.P. Aßhauer, B. Wemheuer, R. Daniel, P. Meinicke (2015) Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data, Bioinformatics (2015) 31 (17): 2882-2884. doi:10.1093/bioinformatics/btv287.

今天为大家带来Tax4Fun的教程。相比于picrust来说,Tax4Fun最大的优点在于可以给予更新速度极快的SILVA数据库进行功能基因的预测。废话不多说,直接上干货。

分析之前需要下载Tax4Fun包及比对库,此文用的最新的 SILVA 123 在Tax4Fun官网上下载 http://tax4fun.gobics.de/

Tax4Fun支持QIIME和SILVAngs的结果

我们需要下载R包,细菌宏基因组数据库(SILVA Reference data),和QIIME格式的SILVA物种注释数据库(SILVA database for usage with QIIME),当然帮助文档Readme也是必读的。


图1. 重要文件下载位置

安装

以Linux为例,Windows类似,它只是个R包。建议在Rstudio中使用。

cd ~/test/example_PE250/tax4fun # 下载帮助文档 wget http://tax4fun.gobics.de/RPackage/Readme_Tax4Fun.pdf # R包安装 # Linux源码版 wget http://tax4fun.gobics.de/Tax4Fun/Tax4Fun_0.3.1.tar.gz # Windows二进级版 wget http://tax4fun.gobics.de/Tax4Fun/Tax4Fun_0.3.1.zip # 安装依赖关系 install.packages("qiimer") install.packages("Matrix") install.packages("biom") # 安装手动下载包 install.packages("~/test/example_PE250/tax4fun/Tax4Fun_0.3.1.tar.gz", repos = NULL, type = "source") # 数据库下载 wget http://tax4fun.gobics.de/Tax4Fun/ReferenceData/SILVA123.zip wget http://tax4fun.gobics.de/SilvaSSURef_123_NR.zip unzip SilvaSSURef_123_NR.zip # 此文件中有fasta和tax,但格式没有注释且长短不一,值得注意 unzip SILVA123.zip # 此为作者制作的功能KEEG注释

以官方示例数据演示使用

下载官网提供的测试数据:基于HMP测序数据用QIIME生成的OTU表

wget http://tax4fun.gobics.de/QIIME/HMP_0.97_table.txt

格式示例:

# Constructed from biom file #OTU ID buccal.mucosa.SRS015921 stool.SRS016203 denovo0 2.0     0.0     0.0     0.0     0.0     denovo1 0.0     1.0     0.0     0.0     0.0     denovo2 0.0     0.0     1.0     0.0     0.0     denovo3 0.0     0.0     1.0     3.0     0.0

1. 导入OTU表

函数描述
importQIIMEBiomData导入1个或多个Biom格式的OTU表,如QIIME输出文件
importQIIMEData导入文本制表符分割格式OTU表
importSilvaNgsData导入SilvaNgs在线生成的结果,逗号分隔格式

我们以QIIME生成的文本OTU表为例演示

# 导入OTUs表 library(Tax4Fun) setwd("~/test/example_PE250/tax4fun") QIIMESingleData <- importQIIMEData("HMP_0.97_table.txt") # 可选步骤 # 发现OTU表为按Taxonomy合并的结果 otu_table = QIIMESingleData$otuTable # 按列求合,发现为每个样品的原始reads colSums(otu_table) # 输出转换后的OTU表,Tax4Fun真实需要的注释格式是这样的;输出自己看看, write.table("ID\t", file="otu_table_tax.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F) write.table(otu_table, file="otu_table_tax.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE) # 警告信息可忽略

我们可以看到importQIIMEData()函数将OTUs表进行了合并,结果如下:

ID      buccal.mucosa.SRS015921 stool.SRS016203 stool.SRS015133 stool.SRS013687 stool.SRS011586 stool.SRS015190 stool.SRS016095 Archaea;Euryarchaeota;Methanobacteria;Methanobacteriales;Methanobacteriaceae;Methanobrevibacter;Methanobrevibacter smithii     Archaea;Thaumarchaeota;Miscellaneous Crenarchaeotic Group;      0       1       0       0       0       0       0       0       Bacteria;Actinobacteria;Acidimicrobiia;Acidimicrobiales;OCS155 marine group;    0       0       0       0       0       0

2. 预测KO功能表

函数Tax4Fun使用格式
Tax4Fun(Tax4FunInput, folderReferenceData, fctProfiling = TRUE,
 refProfile = “UProC”, shortReadMode = TRUE, normCopyNo = TRUE)

参数描述
Tax4FunInput导入的OTU表变量
folderReferenceData数据库位置
fctProfiling默认计算功能谱,FLASE时计算代谢谱
refProfile功能功能谱的计算方式,默认UProC,可选PAUDA
shortReadMode默认基于100bp reads,FALSE时为400bp
normCopyNo校正16S基因拷贝数
# 根据Tax4Fun提供的SILVA123最新数据库进行预测,要求此数据的压缩包拉于此工作目录,这个命令得出来的是KO号的各种酶的基因丰度 Tax4FunOutput <- Tax4Fun(QIIMESingleData, "SILVA123", fctProfiling = TRUE, refProfile = "UProC", shortReadMode = TRUE, normCopyNo = TRUE) # 提取KO表,生成6508个KO相关的通路 KO_table = t(Tax4FunOutput$Tax4FunProfile) # 所有样品标准化为1, 没有原始数据,可以用anova和Limma,无法使用edgeR和DESeq2 colSums(KO_table) # 输出KO表,表头写个制表符,用于对齐表头 write.table("ID\t", file="KO_table.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F) write.table(KO_table, file="KO_table.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)

结果为标准化的KO表,如下图:

ID      buccal.mucosa.SRS015921 stool.SRS016203 stool.SRS015133 stool.SRS013687 stool.SRS011586 stool.SRS015190 stool. K00001; alcohol dehydrogenase [EC:1.1.1.1]      0.000510048165023247    0.000308143029580682    0.000401028969794961   K00002; alcohol dehydrogenase (NADP+) [EC:1.1.1.2]      1.56218072210004e-05    4.09601550909485e-05    2.693692447672 K00003; homoserine dehydrogenase [EC:1.1.1.3]   0.000700178764510118    0.000482493099204403    0.000670762594447745

3. 预测代谢部分

上面的KO表是有几千条的,只关注代谢物,做如下分析,只有几百条结果。

# 上步结果中KO名称,只关注代谢通路的变化,调整fctProfiling=FALSE Tax4FunOutput <- Tax4Fun(QIIMESingleData, "SILVA123", fctProfiling = FALSE, refProfile = "UProC", shortReadMode = TRUE, normCopyNo = TRUE) # 表头写个制表符,用于对齐表头 write.table("ID\t", file="KO_table_fct.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F) write.table(KO_table, file="KO_table_fct.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)

此表只有几百行,方便下游分析,示例如下:

ID      buccal.mucosa.SRS015921 stool.SRS016203 stool.SRS015133 stool.SRS013687 stool.SRS011586 stool.SRS015190 ko00010; Glycolysis / Gluconeogenesis   0.0102626891498075      0.00863427711585646     0.00807436361246528     ko00020; Citrate cycle (TCA cycle)      0.00707255408716388     0.00457664886228007     0.00382738314647692     ko00030; Pentose phosphate pathway      0.00977201885189811     0.00743298942053779     0.00675236272234106

如何制作符何tax4fun注释要求的OTU表

基于之前发布的测试数据,你也可以用你的OTUs表和代表性序列。测试数据位于百度网盘,链接:http://pan.baidu.com/s/1hs1PXcw 密码:y33d。请在共享目录中查找。

# 制作符合tax4fun要求的OTU表 # 按tax4fun silva123 blast注释物种 assign_taxonomy.py -i result/rep_seqs.fa \    -r tax4fun/SILVA_123_SSURef_Nr99_tax_silva.fasta \    -t tax4fun/SILVA_123_SSURef_Nr99_tax_silva.tax \    -m blast -o tax4fun/ # 添加taxonomy至OTUs表的biom格式中 biom add-metadata -i result/otu_table.biom \    --observation-metadata-fp tax4fun/rep_seqs_tax_assignments.txt \    -o tax4fun/otu_table_tax.biom \    --sc-separated taxonomy --observation-header OTUID,taxonomy # 转换为文本格式 biom convert -i tax4fun/otu_table_tax.biom -o tax4fun/otu_table_tax.txt --table-type="OTU table" --to-tsv  --header-key taxonomy

其实有了按要求的OTUs表,只需两条命令即可完成分析,其它代码都是数据表提取、统计、校验等,不是必不可少,但是很有意义的(你考试答完题不检查吗?)。

## 使用我们的数据 QIIMESingleData <- importQIIMEData("otu_table_tax.txt") Tax4FunOutput <- Tax4Fun(QIIMESingleData, "SILVA123", fctProfiling = FALSE, refProfile = "UProC", shortReadMode = TRUE, normCopyNo = TRUE)

基于格式注释结果的分析

比如基于greengene注释的物种,或silva采用qiime标准格式

result/rep_seqs_tax_assignments.txt

OTU_325 k__Bacteria;p__Bacteroidetes;c__Flavobacteriia;o__Flavobacteriales;f__Cryomorphaceae;g__;s__    0.880 OTU_324 k__Bacteria;p__Chlorobi;c__SJA-28;o__;f__;g__;s__       1.000 OTU_327 k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Legionellales;f__Coxiellaceae;g__;s__   0.810

我们对之前发布的../result/otu_table_tax.txt文件进行分析

## 使用greengene注释的结果 QIIMESingleData <- importQIIMEData("../result/otu_table_tax.txt") # 转换格式为tax4fun要求 otudf<-QIIMESingleData$otuTable # 把上图中的otuTable提取出来命名为otudf rownames(otudf)<-gsub("[a-z]__","",rownames(otudf)) # 去掉名字中的特殊符号 QIIMESingleData$otuTable<-otudf # 重新把修改好后的otudf 命名为QIIMESingleData 中的 otuTable Tax4FunOutput <- Tax4Fun(QIIMESingleData, "SILVA123", fctProfiling = FALSE, refProfile = "UProC", shortReadMode = TRUE, normCopyNo = TRUE) # 预测代谢物

KO表结果示例:

图2. 代谢表结果示例

Tax4FunOutput <- Tax4Fun(QIIMESingleData, "SILVA123", fctProfiling = TRUE, refProfile = "UProC", shortReadMode = TRUE, normCopyNo = TRUE) # 根据SILVA123最新数据库进行预测,这个命令得出来的是KO号的各种酶的基因丰度

图3. KO表结果示例


有了结果以后大家用write.table输出一下就可以做代谢分析了,是不是很简单,更多功能可以看Tax4Fun包的帮助。


猜你喜欢

写在后面

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


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

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存