功能预测之Tax4Fun
之前在公众号的文章《根据16S预测微生物群落功能最全攻略》阅读人数近3000人，有需求的用户还是非常多的。其中提到了4个软件，之前已经介绍了其中非常有特点的三种，分别为：
PICRUSt预测宏基因组： 优点是高引流行，但缺点是数据太旧，GreenGene数据库13年5月后再无更新。
FAPROTAX预测元素循环：主要基于物种名预测元素循环相关功能。
bugbase预测细菌表型：主要预测表型，如革兰氏阴阳、氧利用等，也基于Greengene数据库。
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.01. 导入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 02. 预测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.0006707625944477453. 预测代谢部分
上面的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包的帮助。
猜你喜欢
扩增子分析：1图表解读 2分析流程 3统计绘图 4美化进化树ggtree
16S预测功能 1KEEG通路PICRUSt 2元素循环FAPROTAX 3表型bugbase 4KO通路Tax4Fun
宏基因组分析：1拼接MEGAHIT 2基因注释Prokka 3Binning 4MetaPhlAn2
系列教程：1生信入门Biostar 2微生物组 3宏基因组
写在后面
宏基因组公众号经过半年发展，原创文章过百，关注人数12000+。为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外100+ PI，900+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”