MOCAT2: a metagenomic assembly, annotation and profiling framework.

UNLABELLED: MOCAT2 is a software pipeline for metagenomic sequence assembly and gene prediction with novel features for taxonomic and functional abundance profiling. The automated generation and efficient annotation of non-redundant reference catalogs by propagating pre-computed assignments from 18 databases covering various functional categories allows for fast and comprehensive functional characterization of metagenomes.
AVAILABILITY AND IMPLEMENTATION: MOCAT2 is implemented in Perl 5 and Python 2.7, designed for 64-bit UNIX systems and offers support for high-performance computer usage via LSF, PBS or SGE queuing systems; source code is freely available under the GPL3 license at http://mocat.embl.de CONTACT: : bork@embl.de SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

Metagenomics has enabled large-scale studies investigating the structure, function and diversity of microbial communities. The computational analysis of samples, often totaling many gigabases of sequence data, usually involves mapping reads to taxonomic and functional reference databases (which may require the de novo assembly of predicted genes), and subsequent abundance profiling. Whereas taxonomic profiling methodology has matured recently (Segata et al., 2013; Sunagawa ), functional profiling still remains challenging due to the difficulties in assigning functions to millions of reads from metagenomes. Moreover, current metagenomic pipelines (Abubucker ; Bose ; Edwards ; Glass ; Huson ; Lingner ; Markowitz ; Meinicke, 2015; Glass ; Huson ; Lingner ; Abubucker ; Edwards ; Bose ; Silva ) for functional annotation and/or profiling mainly implement metabolic pathway or protein domain databases (Segata et al., 2013) such as KEGG (Kanehisa ), SEED (Overbeek ) or Pfam (Finn ). Here, we present metagenomic analysis toolkit version 2 (MOCAT2), which was developed to enable functional profiling of metagenomes based on a much wider range and diversity of functional gene annotations. Its features are compared to existing tools in Supplementary Table S1.

2 The MOCAT2 pipeline

The metagenomic analysis toolkit (MOCAT) (Kultima ) proceeds through the following steps: raw sequence reads are quality-filtered and subsequently assembled into longer contigs, on which open reading frames are predicted (Fig. 1).

Fig. 1.

The MOCAT2 pipeline. Read quality control, assembly and gene prediction represent the original MOCAT pipeline (dark green box). Blue path: Genes are clustered into reference gene catalogs, which are functionally annotated. Orange path: To quantify functional composition, reads are mapped to the annotated gene catalog and summarized over the respective annotation categories. Taxonomic profiles (mOTU, specI and NCBI) are generated by mapping reads to mOTU and reference marker gene (RefMG) catalogs

Its main extensions in MOCAT2 enable comprehensive functional profiling, in addition to the eggNOG database, by integrating 18 publicly available resources that cover diverse functional properties (Table 1). The databases were selected to include large, widely used protein databases, as well as ones targeting specific functional categories (Supplementary Text). Each database has been filtered for relevance, for example from the eukaryote-centered database DrugBank only the genes with bacterial homologs were extracted.

Table 1.

Databases from which functional properties are obtained

	Proteins	Coverage	Precision	Recall	Reference
Protein domains and families
eggNOG	7 449 593	100	100	100	Huerta-Cepas et al. (2015)
Pfam	16 230*	87	90	94	Finn et al. (2014)
Superfamily	15 438*	93	89	94	Gough et al. (2001)
(Metabolic) pathways
KEGG	7 423 864	98	93	93	Kanehisa et al. (2014)
MetaCyc	388 782	100	89	94	Caspi et al. (2014)
SEED	4 247 700	99	94	94	Overbeek et al. (2014)
Antibiotic resistance
ARDB	25 360	89	99	88	Liu and Pop (2009)
CARD	2 820	100	81	93	McArthur et al. (2013)
Resfams	123*	80	94	94	Gibson et al. (2014)
Virulence factors
MvirDB	29 357	100	95	93	Zhou et al. (2007)
PATRIC	2 194 475	93	93	93	Mao et al. (2015)
vFam	29 655	35	99	86	Skewes-Cox et al. (2014)
VFDB	1 627 380	86	89	91	Chen et al. (2012)
Victors	3 329 893	91	92	94	Mao et al. (2015)
Complex carbohydrate metabolism
dbCAN	333*	76	99	99	Yin et al. (2012)
Bacterial drug targets and exotoxins
DBETH	228	100	99	86	Chakraborty et al. (2012)
DrugBank	3 899	99	88	94	Knox et al. (2011)
Mobile genetic elements
ICEberg	13 984	98	79	91	Bi et al. (2012)
Prophages	119 183	95	88	91	Waller et al. (2014)

