Identifying Sequencing Reads from Multiple Sources using ContextMap

The original implementation of ContextMap presented previously [7] focused on refining initial mappings provided by other mapping algorithms. We recently developed a standalone version that can provide these initial mappings itself (see www.bio.ifi.lmu.de/ContextMap for details). Here, the central concept of both ContextMap versions is the so–called read context. This is defined as a set of reads originating from the same stretch of the genome, indicating that these reads were derived from the same transcript or different transcripts of the same gene. These contexts are defined based on the initial mapping and then extended in a subsequent re–alignment step, allowing a high degree of ambiguity both between and within contexts. For this purpose, ContextMap uses a modified version of Bowtie to identify spliced read alignments in a combination of forward and backward alignments. For each read not only the alignment with the minimum number of mismatches but any alignment to any context with at most a maximum number of mismatches is investigated. The unique mapping for the read to only one context is then determined by first finding the best mapping for the read in each context and subsequently finding the best context. For this purpose, a support score is used, taking into account the number of reads mapping within and around the region to which the read is aligned. Until the final step, contexts are treated independently of each other (see Figure S1).

As we show in this article, the advantage of this approach is that it allows investigating many alternative sources of reads in parallel, such as rRNA sequences, which are generally not included in reference genome assemblies of higher eukaryotes, as well as viral and microbial genomes. Contexts are then identified separately for each genome including the optimal context in each genome for each read. The final step is then used to decide for each read which of these contexts in any of the genomes considered results in the best mapping.

The parallel multi–species mapping is implemented by ContextMap in the following way (Figure 1 A). First, independent Bowtie indices are created for different potential read sources. Separate indices are necessary as Bowtie is limited to –1 characters per index. This is relevant as the human genome alone needs 73% of the maximum index size and all microbial genomes from the NCBI database taken together require 134% of the maximum index size. We, thus, generally use one index for rRNA sequences, one for the host genome, e.g. the human reference genome, one for virus genomes and two for microbe genomes. This can be easily adjusted to more indices as soon as the increasing number of sequenced virus and microbe genomes makes this necessary. After performing the initial alignment against all indices, ContextMap is then run without any further changes to define contexts, the optimal mapping for each read in each context it may belong to and finally the optimal and unique mapping for each read to any context.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Mining RNA–seq data for infections and contaminations.

(A) Approach for mapping sequencing reads in parallel to multiple sources of reads using ContextMap. (B) After obtaining unique mappings to the species in the reference set, different questions can be addressed. Random hits to only a small region of the genome can be identified by investigating coverage. Strong similarities in terms of possible read mappings between different species in the reference set can be identified by analyzing confidence and species clusterings. Finally, by analyzing mismatch distributions in terms of the Jensen–Shannon divergence, it can be determined if reads have been mapped to the correct genome or only to a close relative due to missing genome sequences or local genome similarities.

https://doi.org/10.1371/journal.pone.0073071.g001

In contrast to ContextMap, other RNA–seq mapping tools, which predominantly also use Bowtie, cannot be used for this application as they do not support the use of multiple indices required here due to the size and number of reference sequences and provide no way to distinguish between alternative alignments for a read to two different but related genomes with the same number of mismatches. Thus, they can only be applied sequentially by mapping first all reads e.g. against rRNA sequences, then the unmapped reads against the host reference genome, and then one microbe or virus genome one after the other. However, the latter approach also poses problems as it can lead to different results depending on the order in which genomes are mapped to in case of closely related species or strains.

Analysis of Species Hits

The mapping of reads to reference genomes using any algorithm directly implies a set of species potentially contained in the sample. Please note that in the following we use the term species loosely, in particular in the context of misidentification of species, and it may also refer to a particular strain of a species, represented by a specific genome sequence in the reference database. In particular for bacteria, the distinction between strains and species is not clearly defined and species definition remains a difficult topic. The standard approach is now to use genome sequence differences and a cutoff of 95% average nucleotide identity is often used [22]. However, for species/strains for which no genome sequence is available, nucleotide identity to sequenced species/strains cannot be calculated. Thus, it cannot be determined whether they represent a different species or only a different strain of a species with known genome sequence.

