PhyloPythiaS+: a self-training method for the rapid reconstruction of low-ranking taxonomic bins from metagenomes

PhyloPythiaS

Assignment with PPS proceeds in two steps: In the training step, an ensemble of structured output Support Vector Machines (SVMs) (Joachims, Finley & Yu, 2009) for the specified part of the NCBI taxonomy, defined by the taxa being modeled, are trained using the sample-derived training sequences and additional data for these taxa from a customized reference collection of sequenced genomes and draft genomes (Suplemental Information 1, Section 3.3). The list of modeled taxa and sample-derived data are generated with the + component of PPS+. The list of taxa restricts the taxonomic output space that is modeled, i.e., a sequence from a metagenome sample will be assigned to a leaf node taxon or a corresponding higher-ranking taxon of the learned taxonomy.

In the prediction step, the PPS model ensemble identifies the taxon which best matches a query sequence in terms of its k-mer profile and assigns to it the respective taxonomic identifier. By default, sequences of 1 kb or more are classified (PPS+ configuration parameter: minSeqLen).

The + component of PhyloPythiaS+

The input for the + component of PhyloPythiaS+ is the metagenome sample. This step returns a list of clades and sample-derived data for the subsequent PPS training. The + component performs the following steps:

1. Marker gene identification: DNA sequences from the sample are translated in all six reading frames (i.e., also considering reverse complement sequences) to protein sequences. In both the translated and untranslated sequences, regions with similarity to the DNA or protein Hidden Markov model (HMM) profiles of 34 taxonomically informative marker genes (Wu & Eisen, 2008; Stark et al., 2009; Liu et al., 2011; Wu & Scott, 2012; Segata et al., 2012; Sunagawa et al., 2013) are identified (Supplemental Information 1, Section 3.3 and 6.1). The corresponding DNA marker gene sequences from these regions are used for further analysis.

2. Taxonomic marker gene assignment: The marker gene sequences are assigned a taxonomic identifier using the composition-based Naïve Bayes classifier (Schloss et al., 2009) (Supplemental Information 1, Section 6.2).

3. Taxonomic sequence assignment: If a sequence contains multiple marker genes, multiple taxonomic identifiers are identified in Step 2. Then the highest bootstrap confidence score (hcs) returned by the Naïve Bayes classifier (NBC) for one of the markers on the fragment is identified. We use all marker gene assignments with confidence scores larger than hcs ^∗ (1 − candidatePlTopPercentThreshold). The default setting for the configuration parameter candidatePlTopPercentThreshold is 0.1. From the set of taxonomic identifiers, the lowest taxon t is identified for which all other assignments are either to the same taxon t or defined at higher-ranking parental taxa of t. Taxon t is consequently used for the overall fragment assignment. The confidence score for the fragment is set to the smallest confidence score for the set of retained marker gene assignments.

4. (Optional: Taxonomic scaffold assignment): Scaffolding information (i.e., the mapping of contigs to scaffolds) can be used to obtain more training data for the relevant taxa. Assembled contigs can be grouped into scaffolds based on the paired-end information after the assembly. As all contigs of a particular scaffold originate from the same strain, all contigs of the respective scaffold should have the same taxonomic label. Here, we make use of this scaffolding information, such that unassigned contigs of a particular scaffold can be assigned based on the assigned contigs of the respective scaffold. In the first step, the taxonomic identifiers of all assigned contigs for a scaffold are corrected as follows: Let us consider that n taxonomically assigned contigs of a scaffold are placed along a common path from the root r down to a low-ranking clade lc in the reference taxonomy. The unassigned contigs of a scaffold are not among these n contigs. To obtain a consistent assignment for all the contigs of a scaffold and to correct for ‘outlier’ contig assignments to low ranking taxa, contigs are reassigned according to the following: All n assigned contigs of the respective scaffold are reassigned to the lowest taxon c, which lies on the path from r to lc, where c is chosen such that at least (agThreshold ^∗ n) of the contigs are assigned on the path from c to lc. In the second step, unassigned contigs are assigned to the same taxon c, if a sufficient number of contigs have already been assigned. Let us denote the sum of all contig lengths for a scaffold as l and the sum of all assigned contig lengths of the respective scaffold as al. If al/l ≥ assignedPartThreshold, then the unassigned contigs of a scaffold are also assigned to clade c (see the configuration parameters: placeContigsFromTheSameScaffold = True, agThreshold = 0.3, assignedPartThreshold = 0.5).

5. Assignment path truncation: Contigs assigned to a lower-ranking taxon than the specified lowest rank are reassigned to the parental taxon of this lowest rank (configuration parameter: rankIdCut).

