功能预测之Tax4Fun

最新推荐文章于 2023-09-14 08:44:03 发布

刘永鑫Adam

最新推荐文章于 2023-09-14 08:44:03 发布

阅读量7.7k

点赞数 3

本文链接：https://blog.csdn.net/woodcorpse/article/details/106552537

版权

之前在公众号的文章《根据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.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包的帮助。

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外2600+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。