Independent of which mapping algorithm was used to identify species potentially contained in a sample, a number of problems arise that need to be addressed. First, local similarities in the genome of one species not contained in the sample (species A) to a species contained in the sample (species B) may result in reads erroneously mapped to species A and the reporting of this species for the sample. Second, gaps in the reference database may lead to both missing and incorrect hits. If no genome from the species itself or closely related species is contained in the reference database, fast mapping algorithms, including ContextMap, which tolerate only a limited number of sequence differences, will fail to align the corresponding reads. This type of missing species hits is only a minor problem as a slower but more permissive BLAST run applied to unmapped reads may at least detect the infection by identifying more distant relatives of the infecting pathogen. A more severe problem are misalignments in case that genome sequences are only available for closely related species. In this case, reads are incorrectly aligned to these related species, resulting in the identification of wrong species. For instance, in the recent study by Castellarin et al. [13] several Pseudomonas syringae strains, which are plant pathogens, were likely misidentified in samples of colorectal carcinoma.

In the following, several statistics derived from read mappings are described that can be used to address the described problems and confidently identify the species contained in the sample (see Figure 1 B). Coverage and divergence of mismatch distributions can be calculated based on mappings provided by any algorithm. Calculation of species mapping confidence and distances between species relies on the support score calculated by ContextMap for each read mapping, but can be adapted to methods evaluating only the number of mismatches. All methods are available as part of the ContextMap software suite available at http://www.bio.ifi.lmu.de/contextmap.

Data Sets

HPV–18 Expression in HeLa Cells

RNA–seq data of HeLa cells from the study of Guo et al. [25] were mapped using ContextMap against indices for human, viral and microbe genomes and human rRNA. For the initial Bowtie runs a seed of 25 bp was used allowing up to 1 mismatch in the seed, the same settings used by Guo et al. In total, 5 mismatches were allowed, resulting in 11,040,798 (38.4%) and 10,162,289 (34.3%) mapped reads for the mock and miR–155 transfected cells, respectively. This is only 0.4 and 1.8 million reads less than mapped by Guo et al., although they allowed an arbitrary number of mismatches outside the seed, i.e. up to 12 mismatches. Although Guo et al. did not perform alignment against viral or microbial genomes (but also rRNA), only few () of the reads additionally mapped by ContextMap originated from viral or microbial genomes. Most reads additionally aligned by ContextMap were discarded by Guo et al. due to non–unique alignments. Interestingly, million reads originated from rRNA, which illustrates the importance of including rRNA sequences in the mapping process even though poly–A selection was performed.

Tables S1 and S2 show coverage, mapping confidence and compared to the human reference genome for all species with at least 1,000 mapped reads. Figure S2 illustrates coverage for all microbial or virus hits. Here, HPV–18 is the only virus or microbe with a coverage (0.34–0.37), high confidence () and small () in both samples. This confirms previous reports of HPV–18 expression in HeLa cells [12]. In contrast, no reads were mapped to HPV–16, which is not expressed in HeLa cells. Figure 2 A shows the distribution of reads across the HPV–18 genome both in the mock and miR–155 transfected cells. Here, results were highly reproducible between the two samples with peaks in read heights at the same genomic locations. The mapping to genes showed that only the E6, E7 and E1 genes were strongly expressed. In addition, weaker expression by an order of magnitude was observed for L1 as well as for a region covering the end of E1 and the start of E2. However, as no reads were observed for the rest of E2, it is likely not expressed. The same was true for genes E4, E5 and L2. These observations are in accordance with recent results showing that the oncogenes E6 and E7 are essential for continued proliferation in cervical carcinoma [28]. Both genes are transcriptionally repressed by the E2 protein and loss of E2 expression leads to upregulation of E6 and E7 [29]. Thus, loss of E2 expression in HeLa cells as well as high E6 and E7 expression is consistent with their origin from cervical carcinoma cells and ongoing proliferation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Characterization of HPV–18 infection in HeLa cells.

(A) Distribution of reads across the HPV–18 genome for the mock and miR–155 transfected cells. Read numbers are shown in log scale. Expressed genes include E1 as well as E6 and E7, which are required for ongoing proliferation in cervical carcinoma [28]. L1 also appeared to be weakly expressed, however the expression pattern did not exactly correspond to the annotated gene coordinates. While the start of the gene was not expressed, L1 expression was extended to a region downstream of the gene. (B) Coverage as a function of increasing sequencing depth was evaluated by randomly sampling from the miR–155 data set. Coverage is shown as an average of ten repeated samplings for HPV–18 (black) and other species (gray). Sample size is annotated to the HPV–18 data points.

https://doi.org/10.1371/journal.pone.0073071.g002

This shows that our approach is capable of identifying HPV–18 infection in HeLa cells and distinguishing this from spurious matches to other species. However, as less than 1% of reads in our samples originated from HPV–18 (22,105 and 18,491, respectively), the question remains which sequencing depth is necessary for confidently identifying such an infection. To investigate this question, we randomly sampled reads from the miR–155 data set with sample sizes between and (see Figure 2 B). For each sample size, 10 random repetitions were performed and reads were mapped using ContextMap as described. Here, a sequencing depth of as low as 500,000 reads (1.7% of all reads) was sufficient to clearly distinguish the HPV–18 infection from spurious hits to other species. Although only 303 HPV–18 reads were identified on average at this sample size, almost all of these reads (90%) were mapped to distinct genome positions, resulting in a coverage of . Although this coverage is small, it is more than an order of magnitude larger than for any of the other species at this sequencing depth and increases much faster with increasing sequencing depth.

To compare the proposed method against alternative approaches, we performed megablast alignments for the miR–155 data set against all microbial and viral genomes as well as human rRNA sequences and the human mitochondrial genome. Alignments with an E–value were then evaluated using MEGAN4. A megablast comparison against the complete human genome was aborted as output files already reached 10GB after mapping only 28% of reads against 30% of the genome, which would have resulted in an estimated 120GB of output (for an input of only 0.84GB). Since GASiC and GRAMMy could only be run in reasonable time on the –fold smaller in–vitro microbial community data set by restricting them to the 10 species in question, we did not evaluate them here. MEGAN4 results are shown in Figure S3 both with and without the additional alignment against human rRNA and mitochondrial genome sequences. In both cases, HPV–18 is clearly detected although 762 fewer reads (4.1%) are assigned than by ContextMap despite the fact that an arbitrary number of mismatches and gaps are allowed by BLAST.

However, without additionally BLASTing against human rRNA and the mitochondrial genome, reads each are assigned to one bacterial (Rickettsia rickettsii str. Hino) and one viral (Choristoneura occidentalis granulovirus) species. When including human sequences for mapping, most of these are assigned to the inner nodes “cellular organisms” and “root”, reflecting sequence similarities between human rRNA and the Rickettsia genome (lowest common ancestor (LCA) = “cellular organisms”) and the human mitochondrial genome and the Choristoneura genome (LCA = “root”), respectively. These results show the importance of also including the host species into mapping, as otherwise Rickettsia and Choristoneura would be reported erroneously for this sample. Here, MEGAN4 provides no direct way for identifying these hits as suspicious, e.g. by calculating coverage or mismatch distributions, or for resolving the non–uniquely mapped reads assigned to inner nodes. In contrast, ContextMap correctly assigns 90% of the Choristoneura BLAST hits to human rRNA and only 5% to Choristoneura. Furthermore, 88% of the Rickettsia BLAST hits are correctly identified as originating from human RNA (83% from mitochondrial RNA) by ContextMap and only 1% are assigned to Rickettsia. In addition, the few Choristoneura and Rickettsia reads assigned by ContextMap are clearly flagged as misalignments by very high values of ().

The Microbiome of Colorectal Carcinoma

In the second analysis, we focused on the RNA–seq data of colorectal carcinoma and adjacent normal tissue from the study of Castellarin et al. [13]. This data set was interesting as they identified a Fusobacterium to be enriched in colorectal carcinoma cancer. In addition, they reported a number of microbes that are unlikely to occur in colon tissue, e.g. Pseudomonas fluorescens SBW25, which was found at high levels in all samples, and two Pseudomonas syringae strains. P. fluorescens is mostly found in soil and water, whereas P. syringae are plant pathogens. While in the first case occurrence in colon samples might still be possible, e.g. due to contamination, in the latter case it is very unlikely. Although mapping with ContextMap also identified all three Pseudomonas species in all tumor and normal tissue samples, compared to the human reference genome was larger than 0.2 in all three cases (Figure S4), in particular for Pseudomonas syringae pv. syringae where more than half of the reads had at least 3 mismatches (). This indicates that the actual Pseudomonas species contained in the sample is not yet sequenced, resulting in reads from these species mapped to a number of related Pseudomonas species.

Based on these observations, we performed the same analysis for Fusobacterium. Previously, Castellarin et al. identified Fusobacterium nucleatum subsp. nucleatum as overrepresented in the RNA–seq data of the tumor tissues. Subsequent DNA sequencing of a Fusobacterium culture isolated from the tumor samples and mapping of reads against additional Fusobacterium draft genomes from the Human Microbiome Project (HMP), however identified Fusobacterium sp. 3_1_36A2 as a much better match than F. nucleatum. As F. sp. 3_1_36A2 was extracted from the colon of a patient, this makes more sense than F. nucleatum, which was isolated from the human oral cavity and is most commonly found there.

To re–capitulate their analysis, we performed mapping of tumor samples using ContextMap both with and without the human microbiome in addition to the RefSeq genomes. Without the human microbiome index, F. nucleatum was identified in all tumor samples, in particular in samples from patient 1 ( reads, Figure 3). However, the mismatch distribution differed considerably from the mismatch distribution for the human genome ( for patient 1), clearly indicating that F. nucleatum subsp. nucleatum is not contained in the sample but only a related species. Indeed, when performing mapping including the human microbiome, almost all of the reads originally mapped to F. nucleatum are mapped to contigs of other Fusobacteria species, such as sp. 3_1_33, sp. 11_3_2, sp. D11, sp. 7_1, and sp. 21_1A, which were isolated from biopsy tissues from the gastrointestinal tract. Furthermore, the primer sequences used by Castellarin et al. to confirm the presence of Fusobacterium are a better match to these species, with both primers matching with at most 2 mismatches, whereas for Fusobacterium nucleatum one primer has 3 mismatches.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Comparison of mismatch (mm) distributions for Fusobacteria

. Results are shown for species identified by ContextMap with at least 20 reads on the colorectal carcinoma samples for patient 1 using either only the completed microbe genomes as reference set (left) or also the human microbiome draft genome sequences (right). Distributions are compared against the average mismatch distribution for the human genome. Number of reads mapped to each genome and are indicated in parentheses.

https://doi.org/10.1371/journal.pone.0073071.g003

Among the identified Fusobacteria, F. sp. 3_1_33 has the highest number of reads for patient 1 () and smallest (0.073). It is also enriched in the tumor sample compared to the normal tissue, but not as strongly as some other species from the HMP with fewer mapped reads, in particular some E. coli strains (Table S3). Although a comparatively small number of reads (1,372) are still assigned to F. nucleatum even with the inclusion of the human microbiome, the mismatch distribution still diverges strongly from the human reference () and is unusual in that it has a higher number of reads both with zero and with four mismatches. Together with the observations that for F. sp. 3_1_33 is still higher than in the HPV–18 example, a number of other Fusobacteria are also found with substantial read numbers, and most of the F. sp. 3_1_33 reads can be aligned almost equally well to the other gastrointestinal Fusobacteria, this suggests that F. sp. 3_1_33 is also not the actual strain in the sample. However, it appears to be a much closer relative than F. nucleatum. This also shows that should always be analyzed in combination with read numbers as local similarities may allow the mapping of some reads to a wrong species with few mismatches.

To compare our results against other approaches, we extracted all reads for the tumor sample of patient 1 that were not mapped to human sequences (including rRNA) by ContextMap (404,234 reads) and performed both megablast and novoalign (http://www.novocraft.com) alignments for these reads against virus and microbe genomes and the human microbiome. Novoalign was used by Castellarin et al. to align reads to the bacterial and viral genomes after filtering out all reads that could be aligned to human rRNA, cDNA or the reference genome using BWA [30], a fast short read aligner applying a similar strategy as Bowtie. Thus, we effectively recapitulated their analysis in our study, this time also including the human microbiome. Again MEGAN4 was applied to the BLAST output as shown in Figure S5. Almost all () of the Fusobacteria reads could be aligned to more than one Fusobacterium, thus, resulting in an assignment of these reads to their LCA by MEGAN4. In addition, MEGAN4 allows no further analysis as to which of the identified Fusobacteria is the most likely candidate or closest relative of the species or strain contained in the sample.

Novoalign was applied in two modes: one outputting all alignments for a read with the same maximum score and one outputting only unique alignments (Figure 4). To compare against the ContextMap results and calculate compared to the reads mapped to human by ContextMap, we then extracted only those alignments with at most 5 mismatches and no gaps. Please note that read numbers were hardly increased for the Fusobacteria if gaps were allowed. For evaluation of the novoalign mode allowing multiple alignments, we used only one of the best alignments for each read for each genome, but allowed multiple alignments with equal score to different genomes. Here, almost all reads could be aligned equally well to more than one genome with % unique read alignments per genome. Although reads were still aligned to F. nucleatum, only 44 of these were unique and % were aligned equally well to F. sp. 3_1_33. Again, this illustrates the problems similarities between sequenced genomes present for mapping algorithms that are based only on the individual read alignments. Without taking into account alignments of other reads, they may only either completely exclude or include non–unique alignments. In this application, a restriction to unique alignments would vastly underestimate Fusobacterium expression in the sample, whereas the inclusion of non–unique mappings would result in the reporting of essentially all of the identified Fusobacterium species. In this case, evaluation of is not meaningful as due to the multiple mappings the sets of reads assigned to each species and, consequently, the calculated mismatch distributions are very similar.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Number of reads and mismatch distributions for the novoalign mapping on the Fusobacteria.

Results are shown for species identified by aligning with novoalign against viral and microbial genomes and the human microbiome for the patient 1 colorectal carcinoma sample. Only reads were used that were not mapped to human sequences by ContextMap. Mismatch distributions are compared against the average mismatch distribution for the human genome derived from the ContextMap mapping. Number of reads mapped to each genome and are indicated in parentheses. The left–hand side shows results if multiple read alignments with the same maximum score to different species are allowed. The right–hand side shows the results for unique alignments only.

https://doi.org/10.1371/journal.pone.0073071.g004

Meta–transcriptomics for an in–vitro Simulated Microbial Community

For the final analysis, we analyzed DNA sequencing data for an in–vitro simulated microbial community by Morgan et al. [26] and compared our results against several state–of–the–art metagenomics tools, in particular MEGAN4 [19], GRAMMy [20], and GASiC [21]. This data set was selected as the species contained in the samples were known. Furthermore, it constituted a challenging application due to strong similarities of the genomes of the microbial strains contained in the sample to other sequenced genomes. One example for this is Halobacterium sp. NRC–1, whose genome is almost identical to the Halobacterium salinarum R1 genome [31]. They differ only by 4 base changes, 5 single–nucleotide indels and 3 longer indels between 133 and 10,007 bp long.

We investigated the performance of ContextMap on this data set using a reference containing the yeast genome and all microbial and viral genomes downloaded from NCBI (see methods) and allowing 5 mismatches. To compare our results against BLAST as well as MEGAN4 and GRAMMy, which use BLAST alignments as input, we performed megablast searches against the same genomes and extracted all alignments with the maximum score for each read, using only alignments without gaps and at most 5 mismatches. Here, 12% of reads could be aligned equally well to at least two different RefSeq entries using BLAST. In addition, we applied GASiC to all genomes from the same genus as any of the species contained in the sample (122 RefSeq entries, 92 taxa). The same restriction was applied to GRAMMy as both methods already took more than 7 CPU hours on this smaller set compared to min for ContextMap on all microbes and viruses (Table S4).

Table S5 lists the microbe and virus hits identified by ContextMap with a coverage and at least 20 reads. Here, ContextMap identified all of the microbial species contained in the sample, but also several related strains and prophages. As Myxococcus xanthus had the highest number of mapped reads, we used it as reference for calculating . Interestingly, all microbes that are contained in the sample had a higher mapping confidence than all other hits despite low numbers of reads for some of them.

For five species, identification is straightforward based on this list. A. cellulolyticus, S. amazonensis, L. brevis, and M. xanthus are characterized by high mapping confidence (), low () and high number of reads and coverage. For P. pentosaceus confidence is also high and still relatively low (0.064), but coverage is quite small (). However, as 90% of the reads map to distinct positions, it is clearly a correct hit and the low coverage is likely due to low abundance of P. pentosaceus in the simulated community. In the clustering of species hits, these five species also form distinct clusters with no similarities to any of the other species hits (Figure 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Hierarchical clustering (average linkage) of microbes and viruses

. Results are shown for hits with a coverage and at least 20 mapped reads as determined by ContextMap. Microbes actually contained in the sample are indicated in red and by an asterisk and the three clusters discussed in the text are marked by rectangles. In addition, number of reads, confidence and are indicated next to the microbe names.

https://doi.org/10.1371/journal.pone.0073071.g005

For the remaining hits the situation is less clear–cut. Here, clustering identified three large groups among these: (1) a Halobacterium cluster, (2) a Lactococcus cluster, and (3) a Lactobacillus cluster. In the first case, H. sp. NRC–1 clusters tightly with H. salinarum R1 and the plasmids also cluster together, reflecting the small number of sequence differences between these. In all of these cases, is relatively small (), with the only exception being H. sp. NRC–1 for which all 50 reads have zero mismatches, i.e. fewer mismatches on average (Figure S6). Additionally, confidence for NRC–1 (0.37) is much higher than for R1 (0.07). Thus, although significantly fewer reads are mapped to NRC–1, all other mapping statistics support the presence of the NRC–1 strain rather than the R1 strain. However, since both strains are almost identical, is still very low for the R1 strain (0.02).

In the second cluster, read numbers, confidence and clearly indicate the presence of L. lactis subsp. lactis Il1403. Although (0.07) is somewhat increased compared to the microbes identified unambiguously, it is not yet large enough to question this hit. For all other Lactococcus lactis strains, in particular cremoris SK11, the number of mismatches is significantly increased, indicating that these are not contained in the sample. This is surprising as SK11 was part of the community. The reason for this is that 99% of the reads potentially mapping to SK11 can be aligned equally well to other species, in particular to L. lactis subsp. lactis Il1403. As the latter is more abundant, it ends up with most of the reads, apart from those with too many mismatches. Finally, analysis of the last cluster confirms the presence of L. casei ATCC 334 as it is characterized by high coverage and confidence and sufficiently low (0.06). The other strains in the cluster, in particular the BL23 and Zhang L. casei strain, can be clearly excluded due to high () and low coverage () and confidence ().

In summary, these results show that ContextMap can be used to correctly identify all species in the community including the strain, with the exception of cremoris SK11. However, analysis of results for MEGAN4 (Figure S7), GASiC (Table S6) and GRAMMy (Table S7) shows that none of these identify cremoris SK11, at least not with more confidence than for other species/strains not contained in the community. MEGAN4 assigns almost all of the cremoris SK11 reads to the LCA of the cremoris and lactis subspecies. GASiC assigns a p–value of 1, i.e. considers it an insignificant hit. Finally, GRAMMy, which only estimates relative abundances but performs no read mapping, assigns an abundance of %, less than assigned to L. lactis subsp. lactis KF147 (0.17%), which is not part of the community.

Apart from cremoris SK11, GASiC fails to identify H. sp. NRC–1 and P. pentosaceus but otherwise predicts only microbes contained in the community. Thus, GASiC is the most restrictive of the analyzed approaches. MEGAN4 does not really resolve multiple mappings but assigns reads with multiple mappings to the LCA of these microbes. Based on the number of reads assigned uniquely to any of the children of such an LCA, the correct microbes can then be predicted. In this example, the predictions would be correct with the exception of SK11 and H. salinarum, where the NRC–1 strain cannot be properly distinguished. Nevertheless, even when combining read numbers for the microbes and the LCA, ContextMap generally identifies 2.5–6.5% more reads per microbe (including plasmids). Furthermore, assignments to inner nodes of the phylogenetic tree by MEGAN4 do not allow calculation of mismatch distributions as corresponding genome sequences are not known. Even for the leaves of the taxonomic tree, additional statistics of alignment quality or coverage are not directly accessible and can only be obtained by extracting the assigned reads and performing this analysis using additional scripts.

GRAMMy correctly identifies 7 of the 9 microbes with an estimated abundance of %, but also assigns a very low abundance to P. pentosaceus (0.4%). Remarkably, the relative frequency estimated by GRAMMy and the coverage calculated by ContextMap are highly correlated (correlation coefficient 0.995), in particular for microbes with high coverage (Figure 6). This indicates that coverage as determined by ContextMap provides a reliable estimation of the relative frequencies identified by GRAMMy. As ContextMap is much faster than GRAMMy, it can thus be used to replace GRAMMy for applications where GRAMMy is too inefficient.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Comparison of abundance calculated by GRAMMy and coverage determined by ContextMap on the microbial community data set.

Results are shown for all taxa identified by GRAMMy with a relative abundance of at least 0.1%. The green line indicates a linear fit to the data.

https://doi.org/10.1371/journal.pone.0073071.g006

We also evaluated a number of other metagenomics tools for binning/classifying sequencing reads or identifying relative abundance of species. This includes alignment–based approaches (MG–RAST [32], MetaPhyler [33], SOrt–ITEMS [34], MARTA [35], MLTreeMap [36]), composition–based approaches (PhyloPhytiaS [37], ClaMS [38], Phymm [39]) and a hybrid approach (PhymmBL [39]). Here, comparison of the results was difficult as several approaches only perform classification at the genus– (MetaPhyler) or species–level (MG–RAST, PhyloPythiaS, MARTA), but do not identify individual strains. Thus, we could not evaluate their performance in distinguishing the Halobacterium and Lactococcus lactis strains. Furthermore, only MG–RAST, MetaPhyler, Phymm and PhymmBL were developed for NGS reads as short as 100 bp, while the other tools require longer reads. Thus, the meaningfulness of the comparison against these other approaches is limited. Results for all tools are shown and discussed in Tables S8 (MG–RAST), S9 (MetaPhyler), S10 (SOrt–ITEMS), S11 (MARTA), S12 (MLTreeMap), S13 (PhyloPhytiaS), S14 (ClaMS) and S15 (Phymm and PhymmBL). In summary, although the correct species or at least genera were usually identified, performance at the level of the individual strains was usually poor as often wrong strains were ranked higher than strains contained in the community. A particular poor performance was observed for the composition–based approaches PhyloPhytiaS and ClaMS, which likely suffered from the short sequencing read length.

The analysis of runtime and memory requirements on this data set (Table S4) showed that ContextMap was both faster than almost all other approaches (apart from MetaPhyler) and required less or a comparable amount of memory (with the exception of MLTreeMap and GRAMMy if memory requirements of the BLAST run to provide the input for GRAMMy are not counted). The comparison for the other two data sets was less informative as ContextMap was either applied on more reads as in the case of the colorectal carcinoma data set or more reference sequences as in the case of the HeLa cell RNA–seq data. In the first case, ContextMap took sec per read on the complete 5,343,842 read set for patient 1, whereas BLAST took sec per read on the smaller 404,234 read set without human reads and novoalign only sec per read. Thus, ContextMap was much faster than BLAST but slower than novoalign. In contrast, ContextMap required much less memory with G for the complete 5 million read set compared to the G required by novoalign for only 400 k reads and G required by BLAST. This large memory requirement (and also runtime) of BLAST for this relatively small data set was rather remarkable, in particular in comparison to the similar–sized microbe data set and the much larger HeLa data set, which both only required G. The reason for this is the much larger number of possible hits per read found by BLAST in the colorectal carcinoma data set ( hits per read) as compared to the microbe ( hits per read) and HeLa data set ( hits per read). It should be noted that in the latter case we did not perform simultaneous mapping with BLAST against the complete human genome as the estimated output size was too large. Thus, these results suggest that a complete run of BLAST on the same references as used for the ContextMap run would have resulted in substantially increased runtime and memory requirements compared to ContextMap.