qiime2+biom+qiime1获得16S物种丰度

我们知道,不管是16S等扩增子测序,还是宏基因组,最后最重要的结果,就是物种的丰度情况了,qiime2给出的16S丰度结果是一个计数,对于许多软件来说这是可用的,那么如果我们想获得一个直接的百分比数据应该怎样做呢?

当然,有许多方法可以实现,比如用shell, R, python脚本,或者再简单粗暴点,excel解决,透视表,公式,宏等。自己造轮子总是觉得不怎么踏实,出错咋办。那么,现成的软件有哪些呢,在这里,我抛砖引玉,提出一个曲线救国的方法,使用qiime2的前任qiime1解决,稍微做几步处理即可。如果你有更好的方法,欢迎交流和推荐,我们共同学习的进步!

这里我就从qiime2得出的结果直接开始,参考了生信菜鸟团大神的推文,这个大神的教程以全面著称,推荐学习!

1.导出物种分类信息和置信度

获得taxonomy.tsv,这个文件,其实把qza文件重命名为zip解压,或者qzv可视化文件导出,得到的文件也应该是一样的。

qiime tools export \
   --input-path taxa/taxonmony.qza --output-path taxa

文件是类似这样一个:

Feature ID Taxon Confidence
#q2:types categorical categorical
OTU_1 k__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Actinomycetales; f__; g__; s__ 0.8316610949745203

2.导出 BIOM 表,并加入将物种分类注释信息:

这一步就是处理下表头,让他兼容biom格式。注意,这个sed在mac下命令不能用,暂不确定是什么原因,我是用docker-ubuntu解决的。

#处理表头
sed -i -e '1 s/Feature/#Feature/' -e '1 s/Taxon/taxonomy/' taxa/taxonomy.tsv
#导出otu(feature)表
qiime tools export \
   --input-path deblur_output/table_final.qza \
   --output-path table_exported
#添加物种注释信息
biom add-metadata \
   -i deblur_output_exported/feature-table.biom \
   -o deblur_output_exported/feature-table_w_tax.biom \
   --observation-metadata-fp taxa/taxonomy.tsv \
   --sc-separated taxonomy
#biom转换成txt格式
biom convert \
   -i deblur_output_exported/feature-table_w_tax.biom \
   -o deblur_output_exported/feature-table_w_tax.txt \
   --to-tsv \
   --header-key taxonomy

3.qiime1获利各级分类结果

其实只需要biom格式就好了,唯一不足的是没有把上几级别的分类去除,比如属级别,还包括门纲目科,还不是usearch那种直接就是这个分类的结果。但是根据我的经验,usearch的物种注释结果可能不如qiime2的分类效果好,所以怎样结合这两个方法是个需要解决的问题。

#结果按门、纲、目、科、属五个级别进行分类汇总,对应结果的L2-L6
summarize_taxa.py -i result/otu_table3.biom -o result/sum_taxa # summary each level percentage

好的,我的分享就到这里,期待大家有更好的解决方案。

发表评论