Base-Calling of Automated Sequencer Traces Using Phred. II. Error Probabilities

Overview of Error Probability Issues

Obvious requirements for the error probabilities are that they must be predictive, that is, assigned by use of a method that does not require knowledge of the true sequence; and they must be valid in the sense that they correspond to observed error rates, that is, the set of bases assigned a particular error probability p should have an actual error rate equal to p.

Given these constraints, there are still many possible ways of assigning error probabilities to base-calls. For example, if a set of 1,000,000 base-calls contains 10,000 errors, then a method that assigns an error probability of 0.01 to each base-call is valid; but so is another method that identifies two sub-sets of 500,000 base-calls each, with the first set containing 9000 errors and the second set 1000 errors, and assigns an error probability of 0.018 to every base-call in the first set and 0.002 to every base-call in the second set. While the second method is no more valid than the first, it clearly does a better job at discriminating the less accurate base-calls from the more accurate.

We can formalize the notion of discrimination ability as follows. Intuitively, a method with high discrimination ability should spread out the error probabilities (or quality values) as much as possible. By that criterion, one natural measure of discrimination power is simply the variance of the quality value distribution. However, in practice, it is the most accurate read data that are most important, as those are the data used to derive the consensus sequence, and consequently we prefer to judge different methods by how well they perform at identifying a subset of the bases with a very low error probability. Specifically, given a set B of base-calls and a valid method of assigning an error probability e(b) to each base-callb, then for any subset B of B, the expected number of errors in B is Σ_b∈B e(b)(since the error probability for a base-call is also the expected number of errors for that call); and the expected error rate forB is the expected number of errors in B, divided by the number of base-calls in B. It is easy to see that for any given error rate, r, there is a unique largest set of base-calls, B _r, having the properties that (1) the expected error rate of B_r is ≤r, and (2) whenever B_r includes a base-call b, it includes all other base-calls whose error probabilities are ≤ e(b). The discrimination power at the error rate r is then defined to be that is, the number of base-calls in B_r divided by the total number of base calls. P_r measures the effectiveness of the error probability assignments at extracting a subset of bases having a low error rate r. (There is a precise analogy here to the notions of validity and power of a statistical test.)

Our goal in this study was to develop error probabilities that are valid and have high discrimination power at small values of r(r ≤ 0.01). For this purpose, we used parameters computed from the processed trace from which the base-calls were made, focusing most on those parameters that appear to play a role in intuitive human assessments of data quality as data quality should be predictive of error probabilites. A number of different sets of such parameters were tested, using an algorithm (described below in Error Probability Calibration) that, given a set of parameters and a training set of reads for which it is known which base-calls are correct and which are errors, finds a way of associating parameter values to error probabilities that has (near) maximum discrimination power for smallr.

Trace Parameters

Errors by the basecaller often are attributable to misinterpretation of peaks in a region of the trace, and, as a result, indications of the error may be present in the vicinity of the erroneous peak but not at the peak itself. Consequently the parameters most effective at detecting errors tend to be those that consider a window of the trace that includes several peaks flanking the one whose base-call is being assessed. The parameters that turned out to be the most powerful in our tests were all of this type. A by-product of using parameters computed from such a window is a smoothing of the parameter values (and hence of the error probabilities) from base to base.

The following four parameters were found to be particularly effective at discriminating errors from correct base-calls. In each case, smaller parameter values correspond to higher quality (more accurate sequence).

1.

Peak spacing. The ratio of the largest peak-to-peak spacing, in a window of seven peaks centered on the current one, to the smallest peak-to-peak spacing. The minimum possible value of one corresponds to evenly spaced peaks.

2.

Uncalled/called ratio. The ratio of the amplitude of the largest uncalled peak, in a window of seven peaks around the current one, to the smallest called peak; if there is no uncalled peak, the largest of the three uncalled trace array values at the location of the called base peak is used instead. [An uncalled peak is a peak in the signal that was not assigned to a predicted location by phred(Ewing et al. 1998) and thus does not result in a base call.] If the called base is an N, phred assigns a large value of 100.0. Note that this is not what is sometimes called the signal to noise ratio, as uncalled peaks may be true peaks missed by the base-calling program rather than noise in the conventional sense. The minimum parameter value is 0 for traces with no uncalled peaks.

3.

Same as 2, but using a window of three peaks.

4.

Peak resolution. The number of bases between the current base and the nearest unresolved base, times −1 (to force the parameter to have the right direction). (A base is unresolved if it is called asN or if for at least one of its neighboring bases, there is no point between the two corresponding peaks at which the signal is less than the signal at each peak). The minimum possible parameter value is half the number of bases in the trace, times −1, and the maximum value is 0.

Error Probability Calibration

