WebCARMA: a web application for the functional and taxonomic classification of unassembled metagenomic reads
© Gerlach et al. 2009
Received: 15 September 2009
Accepted: 18 December 2009
Published: 18 December 2009
Skip to main content
© Gerlach et al. 2009
Received: 15 September 2009
Accepted: 18 December 2009
Published: 18 December 2009
Metagenomics is a new field of research on natural microbial communities. High-throughput sequencing techniques like 454 or Solexa-Illumina promise new possibilities as they are able to produce huge amounts of data in much shorter time and with less efforts and costs than the traditional Sanger technique. But the data produced comes in even shorter reads (35-100 basepairs with Illumina, 100-500 basepairs with 454-sequencing). CARMA is a new software pipeline for the characterisation of species composition and the genetic potential of microbial samples using short, unassembled reads.
In this paper, we introduce WebCARMA, a refined version of CARMA available as a web application for the taxonomic and functional classification of unassembled (ultra-)short reads from metagenomic communities. In addition, we have analysed the applicability of ultra-short reads in metagenomics.
We show that unassembled reads as short as 35 bp can be used for the taxonomic classification of a metagenome. The web application is freely available at http://webcarma.cebitec.uni-bielefeld.de.
Metagenomics is a new field of research on natural microbial communities strongly enhanced by the new high-throughput sequencing (HTS) technologies like Roche's 454-sequencing, ABI's SOLiD or Illumina's Genome Analyzer. In traditional genomics, the sequencing of new microbial genomes was limited to those microbes that can be cultivated in a monoculture, which constitute less than 1 percent of all microbes. To explore the other 99 percent, metagenomics makes it possible to retrieve information about the functional and taxonomic composition of the community they originate from .
With decreasing sequencing costs, shotgun metagenomics, i.e. sequencing of random genomic DNA fragments from natural microbial communities, has become feasible [2–4]. Comparing such metagenomic sequences with sequences of known function makes it possible to analyse the biological diversity and the underlying metabolic pathways in microbial communities.
To infer the taxonomic origin of metagenomic reads, two kinds of methods, composition-based and comparison-based, can be distinguished. The composition-based methods extract sequence features, like GC content or k-mer frequencies, and compare them with features of reference sequences with known taxonomic classifications [5–7]. In detail, different techniques, like the calculation of correlation coefficients between oligonucleotide patterns , self-organizing maps (SOMs) , or support vector machines (SVMs)  can be used to classify the metagenomic fragments. The comparison-based methods rely on homology information obtained by database searches, e.g. using search tools like BLAST[11, 12]. An additional challenge is the fact that the new sequencing techniques produce short (100-500 bp with 454) and ultra-short (35 bp with SOLiD, 35-100 bp with Illumina, 30-35 bp with Helicos ) reads. New bioinformatic tools have to be developed that can cope with both, the huge amount of data and the short read lengths. Especially the short read lengths have been considered the main bottleneck for the usage of ultra-short reads in metagenomics .
In 2008, we developed CARMA , a new software pipeline for the characterisation of species composition and the genetic potential of microbial samples using short reads. In contrast to the traditional 16S-rRNA approach for taxonomical classification [16–18], CARMA uses reads that encode for known proteins. By assigning the taxonomic origins to each read, a profile is constructed which characterises the taxonomic composition of the corresponding community. Krause et al.  showed on a synthetic metagenome that even with reads as short as around 100 bp, high accuracy predictions with an average false positive rate of 0.1 to 2.5 percent are possible. The CARMA pipeline has already been successfully applied to 454-sequenced communities including the characterisation of a plasmid sample isolated from a wastewater treatment plant and other communities [19–22].
Here we introduce WebCARMA, a refined version of CARMA available as a web application for the taxonomic and functional classification of unassembled (ultra-)short reads from metagenomic communities. The web application is freely available at http://webcarma.cebitec.uni-bielefeld.de. Furthermore, we examined the applicability of CARMA for simulated ultra-short reads (≥ 35 bp) by comparing the results with earlier ones obtained by short reads (≈230 bp) using samples from a biogas plant .
WebCARMA is based on a new CARMA version 2.1 with some improvements over the last published version 1.2 . We first review the basic concepts behind CARMA, before we present our new WebCARMA.
Next-generation sequencing technologies in metagenomics produce millions of DNA reads, from which CARMA detects those that encode for known proteins. These reads are called environmental gene tags (EGTs) . In a second step they are assigned to taxa from six taxonomical ranks: superkingdom, phylum, class, order, genus and species. The set of classified EGTs provides a taxonomical profile for the microbial community.
In detail, BLASTX [24, 25] is used to search within the set of reads for candidate EGTs that encode for protein sequences contained in the Pfam database . A rather relaxed E-value of 10 and frameshift option -w 15 are used. Each read that has a match to a protein family member is translated according to BLASTX reading frame and frameshift predictions. The final determination of EGTs is done by matching the candidate EGTs against their matching protein families with the corresponding Pfam Hidden Markov Models . For this purpose mmpfam from the HMMER package  is used. Only candidate EGTs with an hmmpfam E-value match of 0.01 or lower are accepted as EGTs.
After the EGTs are identified, they are taxonomically classified. Therefore, each EGT is aligned against the multiple alignment of its family with hmmalign (also contained in the HMMER package). From this new alignment, the pairwise sequence distance is computed for all pairs of sequences, based on the fraction of identical amino acids. This produces a pairwise distance matrix which then is used to compute a phylogenetic tree with the neighbour-joining clustering method . The EGT is then classified depending on its position within this tree. If the EGT is localised within a subtree of family members all sharing the same taxon, then the EGT is classified with the same taxon. For example, if the EGT is localised in a subtree with the three members Bacteria Cyanobacteria Synechococcales Prochlorococcus, Bacteria Cyanobacteria Chroococcales Synechococcus and Bacteria Cyanobacteria Nostocales Nostoc, the EGT is classified as Bacteria Cyanobacteria. For a more formal definition and further details see .
We have reimplemented large parts of CARMA, including a faster construction of the phylogenetic trees by caching the pairwise distances between Pfam family members. The CARMA results now include for each EGT the corresponding hmmpfam E-values and a list of GO-Ids (Gene Ontology Identifiers)  associated with the corresponding Pfam family. The Gene Ontology provides a controlled vocabulary for gene products, distinguishing between their associated biological processes, cellular components and molecular functions, and can therefore be used to create a functional profile of the metagenome. A major modification of CARMA is the usage of the NCBI taxonomy database [31, 32] instead of the Pfam nomenclature. The NCBI taxonomy database currently indexes over 200,000 species , which are classified in a hierarchical tree structure. Each taxon from the taxonomy is represented as a node in the tree with a unique identifier (tax_id) and its taxonomic rank ranging from superkingdom to subspecies. For compatibility with other applications and databases, the output files of CARMA contain, along with the prettyprint taxa classifications, the corresponding tax_id.
The WebCARMA website is built upon an Apache Web Server using Perl and CGI. The CARMA pipeline is executed on the compute cluster of the Bielefeld University Bioinformatics Resource Facility at the Center for Biotechnology (CeBiTec) using Sun Grid Engine http://gridengine.sunsource.net/. In order to use WebCARMA and to upload metagenomic sequences, the users have to register with their e-mail address. After the uploads are finished, CARMA starts with the search for EGTs and the taxonomical classification. By the time the jobs are completed, the user receives an e-mail with a download address pointing to the results.
We provide several pre-computed data formats that allow the user to explore the results in different ways:
The translated EGTs, with additional information about the name of the original metagenomic sequence, reading frame, Pfam ID, HMMER E-value and a list of GO-Ids associated with that corresponding Pfam ID.
A GO-term profile in two variants, as a text data file and visualised as a histogram in PDF format.
The taxonomical classification results as a text data file.
A taxonomic profile, once as a text data file and once visualised by histograms for each taxonomic rank in PDF format.
The profile data files as well as the classification results are provided in TSV-format (Tab Separated Values), which makes it easy to import the data into other programs (e.g. spreadsheet) for different visualisation types or any other further processing.
We provide Perl scripts for download that can easily be used as templates for own data processing pipelines. A manual with further explanations can be found on the WebCARMA site.
WebCARMA produces functional (using GO-Terms) and taxonomic profiles (using NCBI taxonomy) of metagenomic sequences, available in text format and as histograms.
To classify a metagenomic data set as used in the use case study below (600 000 reads and average 230 bp read length) takes about 1 800 CPU hours (each CPU with 2.5 GHz) on our compute cluster. Because the neighbour-joining algorithm has quadratic runtime, cluster jobs which process bigger Pfam families will take much more time to compute than other. Therefore, using 100 CPUs, it takes for this data set about 40 hours until the last job has finished.
In order to evaluate the applicability of short and ultra-short reads in metagenomics, we used a real-world data set and several realistic simulated data sets. The real-world data set consisted of 600 000 reads with an average read length of 230 bp of a microbial community from an agricultural biogas reactor , sequenced with the 454 Genome Sequencer FLX system (Roche Applied Science). Instead of using real ultra-short reads, we decided to simulate ultra-short reads by clipping off suffixes of the 454-reads to get the desired read lengths. We generated nine data sets, each consisting of reads of one of the lengths 35 bp, 40 bp, 50 bp, 60 bp, 70 bp, 100 bp, 150 bp, 200 bp, and 250 bp, respectively.
It is known that some sequencing techniques exhibit correlations between read coverage and GC content [36–38]. By using simulated reads instead of real ultra-short reads we can be sure that any differences we see in the classification results between the data sets are only due to the different read lengths. If there is a bias in the 454 data, then we also have the same bias in the simulated data sets and our comparison should not be much affected by this.
First, we analyse the number and lengths of the EGTs obtained for each data set, then we compare the taxonomic classification results for the different read lengths.
Number of reads in each data set.
In rare cases, it happens that an EGT is one amino acid longer than its read length divided by three. For example, in the set of EGTs produced from the 150 bp-reads data set, the longest EGTs are 51 amino acids long. This can occur when blastx predicts two frameshifts in one read.
For our analysis of the applicability of ultra-short reads with CARMA, we considered seven different taxonomic ranks: superkingdom, phylum, class, order, family, genus and species. A first relative abundance for each taxon is obtained by dividing the absolute number of EGTs that predict this taxon by the total number of EGTs at the same taxonomic rank. The latter do not include EGTs which predict nothing ("unknown"). We consider taxa with a relative abundance below the threshold 0.015 in all data sets, to be false positives. Therefore they are classified as "other".
After applying the threshold we recompute the relative abundances for each taxon, this time subtracting both, "unknown" and "other" from the total number of EGTs at the same taxonomic rank.
With this, we have normalised the relative abundances for the taxa such that they sum up to 1 and therefore ensured comparability between the data sets.
Rate of "Unknown" EGTs. Rate of "Unknown" EGTs that could not be classified further from the complete set of EGTs.
Rate of "Other" EGTs.
Even though the taxonomic predictions on lower taxonomic ranks (order, family, genus and species) are known to be imprecise, we included them in our experiment in order to study the effect of using (ultra-)short reads compared to longer ones at all taxonomic ranks.
The results show that CARMA predicts for all data sets and all taxonomic ranks the same taxa. For higher taxonomic ranks, even the relative abundance levels are similar between the different data sets. Deviations of 35 bp-reads for example can be seen on the level of order, where significantly more of Thermotogales and Haemosporida, and less of Thermoanaerobacterales are predicted. The 40 bp data set does not show these differences. Even more deviations can be found on lower taxonomic ranks, for example species.
Furthermore, as expected, the rate of EGTs which are not classified increases for lower taxonomic ranks for all data sets (Table 2). Interestingly, the rate of unclassified EGTs is smaller for shorter reads than for longer reads. This might be due to the circumstance that shorter EGTs need to have more sequence similarity to the Pfam families than longer EGTs, in order to achieve the same E-value threshold.
CARMA is a software pipeline for the characterisation of species composition and the genetic potential of microbial samples using short reads that encode for known proteins. The species composition can be described by classifying the reads into taxonomic groups of organisms they most likely stem from. CARMA has been successfully evaluated on a synthetic metagenome  and has already been used for the analysis of several microbial communities [19, 20].
Here we have presented our new web application WebCARMA, which makes metagenomic analyses with CARMA easily accessible. WebCARMA provides functional and taxonomic classifications in common data formats as well as basic visualisations of the profiles.
Previous metagenomic analysis relied on reads of length 100 bp or longer. We have shown that ultra-short reads as short as 35 bp can be used for the taxonomic classification of a metagenome. The biogas data set we have used in the analysis is a low complexity data set with only a few prevalent species. Therefore, our results do not necessarily apply to data sets of higher complexity. Still, we think we have shown that ultra-short reads can indeed, in principle, be used for reliable taxonomic classification of a microbial community if the coverage is high enough. We have found most differences between the different data sets in the taxa of higher order, e.g. at the species level, and in general for species with very low abundance.
Metagenomics with CARMA still leaves some room for improvements. Proper statistics to assess the significance of functional and taxonomic predictions based on short reads (Kowalczyk et al.: Significance Tests for Short Read Concentrations, unpublished manuscript) are still missing. The abundance levels of the classification results have to be read with care. Species with larger genomes or more genes than other species will be overrepresented in the taxonomic profiles because more EGTs can be found. Therefore, more accurate results might be achieved by weighting EGTs using additional information like the genome size of the closest known relative.
Project name: WebCARMA
Project home page: http://webcarma.cebitec.uni-bielefeld.de
Operating system(s): Platform independent (web-service); Unix/Linux (stand-alone program)
Programming language: Perl
Other requirements: none
License: GNU GPL
Any restrictions to use by non-academics: none
Upload restriction: Maximum 100 MB of FASTA per month
This is the complete set of figures for the evaluation from the use case study. Only those taxa are visualized, where in at least one of both data sets the predicted relative frequency was 0.015 or higher. Taxa that are not visualised contribute to "Other".
WG was supported by the DFG Graduiertenkolleg 635 Bioinformatik, SJ by the BMBF PathoGeno-Mik Plus network project 0313801N and FT by the BMBF ERA-NET PathoGenoMics project 0313939A.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.