Flexible taxonomic assignment of ambiguous sequencing reads

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][7][8][9][10][11], soil [12][13][14][15][16][17], animal [18][19][20][21][22][23], and plant [24][25][26][27][28][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][32][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][40][41][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.

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 R i 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 R i R, there is a subset M i ⊆ L of leaves whose sequences contain a substring with at most k mismatches to R i ; these leaves are referred to as hits or matches below. The goal is to output, for each R i R with |M i | ≥ 1, one node in T which represents all of M i in a "good" way.
For any R i R, if |M i | = 1 then R i can be trivially assigned to a unique leaf. However, if |M i | ≥ 2 then R i 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 M i for every R i R. 2. Let the user specify a value in the interval [0, 1] for the parameter q.

For each
(c) Else, output all nodes j of T that have the smallest possible penalty score PS i,j with respect to q.
In the following sections we give the formal definition of the penalty score PS i,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 PS i,j values quickly.
Definition of the Penalty Score PS i,j Let T be a (fixed) rooted tree with a leaf set L, let M i ⊆ 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 M i . For every node j in T i , define: • T i,j = The subtree of tree T i rooted at node j. • True positives: TP i,j = Leaves in T i,j that belong to M i . • False positives: FP i,j = Leaves in T i,j that do not belong to M i . • True negatives: TN i,j = Leaves in T i \T i,j that do not belong to M i . • False negatives: FN i,j = Leaves in T i \T i,j that belong to M i . See Figure 1 for an example. Note that for each node j in T i , the leaves of T i are partitioned into four disjoint subsets TP i,j , FP i,j , TN i,j , and FN i,j . The interpretation of this is that in case the node j is selected to be the representive for a read R i whose hits are M i , then each leaf in T i 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.
Finally, we define the penalty score with respect to q for every node j in T i by the following formula: For every node j of T that does not belong to T i , we define PS i,j = ∞. In case |TP i,j | = 0 then we also define PS i,j = ∞.

Different Values of q
Recall that our method for taxonomic assignment assigns each ambiguous read R i to a node j that minimizes the value of PS i,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 R i may have several different optimal assignments, but there always exists an optimal assignment of R i to a leaf in T i since then |FP i,j | (and hence, PS i,j ) will be zero. On the other hand, q = 1 always assigns each R i to the unique LCA of M i because this gives |FN i,j | = 0 and PS i,j = 0. Thus, the special case q = 1 corresponds to the currently commonly used methods for assigning ambiguous reads [39][40][41][42].
Furthermore, we observe that PS i,j is a generalization of the mapping based on the F-measure that we previously introduced in [38]. If the precision of classifying read R i into T j is P i,j = |TP i,j |/(|TP i,j | + |FP i,j |) and the recall is R i,j = |TP i,j |/(|TP i,j | + |FN i,j |), then the F-measure is F i,j = 2P i,j R i,j /(P i,j + R i,j ) = 2|TP i,j |/(|FN i,j | + |FP i,j | + 2|TP i,j |)). This yields:  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 PS i,j
Here, we focus on how to compute the penalty scores PS i,j in Step 3c of our method efficiently. For any tree T, let |T| denote the number of nodes in T. As before, let M i be a set of hits in the reference taxonomy tree T and let T i be the subtree of T that is rooted at the LCA of M i . We may assume that |M i | ≥ 2. Below, we first describe a simple method to obtain PS i,j for every node j in T i in O(|T i |) total time (Theorem 1). Then, we show that if O(|T j |) time preprocessing of T is allowed, we can reduce the time complexity to obtain PS i,j for every node j in T i in O(|M i |) 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 T i , define T i,j as the subtree of tree T i rooted at j. The set of all leaves in T i is denoted by L i , and N i is the set of all leaves in T i that do not belong to M i (hence, L i = M i ∪ N i ). Similarly, the set of all leaves in T i,j that belong to M i is denoted by M i, j , N i,j is the set of all leaves in T i,j that do not belong to M i, j , and L i,j = M i,j ∪ N i, j . Using this notation, we can write the previously defined TP i, j , FP i, j , TN i, j , and FN i,j as: Next, we rewrite the formula for the penalty score in terms of |M i |, |M i, j |, and |N i, j | as: Since M i is given, the value of |M i | is fixed. For any node j in T i , the values of |M i, j | and |N i, j | may be expressed recursively as follows: Hence, the values of PS i,j for all nodes in T i can be obtained by two traversals of T i : a (partial) bottom-up traversal [51,52] to identify the root of the subtree T i of T (start at the leaves belonging to M i and end when a node that is an ancestor of all leaves from M i is reached; which can be determined by storing, for each node, the number of descendent leaves from M i , because the first node in the bottom-up traversal that is an ancestor of all leaves from M i has exactly |M i | descendent leaves from M i ) followed by a top-down traversal to identify the subtree T i of T while computing |M i, j |, |N i, j |, |L i, j |, and PS i,j for all nodes in T i by applying the above relations. There are O(|T i |) nodes in T i , and so we have: Theorem 1. We can find a node j in T that minimizes the value of PS i,j for any 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 |L i | 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 M i ⊆ L of hits. Any node j in T i is called relevant if it is equal to a leaf in M i or equal to the LCA of two or more leaves in M i . We have: Proof. Suppose that j is a node in T i that is not relevant. Let j' be the LCA of the leaves in M i, j . Clearly, j' is relevant and, furthermore, Lemma 2 implies that PS i,j only needs to be computed for nodes in T i that are relevant. Define the topological restriction of T i to M i , denoted by T i || M i , as the tree obtained by deleting from T i all nodes that are not on a path from the root to a leaf in M i along with their incident edges, and then contracting every edge between a node having just one child and its child. Then, the nodes of T i || M i are precisely the relevant nodes in T i .
where |L i ,j| = |L j | 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 PS i,j for any M i ⊆ L in O(|M i |) 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 R i is derived from a randomly selected full-length 16S rRNA gene sequence annotated in the taxonomy, denoted by H i . Because of the errors introduced in R i by the simulation, a search for the most similar full-length sequences to the read produces a set of hits, M i , among which the true hit is usually found. The tree T i rooted at the LCA of M i can therefore include three kinds of leaves: the true hit H i , not true but ambiguous hits M i \{H i }, and false hits When using our q-assignment schema, the sets M i , N i , M i,j , and N i,j are defined with respect to the set of plausible hits but without knowledge of the true hit H i , as shown in Figure 1. Our objective here is to measure how different values of q perform in including H i among the selected leaves. Assignments to the LCA will in most cases include H i but the precision will be very poor if the size of T i is large, while lower values of q can increase the precision at the risk of excluding H i . A common measure of performance for binary classifiers is the area under the ROC curve [56]. For a given read R i and a particular value of q, let us define the true positive rate with respect to the true hit H i as  Figure 2D i equals 2 1 2 ( )/ − FPR H i , and we can define the goodness G q 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: 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 D i differs both in sign and in value depending on whether Hi is in the subtree T i,j or not. We will define D H T i j ∈ , as the distance for read R i when Hi is in the subtree and D H T i j ∉ , when it is not: 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 H i can be found, and those that do not have H i as a descendent. The goodness of an assignment for a value of q can then be rewritten as: Alternatively, and instead of assigning ambiguous reads based on the unique q value that maximizes G q , we can look for the assignment that maximizes the expected distance E(D i ) 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 R i can be assigned to nodes n 1 , ..., n n for each of the q values q 1 , ..., q n . The probability of Hi being among the leaves of the subtree T i,j rooted at n j is p = M i,j /M i , and that of not being included is 1p = (M i -M i,j ) = M i . The expected distance for a read Ri mapped to a node j for a given q is:

⎠ ⎠ ⎟
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 R i 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 H i is among its leaves have expected distance 0, since M i ,j = M i . 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 H i is among its leaves.
Proof. Since the true hit H i is among the leaves of the subtree rooted at the LCA, the TPR H i is 1. The FPR H

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.
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.
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.

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.
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.

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. 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 G q 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 (M i ) and a lower probability of choosing the true hit (M i , j /M i ).
The expected distance metric E(D i ) 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(D i ) than with G q , indicating that a combination of values of q is a better predictor than assigning all reads using a unique value. 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).

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. Percentage of reads assigned at the node selected for each value of q when maximizing E(D i ) in simulations using the full-length 16S rRNA sequence (top) and the V1-V2 hypervariable region (bottom) for reads of length 100, 150, and 200 bp. The column ∑ D i indicates the best sum of distances achieved.

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/.