Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

Open Access Highly Accessed Research article

Flexible taxonomic assignment of ambiguous sequencing reads

José C Clemente12, Jesper Jansson3 and Gabriel Valiente4*

Author Affiliations

1 Center for Information Biology and DNA Databank of Japan, National Institute of Genetics, Yata 1111, Mishima 411-8540, Japan

2 Department of Chemistry and Biochemistry, University of Colorado, Boulder, CO, USA

3 Ochanomizu University, 2-1-1 Otsuka, Bunkyo-ku, Tokyo 112-8610, Japan

4 Algorithms, Bioinformatics, Complexity and Formal Methods Research Group, Technical University of Catalonia, E-08034 Barcelona, Spain

For all author emails, please log on.

BMC Bioinformatics 2011, 12:8  doi:10.1186/1471-2105-12-8

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/12/8


Received:6 July 2010
Accepted:7 January 2011
Published:7 January 2011

© 2011 Clemente et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

To characterize the diversity of bacterial populations in metagenomic studies, sequencing reads need to be accurately assigned to taxonomic units in a given reference taxonomy. Reads that cannot be reliably assigned to a unique leaf in the taxonomy (ambiguous reads) are typically assigned to the lowest common ancestor of the set of species that match it. This introduces a potentially severe error in the estimation of bacteria present in the sample due to false positives, since all species in the subtree rooted at the ancestor are implicitly assigned to the read even though many of them may not match it.

Results

We present a method that maps each read to a node in the taxonomy that minimizes a penalty score while balancing the relevance of precision and recall in the assignment through a parameter q. This mapping can be obtained in time linear in the number of matching sequences, because LCA queries to the reference taxonomy take constant time. When applied to six different metagenomic datasets, our algorithm produces different taxonomic distributions depending on whether coverage or precision is maximized. Including information on the quality of the reads reduces the number of unassigned reads but increases the number of ambiguous reads, stressing the relevance of our method. Finally, two measures of performance are described and results with a set of artificially generated datasets are discussed.

Conclusions

The assignment strategy of sequencing reads introduced in this paper is a versatile and a quick method to study bacterial communities. The bacterial composition of the analyzed samples can vary significantly depending on how ambiguous reads are assigned depending on the value of the q parameter. Validation of our results in an artificial dataset confirm that a combination of values of q produces the most accurate results.

Background

Microbes play a fundamental role as symbionts in the gut of mammals [1], and are strongly correlated with human health [2]. They also control some of the most important environmental processes, such as nitrogen fixation [3], and have been successfully used in the treatment of sewage [4] and to convert waste into usable fuels [5]. The importance of microbes is reflected in the large number of recent studies of bacterial communities in a variety of environments, including aquatic [6-11], soil [12-17], animal [18-23], and plant [24-29] habitats. The use of high-throughput sequencing technologies has greatly benefited the analysis of microbial populations [30], and different methodologies have been developed to characterize the diversity, richness, and similarity of bacterial communities [31-33]. Together with the introduction of new sequencing technologies, several challenges have emerged that need to be overcome to gain a better understanding of the diversity of bacteria that inhabit both the environment and ourselves [34].

Microbial communities are commonly characterized by mapping sequencing reads to a bacterial taxonomy based on the 16S rRNA gene. The effectiveness of this approach is not limited by the length of the read but by the choice of an adequate region of the 16S rRNA gene [35,36]. The structure of microbial communities has a high degree of variability, both in environmental [7] and gut samples [22]. In particular, human microbial communities differ greatly among individuals [23] and depending on the location of the body and time when the sample was taken [21]. Metabolic profiling has nevertheless shown that the functionality of communities is more conserved for particular environments [37], indicating that different species distributions can achieve a similar core functionality. Understanding the correlation between function and distribution of species therefore requires accurate measurements of both variables.

We have previously shown that a large proportion of reads in metagenomic studies can be assigned with equal significance to more than one species in the taxonomy [38]. The assignment of such ambiguous reads to the lowest common ancestor (LCA) of the matched species [39-42] introduces many false positives (leaves in the subtree rooted at the LCA that were not originally matched to the read), and thus we consider other possible nodes below the LCA to assign such reads. Implicitly, assignments at the LCA maximize the coverage but lower the accuracy, and we demonstrated that an assignment based on the F-measure, a combination of precision and recall, produces a significantly different distribution of taxonomic ranks to which reads are assigned.

In the absence of a reference taxonomy, the assignment of ambiguous reads is usually made by either mapping each read to the best BLAST hit in the reference sequences [43] or by using the reference sequences as a template for a multiple alignment of the reads, which defines pairwise similarities that are used to group the reads into clusters of related species [39,42,44]. In the absence of reference sequences, DNA composition can be used to group the reads into clusters of related species [45].

In this paper, we present a new method to assign ambiguous short reads to a node in the reference taxonomy minimizing a penalty score that generalizes our previous mapping based on the F-measure. Our algorithm is both fast and versatile, allowing a fine-grained assignment of reads closer to the LCA or the species depending on the value of a single parameter q. This parameter can be specified to have any value between 0 and 1, where setting q = 0 implies that each ambiguous read has an optimal assignment to a leaf, q = 1 always assigns to the LCA level, and q = 0.5 optimizes a combination of precision and recall. The use of this parameter provides the biologist with an intuitive tool to determine how to assign ambiguous reads, and results on six metagenomic datasets show the usefulness of our approach. The use of information on the quality of the sequencing reads results in a decrease in the number of unassigned reads but increases the number of ambiguous reads, thus making the assignment of such reads even more relevant. A method to validate our assignment algorithm is introduced and results for a set of artificial datasets are presented. Finally we discuss possible causes for the large proportion of ambiguous reads observed in these datasets.

Methods

Materials

We have initially studied six bacterial communities represented by 454 sequencing tags amplified from different 16S rRNA variable regions: a marine environment (V6 region) [7], the human gut (V3 and V6 regions) [20], the gut of lean and obese twins (V2 and V6 regions) [22], chicken gut (V6 region) [46], and rat gut (V4 region) [47]. The samples were between 50 and 329 bp in length and, for each of these communities, our algorithm assigned all the ambiguous sequencing reads at the best possible taxonomic rank, utilizing a reference bacterial taxonomy of 5,165 near-full-length type cultures of high quality [39], ranging in length between 1,202 and 1,780 bp, with a uniform scheme of seven taxonomic ranks (domain, phylum, class, order, family, genus, species). The taxonomy covers the whole spectrum of known bacteria, and the dominant phyla are Proteobacteria (1,925 species), Actinobacteria (1,285), Firmicutes (1,178), Bacteroidetes (355), and Tenericutes (160 species).

To test the effect of sequencing errors in the sequencing reads, we have also studied a dataset obtained from a bacterial community in the Priest Pot Lake [48], which included the quality scores for each read. Two taxonomies were used for this experiment, one based on the 5,165 high-quality type cultures only and the other one using all 322,864 16S rRNA sequences found in the Ribosomal Database Project [39].

