16s分析之Tax4Fun功能预测使用笔记

有一天我在网上发现了这篇文章（当然我是后知后觉）：

看到了这张图片：

穷人总是想法多，处于做功能预测的需求，而且本人做的是土壤，我看到这里，就想尝试一下这种方法：

在readme这里下载

这份完整的教程这里也给出链接：

http://tax4fun.gobics.de/RPackage/Readme_Tax4Fun.pdf

我才用Qiime聚类，通过silva数据库注释：

# seqs.fna为拼接好（单端测序不需要），质控，并且去除嵌合体完成，将所用样本序列合并起来的文件如果没有合并，可以通过cat *.fna > seqs.fna合并

#聚类OTU

pick_otus.py -i seqs.fna -o picked_otus

#挑选丰度最高的OTU为代表序列

pick_rep_set.py -ipicked_otus/seqs_otus.txt -f seqs.fna -o rep_set.fna -m most_abundant

#注释代表序列文件：

assign_taxonomy.py -i ~/Desktop/rep_set.fna-b ~/Desktop/ silvasuuref123nr/SILVA_123_SSURef_Nr99_tax_silva.fasta -t ~/Desktop/silvasuuref123nr/SILVA_123_SSURef_Nr99_tax_silva.tax -m blast -o taxonomy

注意，这里给出参考命令，它使用的115的数据库，但是已经更新到123了，从上面网站下载最新数据库：

这里还需要注意的是：使用Qiime注释非常缓慢，因为默认注释使用单核，可以使用多核。

#这里四核采用如下命令:

parallel_assign_taxonomy_blast.py -i~/Desktop/Shared_Folder/rep_set.fna -b silvasuuref123nr/SILVA_123_SSURef_Nr99_tax_silva.fasta-t silvasuuref123nr/SILVA_123_SSURef_Nr99_tax_silva.tax -o cs_rep_set_taxonomy-O 4 -U start_parallel_jobs.py

# 添加物种信息至OTU表最后一列，命名为taxonomy,我修改了工作目录，请自行修改

biom add-metadata -i otu_table.biom --observation-metadata-fprep_seqs_tax_assignments.txt -o otu_table_tax.biom --sc-separated taxonomy--observation-header OTUID,taxonomy

# 转换biom为txt格式，带有物种注释：

biom convert -i otu_table_tax.biom -ootu_table_tax.txt --to-tsv --header-key taxonomy

#转化txt文件

biom convert -i cs_rep_set_taxonomy/otu_table.biom-o cs_rep_set_taxonomy/otu_table.txt --to-tsv --header-key taxonomy--table-type='OTU table'

最后文件就出来了：

#

#

#

#下一步进入R操作，这里主要还是安装Tax4Fun包和依赖的包，这里我遇到的问题居多，下面截图是readme文件的参考，请自行查看，但是这里我给出我自己的安装过程：

从这里下载包，和参考数据库：

通过本地查找包安装：

找出自己包存放的位置：

当然会提示缺少很多依赖的包，这里我就不一一列举的，因为每个人缺少的包都会不太一样，我建议从 http://www.bioconductor.org/packages/release/bioc/html/BiocGenerics.html下载，然后在进行本地安装，成功率高一些：

选择window版本

最终经过我一个小时的折腾，终于将全部依赖的包，安装完成；

#载入包

library("Tax4Fun")

#进入工作目录

setwd("E:/Shared_Folder/HG_usearch_HG/16s功能预测_usearch/tax4fun")

#importQIIMEBiomData(file)开来不能直接用qiime出来的biom文件和txt文件，格式都有问题用官网例子文件替换格式之后，就没什么问题了

#导入文件

file

QIIMESingleData

# 根据Tax4Fun提供的SILVA123最新数据库进行预测，要求此数据的压缩包拉于此工作目录，这个命令得出来的是KO号的各种酶的基因丰度

Tax4FunOutput

#导出基因丰度表

write.table(KO_table,file="KO_table.txt",quote = FALSE,row.names = F,

col.names = T,sep = "\t")

#这个输出的是功能丰度表

index=print(Tax4FunOutput$Tax4FunProfile)

#转置

index=t(index)

write.table(index, file=" Tax4FunE.txt",quote= FALSE,row.names = F,

col.names = T,sep = "\t")

基因丰度表格：

通路表格：

到这里似乎就得到我们要的结果了，值得注意的是这份表格是相对丰度表格，差异分析就不能使用EDGR包了，这时候我们的limma包就上线了。

微生信生物

学习永无止境，分享永不停息！