Archives

16S扩增子组成谱预测功能谱:本地快速执行otutable_norm和otutable_melt以及otutable_kann

当前可以使用16S扩增子组成谱预测功能谱的工具主要有:PICRUStTax4funPiphilin, 本文主要使用PICRUSt的预测文件,进行重新实现。

PICRUSt 采用了一种ASR(ancestral state reconstruction, 可参考1985年的一篇文章:Phylogenies and the Comparative Method)方法重构物种的16S拷贝数和基因组的基因功能信息(使用KEGG、COG描述),基于系统进化的逻辑是同一个clade分支下的基因组(物种)具有相似的形状,这样就可以使用已测序基因组预测未测序基因组的功能空间,主要依赖参数为系统进化树分支距离。

PICRUSt

图来自PICRUSt文章

原理图分拆解读:

1. 使用Greengene 16S序列构建系统进化树,去除未测序基因组;
2. 预测所有祖先状态的性状(16S拷贝数、基因功能);
3. 预测未知基因组节点的性状(16S拷贝数、基因功能);
4. 标准化,使用预测的拷贝数对OTU丰度进行标准化;
5. 基因功能空间预测(KO、COG)丰度表;

其实关键算法在1、2、3部分。 这篇文章我们直接使用PICRUSt根据 Greengene 13.5 进行的预测结果进行一些数据分析。

需要的文件为双端合并后的序列数据(Usearch sample=’xxx’ 注释文件),我们会在下面一步一步演示如何操作。 准备工作:

下载PICRUSt预计算拷贝数和KO表: 16S_13_5_precalculated.tab.gzko_13_5_precalculated.tab.gz 下载Greengene 97OTU代表序列: 97_otus.fasta.gz

如何使用这些程序,请移步: 如何使用biostack 提供的命令行工具

1. 构建OTU表

虽然有文章讨论过使用closed- 和 open-reference 构建 OTUs 有多不准确,如下:

Accuracy of microbial community diversity estimated by closed- and open-reference OTUs.
Edgar RC.
PeerJ. 2017 Oct 4;5:e3889. doi: 10.7717/peerj.3889. eCollection 2017.
PMID: 29018622

如果有兴趣看看问题在哪里, 可以移步:closed_ref_problems, 这不是这篇文章的主题,因为我们使用基于Greengene 13.5 的PICRUSt预先计算的结果,所以只能选择使用Greengene作为参考构建OTU表。

这里使用Usearchclosed_ref 构建closed OTU表;

closed_ref 命令行接口:

usearch -closed_ref reads.fastq -db gg97.fa -otutabout otutab.txt -strand plus -tabbedout closed.txt

下面简单解读:

reads.fastq: 为具有 'sample=xxx' 注释的16S扩增子序列集合。
-db        : 为greengene 97%聚类 OTU代表序列。
-strand    : plus或者both,搜索查询序列的条件,只都搜索正链还是正负链都搜索。
-otutabout : OTU表输出;
-tabbedout : closed_ref 比对report文件;
-id        : 最低相似度要求,默认参数97%;

下载Greengene 97% OTU代表序列文件:

wget https://github.com/biocore/qiime-default-reference/blob/master/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta.gz

执行构建OTU表:

gunzip 97_otus.fasta.gz
usearch -closed_ref reads.fastq -db 97_otus.fasta -otutabout otutab.txt \
        -strand plus -tabbedout closed.txt -threads 20

2. 16S OTU表 标准化

上一步我们已经获得了OTU表信息,下一步就是使用Greengene的16S 拷贝数进行标准化OTU表; 我们使用 `otutable_norm` 程序实现标准化工作

命令行接口:

$ otutable_norm
Usage: otutable_norm <copy_number> <otutab>
version: 0.0.1

otutable_norm 接受两个参数: 拷贝数表和OTU表;

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/16S_13_5_precalculated.tab.gz
otutable_norm  16S_13_5_precalculated.tab.gz  otutab.txt >otutab_norm.txt

3. 合成功能丰度表

这一步使用 otutable_melt 完成合成功能谱的工作。 命令行接口:

$otutable_melt
Usage: otutable_melt [options] <feature-matrix> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_melt 接受两个参数:PICRUSt兼容的特征表和OTU表(QIIME兼容的OTU表),这里使用上一步产生标准化后的OTU表, otutable_melt也接受 Bugbasedefault_traits_precalculated.txt数据;

-r: 四舍五入,转换成整型输出。

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/ko_13_5_precalculated.tab.gz
otutable_melt   ko_13_5_precalculated.tab.gz otutab_norm.txt  >ko_profiling.txt

或者:

otutable_melt   <(pigz -dk -c  -p 20  ko_13_5_precalculated.tab.gz)  otutab_norm.txt  >ko_profiling.txt

现在我们已经可以获取KEGG的入门券了,后面我们可以使用KO关联KEGG的很多重要信息, 比如 Pathway,Module,Enzyme等。

4. otutable_kann

KO(KEGG Ortholog)是KEGG唯一一个通过序列和功能关联的中间点,我们获得KO的丰度信息, 可以通过合并同一个Pathway或Module的KO丰度,计算Pathway和Module的丰度。

实现该功能我们使用 otutable_kann实现;

命令行接口:

$ otutable_kann
Usage: otutable_kann [options] <feature-map> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_kann接受两个参数, KEGG的KO映射文件和OTU表文件(KEGG或COG的丰度表这里统一为otu_table表)。 文件为两列:

第一列: KO Number 第二列: KO 映射表;

该文件可以通过KEGG的API服务下载,比如:

curl http://rest.kegg.jp/link/pathway/ko | \
sed 's/ko://; s/path://' | grep "map"  >ko_pathway.tab
curl http://rest.kegg.jp/link/module/ko |   \
sed 's/ko://; s/md://'  >ko_module.tab

下载注释文件:

curl http://rest.kegg.jp/list/module | sed 's/md://' >module-definition
curl http://rest.kegg.jp/list/pathway | sed 's/path://' >pathway-definition
curl http://rest.kegg.jp/list/ko | sed 's/path://' >ko-definition

执行(以Module为例子):

otutable_kann   ko_module.tab  ko_profiling.txt >module_profiling.txt

如果需要注释信息需要使用 tabtk_definition 添加:

otutable_kann   ko_module.tab  ko_profiling.txt  | \
tabtk_definition -d " " module-definition  - >module_profiling.txt

5. 总结

到这里我们可以很简单高效完成共物种组成谱到功能谱的预测,我们提供的otutable_normotutable_meltotutable_kann 都比较通用,组合起来可以解决很多问题。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-14 7:46 PM

Comments are closed.