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包就上线了。

微生信生物

学习永无止境,分享永不停息!

  • 发表于:
  • 原文链接:http://kuaibao.qq.com/s/20180116G00F4600?refer=cp_1026

相关快讯

扫码关注云+社区