内容预览
对同一批数据进行注释比较
Metawrap与GTDB-tk的物种注释原理
选择适合自己研究目标的物种注释方法
小彩蛋
前言
宏基因组学作为一种获取环境样本中整个生物群落遗传物质的重要工具，在微生物、环境和生态学中受到广泛的应用。通过宏基因组的Binning组装和注释，我们可以快速地获取到样本中微生物的基因组信息（这些基因组我们称之为MAG，Metagenome-Assembled Genome）。
然而，与微生物组16S/ITS面临的问题一样，宏基因组在处理MAG的物种注释上也没有统一的执行标准。对于MAG的物种注释，有两个因素影响比较大，一是注释原理差异，分为两类：1.基于系统发育的方法；2.基于分类器的数据库检索。由这两种方法延伸出来软件在理论上各有优劣，应用场景重合度低。
二是不同软件的物种分类体系不相兼容，例如GTDB-tk、PhyloPhIAn的分类与NCBI的分类体系不兼容，导致研究人员有时候在对下游数据进行合并操作时，遇到许多无法合并的种属。
| 方法 | 优点 | 缺点 | 研究目的 | 适用软件 |
|---|---|---|---|---|
| 基于系统发育 | 能够提供物种进化历史的信息 | 对于缺乏参考序列的物种，无法做到准确的注释；多序列比对需要耗费大量的计算资源 | 系统发育分析 | PhyloPhIAn、PhyloSift |
| 基于分类器的数据库比对 | 快速，易实现 | 受数据质量和覆盖范围影响大；对新物种识别能力差 | 全面分类分析/快速分类标签 | GTDB-tk、PhyloPhIAn、BAT、Kraken |
尽管在系统发育中，宏基因组也起着非常大的作用，不过相比起其在大规模探索环境样本中微生物物种组成的领域，如何进行准确的高通量物种注释更值得探讨。
因此，本文将对比目前较流行的metawrap的物种注释模块和GTDB-tk物种注释，通过解读其分类原理，帮助大家选择适合自己研究目标的注释工具。
对同一批数据进行注释比较
为了直观的展示metawrap物种注释模块与GTDB-tk对MAG注释的差异，我选择以一批人工混合的微生物样本，输入到metawrap流程中，进行MAG的binning组装，在得到MAG后，分别使用metawrap的物种注释模块（后文均以注释模块中的核心软件taxator来表示）和GTDB-tk进行注释分析。（分析流程见下图）
结果：
物种的注释对齐率
在注释完成后，我们以能对齐人工混合清单的注释结果来比较两者的MAG获得率，metawrap的平均MAG获得率（5%），GTDB-tk的平均MAG获得率（26%）。
可注释到的物种种类
在GTDB-tk中，总共注释到130种菌株，在taxator中，总共注释到60种菌株。有趣的是，taxator还注释到11种GTDB没有注释到的结果。
以上结果表明，两款软件在处理同一批MAG输入数据时，存在明显差异。
我在这里提出两个问题：
👉
1.为什么两者之间差距这么大，哪个因素是对物种注释产生关键影响的？
👉
2.为什么同一批样本，明显注释数量更优的GTDB，还存在部分物种仅taxator注释得到，可能是什么因素影响的？
Metawrap与GTDB-tk的物种注释原理
Metawrap的物种注释模块（Taxator-tk），其功能实现主要是脚本：metawrap/bin/metawrap-modules/classify_bins.sh
通过分析这个脚本，我们可以清晰的看到metawrap的物种注释逻辑
注释使用的物种数据库为NCBI的nt库，而进行注释分类的软件为taxator-tk。
简单介绍一下taxator-tk软件，这款软件主要是用于核苷酸/氨基酸序列分类，主要针对复杂环境的单元分类。这款软件在分配分类群时，采用了一种名为查询段分类标记的算法（realignment placement algorithm/RPA）。该算法通过构建系统发育树（a）、计算类群之间的距离（b）、以及将片段分类单元添加到空集合M（c、d）和分配单元到类群树下的具体位置（e），来实现注释单元的分类。
这款软件在实际调教过程中，参数的设置是采用的比较保守的策略，目的是为了提高其在识别相似性物种上的差异，例如两个单元在属层级的遗传距离上没有大于设定的阈值，软件则会放弃向更深一层（种）的分配。因此通常会导致MAG的注释层级无法达到种水平以下，如果经常运行metawrap的小伙伴仔细观察过BIN_CLASSIFICATION目录下的注释文件，就会发现这一现象非常常见。
与之相对的，GTDB-tk整体分配注释单元的逻辑如下
GTDB-tk依靠于其强大的自建数据库，在进行物种注释之前，会先进行FastANI的快速检索，以帮助判断物种种属的区间，之后结合RED值来决策最终的注释。而GTDB-tk之所以能够在种水平下获得更多的注释结果，就在于其使用RED值进行物种注释决策的设定。
在GTDB-tk论文中，作者这样写道：
测试显示，如果存在新的基因组时，GTDB-tk与手动管理差异最为明显。
理论上来说这些注释结果需要手动矫正，定义为新的分类单元，而GTDB-tk会非常容易地过度分类。
因此，通过分析算法的特点，我们可以很容易的回答前一节中的问题1，两款软件的注释结果差异，在很大程度上是受分类算法特性影响。
对于问题2，在发现taxator注释结果中，存在部分GTDB未注释到的情况后，我检查了GTDB的物种注释数据库清单，单从物种清单来看，taxator注释到的物种均在目录上，因此我怀疑存在两种可能性，一是GTDB算法在应对这部分物种序列分配时存在缺陷，二是GTDB进行fastANI初始判断时，参考数据库中的物种序列与本地的菌株存在一定差异，因此在初筛时就被丢弃。简单来说，就是GTDB算法和参考数据库还值得进一步的探索和讨论。
选择适合自己研究目标的物种注释方法
通过评估两款软件注释结果，我们认识到算法、数据库以及一些特定物种对结果的影响。那么，在应用到实际研究中，我们应该如何设计自己的分析管道呢？
我这里分享一点小小的经验，以研究对象的了解情况，分为两大类：
一是对于完全未知物种的情况
我们此时的目的通常是了解环境样本中可能与实验变量中高度相关的对象，这种情况下，确定一个具体的对象能够帮助我们快速的进行下游meta分析。因此，采用类似GTDB-tk、Kraken2、PhyloPhIAn的软件，能够实现对目标样本的大范围扫描的作用。之后再结合宏基因组关联分析，来进行生物标志物的发现工作。
二是对于已知主要种类的情况
此时我们对目标环境样本具有一定的背景知识，相较于完全不了解环境样本的情况，我们可能已经确定部分高相关性的菌株。因此，可能会涉及到具体菌株的功能发现、核心菌株的系统发育分析等，以更高精度的方式定位菌株，帮助实验端分离、验证微生物的功能。
以上两种情况仅从我个人的项目经验出发，在目前的宏基因组研究中，涉及到的领域非常广泛，如果有一些更好的方法论，欢迎大家与我讨论。
小彩蛋
不知道有没有读者发现，我在前言开篇时，提到了GTDB分类体系与NCBI的分类体系存在差异，但是在后续分析中，我似乎能无缝的对比两组软件的差异。这是不是说明两者之间的分类体系差异并没有我们想象的那么大？
实际上，我这里是采用了GTDB提供的一个三方脚本和对应关系表，将GTDB的注释结果转换成了NCBI体系。这里提供一下操作命令和转换关系表的数据下载链接（对应关系表的下载地址：https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/bac120_metadata.tsv），供有需要的朋友使用：
python /Software/GTDBTk-master/scripts/gtdb_to_ncbi_majority_vote.py --gtdbtk_output_dir BIN_CLASSIFICATION/${sampleID} --output_file BIN_CLASSIFICATION/${sampleID}_ncbi --bac120_metadata_file /mnt/raid1/Databases/GTDB/bac120_metadata.tsv
本文由博客一文多发平台 OpenWrite 发布！