Validation of our assignment algorithm was performed using artificial datasets generated with MetaSim [49], as follows. A first dataset was created importing 5,148 16S rRNA sequences from the taxonomy into MetaSim, and creating a dataset of short sequencing reads with the following parameters: 454 error model; 5,000,000 of reads per run; normal DNA clone size distribution with mean 500 and standard deviation 20; and expected read length 100, 150, and 200. Reads of less than 75% the expected length were discarded. A second dataset was created by extracting first 300 base pairs roughly corresponding to hypervariable regions V1 and V2 in each of the 5,148 rRNA sequences, and then importing those subsequences into MetaSim to generate a new set of short reads using the same set of parameters as before.

Taxonomic Assignment Method

In this section, we introduce a new method for accurately assigning a sequencing read Ri to a single node in a fixed reference taxonomy tree T with a leaf set L. We assume that each leaf in L has an associated known sequence. The input is a set R of sequencing reads and a positive integer k. For each Ri R, there is a subset Mi L of leaves whose sequences contain a substring with at most k mismatches to Ri; these leaves are referred to as hits or matches below. The goal is to output, for each Ri R with |Mi| ≥ 1, one node in T which represents all of Mi in a "good" way.

For any Ri R, if |Mi| = 1 then Ri can be trivially assigned to a unique leaf. However, if |Mi| ≥ 2 then Ri is called an ambiguous read and it is not immediately obvious how to optimally assign Ri to one node in T.

For this purpose, we let the user specify a value in the interval [0, 1] for a new parameter q. Intuitively, setting a low value of q means that ambiguous reads will be assigned to nodes near the leaves, while a high value of q means that assignments near the LCA level are preferred. Our approach for taxonomic assignment of reads is as follows.

1. Apply a read mapping tool, such as GEM [50] for instance, to R to compute the set of hits Mi for every Ri R.

2. Let the user specify a value in the interval [0, 1] for the parameter q.

3. For each Ri R:

(a) If |Mi| = 0 then output null.

(b) If |Mi| = 1 then output the leaf in Mi.

(c) Else, output all nodes j of T that have the smallest possible penalty score PSi,j with respect to q.

In the following sections we give the formal definition of the penalty score PSi,j and study some of its properties. Then, we consider how to implement Step 3c of our method efficiently, that is, how to compute the PSi,j values quickly.

Definition of the Penalty Score PSi,j

Let T be a (fixed) rooted tree with a leaf set L, let Mi L, and let q be a real number in the interval 0[1].

We need some additional notation. Let Ti be the subtree of T that is rooted at the LCA of Mi. For every node j in Ti, define:

Ti,j = The subtree of tree Ti rooted at node j.

• True positives: TPi,j = Leaves in Ti,j that belong to Mi.

• False positives: FPi,j = Leaves in Ti,j that do not belong to Mi.

• True negatives: TNi,j = Leaves in Ti\Ti,j that do not belong to Mi.

• False negatives: FNi,j = Leaves in Ti\Ti,j that belong to Mi.

See Figure 1 for an example. Note that for each node j in Ti, the leaves of Ti are partitioned into four disjoint subsets TPi,j, FPi,j, TNi,j, and FNi,j. The interpretation of this is that in case the node j is selected to be the representive for a read Ri whose hits are Mi, then each leaf in Ti will either be a true positive, a false positive, a true negative, or a false negative depending on whether or not it lies in the subtree rooted at j and if it is one of the hits.

thumbnailFigure 1. Sample Taxonomic Assignment of an Ambiguous Read. Assigning read Ri to the jth node of Ti partitions the leaves of Ti into true positives, false positives, true negatives, and false negatives. In this example, the hits are Mi = {s1, s5, s6, s7}. If we let j be the blue circled node, we obtain TPi,j = {s5, s6, s7}, FPi,j = {s8}, TNi,j = {s2, s3, s4}, FNi,j = {s1}. True hit H = s5 (defined in Section "Validation: Performance in ROC Space").

Finally, we define the penalty score with respect to q for every node j in Ti by the following formula:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M1">View MathML</a>

(1)

For every node j of T that does not belong to Ti, we define PSi,j = ∞. In case |TPi,j| = 0 then we also define PSi,j = ∞.

Different Values of q

Recall that our method for taxonomic assignment assigns each ambiguous read Ri to a node j that minimizes the value of PSi,j for the particular value of q specified by the user. We now study how varying the value of q affects the resulting taxonomic assignments.

First, it is easy to see that selecting q = 0 implies that a read Ri may have several different optimal assignments, but there always exists an optimal assignment of Ri to a leaf in Ti since then |FPi,j| (and hence, PSi,j) will be zero. On the other hand, q = 1 always assigns each Ri to the unique LCA of Mi because this gives |FNi,j| = 0 and PSi,j = 0. Thus, the special case q = 1 corresponds to the currently commonly used methods for assigning ambiguous reads [39-42].

Furthermore, we observe that PSi,j is a generalization of the mapping based on the F-measure that we previously introduced in [38]. If the precision of classifying read Ri into Tj is Pi,j = |TPi,j|/(|TPi,j| + |FPi,j|) and the recall is Ri,j = |TPi,j|/(|TPi,j| + |FNi,j|), then the F-measure is Fi,j = 2Pi,jRi,j/(Pi,j + Ri,j) = 2|TPi,j|/(|FNi,j| + |FPi,j| + 2|TPi,j|)). This yields:

Lemma 1. Any node m that minimizes the penalty score for q = 0.5 also maximizes the F-measure.

Proof.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M2">View MathML</a>

To summarize, the parameter q directly influences where in the reference taxonomy ambiguous reads will be assigned to. The user can adjust q to obtain representatives for ambiguous reads at the leaf level (q = 0), the LCA level (q = 1), or somewhere in-between (0 < q < 1). Interestingly, q = 0.5 is equivalent to maximizing the F-measure, which optimizes a combination of precision and recall. The distribution of taxonomic ranks resulting from setting various values of q in [0, 1] is further investigated for some real datasets in the "Results" Section.

Computation of the Penalty Scores PSi,j

Here, we focus on how to compute the penalty scores PSi,j in Step 3c of our method efficiently. For any tree T, let |T| denote the number of nodes in T. As before, let Mi be a set of hits in the reference taxonomy tree T and let Ti be the subtree of T that is rooted at the LCA of Mi. We may assume that |Mi| ≥ 2. Below, we first describe a simple method to obtain PSi,j for every node j in Ti in O(|Ti|) total time (Theorem 1). Then, we show that if O(|Tj|) time preprocessing of T is allowed, we can reduce the time complexity to obtain PSi,j for every node j in Ti in O(|Mi|) total time (Theorem 2). (The preprocessing of T will be done once, before any other computations in our taxonomic assignment method.) This modification gives a significant speedup in case R contains many reads that induce small sets of hits whose LCA are located at high taxonomic ranks in T.

For every node j in Ti, define Ti,j as the subtree of tree Ti rooted at j. The set of all leaves in Ti is denoted by Li, and Ni is the set of all leaves in Ti that do not belong to Mi (hence, Li = Mi Ni). Similarly, the set of all leaves in Ti,j that belong to Mi is denoted by Mi, j, Ni,j is the set of all leaves in Ti,j that do not belong to Mi, j, and Li,j = Mi,j Ni, j. Using this notation, we can write the previously defined TPi, j, FPi, j, TNi, j, and FNi,j as:

• True positives: TPi,j = Mi, j

• False positives: FPi,j = Ni, j

• True negatives: TNi,j = Ni \Ni, j

• False negatives: FNi,j = Mi \Mi, j

Next, we rewrite the formula for the penalty score in terms of |Mi|, |Mi, j|, and |Ni, j| as:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M3">View MathML</a>

(2)

Since Mi is given, the value of |Mi| is fixed. For any node j in Ti, the values of |Mi, j| and |Ni, j| may be expressed recursively as follows:

• If j is a leaf in Ti and j Mi: Then |Mi, j| = 1, |Ni, j| = 0, and |Li, j| = 1.

• If j is a leaf in Ti and j Mi: Then |Mi, j| = 0, |Ni, j| = 1, and |Li, j| = 1.

• If j is an internal node in Ti: Then |Mi,j| = ∑j'|Mi,j'| and |Li,j| = ∑j'|Li,j'|, where j' ranges over the children of j in Ti, and |Ni, j| = |Li, j| - |Mi, j|.

Hence, the values of PSi,j for all nodes in Ti can be obtained by two traversals of Ti: a (partial) bottom-up traversal [51,52] to identify the root of the subtree Ti of T (start at the leaves belonging to Mi and end when a node that is an ancestor of all leaves from Mi is reached; which can be determined by storing, for each node, the number of descendent leaves from Mi, because the first node in the bottom-up traversal that is an ancestor of all leaves from Mi has exactly |Mi| descendent leaves from Mi) followed by a top-down traversal to identify the subtree Ti of T while computing |Mi, j|, |Ni, j|, |Li, j|, and PSi,j for all nodes in Ti by applying the above relations. There are O(|Ti|) nodes in Ti, and so we have:

Theorem 1. We can find a node j in T that minimizes the value of PSi,j for any Mi L in O(|Ti|) time.

Next, we present an alternative method that improves the time complexity stated in Theorem 1 if preprocessing of the reference taxonomy tree T is allowed.