Coverage of each database in percent, e.g., of the 18 202 orthologous groups in KEGG (KO), 17 773 (98%) are covered and thus propagated by the eggNOG database. Coverage, precision and recall are given as percentages.

*Number of hidden Markov models (HMMs), whereby one HMM can hit several proteins and several HMMs can map to one protein.

To avoid the computational burden of mapping reads to multiple databases, predicted genes are first clustered using CD-HIT (Huang ) into a non-redundant gene set, called a reference gene catalog (Qin ). Next, this gene catalog is mapped to the eggNOG database with wide taxonomic coverage of orthologous groups, to which sequence annotations from other databases have been pre-computed so that functional information from multiple databases can be transferred efficiently to the catalog. This indirect annotation methodology not only provides a 10-fold speed up compared to directly mapping to each database separately, but also enables annotations of short genes, which would otherwise be missed (Supplementary Figure Fig. S1). For computational efficiency MOCAT2 uses DIAMOND (Buchfink ) in the annotation step. Combined, these features yield a more than 1400-fold annotation speedup over a conventional BLAST-based annotation pipeline (Supplementary Text). Users can either create and annotate their own gene catalogs de novo, or use pre-computed and pre-annotated reference gene catalogs for the human gut and skin, mouse gut, or the ocean (Li ; Oh ; Sunagawa ; Xiao ).

Finally, to quantify functional composition, reads from each sample are mapped to the annotated gene catalog and summarized over the respective annotation categories (Fig. 1).

MOCAT2 now also offers several approaches for taxonomic profiling, all of which are based on mapping reads to a benchmarked set of single copy marker genes (Fig. 1). Taxonomic abundance estimates are calculated not only for different NCBI taxonomic levels, but also for species clusters defined based on molecular sequence identity (specI; Mende ) and species that currently lack sequenced reference genomes based on metagenomic operational taxonomic units (mOTU; Sunagawa ).

3 Annotation and profiling benchmarks

As complex functional annotation based on 18 databases via indirect propagation of eggNOG annotations is conceptually new, we benchmarked the (indirect) MOCAT2 annotations and functional profiles (Supplementary Table S2 and Supplementary Text).

First, we compared the indirect annotations to the direct ones (generated using the annotation tool of each individual database or recommended pipeline and cutoffs) for >65 million genes from five diverse datasets (precision and recall are listed in Table 1).

Next, using data from (Zeller ) we compared the direct annotations to ones produced by COGNIZER and UProC (Bose ; Meinicke ), two recently developed annotation tools integrating multiple databases. In our tests, MOCAT2 annotations were either similar to, or more accurate, than those of COGNIZER and UProC (Supplementary Table S3).

Finally, the functional abundance profiles obtained using the indirect MOCAT2 annotations were very similar to those obtained using the direct method (Spearman = 0.95; n = 1300).

4 Conclusions

MOCAT2 is a software pipeline for metagenomics using state of the art assembly, annotation as well as taxonomic and functional profiling approaches in this fast moving field. Generating and annotating gene catalogs with precomputed assignments to a large selection of functional databases allows for comprehensive and efficient functional profiling of complex microbial communities. MOCAT2 thus enables such analysis at an extent far beyond what other tools currently offer and is scalable to the anticipated deluge of metagenomic data from diverse sources.