6. Taxa selection for model specification: Any taxon for which at least 100 kb of sample-derived data have been identified can be modeled. Furthermore, species can be modeled if at least 300 kb of reference sequences are available in the reference sequence database, and higher-ranking taxa can be modeled if data for at least three distinct species with this requirement (>300 kb per species) are available. Contigs assigned to taxa for which there are fewer data are subsequently assigned to higher taxonomic ranks for which sufficient data are available to allow their use as sample-derived training data (configuration parameters: minGenomesWgs = 3 or 1, minBpPerSpecies = 300,000, minBpToModel = 100,000).

7. Abundant taxa selection: To reduce the number of taxa to the most relevant ones, the least abundant taxon is removed iteratively. This is defined as the taxon to which the minimum number of bp is assigned. Sequences assigned to this taxon are reassigned to the closest defined taxon at a parental rank. The algorithm ends when the number of leaf taxa is less than or equal to the maximum number of taxa to be modeled (configuration parameter: maxLeafClades = 50; this can be set realistically up to 800).

Balancing training data: The part of the taxonomy that will be modeled with PPS is defined by the taxa identified in the previous step. It has leaf nodes at different ranks above the specified rank cut-off, and internal nodes. Only leaf node taxa and sample-derived training data assigned to leaf node taxa in the preceding steps are specified as input for PPS training. To balance the training data across clades, a maximum of 400 kb of sample-derived training data are selected for each leaf node taxon (configuration parameter: maxSSDfileSize). For this selection, contigs are used in order of decreasing confidence values and then in order of decreasing length. The balancing of training data can be switched off by setting the configuration parameter (maxSSDfileSize) to a large number.

Simultaneous counting of multiple short k-mers

We provide PPS+ with a new custom k-mer counting algorithm that is based on the Rabin Karp string matching algorithm (Karp & Rabin, 1987). The algorithm is highly optimized to count occurrences of short DNA sequences. It is very fast, as it is memory efficient, because it does not need any large helper data structure similar to suffix trees. It explores the locality of reference, uses very fast bit shift operations and is efficiently implemented in C. Its complexity is O(n), where n is the length of the DNA sequence that is being considered. It enumerates k-mers up to hundred times faster than when using suffix trees that were employed in PPS. This made PPS+ overall up to 3x faster than PPS. Because the algorithm allows to simultaneously enumerate k-mers of consecutive lengths in one run, it is at least 2–7x faster than the state-of-the-art software Jellyfish (Marcais & Kingsford, 2011) and 11x faster than KAnalyze (Audano & Vannberg, 2014) in the scenario used in PPS+, i.e., when calculating k-mers of length 4, 5, and 6 for every sequence (Table S1, Supplemental Information 1, Section 2). We also found that the state-of-the-art k-mer counting methods KMC 2 (Deorowicz et al., 2015) and Turtle (Roy, Bhattacharya & Schliep, 2014) are not applicable to our problem setting, as KMC 2 can count only k-mers ≥ 10 and Turtle is prohibitively slow for sequences ≥ 16 kb.

Benchmarks with simulated datasets

We constructed the simulated datasets by assembling simulated reads with an empirical error profile. The details on how the simulated reads were generated and assembled can be found in (Supplemental Information 1, Section 3.1). For the evaluation, precision and recall were calculated (Supplemental Information 1, Section 3.9). Furthermore, these measures were also calculated with a ‘correction,’ to account for the case where the sequences of one taxon were consistently assigned to a different taxon, as for draft genome reconstruction, it is more important that the sequences are assigned consistently than that the taxonomic identifier is correct. To assess the performance of the different methods in assigning the simulated sequence fragments without related reference genomes being available, ‘new strain,’ ‘new species’ and ‘new genus’ scenarios were simulated by removing all sequence data from the taxa of the simulated test dataset at each rank from the reference data. Furthermore, for PPS+, we distinguished whether the reference data were excluded (masked) from the reference sequence (RS) collection or also from the marker gene (MG) collection, since the MG collection included sequences for 15 times more distinct species than the RS collection. There were therefore two different situations to consider (Table 1).

PPS+ showed a substantially improved precision and recall over the PPS generic model, which demonstrated the impact of the improved selection of training data and modeled taxa (Figs. 2A and 2C, S1A–S1D and S3A–S3D). PPS+ almost always had higher precision and recall than MEGAN4 and Kraken, except when almost all test data were included in the reference sequences (Figs. 2A and 2C, S1A–S1C, S1E, S3A–S3C, S3E, S14A). This was even more pronounced when comparing bin quality using the corrected measures (Figs. 2B and 2D, S2A–S2C, S2E, S4A–S4C, S4E, S14A and S14D). When comparing PPS+ to taxator-tk, PPS+ had substantially improved recall, particularly for lower ranks (Figs. 2A and 2C, S1A–S1C, S1F, S3A–S3C, S3F); while taxator-tk outperformed all other methods in terms of precision (Figs. 2A and 2C, S1A–S1F and S3A–S3F). Both methods were similarly precise when analyzing bin recovery, independent of assigning the taxonomic identifiers to the corrected measures (Figs. 2B and 2D, S2A–S2C, S2F, S4A–S4C and S4F). As a strong point of PPS+ , we also observed that it more rarely predicted wrong taxa that were not a part of the sample than the other methods (Fig. S5). For example, for the genus rank in Scenarios 3 and 8, PPS+ assigned sequences to only 2–5 false positive taxa, while taxator-tk identified 20, MEGAN4 37 and PPS 59 false ones. If PPS+ identified wrong taxa, these were usually very closely related to the true taxa.

Benchmarks with real datasets

Comparison to an expert binning based on marker genes

A taxonomic binning generated by PhyloPythia (PP) with expert guidance for sample-derived model construction (Turnbaugh et al., 2010) was compared to the PPS+ assignments. Scaffolds that were unassigned by either method were not considered. The PP expert binning and the PPS+ binning agreed well, down to the order rank (Table 4). For the family and genus ranks, the overlap of both methods dropped to 69.5–74.1%, which may partly be due to changes in the NCBI taxonomy since the generation of the expert binning in 2009. Both PPS+ and PP assignments were highly consistent with the MG assignments made by the + component of PPS+ alone, though only a small number of scaffolds with marker genes could be compared (7–23% for different ranks). While PPS+ had a larger overlap (‘% agreement’) with the MG assignments at the genus rank, PP had a larger overlap (‘% agreement’) with the MG assignments at the family rank. Moreover, we compared the number of taxonomic assignments for individual methods (Fig. 3): PPS+ assigned sequences to low-ranking taxa down to the species level, in agreement with the MG assignments, while PP often assigned the respective sequences only to the parental taxa. This demonstrates that PPS+ can generate a high quality taxonomic binning in a fully automated manner.

Table 4:

Comparison to an expert binning based on marker genes.

Comparison of the taxonomic assignments of PPS+ versus PhyloPythia (PP), with expert guidance for sample-derived model construction (Turnbaugh et al., 2010) for the human gut scaffolds (161, 343 kb) based on marker genes (MG), using the + component of PPS+. The measure ‘% agreement’ represents the percentage of bp assigned by both methods to the same taxonomic identifiers at a given rank, whereas ‘kb agreement’ is the corresponding number of kb assigned by both methods to the same taxonomic identifier. Scaffolds assigned by only one method are not considered in this comparison. Bold numbers correspond to the best values, whereas italic numbers indicate the worst values.

Comparison	Rank	% agreement	kb agreement
PP vs PPS+	Superkingdom	99.6	160,617
MG vs PP	Superkingdom	99.7	38,314
MG vs PPS+	Superkingdom	99.5	38,220
PP vs PPS+	Phylum	95.4	149,213
MG vs PP	Phylum	96.9	17,771
MG vs PPS+	Phylum	98.7	18,065
PP vs PPS+	Class	97.0	145,887
MG vs PP	Class	98.1	17,599
MG vs PPS+	Class	100.0	17,869
PP vs PPS+	Order	98.0	145,373
MG vs PP	Order	98.3	17,494
MG vs PPS+	Order	100.0	17,764
PP vs PPS+	Family	69.5	95,779
MG vs PP	Family	90.7	13,047
MG vs PPS+	Family	83.7	12,013
PP vs PPS+	Genus	74.1	78,686
MG vs PP	Genus	91.6	12,235
MG vs PPS+	Genus	94.9	11,479

DOI: 10.7717/peerj.1603/table-4

Throughput comparison

The throughput of the individual methods for contig assignments of the human gut sample was calculated (Supplemental Information 1, Section 3.3, 3.4 and 5.3). The throughput of Kraken substantially varied between 38.4 Mb/h and 4.2 Gb/h in our experiments, depending on whether its large (∼200 GB) reference database was already loaded in the main memory or not, therefore Kraken is the fastest method in high performance environments. When only the prediction step of PPS+ was considered, PPS+ assigned up to 0.5 Gb/h and was more than 7 times faster than the homology-based methods (Fig. 4). This is relevant when PPS models are reused for the classification of another sample. Moreover, unlike the homology-based tools and Kraken, PPS+ can be run on a standard laptop, as it requires much less main memory (see Supplemental Information 1, Section 3.4 for the hardware configurations used).