We start by explaining how to preprocess T. Fix an arbitrary left-to-right ordering of the nodes in T and perform a left-to-right postorder traversal of T in O(|T|) time while enumerating the nodes from 1 to |T| in accordance with the order in which they are first visited. Associate each node j with its number and, moreover, keep track of the smallest numbered leaf in the subtree rooted at j and denote it by m(j). Subsequently, for any two nodes j, j' in T, it holds that j is a proper ancestor of j' if and only if m(j) ≤ m(j') ≤ j' < j, and this condition can be checked in O(1) time. (The intervals [m(j), j] induced by nodes in T therefore exhibit a nested structure that will be utilized below.) Next, apply the O(|T|)-time preprocessing of [53,54] to T so that the LCA of any pair of specified leaves from T can be obtained in O(1) time, unless the height of T is bounded by a constant (usual taxonomical classifications as kingdom-phylum-class-order-family-genus-species have height 8, and the NCBI taxonomy [55] has a few more levels to account for finer distinctions), in which case the LCA of any pair of specified leaves from T can readily be obtained in O(1) time, without any preprocessing [52]. Lastly, do a O(|T|)-time bottom-up traversal of T to compute and store the number of leaves |Li| in the subtree rooted at each node i in T .

Now, suppose the preprocessing has been taken care of and we are given a set Mi L of hits. Any node j in Ti is called relevant if it is equal to a leaf in Mi or equal to the LCA of two or more leaves in Mi. We have:

Lemma 2. For each node j in Ti, there exists a relevant node j' such that PSi, j' PSi, j.

Proof. Suppose that j is a node in Ti that is not relevant. Let j' be the LCA of the leaves in Mi, j. Clearly, j' is relevant and, furthermore, |Mi, j| = |Mi, j'| while |Ni, j| > |Ni, j'| since Ti, j' is a subtree of Ti, j. It follows that <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M4">View MathML</a>.

Lemma 2 implies that PSi,j only needs to be computed for nodes in Ti that are relevant. Define the topological restriction of Ti to Mi, denoted by Ti || Mi, as the tree obtained by deleting from Ti all nodes that are not on a path from the root to a leaf in Mi along with their incident edges, and then contracting every edge between a node having just one child and its child. Then, the nodes of Ti || Mi are precisely the relevant nodes in Ti. Observe that Ti || Mi contains O(|Mi|) nodes.

To construct Ti || Mi for any specified Mi L in O(|Mi|) time, proceed as follows. Sort the leaves in Mi in O(|Mi|) time according to their left-to-right postordering numbers by a radix sort and write <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M5">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M6">View MathML</a>. For x ∈ {1, 2, ..., |Mi| , perform an O(1)-time LCA query on the pair (ℓx, ℓx+1) and let kx be the answer. The set <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M7">View MathML</a> then gives the set of nodes in Ti || Mi. To obtain the edges of Ti || Mi, first use O(|Mi|) time to perform a radix sort on the set of ordered pairs {(m(j), j): j U} in non-decreasing order for the first coordinate and decreasing order for the second coordinate so that in the resulting ordering, for any j, j'U, it holds that (m(j), j) < (m(j'), j') if and only if either (1) m(j) < m(j'), or (2) m(j) = m(j') and j > j'. Thus, whenever a node j is a proper ancestor of a node j', the pair (m(j), j) appears somewhere before (m(j'), j') in the sorted list. Then, it is straightforward to recover the edges of Ti || Mi in O(|Mi|) time by traversing the sorted list of pairs while using a stack to store all proper ancestors of the currently considered node (recall that it takes O(1) time to check the condition m(j) ≤ m(j') ≤ j' < j for any node j' in the list and any element j on the top of the stack).

Finally, the values of PSi,j for all relevant nodes can be obtained by a bottom-up traversal of Ti || Mi. There are O(|Mi|) relevant nodes, and so we can compute the values of PSi,j for all relevant nodes j in Step 3c according to formula (2) using O(|Mi|) time. To do this, note that if j is an internal node in Ti || Mi then |Mi, j| = ∑j' |Mi, j'|, where j' ranges over the children of j in Ti || Mi, and |Ni,j| = |Li,j| - |Mi,j|, where |Li,j| = |Lj| has been precomputed. In total:

Theorem 2. After O(|T|) time preprocessing, we can find a node j in T that minimizes the value of PSi,j for any Mi L in O(|Mi|) time.

Validation: Performance in ROC Space

Each read found in a metagenomic dataset must have originated from a unique original 16S rRNA sequence but, due to sequencing errors and incomplete taxonomic information, a significant percentage of the reads end up being assigned at higher taxonomic levels. Using an artificial metagenomic dataset, we can know a priori the original sequence and therefore measure how accurately our algorithm classifies ambiguous reads. We used MetaSim to generate an artificial set of sequencing reads R with a 454 error model, and where each read Ri is derived from a randomly selected full-length 16S rRNA gene sequence annotated in the taxonomy, denoted by Hi. Because of the errors introduced in Ri by the simulation, a search for the most similar full-length sequences to the read produces a set of hits, Mi, among which the true hit is usually found. The tree Ti rooted at the LCA of Mi can therefore include three kinds of leaves: the true hit Hi, not true but ambiguous hits Mi \{Hi}, and false hits Li \Mi = Ni.

When using our q-assignment schema, the sets Mi, Ni, Mi,j, and Ni,j are defined with respect to the set of plausible hits but without knowledge of the true hit Hi, as shown in Figure 1. Our objective here is to measure how different values of q perform in including Hi among the selected leaves. Assignments to the LCA will in most cases include Hi but the precision will be very poor if the size of Ti is large, while lower values of q can increase the precision at the risk of excluding Hi.

A common measure of performance for binary classifiers is the area under the ROC curve [56]. For a given read Ri and a particular value of q, let us define the true positive rate with respect to the true hit Hi as <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M8">View MathML</a> when Hi Ti and 0 otherwise (if Hi Ti both <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M9">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M10">View MathML</a> would be empty). Notice that when previously calculating the assignment we used the sets TP, FP, TN, and FN with respect to Mi, while here we calculate <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M11">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M10">View MathML</a> taking into account Hi only. In a similar way we define <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M12">View MathML</a>. However, <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M13">View MathML</a> can only be 1 (when Hi is in the subtree rooted at the node to which Ri was assigned) or 0 (when Hi is not in the subtree) and therefore, plotting TPR versus FPR would result in degenerated ROC curves.

We need to define <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M14">View MathML</a> as sets of leaves. That is, <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M15">View MathML</a> if Hi Ti, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M16">View MathML</a> otherwise. Then, <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M17">View MathML</a> if Hi Ti, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M18">View MathML</a> otherwise.

Let us define <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M19">View MathML</a> as the point in ROC space that represents read Ri. Points above the diagonal FPR = TPR have good predictive power, points below it are poor classification results, and points on the diagonal have no predictive power; that is, they are a random guess. The distance of pi to the diagonal is denoted by Di, and corresponds to the distance between pi and the point of intersection between the diagonal and the perpendicular that goes through pi. As shown in Figure 2Di equals <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M20">View MathML</a>, and we can define the goodness Gq of an assignment for a particular value of q as the sum of distances to the diagonal for all reads, where distances are negative if a point lies below the diagonal:

thumbnailFigure 2. Validation of Results in ROC Space. Distance in ROC space to the diagonal TPRH = FPRH. Points above the diagonal represent good predictions, the larger the distance the better.

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M21">View MathML</a>

(3)

The value of q that maximizes this sum will be the one with highest predictive power. This measure will be called the q metric of performance. Notice that the distance Di differs both in sign and in value depending on whether Hi is in the subtree Ti,j or not. We will define <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M22">View MathML</a> as the distance for read Ri when Hi is in the subtree and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M23">View MathML</a> when it is not:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M24">View MathML</a>

For a given value of q, our assignment produces two distinct subsets of reads: those that have been assigned to a node among whose descendents the true hit Hi can be found, and those that do not have Hi as a descendent. The goodness of an assignment for a value of q can then be rewritten as:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M25">View MathML</a>

Alternatively, and instead of assigning ambiguous reads based on the unique q value that maximizes Gq, we can look for the assignment that maximizes the expected distance E(Di) for each read, so our mapping would use a combination of different q values depending on the particular read being considered. Let us assume read Ri can be assigned to nodes n1, ..., nn for each of the q values q1, ..., qn. The probability of Hi being among the leaves of the subtree Ti,j rooted at nj is p = Mi,j/Mi, and that of not being included is 1 - p = (Mi -Mi,j) = Mi. The expected distance for a read Ri mapped to a node j for a given q is:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M26">View MathML</a>

We will call this measure of performance the expected distance metric. The value of q (and the corresponding node) that maximizes the expected distance for Ri is chosen as the most appropriate if the distance is positive, otherwise the assignment to the LCA is preferred. It can be easily seen that assignments to the LCA when the true hit Hi is among its leaves have expected distance 0, since Mi,j = Mi. Lemma 3 shows that the true distance for assignments to the LCA is always zero and, therefore, such assignments have no predictive power.

Lemma 3. Assignments to the LCA have no predictive power when the true hit Hi is among its leaves.

Proof. Since the true hit Hi is among the leaves of the subtree rooted at the LCA, the <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M13">View MathML</a> is 1. The <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M27">View MathML</a> is <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M28">View MathML</a> and, since assigning to the LCA selects all possible leaves, <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M29">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M30">View MathML</a>. Therefore, <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M13">View MathML</a> is equal to <a onClick="popup('http://www.biomedcentral.com/1471-2105/12/8/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/12/8/mathml/M27">View MathML</a> and the point representing such assignment in ROC space lies on the main diagonal, thus having no predictive power.

Results

A suffix array was constructed for the 5,165 reference sequences in the bacterial taxonomy, and each of the sequencing reads was mapped to these sequences using the GEM-do-index and GEM-mapper tools [50]. GEM-do-index constructs a suffix array from the set of full-length 16S rRNA sequences using the Burrows-Wheeler Transform [57]. Once the sequences have been efficiently indexed, GEM-mapper finds the closest sequences in the suffix array for each of the short sequencing reads in a metagenomic dataset.

Parameters were set to find all matching sequences with up to 2 mismatches, which is about 99% identity for reads of 200 bp. Figure 3 shows the distribution of sequencing reads mapped to more than one sequence in the taxonomy, and Figure 4 shows the distribution of hits per taxonomic rank. Most gut datasets show a distribution of hits that increases with rank up to the class level and then drops (rat, chicken, human, and twins V6 region), while the twins V2 region samples have a disproportionate number of hits at the domain level and the marine dataset does not seem to show a correlation between the number of hits and the taxonomic rank at which reads get assigned.

thumbnailFigure 3. Distribution of Sequencing Reads. Distribution of sequencing reads (number of hits and reads) ambiguously assigned with up to 2 mismatches to two or more of the 5,165 sequences in the reference bacterial taxonomy.

thumbnailFigure 4. Distribution of Hits per Taxonomic Rank. Distribution of the number of hits (species with up to k mismatches) in ambiguous reads per taxonomic rank.

We performed an alternative mapping using the sequencing reads as BLAST [58] queries against the reference 16S rRNA sequences, defining ambiguous hits as those that could be aligned to more than one species with the same e-value (e-value ≤ 0.001). Table 1 presents a summary of the number of reads that get an assignment with each tool.

Table 1. Performance of GEM and BLAST

Although BLAST hits a larger number of sequencing reads, those reads assigned to one or more species using GEM show a more significative BLAST average e-value than reads with no hits (data not shown). GEM provides a more conservative mapping, discarding those reads that get assigned with lower significance. Allowing for more mismatches with GEM results in a higher number of assigned reads with a higher percentage of ambiguous ones, at the cost of a lower average e-value of the assigned reads (data not shown). The results discussed in this paper, using GEM with up to 2 mismatches, should therefore be considered a conservative estimate.

The sequencing reads that matched two or more sequences in the reference bacterial taxonomy were assigned at the taxonomic rank that minimized the penalty score. The distribution of reads assigned at each taxonomic rank is shown in Figure 5 for values of the q parameter ranging from 0 to 1. These results show how ambiguous reads can be assigned at the desired taxonomic rank using different values of q: low values tend to produce a taxonomic assignment at the genus and species rank, while high values produce a taxonomic assignment at the class, order, and family ranks. The marine dataset seems to have a much higher level of ambiguity, as shown by the large proportion of ambiguous reads that get assigned at the order level for q = 1, and by the fact that lower values of q still assign many reads above the species level. The twins dataset is particularly interesting in that depending on the sequenced region, the reads are assigned quite differently. With region V2 there is a large percentage assigned above the genus level with q = 1, and this percentage is significant even for low values of q. The region V6, on the other hand, has most of its ambiguous reads assigned at the genus level when q = 1 (although with a notable subset of reads mapped very high, at the class level), but most of the reads get assigned at the species level quickly as we decrease the value of q.

thumbnailFigure 5. Distribution of Taxonomic Ranks in Metagenomic Datasets. Ambiguous reads assigned in the bacterial taxonomy at each taxonomic rank for q = 0, ..., 1. Color code: domain: purple; phylum: indigo; class: light blue; order: cyan; family: green; genus: yellow; species: red.

Sequencing Error Bias

Sequencing with 454 suffers mainly from indels in homopolymer runs [30], and such errors can have a significant effect on the final count of bacterial species in a metagenomic sample [48,59]. We analyzed the composition of a bacterial community using a sequencing dataset that included quality scores for each read [48]. Two suffix arrays were constructed: one from the 5,165 high-quality sequences, and one with all the 322,864 16S rRNA sequences found in the Ribosomal Database Project. Each read was then mapped to these taxonomies using the GEM-mapper tool with the parameters described above. The mapping was done with and without the quality scores of the reads, using the error model of 454 sequencing provided by the GEM-mapper tool when incorporating the scores. Table 2 shows the distribution of reads unassigned, assigned with one hit, and assigned with two or more hits (ambiguous reads), with and without quality information.

Table 2. Priest Pot Lake sample: Assignment

Among reads with two or more hits, the maximum number of matches was 82 species (both with and without quality information). Out of 26,458 reads without hits when not using quality information (plain FASTA files), 26,045 also have no hits when incorporating such data (FASTQ files), 158 now have one hit, and 255 have two or more hits. The 642 reads with one hit using FASTA are split into 611 also with a unique hit with FASTQ and 31 with two or more hits. Finally, all 1,261 reads with two or more hits with FASTA also have two or more hits with FASTQ. The distribution of taxonomic ranks with the 5,165 species taxonomy with and without quality scores can be seen in Table 3. Notice how with q = 1.0 the presence of a single incorrect species among the hits of a read results in its mapping to the LCA. The use of lower q values protects against such erroneous assignments when most of the hits belong to a particular taxonomic group, providing evidence of the read belonging to a taxonomic rank lower than the LCA of all hits.

Table 3. Priest Pot Lake sample: Taxonomic Distribution

Assignment Performance

To validate our assignment algorithm we generated six artificial metagenomic datasets using MetaSim, with read lengths 100, 150, and 200 bp for sequences extracted from the whole 16S rRNA sequence or from the V1-V2 hypervariable region only. Out of the original 5,000,000 reads, there were 195,580 (100 bp), 36,462 (150 bp), and 7,637 (200 bp) ambiguous reads when using the whole 16S sequence; and 123,486 (100 bp), 13,147 (150 bp), and 100 (200 bp) ambiguous reads when using the V1-V2 region.

Figure 6 shows the distribution of taxonomic ranks for each of the datasets. For 150 and 200 bp, the percentage of ambiguous reads is lower than that observed for experimental datasets, and assignment of these reads to the LCA produces a taxonomic distribution skewed towards lower ranks, with less reads mapped at high taxonomic ranks. Only the dataset generated using 100 bp short reads extracted from complete 16S rRNA gene sequences shows a similar level of ambiguity to that of experimental datasets.

thumbnailFigure 6. Distribution of Taxonomic Ranks in Simulated Datasets. Simulated ambiguous reads assigned in the bacterial taxonomy at each taxonomic rank for q = 0, ..., 1. Datasets were constructed from whole 16S rRNA sequence and from the V1-V2 hypervariable region. Color code: domain: purple; phylum: indigo; class: light blue; order: cyan; family: green; genus: yellow; species: red.

As observed in Table 4 the q metric performs better for higher values of this parameter. As the simulated read length decreases, there is a higher proportion of ambiguous reads, more values of q result in a positive Gq and the best sum of distances improves: 0.028 (200 bp), 0.048 (150 bp), and 0.096 (100 bp) when using the full-length sequence, and 0.05 (200 bp), 0.015 (150 bp), and 0.036 (100 bp) when using only the V1-V2 region. Assignments with q = 0 do not perform particularly well (-0.42, -0.32, and -0.24 for 100, 150, and 250 bp extracted from the full-length 16S, and -0.23, -0.16, and -0.08 for 100, 150, and 250 bp obtained from the V1-V2 region), indicating that in most cases the true hit H would not be in the chosen subtree if we choose a single leaf out of all possible hits. This is more so as the read length decreases and ambiguous reads are mapped at higher taxonomic ranks, resulting in a larger number of possible hits (Mi) and a lower probability of choosing the true hit (Mi,j/Mi).

Table 4. Validation of Results: Gq

The expected distance metric E(Di) maps half of the reads to the LCA, and the rest get evenly distributed among all q values, as seen in Table 5. As with the q metric, more ambiguous datasets produce better results, and the sum of distances gradually increases as well: 0.073 (200 bp), 0.086 (150 bp), and 0.125 (100 bp) when using the full-length sequence, and 0.031 (200 bp), 0.058 (150 bp), and 0.077 (100 bp) when using the V1-V2 region. The best sum of distances is always larger with E(Di) than with Gq, indicating that a combination of values of q is a better predictor than assigning all reads using a unique value.

Table 5. Validation of Results: Assignments.

Table 6 shows precision and recall values for the q metric in the six artificial datasets. Precision decreases and recall increases with higher values of this parameter, as expected. Moreover, reads extracted from full-length sequences tend to have more matches than reads extracted from the V1-V2 region and thus, precision values are higher for reads extracted from the V1-V2 region than for reads extracted from full-length sequences (the latter contain more false positives), and recall values tend to be higher for reads extracted from full-length sequences than for reads extracted from the V1-V2 region (the latter contain more false negatives).

Table 6. Validation of Results: Precision and Recall.

Discussion

Comparison of results between GEM and BLAST show that, although BLAST can map more reads, there is still a large number of reads that either cannot be assigned to any species in the taxonomy, or can be assigned to more than one species. GEM provides a more conservative assignment, with assigned reads having a more significant e-value in BLAST. The combined speed of GEM in assigning species to each read and of our algorithm in mapping ambiguous reads to the taxonomic rank minimizing the penalty score provides a useful tool to quickly test hypotheses about microbial communities.

Ambiguous reads are assigned to different taxonomic ranks depending on the value of q, as shown in Figure 5 Figure 6 and Table 3. As the value of q increases, more reads get assigned at higher taxonomic ranks, with clearly different distributions between the extreme values q = 0 and q = 1. Ambiguous and unassigned reads could either belong to species not present in the taxonomy or be artifacts due to errors in the experimental process. Bacterial taxonomies are biased towards cultivable species, but the human gut and oceanic environments are known to harbor many rare, uncultivable bacteria [7,20]. Even if a large number of reads come from unknown species, most of them have a small number of hits (as seen in Figure 3) and would only introduce a moderate amount of error. Reads with a large number of hits, such as some of the reads coming from the twins V2 region dataset, might be chimeric sequences product of the PCR amplification step [60], or reads belonging to species from yet to be identified taxonomic groups. Reads not uniquely mapped can also be caused by sequencing errors, most commonly homopolymer indels when using 454 sequencing [61]. We would expect to observe differences in the distributions of homopolymers between the sets of ambiguous and unambiguous reads, but we could not find a significantly higher number of homopolymers for any of these sets across all our datasets (data not shown). Analysis of the distributions of homopolymer lengths versus the number of homopolymers per read was also inconclusive, and we could not clearly differentiate the distributions in ambiguous and unambiguous reads. A BLAST search using reads with no hits as both query and database showed that the vast majority can be aligned with high significance (often perfect matches) to other no-hit reads, indicating that either the same type of sequencing error occurred frequently or that these reads cannot be mapped because the corresponding 16S rRNA sequence is not in the taxonomy. On the other hand, incorporating read quality information reduces the number of unassigned reads but increases the number of ambiguous reads more significantly than that of uniquely assigned reads (see Table 3), stressing the importance of an assignment of ambiguous reads that minimizes bias in the estimation of diversity.

Artificial datasets produced through simulations contained a lower percentage of ambiguous reads when compared to experimental datasets for similar read lengths, and most ambiguous reads were assigned at lower taxonomic ranks. The simulations using the V1-V2 hypervariable region to extract short reads were expected to mimic experimental conditions more closely, but they showed even less ambiguity than simulations using the full-length 16S rRNA sequence to obtain the short reads. The results described in Section "Assignment Performance", therefore, represent a low estimate of performance, and datasets with a higher proportion of ambiguous reads would further benefit from our assignment algorithm. Tables 4 and 5, in fact, show how simulations with a higher proportion of ambiguous reads benefit more from our taxonomic mapping algorithm. Although the sequencing error models utilized in the simulations probably differ slightly from the actual errors, we believe the increased ambiguity observed in the experimental datasets is due to a combination of several factors: sequencing errors, PCR artifacts, and incomplete taxonomic information for some of the species present in the samples. Selecting a unique q value to assign ambiguous reads based on these simulations might therefore produce biased results and, until more realistic simulations can be performed, we suggest the use of the estimated distance metric instead, which does not require an estimation of an optimal value for q. It should be noticed that, both in the q metric and in the expected distance metric, the distances are relevant to determine whether the assignment to the LCA can be outperformed, but the true measure of significance is given by the observed differences in the distribution of taxonomic ranks, as seen in Figure 5 and Figure 6.

Conclusions

In this paper, we have introduced a new method for the taxonomic assignment of ambiguous sequencing reads based on a generalization of the F-measure that minimizes a penalty score combining the precision and recall of the mapping. There is, to the best of our knowledge, no other taxonomic assignment method concerning precision and recall, apart of the assignment to the LCA. By using a suffix array representation of the sequences in the leaves of the taxonomy and preprocessing the taxonomy for fast search of relevant nodes, our assignment algorithm can work in time linear in the number of sequences matching a read. Our algorithm can analyze large metagenomic datasets even on a small PC. For instance, on a MacBook Pro with 8GB of memory, the analysis of the Priest Pot Lake dataset takes approximately 30 minutes for GEM to analyze (up to 7 mismatches), and another 30 minutes to assign ambiguous reads. The use of a single parameter to control whether precision or recall should be prioritized results in assignments with clearly different distributions of taxonomic ranks. The assignment strategy of sequencing reads introduced in this paper is therefore both a versatile and a quick method to study bacterial communities.

The study of six different datasets of environmental and gut samples shows that the composition and diversity of bacterial species observed in a sample can vary significantly depending on whether ambiguous or unambiguous reads are used, and on the particular value of the q parameter. Results with a dataset where read quality information is provided shows that the number of ambiguous reads increases when such information is used, making our algorithm more relevant. Validation of the assignment schema in an artificial dataset shows that a combination of different q values produces the most accurate results. The fact that a unique set of sequencing reads can produce very different distributions depending on how the large number of ambiguous reads are assigned has deep implications for metagenomic studies in general, and in particular for those trying to correlate bacterial composition with disease states. A more accurate characterization of these reads can therefore provide a better understanding of the microbial diversity around and within us.

Availability

The software and data sets are available under the GNU GPL at the supplementary material web page http://www.lsi.upc.edu/~valiente/tango/. webcite

Authors' contributions

All authors conceived the method, prepared the manuscript, contributed to the discussion, and have approved the final manuscript. GV implemented the software. JC also implemented part of the software.

Acknowledgements

JJ was supported by the Special Coordination Funds for Promoting Science and Technology (Japan). The authors would like to thank Daniel L. Hartl, Chaysavanh Manichanh, Paolo Ribeca, Roderic Guigó, and everybody at the Center for Informational Biology (Ochanomizu University) and the Laboratory for DNA Data Analysis (National Institute of Genetics) for helpful discussions.

References

  1. Ley RE, Hamady M, Lozupone C, Turnbaugh PJ, Ramey RR, Bircher JS, Schlegel ML, Tucker TA, Schrenzel MD, Knight R, Gordon JI: Evolution of mammals and their gut microbes.

    Science 2008, 320(5883):1647-1651. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. Dethlefsen L, McFall-Ngai M, Relman DA: An ecological and evolutionary perspective on human-microbe mutualism and disease.

    Nature 2007, 449(7164):811-818. PubMed Abstract | Publisher Full Text OpenURL

  3. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molecular Biology of the Cell. 5th edition. New York, USA: Garland Science; 2008. OpenURL

  4. Gray NF: Biology of Wastewater Treatment. 2nd edition. London, UK: Imperial College Press; 2004. OpenURL

  5. Jeffries T, Jin YS: Metabolic engineering for improved fermentation of pentoses by yeasts.

    Appl Microbiol Biotechnol 2004, 63(5):495-509. PubMed Abstract | Publisher Full Text OpenURL

  6. Venter J, Remington K, Heidelberg J, Halpern A, Rusch D, Eisen J, Wu D, Paulsen I, Nelson K, Nelson W, Fouts D, Levy S, Knap A, Lomas M, Nealson K, White O, Peterson J, Hoffman J, Parsons R, Baden-Tillson H, Pfannkoch C, Rogers Y, Smith H: Environmental genome shotgun sequencing of the Sargasso Sea.

    Science 2004, 304(5667):66-74. PubMed Abstract | Publisher Full Text OpenURL

  7. Sogin ML, Morrison HG, Huber JA, Welch DM, Huse SM, Neal PR, Arrieta JM, Herndl GJ: Microbial diversity in the deep sea and the underexplored "rare biosphere".

    Proc Natl Acad Sci USA 2006, 103(32):12115-12120. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Humbert JF, Dorigo U, Cecchi P, Berre BL, Debroas D, Bouvy M: Comparison of the structure and composition of bacterial communities from temperate and tropical freshwater ecosystems.

    Environ Microbiol 2009, 11(9):2339-2350. PubMed Abstract | Publisher Full Text OpenURL

  9. Pašić L, Rodriguez-Mueller B, Martin-Cuadrado AB, Mira A, Rohwer F, Rodriguez-Valera F: Metagenomic islands of hyperhalophiles: the case of Salinibacter ruber.

    BMC Genomics 2009, 10:570. PubMed Abstract | PubMed Central Full Text OpenURL

  10. Kirchman DL, Cottrell MT, Lovejoy C: The structure of bacterial communities in the western Arctic Ocean as revealed by pyrosequencing of 16S rRNA genes.

    Environ Microbiol 2010, 12(5):1132-1143. PubMed Abstract | Publisher Full Text OpenURL

  11. Revetta RP, Pemberton A, Lamendella R, Iker B, Domingo JWS: Identification of bacterial populations in drinking water using 16S rRNA-based sequence analyses.

    Water Res 2010, 44(5):1353-1360. PubMed Abstract | Publisher Full Text OpenURL

  12. Martín HG, Ivanova N, Kunin V, Warnecke F, Barry KW, McHardy AC, Yeates C, He S, Salamov AA, Szeto E, Dalin E, Putnam NH, Shapiro HJ, Pangilinan JL, Rigoutsos I, Kyrpides NC, Blackall LL, McMahon KD, Hugenholtz P: Metagenomic analysis of two enhanced biological phosphorus removal (EBPR) sludge communities.

    Nat Biotechnol 2006, 24(10):1263-1269. PubMed Abstract | Publisher Full Text OpenURL

  13. Schloss PD, Handelsman J: Toward a Census of Bacteria in Soil.

    PLoS Comput Biol 2006, 2(7):e92. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Tarlera S, Jangid K, Ivester AH, Whitman WB, Williams MA: Microbial community succession and bacterial diversity in soils during 77 000 years of ecosystem development.

    FEMS Microbiol Ecol 2008, 64:129-140. PubMed Abstract | Publisher Full Text OpenURL

  15. Fierer N, Carney KM, Horner-Devine MC, Megonigal JP: The Biogeography of Ammonia-Oxidizing Bacterial Communities in Soil.

    Microb Ecol 2008, 58(2):435-445. Publisher Full Text OpenURL

  16. Hao DC, Ge GB, Yang L: Bacterial diversity of Taxus rhizosphere: culture-independent and culture-dependent approaches.

    FEMS Microbiol Lett 2008, 284(2):204-212. PubMed Abstract | Publisher Full Text OpenURL

  17. Chu H, Fierer N, Lauber CL, Caporaso JG, Knight R, Grogan P: Soil bacterial diversity in the Arctic is not fundamentally different from that found in other biomes.

    Environ Microbiol 2010, 12(11):2998-3006. PubMed Abstract | Publisher Full Text OpenURL

  18. Eckburg PB, Bik EM, Bernstein CN, Purdom E, Dethlefsen L, Sargent M, Gill SR, Nelson KE, Relman DA: Diversity of the human intestinal microbial flora.

    Science 2005, 308(5728):1635-1638. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. Aas JA, Griffen AL, Dardis SR, Lee AM, Olsen I, Dewhirst FE, Leys EJ, Paster BJ: Bacteria of dental caries in primary and permanent teeth in children and young adults.

    J Clin Microbiol 2008, 46(4):1407-1417. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Dethlefsen L, Huse SM, Sogin ML, Relman DA: The Pervasive Effects of an Antibiotic on the Human Gut Microbiota, as Revealed by Deep 16S rRNA Sequencing.

    PLoS Biol 2008, 6(11):e280. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Costello EK, Lauber CL, Hamady M, Fierer N, Gordon JI, Knight R: Bacterial Community Variation in Human Body Habitats Across Space and Time.

    Science 2009, 326(5960):1694-1697. PubMed Abstract | Publisher Full Text OpenURL

  22. Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley RE, Sogin ML, Jones WJ, Roe BA, Affourtit JP, Egholm M, Henrissat B, Heath AC, Knight R, Gordon JI: A core gut microbiome in obese and lean twins.

    Nature 2009, 457(7228):480-484. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  23. Fierer N, Lauber CL, Zhou N, McDonald D, Costello EK, Knight R: Forensic Identification using skin bacterial communities.

    Proc Natl Acad Sci USA 2010, 107(14):6477-6481. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Lambais MR, Crowley DE, Cury JC, Büll RC, Rodrigues RR: Bacterial Diversity in Tree Canopies of the Atlantic Forest.

    Science 2006, 312(5782):1917. PubMed Abstract | Publisher Full Text OpenURL

  25. Leveau JHJ: The magic and menace of metagenomics: prospects for the study of plant growth-promoting rhizobacteria.

    Eur J Plant Pathol 2007, 119(3):279-300. Publisher Full Text OpenURL

  26. Sun L, Qiu F, Zhang X, Dai X, Dong X, Song W: Endophytic Bacterial Diversity in Rice (Oryza sativa L.) Roots Estimated by 16S rDNA Sequence Analysis.

    Microb Ecol 2008, 55(3):415-424. PubMed Abstract | Publisher Full Text OpenURL

  27. Wang HX, Geng ZL, Zeng Y, Shen YM: Enriching plant microbiota for a metagenomic library construction.

    Environ Microbiol 2010, 10(10):2684-2691. Publisher Full Text OpenURL

  28. Adams IP, Glover RH, Monger WA, Mumford R, Jackeviciene E, Navalinskiene M, Samuitiene M, Boonham N: Next-generation sequencing and metagenomic analysis: a universal diagnostic tool in plant virology.

    Mol Plant Pathol 2009, 10(10):2684-2691. OpenURL

  29. Redford AJ, Bowers RM, Knight R, Linhart Y, Fierer N: The ecology of the phyllosphere: geographic and phylogenetic variability in the distribution of bacteria on tree leaves.

    Environ Microbiol 2010, 12(11):2885-2893. PubMed Abstract | Publisher Full Text OpenURL

  30. Shendure J, Ji H: Next-generation DNA sequencing.

    Nat Biotechnol 2008, 26(10):1135-1145. PubMed Abstract | Publisher Full Text OpenURL

  31. Lozupone C, Knight R: UniFrac: A new phylogenetic method for comparing microbial communities.

    Appl Environ Microbiol 2005, 71(12):8228-8235. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  32. Schloss PD, Handelsman J: A statistical toolbox for metagenomics: Assessing functional diversity in microbial communities.

    BMC Bioinformatics 2008, 9:34. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  33. Singleton D, Furlong M, Rathbun S, Whitman W: Quantitative comparisons of 16S rRNA gene sequence libraries from environmental samples.

    Appl Environ Microbiol 2001, 67(9):4374-4376. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. Hamady M, Knight R: Microbial community profiling for human microbiome projects: Tools, techniques, and challenges.

    Genome Res 2009, 19(7):1141-1152. PubMed Abstract | Publisher Full Text OpenURL

  35. Liu Z, Lozupone C, Hamady M, Bushman FD, Knight R: Short pyrosequencing reads suffice for accurate microbial community analysis.

    Nucleic Acids Res 2007, 35(18):e120. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  36. Manichanh C, Chapple CE, Frangeul L, Gloux K, Guigó R, Dore J: A comparison of random sequence reads versus 16S rDNA sequences for estimating the biodiversity of a metagenomic library.

    Nucleic Acids Res 2008, 36(16):5180-5188. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  37. Dinsdale EA, Edwards RA, Hall D, Angly F, Breitbart M, Brulc JM, Furlan M, Desnues C, Haynes M, Li L, McDaniel L, Moran MA, Nelson KE, Nilsson C, Olson R, Paul J, Brito BR, Ruan Y, Swan BK, Stevens R, Valentine DL, Thurber RV, Wegley L, White BA, Rohwe F: Functional metagenomics profiling of nine biomes.

    Nature 2008, 452(7187):629-632. PubMed Abstract | Publisher Full Text OpenURL

  38. Clemente JC, Jansson J, Valiente G: Accurate Taxonomic Assignment of Short Pyrosequencing Reads. In Proc. 15th Pacific Symposium on Biocomputing. Volume 15. World Scientific; 2010::3-9. OpenURL

  39. Cole JR, Wang Q, Cardenas E, Fish J, Chai B, Farris RJ, Kulam-Syed-Mohideen AS, McGarrell DM, Marsh T, Garrity GM, Tiedje JM: The Ribosomal Database Project: Improved Alignments and New Tools for rRNA Analysis.

    Nucleic Acids Res 2009, 37(D):141-145. Publisher Full Text OpenURL

  40. Huson DH, Auch AF, Qi J, Schuster SC: MEGAN Analysis of Metagenomic Data.

    Genome Res 2007, 17(3):377-386. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  41. Liu Z, DeSantis TZ, Andersen GL, Knight R: Accurate Taxonomy Assignments from 16S rRNA Sequences produced by Highly Parallel Pyrosequencers.

    Nucleic Acids Res 2008, 36(18):e120. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  42. Wang Q, Garrity GM, Tiedje JM, Cole JR: Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy.

    Appl Environ Microbiol 2007, 73(16):5261-5267. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  43. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, Fierer N, Peña AG, Goodrich JK, Gordon JI, Huttley GA, Kelley ST, Knights D, Koenig JE, Ley RE, Lozupone CA, McDonald D, Muegge BD, Pirrung M, Reeder J, Sevinsky JR, Turnbaugh PJ, Walters WA, Widmann J, Yatsunenko T, Zaneveld J, Knight R: QIIME allows analysis of high-throughput community sequencing data.

    Nat Methods 2010, 7(5):335-6. PubMed Abstract | Publisher Full Text OpenURL

  44. Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, Lesniewski RA, Oakley BB, Parks DH, Robinson CJ, Sahl JW, Stres B, Thallinger GG, Horn DJV, Weber CF: Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities.

    Appl Environ Microbiol 2009, 75(23):7537-41. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  45. Teeling H, Waldmann J, Lombardot T, Bauer M, Glöckner FO: TETRA: a web-service and a stand-alone program for the analysis and comparison of tetranucleotide usage patterns in DNA sequences.

    BMC Bioinformatics 2004, 5:163. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  46. VAMPS: [http://vamps.mbl.edu/] webcite

    Visualization and Analysis of Microbial Population Structure project. 2009.

    [AGT CKN Bv6--Chicken intestinal microbiota]

    OpenURL

  47. Manichanh C, Reeder J, Gibert P, Varela E, Llopis M, Antolin M, Guigó R, Knight R, Guarner F: Reshaping the gut microbiome with bacterial transplantation and antibiotic intake.

    Genome Res 2010, 20(10):1411-1419. PubMed Abstract | Publisher Full Text OpenURL

  48. Quince C, Lanzén A, Curtis TP, Davenport RJ, Hall N, Head IM, Read LF, Sloan WT: Accurate determination of microbial diversity from 454 pyrosequencing data.

    Nat Methods 2009, 6(9):639-641. PubMed Abstract | Publisher Full Text OpenURL

  49. Richter DC, Ott F, Auch AF, Schmid R, Huson DH: MetaSim - a sequencing simulator for genomics and metagenomics.

    PLoS One 2008, 3(10):e3373. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  50. Ribeca P: [http://gemlibrary.sourceforge.net/] webcite

    GEM--GEnomic Multi-tool. 2009. OpenURL

  51. Valiente G: Algorithms on Trees and Graphs. Berlin, Heidelberg: Springer; 2002. OpenURL

  52. Valiente G: Combinatorial Pattern Matching Algorithms in Computational Biology using Perl and R. Boca Raton, London, New York: Taylor & Francis/CRC Press; 2009. OpenURL

  53. Bender M, Farach-Colton M, Pemmasani G, Skiena S, Sumazin P: Lowest common ancestors in trees and directed acyclic graphs.

    J Algorithms 2005, 57(2):75-94. Publisher Full Text OpenURL

  54. Harel D, Tarjan RE: Fast algorithms for finding nearest common ancestors.

    SIAM J Comput 1984, 13(2):338-355. Publisher Full Text OpenURL

  55. Wheeler DL, Chappey C, Lash AE, Leipe DD, Madden TL, Schuler GD, Tatusova TA, Rapp BA: Database resources of the National Center for Biotechnology Information. [http://www.ncbi.nlm.nih.gov/Taxonomy/] webcite

    Nucleic Acids Res 2000, 28:10-14. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  56. Fawcett T: An introduction to ROC analysis.

    Pattern Recogn Lett 2006, 27(8):861-874. Publisher Full Text OpenURL

  57. Burrows M, Wheeler D: A block sorting lossless data compression algorithm.

    Tech. Rep. 124, Digital Equipment Corporation 1994. OpenURL

  58. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic Local Alignment Search Tool.

    J Mol Biol 1990, 215(3):403-410. PubMed Abstract OpenURL

  59. Kunin V, Engelbrektson A, Ochman H, Hugenholtz P: Wrinkes in the rare biosphere: pyrosequencing errors can lead to artificial inflation of diversity estimates.

    Environ Microbiol 2001, 12:118-123. Publisher Full Text OpenURL

  60. Pääbo S, Irwin DM, Wilson AC: DNA damage promotes jumping between templates during enzymatic amplification.

    J Biol Chem 1990, 265(8):4718-4721. PubMed Abstract | Publisher Full Text OpenURL

  61. Huse SM, Huber JA, Morrison hilaryG, Sogin ML, Welch DM: Accuracy and quality of massively parallel DNA pyrosequencing.

    Genome Biol 2007, 8(7):R143. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL