16S预测宏基因组最强R包-Tax4Fun
以前在公众号的文章《根据16S预测微生物群落功能最全攻略》阅读人数近3000人，有需求的用户仍是很是多的。其中提到了4个软件，以前已经介绍了其中很是有特色的三种，分别为：
- PICRUSt预测宏基因组： 优势是高引流行，但缺点是数据太旧，GreenGene数据库13年5月后再无更新。
- FAPROTAX预测元素循环：主要基于物种名预测元素循环相关功能。
- bugbase预测细菌表型：主要预测表型，如革兰氏阴阳、氧利用等，也基于Greengene数据库。php
Tax4Fun简介
Tax4Fun 2015年发表于Bioinformatics，己被引用68次。ios
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.web
今天为你们带来Tax4Fun的教程。相比于picrust来讲，Tax4Fun最大的优势在于能够给予更新速度极快的SILVA数据库进行功能基因的预测。废话很少说，直接上干货。数据库
分析以前须要下载Tax4Fun包及比对库，此文用的最新的 SILVA 123 在Tax4Fun官网上下载 http://tax4fun.gobics.de/ 网络
Tax4Fun支持QIIME和SILVAngs的结果app
咱们须要下载R包，细菌宏基因组数据库(SILVA Reference data)，和QIIME格式的SILVA物种注释数据库(SILVA database for usage with QIIME)，固然帮助文档Readme也是必读的。ide
图1. 重要文件下载位置svg
安装
以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)
# 提取KO表，生成275个代谢相关的通路
KO_table = t(Tax4FunOutput$Tax4FunProfile)
# 全部样品标准化为1, 没有原始数据，能够用anova和Limma，没法使用edgeR和DESeq2
colSums(KO_table)
# 表头写个制表符，用于对齐表头
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)
基于QIIME结果格式数据库注释结果的分析
好比基于greengene注释的物种，或采用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表结果示例
猜你喜欢
- Top热文：1提DNA发Nature 2论文图表规范 3肠道菌群 4生信宝典 5实验vs分析 6咱们的将来在哪里
- 一文读懂：1微生物组 2进化树 3网络图 4图片配色
- 必备技能：1提问 2搜索 3Endnote 4在线绘图
- 文献阅读 1geenmedical 2热心肠 3SemanticScholar
- 扩增子分析：1图表解读 2分析流程 3统计绘图 4美化进化树ggtree
- 16S预测功能 1KEEG通路PICRUSt 2元素循环FAPROTAX 3表型bugbase 4KO通路Tax4Fun
- 宏基因组分析：1拼接MEGAHIT 2基因注释Prokka 3Binning 4MetaPhlAn2
- 科研经验：1云笔记 2云协做 3公众号 4微生物组入门 5成为团队VIP
- 系列教程：1生信入门Biostar 2微生物组 3宏基因组
- 生物科普 1生命大跃进 2细胞的暗战 3人体奥秘 4人体内旅行 5神奇细胞
写在后面
宏基因组公众号通过半年发展，原创文章过百，关注人数12000+。为鼓励读者交流、快速解决科研困难，咱们创建了“宏基因组”专业讨论群，目前己有国内外100+ PI，900+ 一线科研人员加入。参与讨论，得到专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA