Analysis and comparison of very large metagenomes with fast clustering and functional annotation
© Li. 2009
Received: 16 March 2009
Accepted: 28 October 2009
Published: 28 October 2009
Skip to main content
© Li. 2009
Received: 16 March 2009
Accepted: 28 October 2009
Published: 28 October 2009
The remarkable advance of metagenomics presents significant new challenges in data analysis. Metagenomic datasets (metagenomes) are large collections of sequencing reads from anonymous species within particular environments. Computational analyses for very large metagenomes are extremely time-consuming, and there are often many novel sequences in these metagenomes that are not fully utilized. The number of available metagenomes is rapidly increasing, so fast and efficient metagenome comparison methods are in great demand.
The new metagenomic data analysis method Rapid Analysis of Multiple Metagenomes with a Clustering and Annotation Pipeline (RAMMCAP) was developed using an ultra-fast sequence clustering algorithm, fast protein family annotation tools, and a novel statistical metagenome comparison method that employs a unique graphic interface. RAMMCAP processes extremely large datasets with only moderate computational effort. It identifies raw read clusters and protein clusters that may include novel gene families, and compares metagenomes using clusters or functional annotations calculated by RAMMCAP. In this study, RAMMCAP was applied to the two largest available metagenomic collections, the "Global Ocean Sampling" and the "Metagenomic Profiling of Nine Biomes".
RAMMCAP is a very fast method that can cluster and annotate one million metagenomic reads in only hundreds of CPU hours. It is available from http://tools.camera.calit2.net/camera/rammcap/.
The emerging field of metagenomics enables a more comprehensive understanding of environmental microbial communities [1–9]. However, metagenomic data consists of enormous numbers of fragmented sequences that challenge data analysis methodologically and computationally. To address these challenges, new methods and resources have been developed, such as simulated datasets, IMG/M, CAMERA, MG-RAST, taxonomy tools[14, 15], statistical comparison, functional diversity analysis, binning [18–20] and so on.
The Rapid Analysis of Multiple Metagenomes with a Clustering and Annotation Pipeline (RAMMCAP) presented herein aims to address the particular computational challenges imposed by the huge size and great diversity of metagenomic data. The primary goal is to significantly reduce the computational effort in sequence comparison, as large-scale comparison of metagenomic sequences has become extremely time-consuming. For example, the protein analysis of the Global Ocean Sampling (GOS) study took more than one million CPU hours.
Metagenomic datasets may contain many novel genes that don't show any homology to existing genes. For example, only ~10% of the sequences in the "Metagenomic Profiling of Nine Biomes" (BIOME) study  match known functional genes. Novel genes in metagenomic datasets have not been used in many studies with homology-based gene prediction and analysis, so the second goal of RAMMCAP is to explore whole datasets and make use of the novel sequences. Because the ab initio gene finding approaches developed for complete genomes work poorly with fragmented DNA sequences, recently, several new gene prediction methods were developed for short DNA sequences with high sensitivity and specificity, such as Metagene, MetageneAnnotator, and Neural Networks. In RAMMCAP, ORFs are called with either Metagene or simple six reading frame translation; both methods can identify novel genes.
Since more and more metagenomes will be available in the future, the third goal of RAMMCAP is to provide a new way to compare metagenomes from various environmental conditions and to identity and visualize the statistically significant differences between metagenomes.
In this paper, RAMMCAP was implemented and applied to the two largest metagenomic collections. The first set, GOS [1, 2], features 7.7 million ~800 base Sanger reads from 44 samples. A second, the Biomes  set, has 14.6 million ~100 base 454 reads from 45 microbiomes and 42 viromes samples. With moderate computational effort, RAMMCAP can quickly analyzed these huge datasets and obtained many novel results that could not be achieved by other existing methods.
ORFs are collected from sequence reads with ORF_finder, a ORF calling program implemented here by six reading frame translation in a similar way as the GOS study. Within each reading frame, an ORF starts at the beginning of a read or the first ATG after a previous stop codon; it ends at the first stop codon or the end of that read. The minimal length of ORFs can be specified by users. ORFs can also be called from sequence reads with program Metagene. Since these sequence reads are short, a predicted ORF maybe a portion of a complete ORF. An ORF may also be a translation from a non-coding frame: such an ORF is called a spurious ORF, as defined in the original GOS study . The GOS study also introduced a spurious ORF detection method using nonsynonymous to synonymous substitution test, which is available along with a recent GOS clustering study . This method is not integrated within RAMMCAP, but it can be used independently to identify the spurious ORFs predicted here.
ORFs are first clustered at 90-95% identity to identify the non-redundant sequences, which are further clustered to families (ORF clusters) at a conservative threshold, so that each cluster contains sequences of the same or similar function. A 30% sequence identity indicates significant similarities for full-length proteins. Since ORFs from Sanger reads are long enough, so they are clustered at 30% identity over 80% of ORF length. ORFs from 454 reads are much shorter; they are clustered at 60% identity over 80% of ORF length. ORF clusters are used for functional studies. The size of an ORF cluster is the number of its non-redundant sequences. For one million ORFs, it takes a few CPU hours to cluster at 60% identity and ~100 CPU hours at 30% identity.
The clustering method in RAMMCAP is quite different from the clustering method in the GOS study  and its incremental update , which generated core clusters by all-against-all BLAST search and then merged core clusters into final clusters using sequence profile methods. The final clusters in the GOS study are large and contain sequences of very remote similarities, whereas the clustering method employed here only tries to generate very conservative clusters.
ORFs are annotated from Pfam and Tigrfam with HMMER (accelerated with Hammerhead), and from COG with RPS-BLAST. Hits must be with e-values ≤ 0.001, and meet the significant scores in case of HMMER searches. This annotation process only takes ~100 CPU hours for one million ORFs.
Optionally, ORF annotation may be performed quickly by running only the representative sequence of each ORF cluster and then transferring the annotation to other members in that cluster. But transferred annotation may be wrong in some cases, for example, where the target ORF has fewer domains than the source ORF (the representative). Since the annotation process is very fast, it is preferable to run all the ORFs for more accurate annotation.
In many metagenomic projects, multiple samples from different environmental conditions were studied. This manuscript describes a novel way to compare metagenomes and visualize their differences. Since sequences from multiple metagenomes are clustered into families or classified into reference families (Pfam, Tigrfam, or COG), metagenomes can be compared by their occurrence profiles across all the families or selected families of interests.
Here, an occurrence profile coefficient, rAB = NA ∩ B/NA ∪ B, is defined as the similarity measure between two metagenomes A and B. NA ∩ B is the number of families that are found in both A and B above a minimal occurrence cutoff (defined later) without significant difference. NA ∪ B is NA ∩ B plus the number of families that occur in one metagenome significantly higher than in another metagenome. The value of rAB is between 0-1, with 0 representing no overlap and 1 indicating a perfect match between A and B.
Let N A and N B be the number of sequences in A and B, let H A and H B be the number of sequences that occur in family H. One question is whether the difference between H A and H B is statistically significant. Rodriguez-Brito et al introduced a method to test such statistical significance of differences between two metagenomes . Rodriguez-Brito's method used a large amount (in order of 105) of simulations by randomly picking a certain number of sequences from A, B, and A+B to generate distributions and analyze it, so it is very time-consuming.
In this manuscript, H A is considered significantly higher than H B if both (1) the z score satisfied a user defined confidence level such as 0.95, and (2) P A ≥ f × P B , where f (f>1) is also a user defined parameter, called significant factor.
Since 2N A >>z 2, the minimal cutoff of H A is 4 at 0.95 confidence level (z = 1.96), and 7 at 0.99 confidence level (z = 2.58).
The occurrence profile coefficients are calculated for all the metagenome pairs, and the output matrix is used to hierarchically cluster the metagenomes.
Sensitivity and specificity of Metagene and ORF_finder on simulated metagenomic datasets
ORF_finder, cut off = 30aa
ORF_finder, cut off = 40aa
ORF_finder, cut off = 50aa
ORF_finder, cut off = 60aa
GOS ORFs have more than 1000 clusters that contain ≥ 1000 non-redundant sequences; BIOME ORFs have less than 100 such clusters. About 70% GOS ORFs, 46% BIOME reads, and a similar percent BIOME ORFs are found in non-singleton clusters. Within the BIOME datasets (Figure 3c, d), the microbiomes have more large clusters compared to viromes; suggesting that microbial sequences are more conserved than viral sequences.
Clustering analysis is a powerful tool to recover protein families and to discover the novel ones; it helps to recognize spurious ORFs. Clustering tends to put real ORFs into large clusters and leaves spurious ones in small or singleton clusters because spurious ORFs have more random features. If the ORFs with Pfam, Tigrfam, or COG matches are considered true ORFs, then 93% of these true GOS ORFs are found in clusters of size ≥ 10, which is only 1.3% total GOS clusters; here cluster size is the number of non-redundant sequences in a cluster. Further, 28% of the true BIOME ORFs are in 1.0% of top clusters of size ≥ 5. Many large clusters without any homology to known proteins are found, which may shed light on novel families of environment specific functions.
ORFs called from metagenomic reads are short and fragmented. In addition, errors such as frame shifting and wrong gene boundary may occur due to sequencing errors. Therefore a conservative threshold is used in producing ORF clusters to ensure that a cluster contains sequences of the same or similar function. The quality of clustering was evaluated with Pfam, the manually curated classification of protein families. The domain sequences in Pfam models (release 22.0) were extracted from the alignments and were clustered at 30% identity using the clustering protocol in RAMMCAP.
The CPU time for clustering GOS ORFs, BIOME reads, and BIOME ORFs were 9600, 125, and 513 hours, respectively. GOS ORFs cost relatively more, but still two orders of magnitude less than the original GOS study. The annotation for GOS ORFs and BIOME ORFs took 3800 and 1560 hours. Through clustering analysis, many novel families can be identified and can be used in metagenome comparison. The RAMMCAP software and pre-calculated results are available at http://tools.camera.calit2.net/camera/rammcap/.
This work was supported by Gordon and Betty Moore Foundation through the CAMERA project http://camera.calit2.net and NIH grant 1R01RR025030 from National Center for Research Resources. The CAMERA team provided much technical supports, especially Dr. Ying Huang who gave advices on statistical method. Forest Rohwer from San Diego State University provided the early access to the Nine Biomes data.
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.