Given a set of parameters that can be computed from the trace for each base-call, a set of threshold values for those parameters is said to be optimal for a particular error rate r, if the setB consisting of all bases whose parameter values are less than the threshold values has an error rate ≤r, and no other set of thresholds yields a larger set with error rate ≤r.We describe below a simple greedy algorithm that finds a nearly optimal set of thresholds for small error rates. With respect to linear discriminant analysis (which has similar goals), our method has the advantages of not assuming multivariate normality (which is very far from being true in our case) or other distributional properties for the parameters, of allowing parameters that take on non-numerical values, and of not requiring that parameters be transformed to normality or that outliers receive any special treatment. It is compute-intensive, however, and, in the form described, can only be used with a relatively small number of parameters simultaneously. The only significant assumption about the parameters is that their values should be ordered such that small values tend to correspond to more accurate base-calls and large values to less accurate base-calls.

The algorithm produces a lookup table consisting of a set of lines, each line containing a set of parameter thresholds, together with the error probability and quality value corresponding to those thresholds (there can be multiple lines having the same quality value). Although we have only applied the algorithm to generate a single error probability for each base (combining all types of error—substitution, deletion, and insertion), it can easily be adapted to produce separate probabilities for each error type.

The basic idea of the algorithm is as follows. One starts with a small number of parameters (in our case, four). First, a finite number of threshold values for each parameter is selected; we allow 50 different thresholds, a number large enough to avoid significantly sacrificing resolution but small enough that it is computationally feasible to consider each possible threshold for each parameter in subsequent steps. Then each 4-tuple of parameter thresholds (one for each parameter) is considered in turn, and the empirical error rate is computed for the set of bases defined by those thresholds. The 4-tuple for which the empirical error rate is smallest is selected (in the event of ties, the largest set having a given error rate is taken) and defines the first line of the lookup table. The set of bases defined by these thresholds is then eliminated, and the process is repeated using the remaining bases. Iteration of this procedure produces the desired lookup table.

Note that by construction the first set of thresholds found by this procedure is nearly optimal for its error rate, in the sense defined above (it may fail to be strictly optimal because only a finite set of thresholds are considered for each parameter). Subsequent sets of thresholds also tend to be nearly optimal, but in this case the effect of eliminating bases at earlier steps is another factor which may cause strict optimality to fail. This effect tends to be small in the early steps when error rates are small, but may become more significant later.

We now give a more precise description of the algorithm, in the case where four different parameters r, s, t, and u are used. In the following, if p denotes a parameter, thenp(b) indicates the value of the parameter for a particular base-call b. First, for each parameter pwe find valuesp ₀ < p ₁ < … .  < p ₄₉ < p ₅₀such that all p(b) lie between p ₀and p ₅₀, and the number of bases bsatisfying p _i-1 < p(b) ≤p_i is approximately the same for all i.

We define a cut to be a 4-tuple (i, j, k, m) where i, j, k, and m range between 1 and 50. There are 50⁴, or 6.25 million cuts in all; although this is large relative to the number of data points, we will avoid overfitting by imposing a monotonicity condition. Each cut has a corresponding set of parameter thresholds r_i, s_j, t_k,and u_m, defined as above. For each cut (i, j, k, m), define err _(i,j,k,m) to be the total number of erroneous base calls b below the cut, that is, satisfying r(b) ≤ r_i, s(b) ≤s_j, t(b) ≤ t_k, andu(b) ≤ u_m. Similarly, letcorr _(i,j,k,m) be the total number of correct base-calls below the cut.

The error rate below the cut,e _(i,j,k,m), is defined by and the corresponding quality value by Here 1.0 is a small-sample correction added to ensure that both the numerator and denominator are positive; it produces a slight upward bias in e, that is most pronounced for very small e.We round q _(i,j,k,m) to the nearest integer. The following two steps are now iterated to create the lookup table.

1.

Find the cut (i,j,k,m) for whichq_(i,j,k,m) is largest. In the event of ties, take the one for which corr_(i,j,k,m) +err_(i,j,k,m) is largest; if there is more than one of these, take the one for which the sum of the indices is highest. Output q_(i,j,k,m), e_(i,j,k,m) and the parameter valuesr_i, s_j, t_k, andu_m. Delete (i, j, k, m) from the list of cuts.

2.

For each remaining cut (i′, j′,k′, m′), adjust the countserr_{(i′,j′,k′,m′)}andcorr_{(i′,j′,k′,m′)}by deleting bases below the removed cut (i, j, k, m), and recomputee_{(i′,j′,k′,m′)}andq_{(i′,j′,k′,m′)}using the new values. If err and corr are 0 for all remaining cuts, stop. Otherwise go to step 1.

Note that, by construction, the error probabilitye output in step 1 is the (small-sample corrected) error rate for the set of bases defined by the property that their parameter values are less than or equal to the parameter thresholds output in this step, but not less than or equal to the thresholds output in any previous line. As a result, given the parameter values for a base, one can find an appropriate marginal error probability for that base by looking through the table until the first line is found in which the parameter thresholds equal or exceed the parameter values of the base in question and then reading the error probability on that line.

Cosmid Sets

