INTRODUCTION

Top Abstract Introduction Methods Results & Discussion References

Overview of Sequence Data Processing

At present, nearly all DNA sequencing is done using the enzymatic dideoxy chain-termination method of Sanger (Sanger et al. 1977). One starts with a purified DNA template of interest (usually in single-stranded form) and an oligonucleotide primer complementary to a specific site on the template strand. For each of the 4 bases (A, C, G, T), a reaction is carried out in which DNA polymerase synthesizes a population of labeled single-stranded fragments of varying lengths, each of which is complementary to a segment of the template strand and extends from the primer to an occurrence of that base. These fragments are then separated according to length by gel electrophoresis, whereupon their relative sizes, together with the identity of the final base of each fragment, allow the base sequence of the template to be inferred.

In automated sequencing (Smith et al. 1986), the fragments are labeled with fluorescent dyes attached either to the primer (dye primer chemistry) or to the dideoxy chain-terminating nucleotide (dye terminator chemistry) (Prober et al. 1987). Typically a different dye is used for each of the four reactions, so that they can be combined and run in a single gel lane (in the case of dye-terminator chemistry, this also allows all four reactions to be carried out in a single tube). Multiple templates (36 or more at a time) are analyzed in separate lanes on the same gel. At the bottom of the gel, a laser excites the fluorescent dyes in the fragments as they pass, and detectors collect the emission intensities at four different wavelengths. The laser and detectors scan the bottom of the gel continuously during electrophoresis in order to build a gel image in which each lane has a ladder-like pattern of bands of four different colors, each band corresponding to the fragments of a particular length.

Computer analysis is then used to convert the gel image to an inferred base sequence (or read) for each template. Typically this analysis consists of four distinct steps: lane tracking, in which the gel lane boundaries are identified; lane profiling, in which each of the four signals is summed across the lane width to create a profile, or trace, consisting of a set of four arrays indicating signal intensities at several thousand uniformly spaced time points during the gel run; trace processing, in which signal processing methods are used to deconvolve and smooth the signal estimates, reduce noise, and correct for dye effects on fragment mobility and for long-range electrophoretic trends; and base-calling, in which the processed trace is translated into a sequence of bases. This paper focuses only on the last step.

The processed traces are usually displayed in the form of chromatograms consisting of four curves of different colors, each curve representing the signal for one of the four bases and drawn left to right in the direction of increasing time to detection (increasing fragment size). An idealized trace would consist of evenly spaced, nonoverlapping peaks, each corresponding to the labeled fragments that terminate at a particular base in the sequenced strand. Real traces deviate from this ideal for a variety of reasons having to do with imperfections of the sequencing reactions, of gel electrophoresis, and of trace processing. Because of anomalous migration of very short fragments (caused by relatively greater effects of the dye and specific base sequence on mobility) and unreacted dye-primer or dye-terminator molecules, the first 50 or so peaks of a trace are noisy and unevenly spaced. Toward the end of the trace, the peaks become progressively less evenly spaced as a result of less accurate trace processing, less well resolved as diffusion effects increase and the relative mass difference between successive fragments decreases, and more difficult to distinguish from noise as the number of labeled fragment molecules of a given size decreases. In particular, poorly resolved peaks for the same base may yield a single broad, often lumpy peak.

In better resolved regions of the trace, the most commonly seen electrophoretic anomalies are compressions (Sanger and Coulson 1975; Sanger et al. 1977), which occur when bases near the end of a single-stranded fragment bind to a complementary upstream region, creating a hairpin-like structure that migrates through the gel more rapidly than expected from its length, thus causing a peak to be shifted left of its expected position. This can result in one peak being beneath another or in two successive peaks for the same base being merged into one. GC-rich sequence tends to be more compression prone than AT-rich sequence because of the greater likelihood of stable hairpins. Dye-terminator chemistry appears to resolve most compressions (Lee et al. 1992), possibly because the dye on the terminal nucleotide interferes with base-pairing or because of the use of deoxyinosine triphosphate in place of deoxyguanosine triphosphate; but this chemistry has its own data quality problems caused by reduced polymerase affinity for the dye-labeled terminal nucleotide, one particular problem being a substantial decrease in the signal for a G base following an A (Lee et al. 1992; Parker et al. 1996).

Other frequently seen problems include weak or variable signal strength and noise peaks not corresponding to a base. These can result from variation in the efficiency of chain elongation or termination by the polymerase because of sequence context effects or inefficient incorporation of the dideoxy nucleotide, from impure or degraded reagents, or from laboratory protocol errors. Secondary structure in the template can produce polymerase stops, resulting in a pileup of peaks at one point in the trace. The signal downstream of a run of mononucleotide or dinucleotide repeats frequently is degraded, possibly because of strand slippage during copying by the polymerase.

The goal of base-calling software is to produce a sequence as accurate as possible in the face of the above data problems. Some of the earliest base-calling software was part of the processing software installed on the first ABI sequencing machines (Connell et al. 1987). That software achieves impressive accuracy and remains the standard against which other methods must be judged. Although full algorithmic details have not been published, according to a recent ABI description of its base-calling software (ABI 1996), the program uses mobility curves to predict the peak spacing, and identifies the most likely peak in intervals of the nominal peak spacing, assigning an N in the absence of a good choice. Subsequently, it adds and removes bases using a criterion involving the uniformity of peak spacing.

The advent of high volume sequencing has prompted development of other programs (Giddings et al. 1993; Golden et al. 1993, 1995; Berno 1996). These all perform multiple gel image processing steps including base-calling and have the merit of allowing efficient centralized processing of the data on a computer independent of the sequencing machine; however, none of them report equaling the accuracy of the ABI software.

This paper and the following one describe a computer program phred that performs base-calling and assigns an error probability to each called base. Accuracy is shown to exceed that of the ABI base-caller.

	METHODS

Top Abstract Introduction Methods Results & Discussion References

Base-Calling Algorithm

Overview

Phase 1: Locating Predicted Peaks

1.	Find the predicted peak peakS and period S associated to the starting region S, using the above procedure (in this case the set of permitted periods consists of all possible periods, and the fast Fourier transform is computed to find the optimal period S).
2.	Set the region R to S, and the direction d to "rightward".
3.	Shift R in the direction d by an amount R, and relabel this new region as R. If R has well-defined spacing and d is "rightward," take the set of permitted periods to be {R  0.03, R, R + 0.03}; if R has well-defined spacing and d is "leftward," take the set of permitted periods to be {R  0.25, R  0.20, R  0.15,...,R + 0.25}; if R does not have well-defined spacing, take this set to be {R}.
4.	Find peakR and R. If the end of the trace has been reached, go to step 5; otherwise go to step 3.
5.	If d is "leftward," stop. Otherwise set R to S and d to "leftward," and go to step 3.

Phase 2: Locating Observed Peaks

Phase 3: Matching Observed and Predicted Peaks

Phase 4: Finding Missed Peaks

Accuracy Assessment Methods

In principle, the accuracy of a basecalling program is easily measured by aligning read sequences produced by that program to the correct sequence and tabulating discrepancies. There are a number of subtleties in such an analysis, however, and it is easy to get misleading results. The details of the alignment algorithm can significantly affect how errors are classified. Moreover, it is important to allow for variability of error rates within and between reads, since otherwise the most error prone parts of the reads, or a small subset of reads with a very high error rate, can unduly influence the results. It is also worth recognizing that not all regions of a read are equally important; the higher quality part tends to be the most useful in practice, and the error rates in that part are therefore of most relevance.

We assessed the accuracy of the phred and ABI base-callers in several large sets of reads from cosmid clones sequenced in three laboratories (Table 1). These included (set 1) 9 mammalian (human and mouse) cosmids, sequenced by L. Rowen in L. Hood's laboratory (University of Washington); (set 2) 9 C. elegans cosmids from the Washington University Genome Sequencing Center (R. Waterston); and (set 3) 36 human chromosome 7 cosmids, from the University of Washington Genome Sequencing Center (M. Olson). In each case the sequencing strategy consisted of a moderate to high depth (6×  10×) shotgun phase in which dye primer reads were obtained from M13 subclones, followed by a finishing phase in which additional reads (mostly dye terminator) were used to close gaps and resolve low quality regions. Trace processing was performed using ABI analysis software.

The sets differed somewhat in the machine running conditions that were used and in the sequencing polymerase (see footnote to Table 1); in particular we separate the analyses of the cosmid sets one and two from the analyses of cosmid set three because set three had average read lengths several hundred bases longer than the other two sets, which would complicate interpretation of combined results. Dye primer reads were analyzed separately from dye terminator reads.

For each cosmid we created two FASTA formatted files, one containing the ABI base-called reads, and the other containing the phred base-called reads (as generated from the ABI-processed trace files). Each of these was screened to mask out vector sequence, and the screened reads were then aligned to the finished cosmid sequence and discrepancies tabulated.

Alignment of the reads to the finished sequence was performed by use of a restricted Smith-Waterman (Smith and Waterman 1981) algorithm as implemented in the program cross_match (P. Green, in prep.), which searches bands in the Smith-Waterman matrix surrounding word matches between the query and target sequences. It is important to recognize that the details of the alignment (its precise extent and the positions and types of the discrepancies) depend to some degree on the parameters used in the Smith-Waterman algorithm, which in our case were +1 for a match, 2 for a mismatch, 4 for a gap-initiation penalty, and 3 for a gap-extension penalty. The error rate in sequence reads tends to increase progressively toward the end of the read, and the local alignment stops when the error rate starts becoming high enough to produce a negative score. With our parameter values, this occurs where the error rate starts to exceed roughly 30% (depending on the distribution and type of the discrepancies). In regions containing multiple errors, there may be several possible ways of aligning the read against the cosmid sequence. For example, with compressions (where typically a set of peaks are shifted leftward), frequently the base-caller will get the number of bases correct, but it will omit one base (the first shifted peak in the group) and then compensate by inserting an additional base (or N) following the last shifted peak. Such a miscalled region may be aligned against the true sequence in either of two ways: by having a single deletion and a single insertion flanking several correctly matching bases; or by having no indels but multiple substitution discrepancies. Which of these has the higher score will depend on the relative sizes of the indel and mismatch penalties, as well as on the precise sequence of the shifted bases. Our parameter values tend to modestly favor alignments with substitution discrepancies rather than indels.

The relative performances of the ABI and phred base-callers were assessed by use of a computer program written to compare the cross_match outputs for the two sets of reads for each cosmid. The number of aligned bases for each method (i.e., the number of bases in the Smith-Waterman alignment, summed over all reads) was first determined. For meaningful comparison of the error rates of two different base-calling methods it is important to restrict attention to the part of each trace for which errors can be counted for both methods simultaneously (if one instead counts errors for each method in the entire part of the read alignable for that method, the more accurate method may actually appear less accurate, because in general it will have a longer alignable read length and the additional aligned bases will be in a lower quality part of the trace that has an intrinsically higher error rate). We did this by considering for each trace the segment of the cosmid sequence that is aligned to the ABI-called read by cross_match, the segment of the cosmid that is aligned to the phred-called read, and then taking the intersection of these two segments (this intersection is called the set of jointly alignable bases for the trace). Discrepancies inside this segment were tabulated; those outside it were ignored. Discrepancies were classified as substitutions (not counting Ns), Ns, insertions, and deletions, and tabulated in 100-base intervals with respect to position in the read sequence. The error rate for a particular error type in a particular 100-base interval in a data set is computed as the total number of errors of that type in that interval, divided by the total number of jointly alignable bases in that interval in the data set.

Table 2 shows the number of jointly alignable bases in each interval of 100 read bases in the cosmid sets. From positions 100 to 199, the number of aligned bases is nearly equal to the number of aligned reads multiplied by 100, reflecting the fact that most reads have alignable lengths of 200 or more. The number of aligned bases in positions 1 to 99 is lower, in part because the first few bases of the read are from the sequencing vector and in part because lower trace quality for the first 50 or so peaks sometimes results in enough basecalling errors to prematurely truncate the alignment. The number of aligned bases decreases slowly from position 200 up to about position 400 in cosmid sets one and two, and up to about position 800 in the data sets of cosmid set three, reflecting the fact that most reads have relatively high data quality in these regions. Thereafter, the numbers of aligned bases drop rapidly owing to degradation of trace quality toward the end of the read.

	RESULTS AND DISCUSSION

Top Abstract Introduction Methods Results & Discussion References

We compared the accuracy of the ABI and phred base-callers in three sets of cosmids (see Table 1, and Accuracy Assessment Methods). Although there is some variability within each set, on average cosmid sets 1 and 2 may be considered to be of average quality (Taq polymerase, short gels), whereas cosmid set 3 is of high quality (TaqFS polymerase, long gels). Sets 1 and 2 had very similar error profiles and were combined for the analyses described below. Dye primer and dye terminator reads from each set were analyzed separately.

Figures 1-4 summarize the error rates by error type and read position in the phred and ABI reads for each data set. (To permit meaningful comparison, we count errors only in the part of each trace for which both the ABI and phred base-calls are accurate enough to be jointly alignable; see Accuracy Assessment Methods.) Consistent with our quality designation above, set three has fewer errors at corresponding read positions than do sets 1 and 2, and the dye primer reads in each set have fewer errors at corresponding positions than the dye terminator reads. In general the error rates for both ABI and phred are quite low in the middle (high quality) parts of the reads, rising sharply later in the reads. There is an inflection point (reduction in the rate of increase) in many of the overall error rate plots later in the reads that corresponds closely with the point at which the number of alignable bases drops (Table 2). This appears to be attributable to a subpopulation of reads that have significantly higher quality, and thus longer alignable lengths and lower error rates, than the bulk of the reads.

View larger version (47K): [in this window] [in a new window]   Figure 1   Phred (gray bars) and ABI (black bars) error rates for dye primer data in cosmid sets 1 and 2.

View larger version (43K): [in this window] [in a new window]   Figure 2   Phred (gray bars) and ABI (black bars) error rates for dye terminator data in cosmid sets 1 and 2.

View larger version (54K): [in this window] [in a new window]   Figure 3   Phred (gray bars) and ABI (black bars) error rates for dye primer data in cosmid set 3.

View larger version (54K): [in this window] [in a new window]   Figure 4   Phred (gray bars) and ABI (black bars) error rates for dye terminator data in cosmid set 3.

Phred's total error rates (lower right panel in Figs. 1-4) are lower than ABI's in all parts of the read, in each data set. Overall, for the dye primer data, phred has ~41% fewer errors than ABI in sets 1 and 2, and 52% fewer errors in set 3; for the dye terminator data, phred has 39% fewer errors in all three sets. It is interesting that the relative improvement is greatest in the highest quality data set (set 3 dye-primer reads). Table 1 shows that the phred-called reads on average had about 5% longer alignable read length than the ABI reads, which is consistent with phred's greater accuracy throughout the read. To ensure that these results are not attributable to a small subset of reads with high error rates, we repeated the analyses after excluding the 5% of ABI reads and the 5% of phred reads having the highest error rates. The same relative differential (40%-50% fewer errors with phred) was seen on the pruned data sets.

It is instructive to consider the different error types. Insertion or deletion errors (indels) are more serious than substitution errors for many purposes [e.g., in expressed sequence tag (EST) data, where they may change the reading frame thus making homology detection more difficult]. The combined indel rate for phred is lower throughout the read length in each data set, except past base 1000 in set 3 where ABI does slightly better. In the regions of higher trace quality (i.e., before the sharp dropoff in numbers of aligned bases) phred's insertion and deletion rates are both lower than ABI's. After this point, the rates for both ABI and phred increase substantially but with different tendencies. The phred deletion rates increase somewhat more quickly than ABI's and are substantially higher than ABI's towards the end of the trace; conversely the ABI insertion rates increase more rapidly than phred's. It appears that the ABI base-caller is more likely to add a peak than to remove one, while phred avoids high insertion rates with its tightly constrained peak prediction algorithm. The elevated phred deletion rates later in the trace are probably attributable to the fact that phred freezes the predicted peak spacing value where it can no longer reliably determine it, which tends to underestimate the number of predicted peaks when the true spacing decreases.

Phred makes roughly the same numbers of substitution errors as ABI, except for the terminator reads of cosmid sets one and two where the phred substitution rate is about twice the ABI rate in the high quality part of the trace. However, ABI calls substantially more N's than phred throughout, in part because phred only calls an N when it cannot find an observed peak to associate to a predicted peak location, while the ABI software assigns an N when either it cannot find a peak or the call is uncertain because of noise peaks. The combined substitution +N rate is substantially higher for ABI than for phred in all data sets.

Because we consider any failure to predict the correct base as an error, N's are errors. This is the only sensible policy, as otherwise a base-caller could inflate its accuracy simply by calling N's in any low-quality region; for example, phred could attain an error rate of <1 in 10,000 by designating every base with quality value <40 as an N (Ewing and Green 1998). Assigning an N to an ambiguous peak has some merit in the absence of quality values because it helps to indicate low trace quality, but is not useful when quality values are available. Nonetheless in fairness to the ABI software one must acknowledge that its performance under our criteria could be improved by modifying it to always guess a base. If we assume (optimistically) that by doing so the ABI software could convert 75% of its N's to correct base-calls, with the remaining 25% becoming substitution errors, then the ABI substitution error rates would approximately equal the phred rates in the high quality regions of the traces, but would exceed the phred rates later in the traces in all cases except the terminator data for cosmid sets one and two. Total error rates would still be less for phred than for ABI.

Sensitivity to Lane-Processing Method

We believe that phred's observed peak detector and its methods for aligning observed peaks to predicted peaks will work well regardless of the trace source. This is not true of the peak prediction algorithm, however, which was tuned on ABI-processed traces and works well only when the trace processing includes a step to render the peak spacing relatively uniform; without such peak spacing normalization, phred may predict peaks poorly in some regions of the trace, resulting in increased numbers of indels. Other specifics of the lane processing may also affect phred performance. The ABI mobility correction appears to be less accurate toward the beginning and end of the trace, affecting basecalling accuracy there. In addition, it appears that the ABI trace processing broadens peaks excessively, which has the effect of reducing resolution and increasing the numbers of deletion errors.

We have recently developed a trace-processing program, plan, that reads raw trace data and processes the trace arrays. Preliminary tests indicate that phred's performance on traces processed using plan is comparable to, or slightly better than, that using ABI processed traces.

Further Improvements

What further accuracy improvements are possible? The high-quality part of the trace presents different issues from the later, lower quality part. For most applications where an accurate finished sequence is required, the high-quality part of the trace is by far the most important, since it is the only part usable in deriving the final sequence. Raw data quality in the later parts of the trace is such that, although improved basecalling accuracy is certainly possible, it would be extremely difficult to reduce the error rate to the level necessary to infer highly accurate sequence.

In the high-quality part of the trace the major sources of error at present are compressions in dye primer data, and highly variable peak heights or missing peaks in dye terminator data. Both of these are sequence-context dependent, and classifying the specific sequence motifs that are most error prone is a promising approach to improving accuracy. It is particularly important to improve phred's ability to identify and resolve CC and GG compressions, which are often difficult even for skilled human finishers to detect by eye, and we have recently made some progress in this regard using rules for hairpin-prone sequences.

Improved accuracy in the lower quality part of the trace would be useful in single read applications (e.g., EST sequencing), and to improve repeat discrimination and assist in making joins in the early phases of shotgun sequencing projects (which helps make collection of additional data more efficient by delineating the sizes of regions where more data are needed). The major gains to be made here are in improved peak spacing estimation and in processing to improve peak resolution.

Program Availability

C source code is available at no charge to academic researchers for research purposes, and by commercial license from the University of Washington to other users; contact Brent Ewing at bge{at}u.washington.edu.