BACKGROUND
microbial diversity within the human body has recently been quantified through 16s rrna surveys  <cit>  and metagenomic methods.
the latter provide a detailed view of the genomic composition and functional potential of human-associated microbial communities through shotgun sequencing  <cit> .
however this level of resolution comes with a high price-tag - billions of base-pairs need to be sequenced to ensure a sufficient level of sampling of complex communities  <cit>  such as those found in the human gastrointestinal tract.
16s rrna surveys provide limited insight into the composition of the commensal microbiome.
however, due to substantially lower costs, such studies are currently the only practical approach for studying large numbers of samples .
in this paper, we explore the limitations of unsupervised clustering methods used to analyze 16s rrna data, particularly the large impact of small changes in the parameters of the analysis process.
we specifically focus on the most common strategy - the clustering of 16s rrna sequences into a collection of operational taxonomic units  or phylotypes on the basis of sequence similarity.
an evaluation of taxonomic classification through database searches  <cit>  or other fully supervised classification methods  <cit>  are beyond the scope of this paper.
fully supervised approaches are inherently limited due to the current undersampling of the global microbial population, only allowing accurate classification of a fraction of sequences .
as an alternative, the unsupervised clustering of 16s sequences allows for detection of species-like units in complex unknown bacterial environments, even if a precise taxonomic identity cannot be assigned to these units.
the otu clustering process typically begins by constructing a multiple alignment  of the 16s rrna sequences.
the msa is then used to estimate pairwise distances between individual sequences, expressed as the fraction of nucleotides that have changed as the sequences have evolved from their most recent common ancestor.
to accurately reflect evolutionary processes, the distances inferred from the msa are corrected using one of several models of evolution  <cit> .
the distances are provided as input to a hierarchical clustering algorithm .
sub-clusters or otus are defined by applying a distance threshold, selected to roughly approximate a specific taxonomic level: thresholds between 1-3% are typically used to approximate individual species, 5% for individual genera, 15% for classes, etc.
<cit> .
note that alternative approaches to otu creation exist   <cit> , or use databases  <cit> ), but the general technique we describe above has been widely employed across bacterial diversity studies.
the choice of msa, distance correction, clustering algorithm, and distance threshold vary considerably between studies, and to our knowledge, there has been no rigorous evaluation of the impact of methodological choices on the ecological conclusions of the analysis process.
previously, simulated datasets have been successfully used to evaluate methods for the assembly and binning of metagenomic data  <cit> .
in this study, we rely on simulated datasets to provide a comprehensive assessment of the extent to which individual parameters in the otu clustering process affect the estimated diversity and composition of a microbial environment.
we evaluate methodological choices in terms of how well the clustering of sequences into a set of otus matches the clustering imposed by the known  membership of the sequences to individual bacterial species.
the otu clusterings are compared using a mathematically robust metric - the variation of information   <cit>  - an information-theoretic measure of the amount information lost or gained by changing from one clustering to another.
through this evaluation framework, we will demonstrate that otus are highly sensitive to small changes in the clustering methodology and reveal a surprising observation that reducing the stringency of clustering distance thresholds tends to produce more accurate species-level representations of a community.
we will further assess the impact of otu variability on common ecological measures of diversity and provide an example of how semi-supervised clustering could produce more robust otu structures by accounting for varying evolutionary rates across the microbial phylogeny.
methods
creation of simulated datasets
our simulated dataset comprises  <dig> 16s rrna sequences from the rdp database   <cit> , that satisfy the following properties:  at least  <dig> bp long;  can be aligned by nast  <cit>  requiring a match to the template alignment of at least 75% identity;  have full taxonomic identification at all levels from phylum to species according to the rdp taxonomy ; and  the taxonomic identity of each of the sequences is confirmed by the rdp naïve bayesian classifier  <cit>  at ≥ 95% confidence, greengenes simrank  <cit> , and through blastn  <cit>  searches against a reduced rdp database .
thus, while it is impossible to verify the correct taxonomic membership of all these sequences, we can guarantee that their annotation is consistent across multiple databases and classification procedures.
these sequences were largely obtained from bacterial isolates  and had unambiguous taxonomic assignment at the species level.
our simulated environment spans  <dig> species,  <dig> genera,  <dig> families,  <dig> orders,  <dig> classes, and seven phyla including proteobacteria, bacteroidetes, firmicutes, and actinobacteria.
alpha-, beta-, and gammaproteobacteria make up 66% of the sequences in roughly equal proportions.
a similar class distribution has been reported for microbial communities found in the phyllosphere of the atlantic rainforest  <cit> .
it is important to note that, by choice, these sequences are high quality and belong to relatively well-characterized taxonomic groups.
therefore any results obtained on this highly curated dataset represent an upper bound on the performance that can be achieved when analyzing noisy data from environmental surveys.
parameter evaluation: multiple sequence alignments, distance corrections, clustering methods and distance thresholds
all  <dig> sequences were aligned using muscle  <cit> , clustalw  <cit> , and nast  <cit>  using default parameters.
clustalw was run with the "fast" option for pairwise alignments, a heuristic setting that dramatically improves running time  at the cost of a lower quality alignment.
in the nast alignment, all columns containing only gaps were removed, and each msa was trimmed so that every sequence spanned the entire alignment.
the trimmed msas covered the range of v <dig>  v <dig>  and v <dig> hypervariable regions within the e. coli o157:h <dig> str.
tw <dig> 16s rrna gene.
distance matrices were constructed with dnadist from the phylip package  <cit>  using each of the jukes-cantor, kimura 2-parameter, and felsenstein <dig> corrections, keeping the remaining parameters at their default values .
distance matrices were then provided as input to dotur  <cit> , and the resulting data clustered using, in turn, nearest neighbor, average neighbor, and furthest neighbor clustering procedures on each individual matrix.
individual clusters were generated, for each clustering algorithm, by varying a distance threshold parameter ranging from  <dig>  to  <dig>  .
the distance threshold d has a different meaning depending on the particular algorithm: in furthest/average/nearest neighbor, two clusters are merged if the maximum/average/minimum distance between any two elements in the combined cluster is less than or equal to d.
the process described above generated  <dig> sets of clusters/otus ,  <dig> of which are distinct .
comparing clusterings
we employed the variation of information  metric  <cit>  as a measure of similarity between two partitions  of a given set  <cit> .
for this study, the set comprises the  <dig> 16s sequences selected for the artificial environmental sample.
mathematically, a given clustering c, is a partition of a set s into disjoint subsets  where:
if there are m elements in set s, and we let mi be the number of elements in cluster ci, then to compute the variation of information between two clusterings, we first find the probability that a randomly selected sequence is in a particular cluster, that is, .
given this discrete probability distribution, the uncertainty of the random variable i is the entropy associated with clustering c, defined as:
now, suppose we have two clusterings c = {c <dig>  c <dig>  ..., cm}, and d = {d <dig>  d <dig>  ..., dm'}.
then we calculate the joint distribution  describing the similarity of all pairs of clusters between c and d. the mutual information between the clusterings c and d is then defined to be , and finally, the variation of information between c and d is defined as the sum of the individual clustering entropies less  <dig> times the mutual information:
if c and d are identical clusterings, then h = h = i, and the vi =  <dig>  the vi distance is a true metric, satisfying symmetry, non-negativity, and the triangle inequality.
in order to provide a reference set of vi distances for known clusters, we measured the vi between the "true" species clustering  and the annotated phylum, class, order, family, and genus clusterings .
vi-cut method for defining otus
vi-cut is a procedure that finds a clustering from a hierarchical tree decomposition t that optimally matches a subset of known labels, as defined by the variation of information metric  <cit> .
a clustering is defined by choosing a set of nodes in t. each chosen node c corresponds to a single cluster consisting of all the leaves  in the subtree rooted at c. collectively, the chosen nodes correspond to a node-cut k, which induces a non-overlapping clustering ak.
let d represent the partial clustering of labeled sequences such that sequences with the same label are grouped together.
the vi-cut algorithm finds the ak that minimizes the vi distance to d:
although there are exponentially many possible node-cuts in t, an optimal node-cut can be found efficiently using dynamic programming  <cit> .
as described, the vi-cut algorithm creates overly-general clusters when faced with sparsely labeled data.
to correct for this, we used a strategy called forbidden nodes.
specifically, any node n with a corresponding diameter ≥  <dig>  was labeled as "forbidden", and we required that the clusters produced by vi-cut do not contain any forbidden nodes.
this requirement implicitly restricts the maximum diameter of an otu to 7% divergence.
to incorporate forbidden nodes into the vi-cut algorithm, we first ran the standard vi-cut algorithm.
any cluster that contained a forbidden node n was then sub-divided into non-overlapping clusters by identifying subtrees dominated by n that do not contain any forbidden nodes.
assessment of otu methodology factors using anova
to isolate the individual impact of each component in an otu methodology, the  <dig> methodologies resulting in the lowest vi distance from the annotated species clustering were analyzed using multi-way analysis of variance  considering four factors: multiple sequence alignment, evolutionary distance correction, clustering algorithm, and distance threshold.
RESULTS
otu variability
the results of our analysis  reveal a large variation in the level of concordance of the otu clustering with the annotated species-level composition of the environment when varying methodology parameters.
we specifically highlight the parameters with most impact - distance threshold , msa , and clustering strategy  - parameters that accounted for 56%, 33%, and 7% of the total variation, respectively .
distance correction measures did not significantly contribute to the variability , implying that such corrections are roughly equivalent within the short evolutionary time frame defining a microbial species.
we extended this comparison to include the olsen distance correction in arb  <cit> , which we found produced otus virtually identical to those created using the f <dig> correction.
the factor with the largest effect on the quality of the otus was the distance threshold, followed by the msa, and then the clustering algorithm.
the distance correction explained <  <dig> % of the variance and no statistically significant difference could be detected between the corrections .
none of the parameter combinations perfectly captured the annotated species composition .
large variation in the otu content is observed even when we fix the similarity threshold to  <dig>   - the number of otus ranges from  <dig> to  <dig> at this similarity level, depending on the choice of msa or clustering strategy.
surprisingly, the otu clustering closest to the annotated species clustering was obtained using a similarity threshold of  <dig>  - a value larger than the cutoffs usually used to approximate the species-level composition of an environment .
in terms of alignment, methodologies employing clustalw or nast were roughly similar and performed slightly better than those using muscle .
however, applying our validation methodology to ten randomly selected subsets of the original data, we discovered that no msa program consistently outperformed the others .
in terms of clustering strategy, furthest neighbor resulted in the best agreement with the annotated species structure of our environment .
even the combination of analysis parameters with the lowest vi distance  led to an overestimate of the number of species in our sample, resulting in  <dig> otus.
this result highlights a fundamental limitation of hierarchical clustering strategies for 16s rrna analysis - only  <dig> of the  <dig> species present in our sample corresponded to a homogeneous sub-tree within the best hierarchical clustering of the data.
the remaining  <dig> species cannot be correctly clustered irrespective of the similarity threshold chosen.
the results presented above highlight a wide variation in the otu structure as we explore the parameters of the analysis process.
to determine whether such variation is also present in the methodologies used in practice, we compared three analysis methodologies that performed well in our combinatorial search to several methodologies reported in published literature.
the results shown in table  <dig> indicate that the published methodologies can overestimate the diversity of the simulated environment, sometimes by more than 3-fold .
the fragmentation of the resulting otus is particularly striking among the most abundant phylotypes , where sequences belonging to the same species are distributed among multiple otus.
the "vi" column indicates the vi distance of each clustering from the annotated species clustering.
best methods chosen by our validation methodology for each msa  are contrasted with five recently published methodologies.
the "correction" column corresponds to the evolutionary distance correction.
note that for the optimal methods using clustalw and nast alignments, the f <dig> and k2p corrections produced identical otu sets because the distance matrices were very similar, though not identical.
all methods in this table used furthest neighbor  clustering.
the ace, chao <dig>  and shannon diversity estimators are also provided.
nonparametric estimators of richness and diversity
the large variability in the otu estimates produced by different methodologies also had a significant effect on commonly inferred ecological parameters.
the chao <dig>   <cit>  and ace   <cit>  richness estimators and the shannon diversity  index  <cit>  are measures commonly used to approximate the level of diversity present in an environment.
restricting our analysis to methodologies with distance thresholds from  <dig> - <dig> , these three measures were found to be highly sensitive to differences in otu structure.
under the true  species clustering, sace =  <dig>  schao <dig> =  <dig>  and h =  <dig> .
sace and schao <dig> estimates for the computed otu clusterings ranged from  <dig> to  <dig> and  <dig> to  <dig> phylotypes, respectively, while shannon diversity indices ranged from  <dig>  to  <dig>  .
accurate estimates of the diversity of an environment are essential when planning whole-genome shotgun metagenomics sequencing projects, and an environment whose diversity is incorrectly perceived to be too high might not be studied at a metagenomic level.
partial masking of msas
to improve phylogenetic analyses, researchers often remove hypervariable segments of msas either manually or using a filter such as lanemask  <cit> .
we explored the impact of this approach on otu clustering.
specifically, we used the greengenes lanemask filter, which reduces a nast alignment to  <dig> highly conserved columns.
the results are surprising - on average, employing lanemask resulted in a worse approximation of the true species composition than the unmasked alignment .
this suggests that the use of a generic mask might be inappropriate when analyzing highly similar sequences with the goal of constructing otus, though its use may still be relevant to phylogenetic analyses of divergent sequences.
semi-supervised clustering alternatives
our study has so far made the assumption that one of the primary goals of a 16s analysis pipeline is to estimate the composition of an environment at a pre-specified taxonomic level .
as demonstrated by our results, the otu methodologies proposed in the literature fail to achieve this goal, generally overestimating the number of species.
even by systematically evaluating various settings for the parameters of the analysis process, we could not obtain perfect concordance between the otu structure and the species composition of the environment.
this is in part due to the fact that the concept of "species" is born out of gross morphological and phenotypic traits of microorganisms, and therefore cannot be precisely mapped to fine-scale molecular measurements.
furthermore, the rate of evolution varies across the tree of life, making it unrealistic to rely on a single distance threshold.
as an alternative, we investigated the use of a semi-supervised clustering method to adaptively select a set of local distance thresholds that lead to otus that better fit the species composition of the environment.
specifically, we employed vi-cut  <cit> , a clustering approach that identifies a cut within a hierarchical clustering tree that maximizes the fit with a labeled subset of the sequences.
in the case of 16s analysis, vi-cut constructs a set of otus that optimally matches  the species structure of an environment as inferred from a small subset of sequences that have known taxonomic assignments .
we applied vi-cut to our data by simulating partial taxonomic knowledge of the dataset.
for each msa and the optimal distance correction , we randomly selected 10% of the sequences and provided vi-cut with their true labels.
to assess the variability in the algorithm's results, we repeated this procedure  <dig> times.
as seen in figure 5a, vi-cut outperforms methodologies that employ a single distance threshold, irrespective of the msa employed or the random selection of labeled sequences.
the need for an adaptive threshold  is highlighted in figure 5b - the diameter of clusters corresponding to a single species in our data varies considerably among our sequences  and the semi-supervised learning algorithm implemented in vi-cut is able to closely approximate the true distribution of distance thresholds.
note that perfect concordance between otus and species cannot be achieved even with the best hierarchical clustering tree constructed from our data.
it is an open question whether other clustering approaches could perform better in this context.
consistency of methods across multiple datasets
to investigate the consistent improvement of the vi-cut methodology over other methods, we created ten additional 16s environmental samples - each sample containing  <dig> randomly selected sequences from the original dataset.
we repeated our comparison of vi-cut to other methods for these ten simulated samples.
examining the results across each msa, we found that vi-cut consistently produced the best species-level approximation compared to standard methodologies .
importantly, we point out that no msa program was consistently optimal across all datasets.
methods are ranked by their mean vi-distance over  <dig> simulated datasets.
we constrained the results to commonly accepted methods using furthest neighbor clustering and distance thresholds less than  <dig> .
a distance threshold of  <dig>  is consistently among the worst performing methodologies.
vi-cut consistently results in the best clustering for each msa.
discussion and
CONCLUSIONS
in this study we have shown how small differences in otu methodologies can lead to significant variability in the resulting otu structure, thereby affecting estimates of microbial diversity and ecological conclusions.
our results indicate that the most important factor in an otu methodology is the distance threshold imposed during clustering, and that small changes in this parameter lead to a substantial variance in the estimated diversity of a community, making it difficult or even impossible to directly compare the results of studies utilizing different thresholds.
commonly employed thresholds of  <dig> - <dig>   fail to capture the underlying species composition of an environment and are frequently too stringent, producing inflated estimates of diversity.
this result is in large part due to the fact that the overly-general biological definition of a species cannot be directly mapped at the fine-scale resolution provided by molecular-level observations.
in addition, diversity estimates are highly sensitive to the abundance of rare members of a community and, thus, can be easily confused by the "noise" caused by sequencing errors, transient organisms, or "naked" dna not originating from one of the members of the community.
for all of these reasons, aggregate diversity measures do not necessarily correlate with the biological functions performed by members of a community, and must be augmented by additional, more specific, measurements of the community structure .
while the community simulated in our study was dominated by proteobacteria, a significant proportion of the data belongs to other microbial phyla, thus the general results of our study likely hold in most other datasets.
furthermore, the organisms present in our samples are sufficiently distant at the 16s rrna level to allow unambiguous clustering once a distance threshold was selected, i.e.
the results we obtained reflect actual characteristics of 16s rrna data rather than artifacts of the analysis process.
our simulated samples were composed of a small number of high-quality sequences with exceptional length  - a relatively simple challenge compared to current datasets generated by pyrosequencing .
16s surveys employing pyrosequencing technologies generate considerably larger datasets , comprising reads of shorter lengths .
the results we describe for long reads will only be exacerbated in the context of short-read data.
further, these data present several new computational challenges for otu clustering including the need for faster alignment algorithms, the computation and storage of evolutionary distances , and the selection of distance thresholds appropriate for sequences with significantly less phylogenetic information than the sanger-based reads used in our study.
the analysis of the new sequencing data require the development of robust methods for otu clustering, as well as rigorous validation of analysis methods.
we finally demonstrated that a semi-supervised clustering approach  can significantly improve analysis quality, highlighting the potential for semi-supervised clustering approaches to fill the gap between the two extreme classes of approaches commonly used in analyzing 16s rrna datasets: fully-supervised database searches that rely on expensive highly-curated datasets; and fully unsupervised otu clustering procedures that use ad hoc similarity thresholds in an attempt to match poorly-defined taxonomic labels.
semi-supervised procedures such as vi-cut can generate accurate clusterings while only requiring high-quality labels for only a small subset of the sequences.
finally, it is important to observe that clustering 16s rrna sequences into a set of otus is a valuable analysis tool even if the resulting otus do not correlate with pre-defined taxonomic entities.
however, the ad hoc choice of analysis parameters, in particular the selection of different distance thresholds, complicates cross-study comparisons and even basic descriptions of diversity.
as 16s rrna surveys are increasingly applied in a clinical setting  <cit>  , it is critical to accurately measure taxonomic diversity and identify individual species.
our results highlight the need for standardizing 16s rrna analysis methods, or in the very least, reporting results obtained with multiple distance thresholds or clustering algorithms.
the data used in this study have been deposited in the fames online database http://fames.jgi-psf.org - a repository for metagenomic analysis benchmarks  <cit> .
software used for semi-supervised clustering analysis is also available at http://www.cbcb.umd.edu/vicut/.
abbreviations
otu: operational taxonomic unit; rrna: ribosomal ribonucleic acid; vi: variation of information; msa: multiple sequence alignment; rdp: ribosomal database project; anova: analysis of variance;
authors' contributions
jrw and sn ran the experiments and analyzed the data.
m-rg provided computational expertise.
jw, sn, nn, m-rg, ck and mp wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
variation of information distances of high-level taxonomic clusterings from the annotated species clustering.
to give the reader some intuition about the vi distance metric, we computed vi distances between the annotated species-level clustering and other clusterings based on phylum, class, order, family, and genus annotations.
this file contains a table of these reference distances.
click here for file
acknowledgements
we thank pawel gajer and jacques ravel for helpful conversations.
this work was supported in part by the bill and melinda gates foundation , the nih , and the nsf .
ck and sn partially supported by nsf grant  <dig> to ck.

 microbial diversity within the human body has recently been quantified through 16s rrna surveys and metagenomic methods. the latter provide a detailed view of the genomic composition and functional potential of human-associated microbial communities through shotgun sequencing. however this level of resolution comes with a high price-tag - billions of base-pairs need to be sequenced to ensure a sufficient level of sampling of complex communities such as those found in the human gastrointestinal tract. 16s rrna surveys provide limited insight into the composition of the commensal microbiome. however, due to substantially lower costs, such studies are currently the only practical approach for studying large numbers of samples. in this paper, we explore the limitations of unsupervised clustering methods used to analyze 16s rrna data, particularly the large impact of small changes in the parameters of the analysis process. we specifically focus on the most common strategy - the clustering of 16s rrna sequences into a collection of operational taxonomic units or phylotypes on the basis of sequence similarity. an evaluation of taxonomic classification through database searches or other fully supervised classification methods are beyond the scope of this paper. fully supervised approaches are inherently limited due to the current undersampling of the global microbial population, only allowing accurate classification of a fraction of sequences. as an alternative, the unsupervised clustering of 16s sequences allows for detection of species-like units in complex unknown bacterial environments, even if a precise taxonomic identity cannot be assigned to these units. the otu clustering process typically begins by constructing a multiple alignment of the 16s rrna sequences. the msa is then used to estimate pairwise distances between individual sequences, expressed as the fraction of nucleotides that have changed as the sequences have evolved from their most recent common ancestor. to accurately reflect evolutionary processes, the distances inferred from the msa are corrected using one of several models of evolution. the distances are provided as input to a hierarchical clustering algorithm. sub-clusters or otus are defined by applying a distance threshold, selected to roughly approximate a specific taxonomic level: thresholds between 1-3% are typically used to approximate individual species, 5% for individual genera, 15% for classes, etc.. note that alternative approaches to otu creation exist, or use databases ), but the general technique we describe above has been widely employed across bacterial diversity studies. the choice of msa, distance correction, clustering algorithm, and distance threshold vary considerably between studies, and to our knowledge, there has been no rigorous evaluation of the impact of methodological choices on the ecological conclusions of the analysis process. previously, simulated datasets have been successfully used to evaluate methods for the assembly and binning of metagenomic data. in this study, we rely on simulated datasets to provide a comprehensive assessment of the extent to which individual parameters in the otu clustering process affect the estimated diversity and composition of a microbial environment. we evaluate methodological choices in terms of how well the clustering of sequences into a set of otus matches the clustering imposed by the known membership of the sequences to individual bacterial species. the otu clusterings are compared using a mathematically robust metric - the variation of information - an information-theoretic measure of the amount information lost or gained by changing from one clustering to another. through this evaluation framework, we will demonstrate that otus are highly sensitive to small changes in the clustering methodology and reveal a surprising observation that reducing the stringency of clustering distance thresholds tends to produce more accurate species-level representations of a community. we will further assess the impact of otu variability on common ecological measures of diversity and provide an example of how semi-supervised clustering could produce more robust otu structures by accounting for varying evolutionary rates across the microbial phylogeny. methods creation of simulated datasets our simulated dataset comprises 16s rrna sequences from the rdp database, that satisfy the following properties: at least bp long; can be aligned by nast requiring a match to the template alignment of at least 75% identity; have full taxonomic identification at all levels from phylum to species according to the rdp taxonomy; and the taxonomic identity of each of the sequences is confirmed by the rdp naïve bayesian classifier at ≥ 95% confidence, greengenes simrank, and through blastn searches against a reduced rdp database. thus, while it is impossible to verify the correct taxonomic membership of all these sequences, we can guarantee that their annotation is consistent across multiple databases and classification procedures. these sequences were largely obtained from bacterial isolates and had unambiguous taxonomic assignment at the species level. our simulated environment spans species, genera, families, orders, classes, and seven phyla including proteobacteria, bacteroidetes, firmicutes, and actinobacteria. alpha-, beta-, and gammaproteobacteria make up 66% of the sequences in roughly equal proportions. a similar class distribution has been reported for microbial communities found in the phyllosphere of the atlantic rainforest. it is important to note that, by choice, these sequences are high quality and belong to relatively well-characterized taxonomic groups. therefore any results obtained on this highly curated dataset represent an upper bound on the performance that can be achieved when analyzing noisy data from environmental surveys. parameter evaluation: multiple sequence alignments, distance corrections, clustering methods and distance thresholds all sequences were aligned using muscle, clustalw, and nast using default parameters. clustalw was run with the "fast" option for pairwise alignments, a heuristic setting that dramatically improves running time at the cost of a lower quality alignment. in the nast alignment, all columns containing only gaps were removed, and each msa was trimmed so that every sequence spanned the entire alignment. the trimmed msas covered the range of v v and v hypervariable regions within the e. coli o157:h str. tw 16s rrna gene. distance matrices were constructed with dnadist from the phylip package using each of the jukes-cantor, kimura 2-parameter, and felsenstein corrections, keeping the remaining parameters at their default values. distance matrices were then provided as input to dotur, and the resulting data clustered using, in turn, nearest neighbor, average neighbor, and furthest neighbor clustering procedures on each individual matrix. individual clusters were generated, for each clustering algorithm, by varying a distance threshold parameter ranging from to. the distance threshold d has a different meaning depending on the particular algorithm: in furthest/average/nearest neighbor, two clusters are merged if the maximum/average/minimum distance between any two elements in the combined cluster is less than or equal to d. the process described above generated sets of clusters/otus, of which are distinct. comparing clusterings we employed the variation of information metric as a measure of similarity between two partitions of a given set. for this study, the set comprises the 16s sequences selected for the artificial environmental sample. mathematically, a given clustering c, is a partition of a set s into disjoint subsets where: if there are m elements in set s, and we let mi be the number of elements in cluster ci, then to compute the variation of information between two clusterings, we first find the probability that a randomly selected sequence is in a particular cluster, that is,. given this discrete probability distribution, the uncertainty of the random variable i is the entropy associated with clustering c, defined as: now, suppose we have two clusterings c = {c c..., cm}, and d = {d d..., dm'}. then we calculate the joint distribution describing the similarity of all pairs of clusters between c and d. the mutual information between the clusterings c and d is then defined to be, and finally, the variation of information between c and d is defined as the sum of the individual clustering entropies less times the mutual information: if c and d are identical clusterings, then h = h = i, and the vi = the vi distance is a true metric, satisfying symmetry, non-negativity, and the triangle inequality. in order to provide a reference set of vi distances for known clusters, we measured the vi between the "true" species clustering and the annotated phylum, class, order, family, and genus clusterings. vi-cut method for defining otus vi-cut is a procedure that finds a clustering from a hierarchical tree decomposition t that optimally matches a subset of known labels, as defined by the variation of information metric. a clustering is defined by choosing a set of nodes in t. each chosen node c corresponds to a single cluster consisting of all the leaves in the subtree rooted at c. collectively, the chosen nodes correspond to a node-cut k, which induces a non-overlapping clustering ak. let d represent the partial clustering of labeled sequences such that sequences with the same label are grouped together. the vi-cut algorithm finds the ak that minimizes the vi distance to d: although there are exponentially many possible node-cuts in t, an optimal node-cut can be found efficiently using dynamic programming. as described, the vi-cut algorithm creates overly-general clusters when faced with sparsely labeled data. to correct for this, we used a strategy called forbidden nodes. specifically, any node n with a corresponding diameter ≥ was labeled as "forbidden", and we required that the clusters produced by vi-cut do not contain any forbidden nodes. this requirement implicitly restricts the maximum diameter of an otu to 7% divergence. to incorporate forbidden nodes into the vi-cut algorithm, we first ran the standard vi-cut algorithm. any cluster that contained a forbidden node n was then sub-divided into non-overlapping clusters by identifying subtrees dominated by n that do not contain any forbidden nodes. assessment of otu methodology factors using anova to isolate the individual impact of each component in an otu methodology, the methodologies resulting in the lowest vi distance from the annotated species clustering were analyzed using multi-way analysis of variance considering four factors: multiple sequence alignment, evolutionary distance correction, clustering algorithm, and distance threshold.  otu variability the results of our analysis reveal a large variation in the level of concordance of the otu clustering with the annotated species-level composition of the environment when varying methodology parameters. we specifically highlight the parameters with most impact - distance threshold, msa, and clustering strategy - parameters that accounted for 56%, 33%, and 7% of the total variation, respectively. distance correction measures did not significantly contribute to the variability, implying that such corrections are roughly equivalent within the short evolutionary time frame defining a microbial species. we extended this comparison to include the olsen distance correction in arb, which we found produced otus virtually identical to those created using the f correction. the factor with the largest effect on the quality of the otus was the distance threshold, followed by the msa, and then the clustering algorithm. the distance correction explained < % of the variance and no statistically significant difference could be detected between the corrections. none of the parameter combinations perfectly captured the annotated species composition. large variation in the otu content is observed even when we fix the similarity threshold to - the number of otus ranges from to at this similarity level, depending on the choice of msa or clustering strategy. surprisingly, the otu clustering closest to the annotated species clustering was obtained using a similarity threshold of - a value larger than the cutoffs usually used to approximate the species-level composition of an environment. in terms of alignment, methodologies employing clustalw or nast were roughly similar and performed slightly better than those using muscle. however, applying our validation methodology to ten randomly selected subsets of the original data, we discovered that no msa program consistently outperformed the others. in terms of clustering strategy, furthest neighbor resulted in the best agreement with the annotated species structure of our environment. even the combination of analysis parameters with the lowest vi distance led to an overestimate of the number of species in our sample, resulting in otus. this result highlights a fundamental limitation of hierarchical clustering strategies for 16s rrna analysis - only of the species present in our sample corresponded to a homogeneous sub-tree within the best hierarchical clustering of the data. the remaining species cannot be correctly clustered irrespective of the similarity threshold chosen. the results presented above highlight a wide variation in the otu structure as we explore the parameters of the analysis process. to determine whether such variation is also present in the methodologies used in practice, we compared three analysis methodologies that performed well in our combinatorial search to several methodologies reported in published literature. the results shown in table indicate that the published methodologies can overestimate the diversity of the simulated environment, sometimes by more than 3-fold. the fragmentation of the resulting otus is particularly striking among the most abundant phylotypes, where sequences belonging to the same species are distributed among multiple otus. the "vi" column indicates the vi distance of each clustering from the annotated species clustering. best methods chosen by our validation methodology for each msa are contrasted with five recently published methodologies. the "correction" column corresponds to the evolutionary distance correction. note that for the optimal methods using clustalw and nast alignments, the f and k2p corrections produced identical otu sets because the distance matrices were very similar, though not identical. all methods in this table used furthest neighbor clustering. the ace, chao and shannon diversity estimators are also provided. nonparametric estimators of richness and diversity the large variability in the otu estimates produced by different methodologies also had a significant effect on commonly inferred ecological parameters. the chao and ace richness estimators and the shannon diversity index are measures commonly used to approximate the level of diversity present in an environment. restricting our analysis to methodologies with distance thresholds from -, these three measures were found to be highly sensitive to differences in otu structure. under the true species clustering, sace = schao = and h =. sace and schao estimates for the computed otu clusterings ranged from to and to phylotypes, respectively, while shannon diversity indices ranged from to. accurate estimates of the diversity of an environment are essential when planning whole-genome shotgun metagenomics sequencing projects, and an environment whose diversity is incorrectly perceived to be too high might not be studied at a metagenomic level. partial masking of msas to improve phylogenetic analyses, researchers often remove hypervariable segments of msas either manually or using a filter such as lanemask. we explored the impact of this approach on otu clustering. specifically, we used the greengenes lanemask filter, which reduces a nast alignment to highly conserved columns. the results are surprising - on average, employing lanemask resulted in a worse approximation of the true species composition than the unmasked alignment. this suggests that the use of a generic mask might be inappropriate when analyzing highly similar sequences with the goal of constructing otus, though its use may still be relevant to phylogenetic analyses of divergent sequences. semi-supervised clustering alternatives our study has so far made the assumption that one of the primary goals of a 16s analysis pipeline is to estimate the composition of an environment at a pre-specified taxonomic level. as demonstrated by our results, the otu methodologies proposed in the literature fail to achieve this goal, generally overestimating the number of species. even by systematically evaluating various settings for the parameters of the analysis process, we could not obtain perfect concordance between the otu structure and the species composition of the environment. this is in part due to the fact that the concept of "species" is born out of gross morphological and phenotypic traits of microorganisms, and therefore cannot be precisely mapped to fine-scale molecular measurements. furthermore, the rate of evolution varies across the tree of life, making it unrealistic to rely on a single distance threshold. as an alternative, we investigated the use of a semi-supervised clustering method to adaptively select a set of local distance thresholds that lead to otus that better fit the species composition of the environment. specifically, we employed vi-cut, a clustering approach that identifies a cut within a hierarchical clustering tree that maximizes the fit with a labeled subset of the sequences. in the case of 16s analysis, vi-cut constructs a set of otus that optimally matches the species structure of an environment as inferred from a small subset of sequences that have known taxonomic assignments. we applied vi-cut to our data by simulating partial taxonomic knowledge of the dataset. for each msa and the optimal distance correction, we randomly selected 10% of the sequences and provided vi-cut with their true labels. to assess the variability in the algorithm's results, we repeated this procedure times. as seen in figure 5a, vi-cut outperforms methodologies that employ a single distance threshold, irrespective of the msa employed or the random selection of labeled sequences. the need for an adaptive threshold is highlighted in figure 5b - the diameter of clusters corresponding to a single species in our data varies considerably among our sequences and the semi-supervised learning algorithm implemented in vi-cut is able to closely approximate the true distribution of distance thresholds. note that perfect concordance between otus and species cannot be achieved even with the best hierarchical clustering tree constructed from our data. it is an open question whether other clustering approaches could perform better in this context. consistency of methods across multiple datasets to investigate the consistent improvement of the vi-cut methodology over other methods, we created ten additional 16s environmental samples - each sample containing randomly selected sequences from the original dataset. we repeated our comparison of vi-cut to other methods for these ten simulated samples. examining the results across each msa, we found that vi-cut consistently produced the best species-level approximation compared to standard methodologies. importantly, we point out that no msa program was consistently optimal across all datasets. methods are ranked by their mean vi-distance over simulated datasets. we constrained the results to commonly accepted methods using furthest neighbor clustering and distance thresholds less than. a distance threshold of is consistently among the worst performing methodologies. vi-cut consistently results in the best clustering for each msa. discussion and  in this study we have shown how small differences in otu methodologies can lead to significant variability in the resulting otu structure, thereby affecting estimates of microbial diversity and ecological conclusions. our results indicate that the most important factor in an otu methodology is the distance threshold imposed during clustering, and that small changes in this parameter lead to a substantial variance in the estimated diversity of a community, making it difficult or even impossible to directly compare the results of studies utilizing different thresholds. commonly employed thresholds of - fail to capture the underlying species composition of an environment and are frequently too stringent, producing inflated estimates of diversity. this result is in large part due to the fact that the overly-general biological definition of a species cannot be directly mapped at the fine-scale resolution provided by molecular-level observations. in addition, diversity estimates are highly sensitive to the abundance of rare members of a community and, thus, can be easily confused by the "noise" caused by sequencing errors, transient organisms, or "naked" dna not originating from one of the members of the community. for all of these reasons, aggregate diversity measures do not necessarily correlate with the biological functions performed by members of a community, and must be augmented by additional, more specific, measurements of the community structure. while the community simulated in our study was dominated by proteobacteria, a significant proportion of the data belongs to other microbial phyla, thus the general results of our study likely hold in most other datasets. furthermore, the organisms present in our samples are sufficiently distant at the 16s rrna level to allow unambiguous clustering once a distance threshold was selected, i.e. the results we obtained reflect actual characteristics of 16s rrna data rather than artifacts of the analysis process. our simulated samples were composed of a small number of high-quality sequences with exceptional length - a relatively simple challenge compared to current datasets generated by pyrosequencing. 16s surveys employing pyrosequencing technologies generate considerably larger datasets, comprising reads of shorter lengths. the results we describe for long reads will only be exacerbated in the context of short-read data. further, these data present several new computational challenges for otu clustering including the need for faster alignment algorithms, the computation and storage of evolutionary distances, and the selection of distance thresholds appropriate for sequences with significantly less phylogenetic information than the sanger-based reads used in our study. the analysis of the new sequencing data require the development of robust methods for otu clustering, as well as rigorous validation of analysis methods. we finally demonstrated that a semi-supervised clustering approach can significantly improve analysis quality, highlighting the potential for semi-supervised clustering approaches to fill the gap between the two extreme classes of approaches commonly used in analyzing 16s rrna datasets: fully-supervised database searches that rely on expensive highly-curated datasets; and fully unsupervised otu clustering procedures that use ad hoc similarity thresholds in an attempt to match poorly-defined taxonomic labels. semi-supervised procedures such as vi-cut can generate accurate clusterings while only requiring high-quality labels for only a small subset of the sequences. finally, it is important to observe that clustering 16s rrna sequences into a set of otus is a valuable analysis tool even if the resulting otus do not correlate with pre-defined taxonomic entities. however, the ad hoc choice of analysis parameters, in particular the selection of different distance thresholds, complicates cross-study comparisons and even basic descriptions of diversity. as 16s rrna surveys are increasingly applied in a clinical setting, it is critical to accurately measure taxonomic diversity and identify individual species. our results highlight the need for standardizing 16s rrna analysis methods, or in the very least, reporting results obtained with multiple distance thresholds or clustering algorithms. the data used in this study have been deposited in the fames online database http://fames.jgi-psf.org - a repository for metagenomic analysis benchmarks. software used for semi-supervised clustering analysis is also available at http://www.cbcb.umd.edu/vicut/. abbreviations otu: operational taxonomic unit; rrna: ribosomal ribonucleic acid; vi: variation of information; msa: multiple sequence alignment; rdp: ribosomal database project; anova: analysis of variance; authors' contributions jrw and sn ran the experiments and analyzed the data. m-rg provided computational expertise. jw, sn, nn, m-rg, ck and mp wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 variation of information distances of high-level taxonomic clusterings from the annotated species clustering. to give the reader some intuition about the vi distance metric, we computed vi distances between the annotated species-level clustering and other clusterings based on phylum, class, order, family, and genus annotations. this file contains a table of these reference distances. click here for file acknowledgements we thank pawel gajer and jacques ravel for helpful conversations. this work was supported in part by the bill and melinda gates foundation, the nih, and the nsf. ck and sn partially supported by nsf grant to ck.
here we provide a critical assessment of current analysis methodologies that cluster sequences into operational taxonomic units and demonstrate that small changes in algorithm parameters can lead to significantly varying results. molecular studies of microbial diversity have provided many insights into the bacterial communities inhabiting the human body and the environment. to date, however, there exists significant variability in analysis methods employed in these studies. our results highlight the need for systematic and open evaluation of data analysis methodologies, especially as targeted 16s rrna diversity studies are increasingly relying on high-throughput sequencing technologies. all data and results from our study are available through the jgi fames website http://fames.jgi-psf.org/. a common first step in such studies is a survey of conserved marker genes to characterize the taxonomic composition and diversity of these communities.
 molecular studies of microbial diversity have provided many insights into the bacterial communities inhabiting the human body and the environment. a common first step in such studies is a survey of conserved marker genes to characterize the taxonomic composition and diversity of these communities. to date, however, there exists significant variability in analysis methods employed in these studies.  here we provide a critical assessment of current analysis methodologies that cluster sequences into operational taxonomic units and demonstrate that small changes in algorithm parameters can lead to significantly varying results. our analysis provides strong evidence that the species-level diversity estimates produced using common otu methodologies are inflated due to overly stringent parameter choices. we further describe an example of how semi-supervised clustering can produce otus that are more robust to changes in algorithm parameters.  our results highlight the need for systematic and open evaluation of data analysis methodologies, especially as targeted 16s rrna diversity studies are increasingly relying on high-throughput sequencing technologies. all data and results from our study are available through the jgi fames website http://fames.jgi-psf.org/. 