For these studies we used the ABI-processed trace data from four sets of cosmids (Table 1): two training sets consisting of 9 mammalian cosmids from L. Rowen (L. Hood’s laboratory, University of Washington), and 9 Caenorhabditis eleganscosmids from the Washington University Genome Sequencing Center (R. Waterston), and two test sets consisting of 22 C. eleganscosmids from the Washington University Genome Sequencing Center, and 36 human chromosome 7 cosmids from the University of Washington Genome Sequencing Center (M. Olson). Phred basecalls (Ewing et al. 1998) were generated for each trace, and reads were then screened for sequencing vector and aligned against the finished cosmid sequence using cross_match as described previously (Ewing et al. 1998). Each read base in the alignable part of the read was classified as correct (matching the final sequence) or erroneous (discrepant with the final sequence); deletion errors (i.e. cases where one or more bases in the cosmid sequence are missing from the read) were assigned randomly to one of the two read bases adjacent to the deletion. Unaligned bases are ignored in the following.

Quality Assignment

Application of the error probability calibration algorithm, with the parameters described above, to a single combined training set (containing all reads in cosmid sets 1 and 2) produced a lookup table with 2011 lines. This table is used to assign quality values to the base-calls from a particular read, as follows. For each base call,phred computes the four parameter values, and then searches the lookup table line by line, in order, until it finds a line in which each of the four parameter values is at least as large as the corresponding parameter value for the base-call. The quality value associated to that line is then assigned to the base. If no such line is found, the base-call is assigned a quality of 0.

This procedure was used to assign a quality value to each base-call in the two test sets of cosmids (Table 1). Initial examination of errors and quality values in the test data sets showed more errors than expected among bases with predicted quality values of 40 and above. On inspection, we found five kinds of spurious contribution to the error counts:

1.

Chimeras. Chimeric reads contain two or more noncontiguous segments that are incorrectly juxtaposed, as a result of a cloning artifact (a chimeric subclone) or a gel lane tracking error. In such cases, cross_match finds an alignment involving one of the two pieces, but occasionally the alignment extends spuriously for a few bases into the adjacent piece because of fortuitously matching nucleotides and includes spurious high quality discrepancies in this extension.

2.

Contaminant or Misassembled Reads. Most data sets have a small number of contaminant reads arising either from sample mistracking or from low levels of contamination of the subclone libraries. If the contaminating read contains a repeated sequence, it may spuriously appear to match the cosmid sequence with multiple high-quality discrepancies. Similarly a noncontaminant read lying in a near-perfect repeat in the cosmid may be aligned against the wrong copy of the repeat by cross_match if errors in the low-quality part of the read result in a higher Smith–Waterman score against the wrong copy. Such cases are generally revealed by examining all matches of the read to the same and other cosmids.

3.

Subclone Mutation. Some errors in high-quality trace regions appeared as a single inserted or deleted base in a long mononucleotide run, or as a single deleted unit of a microsatellite repeat. In several such cases another read obtained from the same subclone showed the same discrepancy. These examples appear to represent spontaneous subclone mutations. Another type of subclone mutation seen was a larger deletion, apparently mediated by recombination involving an imperfect direct repeat in the subclone. In this case the alignment found by cross_matchincluded the part of the read lying to one side of the deletion, but extended spuriously into the region on the other side of the deletion (often for a significant distance) because of the direct repeat, resulting in multiple high-quality discrepancies.

4.

Unremoved Vector. Occasionally multiple errors in the sequencing vector part of a read preventedcross_match from finding the match between it and the vector sequence, so that it remained unmasked; and the alignment of the insert against the cosmid sequence found bycross_match continued into this unmasked sequence because of a few fortuitous base matches at the end of the vector sequence, and the extension included spurious high-quality discrepancies.

5.

Misalignment. In some cases phred deleted the first base in a mononucleotide run of three or more bases, when imperfect mobility correction shifted the run close to the preceding peak, but (owing to the details of the Smith–Waterman algorithm) the alignment constructed by cross_match instead indicated the site of the deletion to be the last base in the run. The called bases on each side of the deleted peak received low-quality values because of bad peak spacing and signal to noise, but the last base in the run received a high-quality value.

To identify and remove as many of the above cases as possible, we examined each read having two or more discrepancies with quality values of 40 or more, using consed (Gordon et al. 1998). Reads that appeared to be unambiguous instances of one of the first four types above were eliminated; in all, 10 reads (from a total of 26091) were removed from cosmid set 3, and 35 (from a total of 27184) were removed from cosmid set 4. It is likely that a small number of spurious high-quality discrepancies of the above type occurred only once in a read and thus escaped detection by this procedure and that others occurring at quality values below 40 were also not detected. Instances of high quality discrepancies occurring in mononucleotide runs of two or more bases were also inspected and (in cases where the wrong base was flagged) the quality value altered to that of the base actually deleted by phred. This resulted in changing 10 quality values in cosmid set 3, and 30 in set 4. Following these corrections to the data set, the observed quality values were recalculated.