MUSiCC: a marker genes based framework for metagenomic normalization and accurate profiling of gene abundances in the microbiome
© Manor and Borenstein; licensee BioMed Central. 2015
Received: 13 November 2014
Accepted: 10 February 2015
Published: 25 March 2015
Functional metagenomic analyses commonly involve a normalization step, where measured levels of genes or pathways are converted into relative abundances. Here, we demonstrate that this normalization scheme introduces marked biases both across and within human microbiome samples, and identify sample- and gene-specific properties that contribute to these biases. We introduce an alternative normalization paradigm, MUSiCC, which combines universal single-copy genes with machine learning methods to correct these biases and to obtain an accurate and biologically meaningful measure of gene abundances. Finally, we demonstrate that MUSiCC significantly improves downstream discovery of functional shifts in the microbiome.
MUSiCC is available at http://elbo.gs.washington.edu/software.html.
The study of naturally occurring microbial communities through shotgun metagenomic assays has become a routine procedure in recent years [1-6]. Such assays are used, for example, to catalog the collection of genes in the metagenome, to estimate their abundances, and ultimately, to characterize the functional capacity of the community [1,3,7,8]. This process involves two steps. First, genomic DNA is extracted from the sample and sequenced using next-generation technologies. Next, sequenced reads are aligned to a database of reference genes or genomes, and the number of reads that map to each gene is used as a proxy for its abundance in the sample [7,9,10]. Clearly, however, the resulting read counts are highly dependent on the sequencing depth in each sample, and some normalization method is required to allow comparison across samples. This is most commonly achieved by a simple compositional normalization process, whereby the obtained abundance value associated with each gene is divided by the sum of abundance values for all genes identified in the sample (for example, [2,11]). The resulting normalized value therefore represents a measure of relative abundance and is used in subsequent comparative analyses of the samples.
This normalization scheme, however, while extremely prevalent, has several fundamental weaknesses that may influence downstream analysis and ultimately impact the identification of functional shifts across samples. First, the resulting relative abundance values are unitless and do not necessarily represent a meaningful biological quantity. Second, in this normalization scheme, the scaled abundance of each gene crucially depends on the measured abundances of all other genes. As many different sample-specific factors could affect these quantities, abundance values could be disproportionately scaled in different samples, dramatically biasing any downstream comparative analysis. Compositional normalization is also associated with several statistical drawbacks and may give rise to misleading patterns [4,12]. For example, as a marked increase in the abundance of one element decreases the apparent relative abundance of other invariant elements, this normalization scheme tends to induce spurious correlations between various elements in the sample. As a result, comparative analyses of these values across samples may be hard to interpret. These drawbacks call for an alternative normalization procedure, one that can produce accurate and easy to interpret abundance measures that can be reliably compared across samples.
Notably, a few previous metagenomics-based studies have already highlighted the challenges involved in compositional normalization. Specifically, studies of species composition have previously demonstrated that compositional normalization of taxonomic data could both mask true correlations between pairs of taxa and introduce false correlations [13-16]. Other studies of oceanic communities have further emphasized the biases introduced by compositional normalization of environmental metagenomic samples, specifically highlighting the potential contribution of the average genome size in each sample to these biases [17-19]. To date, however, the impact of compositional normalization on functional metagenomic studies of the human microbiome has never been shown or characterized, nor have the various sample-specific properties that may contribute to inaccuracies in abundance measures. Furthermore, previous studies of environmental metagenomes that aimed specifically to address genome-size induced bias still failed to provide biologically meaningful and interpretable measures of gene abundance.
Finally, even within each sample, various gene-specific properties may bias measured abundances. Compositional normalization, or for that matter, any normalization scheme that applies an identical processing protocol to all genes, will inevitably fail to account for such errors. Indeed, to date, no attempts to characterize or correct within-sample biases in genes’ abundances have been introduced, potentially neglecting important factors that may further contribute to inaccuracies in gene abundance measures.
In this study, we analyze samples from the Human Microbiome Project (HMP) [2,7,11,20], as well as samples from two additional independent studies of the human gut microbiome [4,11,21], and demonstrate that compositional normalization has a clear and measurable effect on the obtained metagenomic functional profiles. We specifically show that this normalization protocol induces spurious inter-sample variation in the calculated abundances of various genes across samples from the same body site, across samples from different body sites, and across samples from different studies. We identify three sample-specific properties that play a key role in generating this spurious variation, including the average genome size, species richness, and mappability of genomes in the sample, and suggest a simple and more biologically meaningful normalization method that aims to quantify the typical genomic copy number of each gene in the sample. We additionally show that gene-specific properties, such as sequence conservation and nucleotide content, further induce spurious intra-sample variation in the measured gene abundances, and provide an additional machine learning-based method to correct these inaccuracies. Finally, we demonstrate that our correction paradigm indeed provides improved abundance estimates and has clear and significant benefits for downstream comparative analyses. Specifically, we show that our method aids in the discovery of differentially abundant genes (and corrects false discovery of invariant genes), markedly increases the power to detect disease-associated pathways, and supports cross-study comparative analyses. Combined, these benefits allow us to move towards a more rigorous and unbiased estimate of the average genomic copy number of each gene in human microbiome samples, facilitating an accurate identification of functional shifts that may be associated with disease.
Results and discussion
Spurious inter-sample variation in HMP samples and its determinants
Correlations and controlled correlations between the abundance of USiCGs and various sample-specific properties in the stool samples
Correlation with USiCGs abundance
Controlled correlation with USiCGs abundance a
-0.82 (P <10-14)
-0.46 (P <10-3)
-0.85 (P <10-17)
-0.58 (P <10-6)
0.44 (P <10-3)
0.37 (P <0.003)
Correlations between sample-specific properties
0.8 (P <10-300)
-0.25 (P <0.04)
-0.27 (P <0.02)
Spurious inter-sample variation across different studies
We further wished to examine whether the same sample-specific properties that we identified above as contributing to this bias in HMP samples play a role also in the observed spurious variation across the different studies. However, as full taxonomic data (in the form of OTU counts) was only available for HMP samples, a direct comparison was not feasible. Instead, we therefore searched for a sample-specific property that was available for all three studies and that could serve as a statistical proxy for the OTU-based properties examined above. Using the HMP dataset, we found that the total number of assembled contigs in a sample was highly correlated with the average genomic size, species richness, and average mappability in the sample (Pearson correlation R = -0.52, R = 0.55, and R = -0.53, respectively). Since all three studies performed de novo contig assembly, this measure could be used in lieu of the more explicit sample-specific properties utilized above. Comparing the number of contigs per sample in the various studies with the median relative abundance of USiCGs, both between and within the different studies, we indeed found a similar trend as observed above. Specifically, the study with the highest number of assembled contigs (HMP; median value of 119,231 contigs per sample) was also the one with the highest median relative abundance of USiCGs (8.65 × 10-4), whereas the study with the fewest contigs (T2D; median value of 27,274 contigs per sample) was the one with the lowest median relative abundance of USiCGs (2.63 × 10-4). Similarly, a modest to strong correlation was found between the number of assembled contigs and USiCGs’ abundance across the samples within each study (R = 0.6, R = 0.11, and R = 0.28 for HMP , IBD , and T2D , respectively).
These findings suggest that spurious variation in the relative abundance of USiCGs is also evident across multiple independent studies and may be accounted for by similar sample-specific properties. Importantly, this implies that compositional normalization may hinder future attempts to pool datasets from multiple studies in order to increase the power to discover functional variation in the human microbiome (and see also our cross-study comparative analysis below).
Inter-MUSiCC: correcting spurious inter-sample variation
To correct the spurious inter-sample variation demonstrated above, we present here a fundamentally different normalization scheme. Clearly, with the identification of sample-specific properties that contribute to spurious variation, one potential approach for addressing this challenge would be to carefully characterize how each sample-specific property affects normalized abundances and apply some mathematical formulation to reverse this effect. Previous attempts to specifically correct the impact of average genome length indeed took this approach [18,19]. However, considering the multiple sample-specific properties identified in the previous section and their complex intertwined contribution to spurious variation, we took a more direct approach that corrects spurious sample-specific variation, regardless of its determinants, and provides a biologically intuitive and meaningful measure for gene abundances in a metagenome (and see also our analysis of simulated samples below). Specifically, our method aims to describe the abundance of each gene in the microbiome, not as its relative abundance in the sample, but rather as the average copy number of this gene (with respect to some copy number proxy as described below) across all microbial cells in the sampled community. The functional profile obtained by this normalization method for a given sample can accordingly be conceived as a description of the gene content of a ‘typical microbe’ from that sample. This profile therefore has a clear biological interpretation, and can be reliably and meaningfully compared across samples.
In practice, this normalization paradigm is achieved by using the same set of USiCGs described above as the yardstick for estimating the average copy number of all other genes. Technically, the measured abundance of each gene in a sample (after correcting for gene length) is divided not by the sum of all abundances, but rather by the median abundance of the USiCGs in the sample. With this correction, USiCGs in every sample have a consistent median value of one (as expected by their definition), and all other genes are measured in comparison to these USiCGs. We therefore term this normalization method ‘Inter-sample Metagenomic Universal Single-Copy Correction’ (Inter-MUSiCC).
Evaluating the impact of Inter-MUSiCC on functional metagenomic profiles
Inter-MUSiCC, by definition, corrects the spurious inter-sample variation we demonstrated above in USiCGs. For this normalization scheme to be useful, however, it should clearly also yield improved abundance estimates for all other genes. We focus on HMP stool samples, owing to the large number of samples available and the ability to examine also the impact of MUSiCC on multiple studies of the gut microbiome. Focusing on a single body site also allows us to assess the performance of MUSiCC on a more homogenous set of samples and therefore to evaluate MUSiCC’s ability to correct even the relatively fine variation observed in such samples. Notably, however, as the actual average copy number of each gene in these samples is unknown, evaluating the impact of inter-MUSiCC and assessing whether it produces reasonable estimates of average copy numbers in these metagenomes is a challenging task.
To address this challenge, we considered the set of genes that occur in only one OTU per sample (Methods). For such OTU-Specific Genes (OSGs), the relationship between the abundance of each gene and the abundance of the corresponding OTU is not complicated by the presence of other OTUs in the sample, making it easier to evaluate the impact of Inter-MUSiCC on the corrected abundance values. Specifically, as each OSG occurs in only one OTU, clearly, the abundance of the OSG across the various samples should positively correlate with the abundance of the respective OTU. If Inter-MUSiCC indeed provides a more accurate estimation of gene abundance compared to compositional normalization, this correlation between OSGs’ and OTUs’ abundances should improve with the application of Inter-MUSiCC. Since many OTUs observed in each sample are not yet associated with a fully sequenced genome, we used a recently introduced tool to predict the genomic content of each OTU (Methods) . In total, we identified 3,821 OSGs in 993 OTUs across 65 HMP stool samples. For each such OSG, we calculated the correlation between its abundance across the various samples and the abundance of its respective OTU with and without the application of Inter-MUSiCC. We found that Inter-MUSiCC indeed significantly increased the average correlation between the abundances of OSGs and their respective OTUs (P <10-32, paired t-test; Methods), suggesting that the abundance values obtained by Inter-MUSiCC more accurately reflect gene abundances. To further confirm the beneficial impact of Inter-MUSiCC on the abundance of OSGs, we additionally examined the ratio of OSGs’ (and OTUs’) abundances across pairs of samples. Such a ratio-based analysis may cancel out various factors that confound the relationship between OSGs’ and OTUs’ abundances (such as false detection of OSGs), further highlighting the accuracy of calculated OSGs’ abundances. Specifically, as each OSG occurs in only one OTU, the ratio between its abundance in one sample to its abundance in another should be similar to the ratio between the abundances of the respective OTU in these two samples. For example, if a certain OTU exhibits a three-fold increase in abundance between two samples, any OSG that is associated with only this OTU in these two samples is also expected to exhibit a three-fold increase in abundance. To confirm this, for each case where the same OSG appeared in two different samples, we quantified the difference between the fold-change in the OSG abundance between the two samples and the fold-change in the OTU abundance between these samples (Methods). We found that both when using compositional normalization and when using Inter-MUSiCC, OSGs exhibit on average a higher fold-change than expected from the fold-change in their respective OTU abundances. This may reflect subsampling issues, wherein low-abundance OTUs can be accurately measured through 16S sequencing, while the abundance of their genes may still be underestimated by shotgun metagenomic sequencing. Yet, with Inter-MUSiCC, OSGs’ abundance fold-change was significantly closer to the measured fold-change in OTU abundance (P <10-315, paired t-test), suggesting that at least some of the bias in the abundances of each gene was corrected.
Finally, since compositional normalization is known specifically to introduce false correlations between normalized elements , we set out to examine the impact of Inter-MUSiCC on pairwise gene correlation. As before, a gold standard for the true correlation structure in the human microbiome is not available. However, since metagenomes represent a large collection of genomes, one could expect that pairs of genes whose occurrences across genomes are highly correlated (that is, they tend to either both be present in a genome or both be absent), would also be correlated in their abundance across metagenomes. We therefore compared for each gene the correlation of its occurrence and the occurrence of every other gene across fully sequenced genomes in the KEGG database [1,3], with the correlation of its abundance and the abundance of every other gene in HMP stool metagenomic samples (Methods). We found that with Inter-MUSiCC the similarity between the genomic and metagenomic correlation structure is significantly higher than with compositional normalization (P <10-121, Student’s t-test), again suggesting that Inter-MUSiCC better captures the abundance levels of the various genes in the metagenome.
Spurious intra-sample variation and its determinants
Intra-MUSiCC: correcting spurious intra-sample variation and evaluating its impact on functional metagenomic profiles
Our findings above suggest that spurious intra-sample variation in USiCGs could be corrected by learning a predictive linear model that accounts for the impact of various gene-specific properties on measured gene abundances. Assuming that such spurious intra-sample variation also occurs in other, non-USiCG genes, we can use the USiCGs-based model as a proxy for the effect of gene-specific properties and apply this model to correct the measured abundances of all other genes in the sample. Since this method again relies on the set of USiCGs identified above, we term it ‘Intra-sample Metagenomic Universal Single-Copy Correction’ (Intra-MUSiCC). Notably, this predictive model can be learned once (for example, using a large set of metagenomic samples) and applied to every new sample without modification. In what follows, however, we used a somewhat more sophisticated approach wherein a predictive model is learned in each sample separately (based on the abundances of USiCGs in that sample), cross-validated on unseen USiCGs abundances, and used to correct variation within that sample. This approach assumes that the impact of gene-specific properties on measured abundances may slightly vary from sample to sample and aims to capture the potentially unique effect of gene-specific properties in the sample. We confirmed, however, that the more generic approach, in which the predictive model is learned once and applied to all samples, produced qualitatively the same results as reported below.
In the previous section, we already demonstrated that Intra-MUSiCC corrects much of the spurious intra-sample variation in USiCGs’ abundances (see Figure 5). However, to confirm that Intra-MUSiCC is a useful and applicable correction scheme, we next set out to validate that it indeed yields improved abundance estimates also for non-USiCGs genes, focusing, as before, on the set of HMP stool samples (see above). As was the case for our analysis of Inter-MUSiCC, this task is challenging since the true abundance values of non-USiCGs across samples is not available. We therefore used several different approaches, each focusing on a specific set of genes, to evaluate the performance of Intra-MUSiCC and to confirm that the corrected gene abundance values are indeed more accurate and more informative than the widely-used relative abundance values.
First, we identified a set of 72 genes that did not meet our criteria for USiCGs, but that could still be considered universal single-copy genes under a more relaxed set of requirements (Additional file 6: Table S3; Methods). Specifically, such ‘semi’-USiCGs tend to have a single copy in the vast majority of genomes, but are not as prevalent as the USiCGs described above. While this set is expected to be more variable than USiCGs (due to true variation in their occurrence across genomes), we hypothesized that at least some of the observed intra-sample variation among these genes might be spurious and accounted for by the various gene-specific properties detected above. Indeed, we found that when using the USiCGs-based model (Intra-MUSiCC) to correct the abundance values of these genes across all HMP stool samples, 17% of the variation in these semi-USiCGs was corrected. Notably, this reduction in variation is observed even though these semi-USiCGs were never used in the construction of the model.
Next, we focused on pairs of genes whose occurrence is highly correlated across genomes. Specifically, examining the gene content of all reference genomes in KEGG [1,3], we identified 1,074 pairs of genes whose presence/absence profiles agree across >95% of the genomes in which they appear (Methods). Importantly, we excluded all USiCGs (as well as the semi-USiCGs described above) from this analysis. For each such pair, we computed the average absolute difference in abundance across HMP stool samples, before and after applying Intra-MUSiCC. Clearly, high co-occurrence across genomes does not necessarily imply perfectly correlated abundance across metagenomes. A metagenome may, for example, include a relatively large proportion of exactly those genomes in which the two genes do not co-occur. Yet, it is reasonable to assume that such genomically co-occurring genes will tend to exhibit similar abundances across metagenomes. Accordingly, we hypothesized that the observed differences between the abundances of such gene pairs can be partly accounted for by spurious intra-sample variation that arises from gene-specific bias as opposed to true biological variation. In line with this hypothesis, we indeed found that both the mean and variance of the pairwise differences in the abundance of each of the 1,074 gene-pairs were significantly lower with the application of Intra-MUSiCC (P <10-66 and P <10-46, for mean and variance, respectively; paired t-test for mean, paired F-test for variance).
Finally, we extended the analysis of genomically co-occurring genes, focusing on clusters (rather than pairs) of genes whose occurrence is highly correlated across genomes. We again used the set of KEGG reference genomes to identify 25 gene clusters, each consisting of at least five genes (with a total of 216 genes) whose pairwise presence/absence profiles agree in >95% of genomes (Additional file 7: Table S4; Methods). Again, we assumed that such genomically co-occurring genes are likely to exhibit similar abundances across HMP gut metagenomic samples. We used Intra-MUSiCC to correct the abundance of each gene in the various clusters and compared the corrected abundance measures to the original relative abundance values in each cluster. Importantly, using clusters of genes rather than independent pairs allowed us to examine not only the pairwise difference in abundance but also the proportion of variance explained within each cluster. We found that on average Intra-MUSiCC corrected 8.2% of the variance within each gene cluster (Additional file 7: Table S4). As an example, one such gene cluster contained all eight genes in the F-type ATPase structural module. Although these eight genes are highly consistent in their presence/absence patterns across reference genomes (median pairwise Jaccard similarity of >98%), their abundances across HMP stool samples vary by up to approximately 2.1-fold on average (and >13.4-fold in some samples). Intra-MUSiCC (which notably employs a model learned on USiCGs alone and for which these eight genes can be conceived as unseen data) corrected on average 32% of this variance.
MUSiCC markedly enhances the discovery of functional shifts in the microbiome
Above, we demonstrated that MUSiCC successfully corrects both inter- and intra-sample spurious variation in gene abundances. Clearly, however, the key goal of any metagenomic normalization scheme is to facilitate comparative analysis and to enable the discovery of functional shifts in the metagenome that may be associated with a given host phenotype. Next, we therefore set out to examine whether MUSiCC (that is, the combined application of Inter-MUSiCC and Intra-MUSiCC) has a concrete and significant impact on such comparative analyses and whether it affects the identified set of differentially abundant genes or pathways in various settings.
Crucially, MUSiCC’s impact on the detected set of differentially abundant genes is observed not only when comparing different body sites in healthy individuals, but also when comparing disease cases vs. healthy controls. Specifically, analyzing the two disease-related studies discussed above (namely, the Type 2 diabetes study ; T2D, and the inflammatory bowel disease study ; IBD), a similar pattern was found, with many genes detected as differentially abundant by only one of the two normalization methods. In fact, in these datasets, the set of differentially abundant genes discovered with and without MUSiCC exhibit an even lower agreement, with only 78% and 48% overlap between the methods for T2D and IBD, respectively (Figure 6B and C). As before, the sets of disease associated differentially abundant genes discovered only by MUSiCC were enriched for pathways previously reported to be linked to T2D and IBD, including glycolysis/gluconeogenesis [13,15] and glutathione metabolism , respectively. These findings suggest that MUSiCC not only corrects various biases in metagenomic functional data, but can also lead to significant changes in the identification of disease-associated genes, pointing to novel candidates for microbiome-based intervention (and avoiding false discovery of other, clearly unrelated genes).
Considering the overall increase in the significance level of disease-associated pathways observed above with MUSiCC (Figure 7A and B) and since many disease-related studies can obtain or process only a limited number of samples, we further examined whether MUSiCC can indeed enhance the power to detect disease-associated pathways when data are limited. We therefore repeated the comparative analysis above, focusing on two disease-associated pathways of interest, and measured how significant the obtained association signal was when using only a subset of the available samples. Indeed, we found that with MUSiCC, significance is reached with far fewer samples than with compositional normalization (Figure 7C and D), suggesting that MUSiCC can successfully uncover novel disease-associated pathways that may be masked by noise or by sparse sampling when using standard compositional normalization. Using this disease-associated pathway discovery as the ultimate benchmark for the applicability of the various normalization methods, we additionally tested an alternative approach that was previously used to process a set of oceanic samples (; Methods). We found that this method outperformed the standard compositional normalization approach but was nonetheless far less successful than MUSiCC and still failed to identify many of the relevant pathways as significant (Additional file 9: Figure S5).
MUSiCC significantly reduces spurious variation in simulated bacterial samples
Above, we demonstrated that MUSiCC reduced spurious variation and improved our ability to detect disease-associated pathways in real metagenomic datasets. Such datasets provide a means to assess the full range of factors that potentially impact functional profile measurements. However, since the exact underlying taxonomic and functional compositions in these real datasets are unknown they cannot serve as a gold standard for evaluating our method or for comparison. We therefore wanted to further examine the ability of MUSiCC to reduce spurious variations on a synthetic dataset, where the true abundances of genes and pathways are available. To this end, we generated 20 simulated metagenomic samples, each of which consisted of 500,000 reads generated at random from a mock community of reference genomes that were randomly assigned different relative abundances (see Methods and Additional file 10: Table S5). We then mapped the reads in each sample to the KEGG database and calculated the read count of each KEGG Orthology group (KO) in each sample. We finally normalized and corrected the obtained KO abundances using either MUSiCC or standard compositional normalization.
Next, since many comparative analyses of metagenomes are performed at the pathway level (for example, by summing the abundances of all the KOs associated with a given pathway), we wanted to specifically examine the impact of MUSiCC on spurious pathway abundance variation. To this end, in generating the simulated samples discussed above, we intentionally limited our selection of bacterial genomes to those that contain the entire set of flagellar assembly genes with the same exact copy numbers. Accordingly, any observed variation in the flagellar assembly pathway across the various samples is, by construction, spurious. Indeed (see Figure 9C), the coefficient of variation (CoV; standard deviation over mean) of the corrected abundance of this pathway across samples was markedly lower with MUSiCC (CoV = 0.03) than with compositional normalization (CoV = 0.175). To allow a more rigorous statistical analysis of the performances of each normalization method and to quantify the robustness of the various methods to subsampling, we further used a bootstrapping approach, repeatedly selecting a subset of the simulated samples and calculating the CoV observed in the abundance of this pathway across each subset. Comparing the distribution of CoV obtained for 100 subsets, we confirmed that inter-sample variation is indeed significantly lower with MUSiCC than with compositional normalization (Figure 9C; P <10-72, Student’s t-test for the difference in mean CoV). We additionally found that the distribution of CoV values obtained with MUSiCC across these 100 subsets exhibited lower variance than the distribution of CoV values obtained with compositional normalization (Figure 9C; P <10-29, F-test), suggesting that MUSiCC is also more robust to subsampling.
Finally, since our analysis revealed that variability in the average genome size across real samples is one of the factors contributing to spurious inter-sample variation, we wished to use the simulated samples described above to compare MUSiCC to normalization that is based on the average genome size in each sample. In addition to using the real average genome sizes (which are known for the simulated samples), we utilized two previously introduced methods, namely GAAS  and Raes et al. , for metagenomic-based estimation of the average genome size in each sample (see Methods). While these methods were not developed with the aim of inter-sample normalization in mind, the average genome size estimates they provide may still be applied to correct observed spurious variation. We found that while both real and estimated genome size normalization markedly decreased spurious variation across samples in KO abundances (Additional file 11: Figure S6) and in pathway abundances (Figure 9C), MUSiCC significantly outperformed these genome size based approaches. For example, the mean CoV for the abundance of the flagellar assembly pathway across samples was significantly lower in MUSiCC than in any genome size normalization methods (P <10−50, P <10-50, and P <10-60 , for real average genome size, GAAS, and Raes et al. respectively; Student’s t-test), as was the variance in the CoV across the 100 subsets described above (P <0.02, P <0.03, and P <10-4, for real average genome size, GAAS, and Raes et al. respectively; F-test). These results strengthen our finding that other sample-specific properties impact inter-sample variation and suggest that our USiCGs-based approach indeed captures a wide range of determinants in reducing spurious variation (and see also the discussion below), providing a more reliable characterization of the functional capacity of different metagenomes.
The development of computational and statistical protocols for analyzing high-throughput metagenomic data has been a major line of research in recent years. Such protocols are essential for accurately characterizing the composition of various microbial communities and ultimately for reliably detecting functional shifts that may be associated with disease. These protocols often include a compositional normalization step, wherein read counts are converted into relative abundances, allowing researchers to compare different samples with potentially markedly different coverage. Importantly, this normalization procedure is extremely common and has become the de facto standard in metagenomic assays and especially in human microbiome studies. Above, however, we systematically demonstrated that this normalization scheme introduces marked inter-sample spurious variation both across samples from the same study and across different studies of the human microbiome. We further rigorously characterized various sample-specific properties that give rise to this variation, demonstrating that the average genome size, the species richness, and the average genome mappability in the sample all independently contribute to spurious inter-sample variation. Rather than trying to reverse these complex and intertwined effects, we suggested a simple correction method, using a large set of universal single-copy genes to normalize the abundance of all other genes. This method accordingly not only corrects spurious variation across samples and facilitates accurate downstream comparative analyses, but also provides a more biologically relevant and meaningful abundance measure, estimating the average genomic copy number of every gene in the sample.
We additionally demonstrated that various gene-level properties, such as sequence conservation and nucleotide content, induce measurable and clear spurious intra-sample variation, further biasing calculated gene abundance profiles. We showed that the impact of these gene-specific properties on the measured abundance of universal single-copy genes can be modeled through an L1-regularized linear regression approach, and that the learned model can be used to partially correct such spurious intra-sample variation across all other genes in the metagenome. Notably, correcting intra-sample variability can be done by either using the built-in intra-MUSiCC trained model or by learning a specific model de novo for each new sample. Since our built-in model was learned using the HMP samples, this model is especially appropriate for samples that were obtained through a similar protocol. For samples obtained through different technologies or different processing pipelines that could potentially give rise to different biases, the intra-MUSiCC model should be ideally retrained on each dataset. Importantly, however, model training is extremely fast (see Methods), and since it is based only on the observed abundances of USiCGs (rather than on raw reads), running time is independent of the number of reads or their length.
Correcting inter- and intra-sample biases is clearly crucial for obtaining an accurate measure of gene abundance and for faithfully characterizing the functional profile of the metagenome. Since the true functional profile of the human microbiome is generally unknown, we used multiple approaches to confirm that our correction methods indeed yield improved abundance values and better recover the expected abundance and correlation structure of various sets of genes. Most importantly, however, we demonstrated that our methods not only provide a more accurate gene abundance measure, but also has a significant impact on the discovery of differentially abundant genes and of enriched pathways. Specifically, we showed that our MUSiCC pipeline markedly enhances the identification of disease-associated pathways, offering substantially increased statistical power both in terms of the number of pathways identified and the number of samples required for identification. Perhaps most remarkably, we found that MUSiCC facilitates efforts to pool data from several human microbiome studies, correcting much of the study-specific variation and laying the foundation for future cross-study comparative analyses. Notably, without MUSiCC, such study-specific variation was so dramatic that it masked practically any genuine, disease specific variation in the data.
In addition, we also evaluated MUSiCC using a set of simulated samples in which the underlying copy number of each gene in each sample is known. We demonstrated that MUSiCC significantly reduced spurious inter-sample variation both at the KO level and at the pathway level compared to compositional normalization. We further demonstrated that using the average genome size in each sample (either real or predicted) for normalization is not sufficient for removing inter-sample variation and that MUSiCC significantly outperformed such a normalization approach even in the ideal case of simulated communities. This finding highlights the benefit of a marker gene based normalization scheme and specifically the advantage of using USiGCs as a yardstick for copy number estimation, since it can capture the full range of factors that may introduce inter-sample variation. Moreover, it should be noted that methods for estimating average genome size often employ scenario-specific optimized parameters or require complete alignment, whereas MUSiCC is a parameter-free method that relies solely on KO measured abundances.
Clearly, the development of robust computational and statistical methods for an accurate characterization of gene abundances is an ongoing effort. Specifically, while inter-MUSiCC potentially corrects most of the inter-sample variation introduced by sample-specific properties, correcting gene-specific intra-sample biases is a much more challenging task. For one, the set of gene-specific properties that could affect the measured abundance of a gene may be markedly larger than the set we analyzed. Moreover, the directionality and magnitude of such effects may not be consistent across different groups of genes (note, for example, that a USiCGs-based model corrected >50% of the variation in unseen USiCGs, but <20% of the variation in semi-USiCGs; and see ). Distinguishing gene-specific effects from sampling noise is also challenging. We therefore hope that future work could further elucidate the contribution of additional properties and fine-tune the model relating gene-specific features to functional metagenomic profiling biases. Another opportunity for future extension of our method is going beyond bacterial and archaeal organisms and accounting also for other domains of life. Recent studies, for example, have highlighted the importance of the fungal component in the human microbiome . Since our method was designed primarily with bacterial and archaeal genomes in mind, it may not be optimized to fungi-rich samples and may not accurately correct the abundances of fungal genes. Importantly, however, although the set of USiCGs used in our study was selected based on prokaryotic genomes, it is in fact well represented also in fungi, with 55 of the 76 USiCGs found on average in 58 of the 71 (82%) fungal genomes currently in KEGG, with a median copy number of 1.03. Rigorous processing of fungi-rich samples would thus ideally involve a modified and carefully selected set of marker genes and specifically tailored settings. Similarly, several recent studies have demonstrated that specific body sites, such as the skin, can harbor a substantial viral component [26-28]. Clearly, the human virome is an active and important component of the human microbiome, and as future efforts provide better characterization of this component, our methods could be further refined to account for the viral genes that may be identified across samples. One approach, for example, would be to augment our analysis with a sequence-based preprocessing step to distinguish reads that originated from this viral component. Combined, such extensions will allow researchers to improve various metagenomic processing pipelines and move towards a more accurate estimation of the functional composition of the metagenome. We believe that the MUSiCC framework and the improved estimation of the average genomic copy number of each gene in the metagenome represent an important step in this direction and can ultimately aid in the discovery of microbiome-based therapeutic targets.
Software implementation and distribution
MUSiCC was implemented in Python and is available for download or as a web-based application at: http://elbo.gs.washington.edu/software.html. In addition, it is available both in GitHub (https://github.com/omanor/MUSiCC) and as a pip-installable python package (using: pip install -U MUSiCC). Input files for MUSiCC should be processed for host contamination removal and annotated by the KEGG orthology group database.
Functional metagenomic data from multiple body sites were obtained from the HMP  and were downloaded from the Data Analysis and Coordination Center (DACC) website. Genes were labeled by their KEGG Orthology (KO) Group and relative abundances were normalized for length. Sample and gene (KO) abundance data for the inflammatory bowel disease study  were obtained from a subsequent analysis of these samples performed by . Sample and gene (KO) abundance data for the type 2 diabetes study were obtained from .
PiCRUSt  pre-calculated matrices for operational taxonomic units (OTUs) and their predicted KOs were downloaded from the developer’s github site. OTU abundances mapped to GreenGenes IDs in the HMP body sites were obtained through personal communication with the lead author of PICRUSt.
Detecting Universal Single-Copy Genes (USiCGs)
The list of USiCGs was compiled to include KOs that are both universal and appear in a single-copy in each genome. To determine the level of universality required, we used the list of 31 marker genes from the PhlyoSift pipeline , and examined the number of KEGG genomes in which each of these genes appear. We found that at minimum (after removing one outlier) these genes appear in 91.5% (2,313) of the bacterial and archaeal genomes in KEGG, and therefore considered as universal any gene that appears in at least 91.5% of KEGG genomes. Of these genes, we further selected those that had an average number of copies per genome <1.1 in the genomes they appeared in, resulting in a list of 76 genes.
Estimating average genome size in HMP samples
To estimate the average genome size in each HMP sample, we used PICRUSt pre-calculated files (see above), listing the predicted set of KOs and their copy number in each OTU. We multiplied the KO copy number by the average KO length in the KEGG database [1,3] to obtain an estimate of the effective genome size of each OTU. We then computed the weighted average genome size in each sample, weighting the estimated genome size of each OTU by its relative abundance.
Estimating the species richness and genome mappability in HMP samples
To estimate the species richness of a given sample, we counted the number of OTUs identified in that sample. To estimate the average mappability of short reads in a given sample, we utilized a recently introduced measure (implemented as part of the PICRUSt software ), termed Nearest Sequenced Taxon Index (NSTI), which aims to evaluate the evolutionary distance between each OTU in a sample to its closest reference genome. Specifically, we used PICRUSt pre-calculated NSTI values (see above) of each OTU, and for each sample computed the weighted average of (1-NSTI) across all OTUs in the sample, weighted by their relative abundance.
Controlled correlations between USiCGs’ abundances and sample-specific properties
When calculating the correlation between the abundance of USiCGs and the various sample-specific properties, we used partial correlation analysis to control for the effect of other sample-specific properties. Partial correlations (and the corresponding P values) were calculated using the ‘partialcorr’ function in MATLAB. We additionally examined the correlation between species richness (and genome mappability) and USiCGs’ residuals after correcting for average genome size (see Additional file 4: Figure S3C and D and Additional file 3: Figure S2E and F). Specifically, we used regression analysis with the median abundance of USiCGs in each sample as the response and the average genome size as the regressor, using the ‘robustfit’ function in MATLAB to obtain the residuals. P values for the correlation with the residuals are reported in Additional file 4: Figure S3C and D and in Additional file 3: Figure S2E and F.
Identifying and analyzing OTU-Specific Genes (OSGs)
For each sample, we used the PICRUSt pre-calculated file (see above) that describes the predicted KO content of each OTU , combined with the list of OTUs that appear in the sample, to generate a list of KOs that are predicted to originate from only a single OTU in that sample (termed OTU-Specific Genes, or OSGs). For each OTU and each of its respective OSGs, we computed the Pearson correlation between the relative abundance of the OTU and the abundance of the OSG (using either compositional normalization or Inter-MUSiCC) across HMP stool samples. To be included in this analysis we required that both the OTU and the OSG were present in at least five samples. We then compared the correlation coefficients obtained when using compositional normalization to those obtained when using Inter-MUSiCC. To control for various factors that may confound the relationship between OSG and OTU abundances, we further performed a ratio-based analysis. For a pair of samples, sharing at least one OTU that contains OSGs, the fold-change between the OTU across the two samples was compared to the fold-change of each of its OSGs across the two samples. To avoid transitivity effects between pairs of samples, this process was performed as follows: for each OTU, the samples were ordered by the OTU’s relative abundance, and pairwise comparisons were performed only between consecutive samples in this ordered list. This was repeated for all OTUs. The distribution of the ratios between OSG and OTU fold-changes obtained when using compositional normalization was compared to the distribution obtained when Inter-MUSiCC was applied. The Statistical significance of the reduction in the mean and variance of this distribution were computed using the ‘ttest’ and ‘vartest2’ functions in MATLAB.
Examining pairwise gene correlation structure with compositional normalization and with Inter-MUSiCC
We first downloaded all fully sequenced and annotated genomes from KEGG (downloaded on 15 July 2013). For each gene (KO), we computed the Jaccard similarity between its presence/absence pattern across these genomes and the presence/absence pattern of all other genes. For each gene, we then also computed the Pearson correlation between its abundance across HMP stool samples and the abundances of all other genes. Given these two metrics, finally, for each gene, we computed the Pearson correlation between its correlation with all other genes across genomes and its correlation with all other genes in metagenomes. Distances and statistical significance were computed using the ‘pdist’ and ‘ttest’ functions in MATLAB.
Gene length and species statistics properties were downloaded from KEGG. Gene length was calculated as the average length of all genes labeled with the associated KO. Conservation and alignment properties were calculated by first downloading the gene sequences associated with each KO from KEGG. Next, MAFFT  was run on the set of genes for each KO, and statistics were calculated from the MAFFT multiple alignment output. Nucleotide content properties (for example, mean GC%) were calculated from the set of gene sequences assigned to each KO in KEGG. KO recall and precision properties were obtained from a large-scale study of short read annotation  and describe the average recall and precision (across multiple genomes) in annotating simulated short shotgun reads originating from each KO.
Regularized linear model linking gene-specific properties to intra-sample variation
For a given sample, we first defined the observed response as the fold-change between the abundance of each USiCG and the mean abundance of all USiCGs in the sample. The model covariates were defined for all samples as the standardized values of the various USiCGs’ gene-specific properties (Additional file 5: Table S2). When learning the model for a specific sample, we used a strict five-fold cross-validation (CV) scheme to learn an Elastic-Net regularized linear model [33-35] that predicts the fold-change of each USiCG in the sample. Importantly, in each CV partition to training and test sets, we first learned the penalty parameter and model weights using solely the training data (by using an internal CV scheme with the MATLAB version of glmnet ), and evaluated the performance of our learned model by quantifying the fraction of variation the model explains when predicting the fold-change of USiCGs held-out as test data. Running time for the learning step on a typical sample (approximately 13,000 KOs) took less than 1s on a single core processor.
Identifying Semi-Universal Single-Copy genes (semi-USiCGs)
Semi-USiCGs were selected as genes that were present in at least 2,148 (85%) of bacterial and archaeal genomes in KEGG with an average number of copies per genome <1.1, but did not reach the USiCGs threshold of 2,313 (91.5%) genomes, resulting in a list of 72 such genes.
Identifying highly-correlated genes across genomes
For all KO pairs, we computed the Jaccard similarity across all KEGG genomes. Only the KO presence/absence pattern was considered, ignoring KO copy number. KO pairs with Jaccard similarity >0.95 (that is, they agree on at least 95% of the genomes they appear in) were selected, resulting in 1,074 such pairs. Clusters were defined as sets of five or more genes, all with pairwise Jaccard similarly >0.95.
Identifying differentially abundant genes
For a given gene and two sets of metagenomic samples (for example, stool vs. tongue or IBD cases vs. controls), we compared the distribution of abundances between the two sets using the Wilcoxon rank-sum test. Genes that passed the Bonferroni correction (for comparing HMP body sites) or the FDR correction (for T2D and IBD cases vs. controls) with corrected P values <0.05 were defined as differentially abundant.
Identifying disease-associated pathways
For a given pathway and two sets of metagenomic samples (for example, disease cases and controls), we first computed the pathway abundance in each sample as the sum of the abundances of genes associated with that pathway. Next, we computed an association P value for the pathway by comparing the distribution of pathway abundance values in the two sets using the Wilcoxon rank-sum test. Pathways that had a higher median in disease samples and passed FDR correction (<0.05) were defined as disease-associated pathways.
Evaluating an alternative normalization approach
To compare the ability of MUSiCC to discover disease-associated pathways with an alternative normalization approach previously used to process a set of oceanic samples , we followed the scheme applied in  to correct the gene abundances in the IBD dataset. Specifically, we first estimated the average genome size within each sample, using the set of eight universal genes suggested in . For each gene we estimated the average genome size, G, as G = (g + L-2 m)/f, where f denotes the relative abundance of the gene in the sample, L denotes the sequence read length in the data (set to 75 bp for the IBD dataset), and m denotes a minimum overlap parameter (set to 90 as in ). Next, we averaged the eight estimations and corrected gene abundances by multiplying the relative abundance of each gene by this estimated average genome size of the sample. Finally, we used these corrected abundances and applied the same pathway-level comparative analysis to identify IBD-associated pathways. Since this normalization approach requires raw read counts, we could not apply it the other datasets analyzed in our study.
Measuring T2D-associated pathway recovery when using T2D cases with HMP or IBD controls
First, we identified the original T2D-associated pathways as described above. Next, we repeated the same process but as controls used either HMP stool samples or healthy samples from the IBD study. To compute the number of recovered pathways, we counted the number of pathways that were discovered both in the original setting and in the cross-study setting. In order to prevent a situation where high recovery stems from the identification of many pathways, we limited the number of discovered pathways in the cross-study setting to be the number of pathways discovered in the original setting.
Simulating microbial samples
To evaluate the performance of MUSiCC on a dataset in which the underlying KO and pathway abundances are known, we generated a dataset of synthetic metagenomic samples, following the procedure described in . Specifically, we generated 20 simulated metagenomic samples, each of which consisted of 500,000 101 bp reads generated at random from a collection of reference genomes that were randomly assigned different relative abundances (up to 100-fold) in each sample. To facilitate analysis of pathway level variation, we limited the set of genomes used to 21 bacterial genomes from KEGG that contained the entire set of KOs associated with the flagellar assembly pathway (Additional file 10: Table S5), and each sample harbored 10 genomes randomly selected from this set. We then mapped the simulated reads to known KOs and calculated the read count of each KO in each sample, as previously described . The underlying true average copy number for each KO was calculated based on the genomes included in each sample and weighted by their relative abundances.
Using average genome size estimation methods for sample normalization
where L is the read length (101 bp), and a, b, and c are parameters optimized in the original study (a = 21.2, b = 4230, c = 0.733) as described in Raes et al. . The effective genome sizes obtained were highly correlated with the real average genome sizes in each sample (r = 0.99; Pearson correlation). Using other sets of parameters described in the original study did not improve this correlation. Finally, we used these calculated effective genome sizes to normalize the measured KO abundances in each simulated sample.
Second, to apply the method introduced in , we installed the GAAS software (version 0.17). Since GAAS requires the complete BLAST alignment files, we ran BLAST for each simulated sample against the 21 genomes in our simulation (this represents both a best case scenario and is feasible computationally), using the same parameters as described in the original study: blastn -query < sample fasta file > -db < simulation genomes nucleotide database > -out < GAAS blast results file > -outfmt 6 -evalue 0.001 -gapopen 5 -gapextend 2 -word_size 11 -penalty −3 -reward 1. Next, we ran GAAS on the alignment results with the command: gaas -f < sample fasta file > -d < GAAS nucleotide database > -m < GAAS blast results file > -gt 0 -gp 0 -gs 0 -j 1, and obtained GAAS estimated average genome size in each sample. Overall, the median error in average gnome size estimation was <1%, comparable to the accuracy reported in the original study. These GAAS average genome size estimations were again used to normalize the measured KO abundances.
We thank Rogan Carr for valuable advice and for data on the mappability of genes used in this work. We thank Morgan Langille for information about the PICRUSt pipeline. We are grateful to all members of the Borenstein lab for helpful discussions. We thank the reviewers for their in-depth review and their suggestions for the improvement of our manuscript and software. This work was supported in part by NIH P30 DK089507 and by New Innovator Award DP2 AT007802-01 to EB.
- Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.View ArticlePubMed CentralPubMedGoogle Scholar
- Consortium THMP. Structure, function and diversity of the healthy human microbiome. Nature. 2013;486:207–14.Google Scholar
- Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40:D109–14.View ArticlePubMed CentralPubMedGoogle Scholar
- Qin J, Li R, Jeroen R, Arumugam M, Burgdorf KS, Manichanh C, et al. A human gut microbial gene catalogue established by metagenomic sequencing. Nature. 2010;464:59–65.View ArticlePubMed CentralPubMedGoogle Scholar
- Yatsunenko T, Rey FE, Manary MJ, Trehan I, Dominguez-Bello MG, Contreras M, et al. Human gut microbiome viewed across age and geography. Nature. 2012;486:222–7.PubMed CentralPubMedGoogle Scholar
- David LA, Maurice CF, Carmody RN, Gootenberg DB, Button JE, Wolfe BE, et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature. 2013;505:559–63.View ArticlePubMed CentralPubMedGoogle Scholar
- Abubucker S, Segata N, Goll J, Schubert AM, Izard J, Cantarel BL, et al. Metabolic reconstruction for metagenomic data and its application to the human Microbiome. PLoS Comput Biol. 2012;8:e1002358.View ArticlePubMed CentralPubMedGoogle Scholar
- Shafquat A, Joice R, Simmons SL, Huttenhower C. Functional and phylogenetic assembly of microbial communities in the human microbiome. Trends Microbiol. 2014;22:261–6.View ArticlePubMedGoogle Scholar
- Cani PD, Amar J, Iglesias MA, Poggi M, Knauf C, Bastelica D, et al. Metabolic endotoxemia initiates obesity and insulin resistance. Diabetes. 2007;56:1761–72.View ArticlePubMedGoogle Scholar
- McHardy IH, Goudarzi M, Tong M, Ruegger PM, Schwager E, Weger JR, et al. Integrative analysis of the microbiome and metabolome of the human intestinal mucosal surface reveals exquisite inter-relationships. Microbiome. 2013;1:17.View ArticlePubMed CentralPubMedGoogle Scholar
- Qin J, Li Y, Cai Z, Li S, Zhu J, Zhang F, et al. A metagenome-wide association study of gut microbiota in type 2 diabetes. Nature. 2013;490:55–60.View ArticleGoogle Scholar
- Aitchison J. The statistical analysis of compositional data. London: Chapman & Hall Ltd.; 1986.View ArticleGoogle Scholar
- Hattersley AT, Turner RC, Patel P, O’Rahilly S. Linkage of type 2 diabetes to the glucokinase gene. Lancet. 1992;339:1307–10.View ArticlePubMedGoogle Scholar
- Friedman J, Alm EJ. Inferring correlation networks from genomic survey data. PLoS Comput Biol. 2012;8:e1002687.View ArticlePubMed CentralPubMedGoogle Scholar
- Bouché C, Serdy S, Kahn CR, Goldfine AB. The cellular fate of glucose and its relevance in type 2 diabetes. Endocr Rev. 2004;25:807–30.View ArticlePubMedGoogle Scholar
- Paulson JN, Stine OC, Bravo HC, Pop M. Differential abundance analysis for microbial marker-gene surveys. Nat Meth. 2013;10:1200–2.View ArticleGoogle Scholar
- Morgan XC, Tickle TL, Sokol H, Gevers D, Devaney KL, Ward DV, et al. Dysfunction of the intestinal microbiome in inflammatory bowel disease and treatment. Genome Biol. 2012;13:R79.View ArticlePubMed CentralPubMedGoogle Scholar
- Beszteri B, Temperton B, Frickenhaus S, Giovannoni SJ. Average genome size: a potential source of bias in comparative metagenomics. ISME J. 2010;4:1075–7.View ArticlePubMedGoogle Scholar
- Frank JA, Sørensen SJ. Quantitative metagenomic analyses based on average genome size normalization. Appl Environ Microbiol. 2011;77:2513–21.View ArticlePubMed CentralPubMedGoogle Scholar
- Le Chatelier E, Nielsen T, Qin J, Prifti E, Hildebrand F, Falony G, et al. Richness of human gut microbiome correlates with metabolic markers. Nature. 2013;500:541–6.View ArticlePubMedGoogle Scholar
- Mathur R, Goyal D, Kim G, Barlow GM, Chua KS, Pimentel M. Methane-producing human subjects have higher serum glucose levels during oral glucose challenge than non-methane producers: a pilot study of the effects of enteric methanogens on glycemic regulation. Res J Endocrinol Metab. 2014;2:2.View ArticleGoogle Scholar
- Angly FE, Willner D, Prieto-Davó A, Edwards RA, Schmieder R, Vega-Thurber R, et al. The GAAS metagenomic tool and its estimations of viral and microbial average genome size in four major biomes. PLoS Comput Biol. 2009;5:e1000593.View ArticlePubMed CentralPubMedGoogle Scholar
- Raes J, Korbel JO, Lercher MJ, Von Mering C, Bork P. Prediction of effective genome size in metagenomic samples. Genome Biol. 2007;8:R10.View ArticlePubMed CentralPubMedGoogle Scholar
- Langille MGI, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol. 2013;31:814–21.View ArticlePubMedGoogle Scholar
- Benjamini Y, Speed TP. Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Res. 2012;40:e72.View ArticlePubMed CentralPubMedGoogle Scholar
- Oh J, Byrd AL, Deming C, Conlan S, Program NCS, Kong HH, et al. Biogeography and individuality shape function in the human skin metagenome. Nature. 2014;514:59–64.View ArticlePubMed CentralPubMedGoogle Scholar
- Wylie KM, Mihindukulasuriya KA, Zhou Y, Sodergren E, Storch GA, Weinstock GM. Metagenomic analysis of double-stranded DNA viruses in healthy adults. BMC Biol. 2014;12:71.View ArticlePubMed CentralPubMedGoogle Scholar
- Minot S, Bryson A, Chehoud C, Wu GD, Lewis JD, Bushman FD. Rapid evolution of the human gut virome. Proc Natl Acad Sci U S A. 2013;110:12450–5.View ArticlePubMed CentralPubMedGoogle Scholar
- Greenblum S, Turnbaugh PJ, Borenstein E. Metagenomic systems biology of the human gut microbiome reveals topological shifts associated with obesity and inflammatory bowel disease. Proc Natl Acad Sci U S A. 2012;109:594–9.View ArticlePubMed CentralPubMedGoogle Scholar
- Darling AE, Jospin G, Lowe E, Matsen IV FA, Bik HM, Eisen JA. PhyloSift: phylogenetic analysis of genomes and metagenomes. Peerl. 2014;2:e243.View ArticleGoogle Scholar
- Katoh K, Misawa K, Kuma K-I, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.View ArticlePubMed CentralPubMedGoogle Scholar
- Carr R, Borenstein E. Comparative analysis of functional metagenomic annotation and the mappability of short reads. PLoS One. 2014;9:e105776.View ArticlePubMed CentralPubMedGoogle Scholar
- Friedman J, Hastie T, Tibshirani R. glmnet: lasso and elastic-net regularized generalized linear models. J Stat Softw. 2010;33:1–22.PubMed CentralPubMedGoogle Scholar
- Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Methodol. 1996;58:267–88.Google Scholar
- Zou H, Hastie T. Regularization and variable selection via the elastic net. J R Stat Soc B Stat Meth. 2005;67:301–20.View ArticleGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.