All Your Base: a fast and accurate probabilistic approach to base calling
© Massingham et al.; licensee BioMed Central Ltd. 2012
Received: 29 July 2011
Accepted: 29 February 2012
Published: 29 February 2012
The accuracy of base calls produced by Illumina sequencers is adversely affected by several processes, with laser cross-talk and cluster phasing being prominent. We introduce an explicit statistical model of the sequencing process that generalizes current models of phasing and cross-talk and forms the basis of a base calling method which improves on the best existing base callers, especially when comparing the number of error-free reads. The novel algorithms implemented in All Your Base (AYB) are comparable in speed to other competitive base-calling methods, do not require training data and are designed to be robust to gross errors, producing sensible results where other techniques struggle. AYB is available at http://www.ebi.ac.uk/goldman-srv/AYB/
There can be little doubt that the vastly increased throughput of Next-Generation Sequencing (NGS) machines has revolutionised DNA sequencing, but the reads produced are both shorter and less accurate than those from capillary sequencing and discoveries from NGS are often verified using traditional sequencing . The challenges to overcome to improve the accuracy and read length of NGS platforms are different from those that were faced by capillary sequencing  and require different strategies to tackle them. In particular, the phasing process -- individual molecules of DNA becoming out-of-step with others in the same cluster -- is complex and ultimately limits the length of reads which can be obtained from cluster-based sequencing-by-synthesis methods . Here we develop an explicit statistical model of the sequencing process, including phasing and other signal-degrading processes. By implementing a base calling algorithm based on this model, our AYB software is able to produce more accurate reads.
Our statistical model is quite generic and so applicable to all sequencing-by-synthesis and similar platforms (sequencing-by-ligation, pyrosequencing; see Metzker  for a comparison) that rely on large numbers of clusters consisting of many homogeneous DNA molecules. We concentrate on the Illumina Genome Analyser (GA-II), both to provide a concrete foundation to aid exposition of the methods and because of local availability of data for testing and comparison. The mechanics of sequencing-by-synthesis on the GA-II platform have been described elsewhere in detail . Here we present an overview of an idealised sequencing process to establish context and terminology for the rest of this paper and then a critique to show how errors arise.
Fragmented single-stranded DNA is washed through the lanes of a slide, where it attaches and is amplified to form a sequence-homogeneous cluster of molecules. Sequencing progresses in steps, referred to as cycles, with each cycle conceptually sequencing one position of DNA. For each cycle, a mixture of Fluorophore-Labeled Nucleotides (FLNs) is washed through the lanes of the slide and attach to the molecules in each cluster; attachment of more than one nucleotide in a given cycle is prevented by the presence of a reversible terminator element on each FLN. A different fluorophore is associated with each of the four nucleotides (A, C, G, T) and so the nucleotide sequence of the DNA fragments can be uniquely identified from the fluorescence. After the attachment process has run to completion, the intensity of fluorescence from each cluster is recorded in four channels, each channel being a combination of illumination with a specific laser and imaging through a specific filter. Clusters are artificially grouped into tiles, regions of the lane consistent over cycles, whose size is constrained by the capacity of the imaging equipment. The terminator elements and fluorophores are then cleaved from the FLNs, setting up each cluster so the FLNs in the next cycle will attach to the next position of sequence.
After processing the images to pick out individual clusters, the output of the sequencing machine is many channel × cycle matrices of intensities, one matrix for each cluster. In principle the bases could be called straight from these intensities but there are several complicating factors  that must be dealt with, cross-talk, phasing and dimming being of particular importance.
Cross-talk is the recording of light from a single fluorophore in multiple channels. This occurs because, although they are chosen to be distinguishable, the fluorophores' emission spectra overlap. There is not a one-to-one correspondence between channels and FLNs and the relationship between the emission of each fluorophore and the intensity observed in each channel needs to be ascertained and corrected for. Phasing refers to the deterioration in relationship between sequencing cycle and sequence position as the cluster loses coherence: on a given cycle, FLNs may be attaching to different positions on different molecules within the cluster. There are many possible explanations for phasing: for example, a FLN might have a defective reversible terminator element leading to the attachment of two FLNs to a molecule on a single cycle, allowing the molecule to get ahead in the sequencing process ('pre-phased'), or the cleaving of the reversible terminator might fail for a cycle so the molecule lags behind when the element is finally removed ('post-phased'). A further possible cause of post-phased molecules is the chemistry not running to completion, resulting in either no FLN being attached that cycle or cleaving failure as previously mentioned. Finally, molecules within a cluster gradually stop contributing to the total signal, possible causes being laser damage to the individual molecules or problems reversing the terminator element, and this leads to a decrease (dimming) in the overall emission observed from each cluster in later cycles of sequencing.
The cross-talk is a consequence of the physics of fluorophore excitation and methods for estimating it have already been developed for dye-terminated capillary electrophoresis sequencing platforms . Phasing and dimming are more specific to NGS methods, the Illumina platform in particular, and have been approached in a variety of ways. The Illumina base calling software (Bustard) assumes a constant rate of post-phasing and pre-phasing for all cycles , as do the Alta-Cyclic caller  and Rolexa , whereas BayesCall  allows the phasing at each position of the sequence to depend on several of the neighbouring bases and Ibis  assumes that all information about the phasing at a given cycle is contained in intensities of the cycles either side. The Seraphim base-caller  assumes that cross-talk and phasing have already been corrected for, empirically estimating the sequencing noise and calling bases using a log-linear model over cycles to model base-specific signal degradation. In contrast our method uses a completely empirical model, generalising both cross-talk and phasing, that allows all aspects of the sequencing process to be determined by the data on a run-by-run, indeed tile-by-tile basis. This model is unique to AYB and incorporates effects such as cycle-wise variation in cross-talk and allows context-specific phasing rates that may account for some of the reported sequence specific errors in GC rich reads [11, 12].
Results and Discussion
In a recent comparison , Ledergerber and Dessimoz compared several different base callers and Ibis  was clearly the most accurate and considerably quicker than any comparable base caller. While Naïve Bayescall  is almost as accurate as Ibis and boasts a much improved base calling time over the original Bayescall algorithm  (37 min. compared to the 2266 min. for 300K reads, respectively, for Ledergerber and Dessimoz's example), it is still c. 24 × slower than Ibis and requires an already trained model to be available (Ledergerber and Dessimoz report 1842 min. for training). Despite its considerable improvement over the original algorithm, the computational requirements of Naïve Bayescall restrict its use and, since it is also strictly dominated by Ibis in both accuracy and speed, we do not include it in the comparisons below. To quantify the performance of AYB, we therefore use Bustard and Ibis as comparators to represent standard practice (Bustard; various versions depending on age of data) and the current state-of-the-art (Ibis; version 1.1.5).
Bustard is the standard base calling software for the Illumina platform and a detailed description of the algorithm is given in . Ibis uses a machine learning approach to base calling, training a Support Vector Machine (SVM) at each cycle on the intensities of the current, previous and next cycles. The SVM is trained on true calls for some subset of the data, and these are estimated by mapping Bustard calls back to the reference; the trained SVMs are specific to given read lengths and cluster densities and may also incorporate artifacts that are specific to given run. If sequencing de novo then correct calls may not be available but an SVM from a similar run could be used, with potential reduced performance, or a small amount of a known genome could have been 'spiked-in' and the SVM trained using these reads. Ibis was trained by mapping the Bustard calls for a training set to the reference genome using the default aligner (a modified version of SOAP ), the training set comprising every tenth tile starting from the fifth for our full lane sets of data (B. pert., BGI and Illumina) and all reads for the reduced sets (ϕX174 L2, ϕX174 L4, ϕX174 L6, Ibis Test and HiSeq; see below for full details of test data sets).
Summary of data sets analysed and reference genomes used for mapping.
B. pertussis Tohama I
H. sapiens GRCh37
H. sapiens GRCh37
H. sapiens GRCh37 + ϕX174
Performance comparison of Bustard, Ibis and AYB on several sets of reads of varying read length and chemistry versions.
Reads mapped, %
Reads perfect, %
Training and base calling time for Ibis and AYB.
ϕX174, 76 cycle and 51 cycle
Two test sets from the bacteriophage ϕX174 were used: nine tiles each from three lanes of 76 cycle data produced by the Sanger Institute and differing in cluster density (named ϕX174 L2, ϕX174 L4, ϕX174 L6), and a decimated run containing 200K clusters of 51 cycle data that is distributed as a test set with Ibis  (Ibis Test). Reads from the 27 tiles from the Sanger Institute were aligned against a SNP-corrected genome (two SNPs corrected), such a correction being possible because of the high coverage produced by these tiles. Reads called from the Ibis Test set were mapped to the genome distributed with it. Given the small number of clusters in the Ibis Test set, it was analysed as whole with AYB rather than on a tile by tile basis.
Bordetella pertussis, 76 cycle paired-end
A second comparison was based on a data set comprising an entire lane (100 tiles) of 76 cycle paired-end reads from the coccobacillus Bordetella pertussis (data sets denoted B. pert./1 and B. pert./2 for the two paired ends, respectively), using the complete genome of the Tohoma I strain as a reference. The tiles from this run showed a large variation in the number of clusters, ranging from 345 to 82,000, with the tiles close to the ends of the flow cell containing fewer clusters. The sequence produced was generally of low quality, with Bustard producing an error rate of 50% for the final five bases of the first end of the read-pairs, which suggests that a problem occurred during the sequencing. Oddities in the cross-talk matrix, the channels corresponding to A and C nucleotides being noticeably brighter than those corresponding to G and T, suggest that there may have been illumination problems with one of the lasers. Because of the problems with this run and the presence of polymorphisms relative to the reference genome, it provides a useful comparison between the base callers when problems occur.
Mapping back to the reference revealed a marked difference in base-caller performance (table 2): AYB produced more than four times as many perfect reads as Bustard (1,133k vs. 265k reads) and 1.5 times as many as Ibis (752k reads) for the second end of the read-pairs. AYB produces 56% more mapped reads than Bustard and 13% more than Ibis. The increased number of perfect and close to perfect reads produced by AYB has real consequences for down-stream analysis, with the genome being covered to an average depth of 29.3 × for AYB, 9.5 × for Bustard and 21.7 × for Ibis. Greater coverage means more confident SNP and variant detection, which in turn leads to improved mapping of reads. For the first end of the read, AYB produces a greater number of perfect reads (176k reads) than Bustard (108k) or Ibis (133k) and 59% and 17% more mapped reads than Bustard or Ibis, respectively.
As well as mapping to a reference genome, the length of contigs produced by de novo assembly is a useful guide to the quality of reads produced and of relevance in cases where a reference is not available. Applying Velvet  to the second end of the B. pertussis paired-end reads (kmer length 31; automatic coverage cut off; default options otherwise) produces a N50 contig length of 6690 bases for the AYB reads; the reads from both Bustard and Ibis produce much shorter contigs on average, with N50 lengths of 2029 and 4473 bases, respectively.
Human NA19240, 45 cycle and 51 cycle paired-end
Much of the pilot data for the 1000 Genomes project  has been archived and is publicly available for reanalysis, allowing for a further comparison between the base-callers and showing that AYB can be usefully applied to improve the analysis of existing data. Two sequencing runs for NA19240 (Yoruban daughter) were reanalysed: ERR000479 (9.6 million 45bp paired-end reads, part of ERA000013 by the Beijing Genomics Institute, referred to as sets BGI/1 and BGI/2) and ERR000610 (14.0 million 51bp paired-end reads, part of ERA000023 by Illumina Inc., referred to as sets Illumina/1 and/2). The accuracy of the base-calling was assessed by mapping to the human reference provided by the 1000 Genomes Consortium, based on GRCh37. This genome has variants relative to the sample sequenced but their presence penalises all callers equally. The raw intensities submitted to the archive have apparently been filtered for quality since an abnormally high proportion of the reads map back to the reference when called with Bustard: 86% for the BGI run and 97% for the Illumina run.
Despite the limited scope for improvement, both AYB and Ibis produce slightly more mappable reads than Bustard (table 2): a 1.97% or 1.92% increase respectively for the BGI data and 0.41% or 0.38% increase for the Illumina data averaged over both ends of each set. Large increases over Bustard are observed for the proportion of reads that match the reference exactly, AYB showing a clear lead over Ibis on both sets of data with 14.65% and 5.73% increases over Bustard for the BGI and Illumina sets respectively, compared to increases of 12.74% and 4.56% for Ibis. The superiority of AYB here is surprising as this is not a situation where it would be expected to do particularly better: the number of cycles, and so phasing, is comfortably small in both cases. The read length, vintage and error-rate of the BGI run is consistent with the older "sticky-T" chemistry (incomplete cleavage of the 'T' FLN, leading to an increased concentration in later cycles) and the improvement seen is typical for similar data.
The per-mapped-base error rate for reads produced from all BGI and Illumina data sets by AYB is lower than that of either Bustard or Ibis for a variety of mapping criteria (Figure 1).
Variant calling of NA19240
Following the guidance for best practice variant calling  using the Genome Analysis Tool Kit , we implemented a SNP calling pipeline including the full recalibration and realignment of reads. For each base caller, the both BGI and Illumina sets of paired-end reads were combined and then variants for the NA19240 data set were called against the human reference GRCh37. After variant recalibration, 399091 of AYB's 438486 SNP predictions passed the truth sensitivity filter (set to 99.0), compared to 371959 of 407110 for Bustard and 386056 of 413198 for Ibis. Base calls produced by AYB show a 7.3% increase in filtered variant predictions over Bustard, compared with a 3.8% increase shown for variant calls produced by Ibis. The median quality of the calls was 29.9, 29.0 and 27.0 for AYB, Bustard and Ibis, respectively, so calling bases using the AYB base caller results in more variant calls with a higher average recalibrated quality than either of the other two callers.
Comparison between SNP predictions made by the three base callers and the 'CG-HQ' high-quality set of variant calls from Complete Genomics.
8.79 × 10-4
8.90 × 10-4
10.68 × 10-4
All three callers predicted a number of SNPs that were not contained in the CG-HQ set: 9.4% of the total for AYB, 8.8% for Bustard and 11.4% for Ibis table 4 (column labelled 'Extra'). Although potentially false positives, these predictions show a good degree of congruence with 40% of the total being common to all three base-callers and 57% being made by at least two base-callers, and so cannot all be discounted as errors. Finally, the median quality of these 'Extra' SNP predictions was lower than those that the CG-HQ set validates ('Match'): 24.1, 24.1 and 20.2 versus 30.0, 29.0 and 27.0 for AYB, Bustard and Ibis, respectively. For comparison, the average quality of SNP predictions at sites where the CG-HQ set predicted a different SNP ('Mismatch') was 24.1, 23.0 and 20.4, respectively. Possible explanations for these discrepancies are variation between the sequence libraries or in the cell lines used .
HiSeq, 101 cycle paired-end
Illumina Inc. made available to us a decimated set of data produced on an Illumina HiSeq machine, comprising 7.8 million paired-end reads (8 lanes with a read length of 101 bases) of human sequence with a ϕX174 spike-in (approximately 0.43% of reads). Ibis was trained separately for each of the 8 lanes, although cross validating revealed little difference in performance since all lanes were from the same sequence library. All statistics reported for Ibis are an average over all 8 models and lanes. The reads were mapped against the reference human genome and that of ϕX174 and, again, the error rates reported are all slightly inflated due to the genomes sequenced having variants relative to the references.
Even on the modern HiSeq chemistry AYB and Ibis improve on the base calls produced by Bustard, although the increase in mapped reads is quite small: 0.32% and 1.18% for the first end for Ibis and AYB, respectively. Larger improvements are seen for the number of perfect reads where both base callers improve on Bustard by several percent (3.75% for Ibis, 6.98% for AYB) so AYB improves on Ibis by almost as much as Ibis improves on Bustard. Greater improvements are seen for AYB on the second end of these reads, with an 11.63% increase over Bustard in the number of perfect reads. We were unable to get Ibis to train on the second end of the reads so results are not available.
For the HiSeq data sets, the per-mapped-base error rate for reads produced by AYB is lower than that of either Bustard or Ibis for a variety of mapping criteria (Figure 1).
Looking only at per-mapped-base error rates does not tell the whole story of base caller accuracy. Some reads are produced from extremely clean intensities whereas others may have been extracted from very noisy data, and base callers assign each base a quality score to indicate their confidence in that call.
Typically, the quality score is a discrete value related to the estimated probability that the call is correct. Given a set of mapped reads, the actual proportion of bases in error can be found for each (estimated) quality value assigned by the base caller; these proportions can be used to calculate empirical quality values to which the estimated values can be compared to assess their accuracy. If the estimated quality values were perfectly calibrated then they would agree, to within sampling error, with the empirical quality values (a linear relationship with unit slope and zero intercept). For model-based methods like AYB, major discrepancy between estimated and empirical qualities is indicative of poor fit of the model to experimental data. Importantly, knowledge of the accuracy of quality scores, often dependent on per-experiment differences, can be used to calibrate them to give better error probability estimates.
where the discrete range of quality values is indexed by a and a base caller assigns n a bases to quality q a , for which the empirical quality is . Notice that the second term of the numerator is negative when the quality is over-estimated and its magnitude increases exponentially with increasing error, whereas the 'bonus' for making a conservative prediction is bounded.
Assessment of quality score accuracy across all data sets for the three base callers under a variety of criteria; see text for details.
Q tot, ×109
The maximum value of assuming perfect calibration, Max, is also shown in table 5 as a measure of the maximum amount of information that could be extracted from the base calls assuming further, probably reference-based, calibration. For these scores, Bustard outperforms Ibis, with the sole exception of the first end of the B. pertussis data set, in complete contrast to their relative performance for the other criteria. While Ibis makes fewer base calling errors than Bustard, it does not do such a good job of separating the bases according to confidence and so the total information content of the calls is lower. AYB outperforms both Bustard and Ibis, particularly on the four sets of data that are similar to that which its calibration table was estimated (ϕX174 L2, ϕX174 L4, ϕX174 L6 and Ibis Test).
The 'Total Quality', Q tot, is a measure of the total information content of the entire set of called bases and is equal to Max multiplied by the number of mapped bases. As it is a sum rather than an average, the Q tot criterion allows base callers to compensate in bulk for producing low quality calls. The values of Q tot are also shown in table 5 and have trends that are broadly similar to those for Max, with Ibis generally having a slightly worse scores than Bustard and AYB performing best for all data sets. A notable difference between the Max and Q tot results is that Ibis has a much better Q tot score than Bustard on both ends of the B. pertussis data, reflecting the increased proportion of mapped reads (see table 2). Ibis has an important advantage over the other base-callers on the RMS and S p criteria as it was always trained on representative data, and so effectively recalibrated using a reference each time. This advantage is evidenced by the poor performance of Bustard and the improvement in AYB's performance when run-specific calibration is used. These results should not be taken as a indication of the superiority of the quality scores or calibration of Ibis since whatever data was used to train Ibis could also have been used to come up with a good calibration for either AYB or Bustard. While the calibration of neither Bustard nor AYB is exceptional, their performance on the Max and Q tot criteria, as well as the scores for AYB on the ϕX174 sets of data and after recalibration for other data sets, suggest that bases from both callers can be recalibrated to produce qualities whose informativeness equals or exceeds those from Ibis.
A particular focus when developing AYB was to make the algorithms robust to problems that might arise during normal use, so it can be used confidently in cases where other base-callers require manual intervention to get the best results. The B. pertussis example was presented as such a case; another is data produced using the TraDIS technique  where the first few cycles of every cluster consist of known identical sequence, causing algorithms that estimate cross-talk from a single early cycle to fail.
The statistical model underlying AYB has several weaknesses that could be addressed in future work. The model assumes that the descriptive parameters of the sequence process are constant across a tile but this is only going to be approximately true in practice: differences in illumination (e.g. mode scrambler problems) and the relative intensities of the two lasers will affect the cross-talk and background noise; also the expected amount of phasing might be affected by fluctuations in the chemistry. The phasing matrix represents an average over many clusters and the actual amount of phasing at a particular cluster is subject to stochastic variation. The fewer molecules contained in the cluster, the further from the average it is likely to deviate and this can lead to counter-intuitive consequences as a small cluster that has, by chance, undergone little phasing may fit the average model as poorly as one that has undergone a lot of phasing -- clusters could then be penalised despite giving clear signal.
Sequence-like errors, for example mutations introduced during sample preparation, short fragments ligating together or adapter sequence, are essentially invisible to any base-caller and render it impossible to call the original sequence accurately. Other sources of error may not appear sequence-like: for example, microscopic particles of dust can get entangled in a cluster and produce bright artefacts for one or more cycles. Since very bright peaks deviate from the average brightness of the read, AYB penalises these calls heavily and they rarely contribute to the higher quality base-calls; they also reduce the quality of the surrounding calls due to over-correction for pre- and post-phasing. Ideally over-bright peaks would be removed prior to analysis and treated as missing data, with the actual intensity and base-call imputed from the remaining three intensities and the effect the position has on the neighbouring cycles through the phasing correction. A similar idea could be used to deal with clusters where intensities are missing (i.e. unrecorded, perhaps due to image registration problems) for some cycles, producing low quality calls rather than as at present where they are treated as a cycle with four exactly zero intensities.
A final issue that AYB fails to account for is that of heterogeneous clusters of sequence, a common cause of which is two clusters merging into each other during the amplification step, since there is an implicit assumption that each cluster only contains fragments from one particular sequence. The intensities from such clusters appear to be extremely noisy, far above the stochastic background, and AYB's criteria to assess model fit are misled since the effects of both constituent sequences need to be be removed to get the residual noise. Failure to do so means that the calls from the strongest sequence get penalised for badly fitting the model; in particular, cycles where the two constituent clusters have the same base appear much brighter than expected given the intensities from other positions and are thus penalised despite the fact we should be more confident about these calls. In principle heterogeneous sequence could be explicitly estimated for each cluster, the relative brightness being used to separate contributions, but this will result in a loss of power in the majority of cases where the cluster is homogeneous and may not result in high-quality calls otherwise.
Despite being noticeably better than those produced by Bustard, the uncalibrated quality values for AYB are worse than might be desired. Improvements to these could be the subject of further work. Run-specific recalibration leads to significant improvement, with AYB outperforming Ibis. As with Bustard, and indeed all other base callers, we recommend that the qualities produced by AYB should be recalibrated whenever possible as part of the analysis pipeline. The superiority of AYB on both the S p Max and Q tot criteria suggest that such recalibration would be fruitful. Where a good quality reference is not available, recalibration based on reads from a spiked-in known genome is a promising approach that could be taken advantage of. Such spike-in data may also help with convergence and improve the estimate of the interaction matrix since it provides a set of reads whose sequence does not have to be estimated.
The speed at which tiles can be analysed is extremely important given the vast amount of data produced by current and future platforms. AYB is much quicker than many competitive base callers, taking only a few minutes to analyse each tile on ordinary computing hardware and of comparable speed to Ibis, the most accurate alternative base caller (table 3).
Even AYB's speed could be prohibitive if computing resources are limited. There are, however, avenues to increase the speed of AYB with possible trade-offs against accuracy. As noted previously, AYB assumes that the sequencing process is constant within a tile and this assumption could be to strengthened to assuming it is constant across multiple tiles or across lanes, a similar assumption to that which Bustard and other base calling programs implicitly make when they train or estimate parameters on a subset of the data. The parameters describing the sequencing process could be estimated from a subset of data and then held fixed so AYB need only perform a base calling step for the majority of the clusters with a considerable reduction in processing time. The major bottleneck for AYB is processing the raw intensities for each cluster, a step that is repeated every iteration and is quadratic in the number of cycles, and speeding up this calculation would greatly accelerate the algorithm. One potential approach would be to assume that the interaction matrix A (see Methods) is sparse, the intensities at one cycle only depending on the sequence at nearby cycles for example.
AYB is more accurate than other methods of base-calling. In comparison with the leading competitor, Ibis , it generally gives similar or improved performance in the number of mapped reads, and in our tests it always performed considerably better in the number of perfect (error-free) reads (table 2) and almost always achieves a lower per-mapped-base error rate (Figure 1). As the yield from sequencing machines increases, speed of analysis becomes important and our base-calling method offers a unique combination of speed and accuracy. Such a combination is ideally suited for use with more recent 'personal' platforms, such as Illumina's MiSeq , which are aimed at smaller institutes and research groups who will be interested in a diverse range of organisms for which a good reference genome is not likely to be available. In addition to its speed and accuracy, AYB has two other desirable properties. First, it does not require training data so calls can be made where a reference sequence is unknown. Second, it uses robust statistical methods to limit undesirable consequences of gross errors in a few clusters.
The AYB base-calling software is written in C and available under the GPL v. 3 licence from http://www.ebi.ac.uk/goldman-srv/AYB/ A set of utilities for extracting and manipulating CIF format intensity data files, under the same licence as AYB, is available from http://www.ebi.ac.uk/goldman-srv/ciftools/
Materials and methods
The two major differences between AYB and other base-callers are its empirical model of the sequencing process, potentially allowing the intensities at a given cycle to depend on the entire sequence rather than just a few neighbouring cycles, and its focus on robust algorithms so that sensible base calls are still made even when problems have occurred during a run. Here we describe the underlying statistical model used by AYB, the method of estimation and the techniques used to make the procedure robust.
The foundation of AYB is a mechanistic model of the sequencing process, relating what is observed at each cycle to the underlying sequence of nucleotides. Clusters are analysed in groups, the natural such group being a tile, with the interaction between cycles assumed to be constant and common to all clusters within each group. Other parameters such as the luminescence and the sequence are specific to each cluster. We first describe a simple model of how the observed intensities might be related to the underlying sequence and then show how AYB generalises it.
Each cluster (indexed by i) is considered to contain homogeneous sequence, represented by the base × position matrix S i whose (b, j) entry is one if the base at the j th position of the sequence is base 'b' or zero otherwise. Each column of S i therefore contains exactly one non-zero entry. The amount of light emitted by a cluster in a given cycle is proportional to the number of FLNs bound to the cluster, which in turn is proportional to the number of molecules in the cluster; this cluster-specific scaling is represented by the scalar λ i , referred to as the luminescence since it also incorporates a factor representing the intensity of light incident on the cluster and implicitly models variation in incident radiation across the slide.
Due to phasing, the molecules within a cluster lose synchronicity with each other and the relationship between position and cycle becomes blurred; the procession from one cycle to the next of an average cluster is described by the position × cycle phasing matrix P. Each column of P corresponds to one cycle and describes the distribution of sequence positions where the FLNs bind, so the (j, k) entry is the relative proportion of FLNs bound to position j of the sequence on cycle k of the sequencing process. As sequencing progresses, the signal decreases as molecules randomly become inactive and stop contributing (dimming) and this is incorporated into P by scaling its columns so each sums to the proportion of molecules in the cluster expected to be still active. An ideal P would have ones down its leading diagonal with all other elements being zero; a good P will be dominated by its diagonal and each column sum will be close to one. Elements of P are non-negative, and its column sums are ≤ 1.
Finally, the emissions from each cluster are observed via the four channels and the cross-talk, the relationship between fluorophore emission and what is observed in each channel, is represented by a channel × base (4 × 4) matrix M. Column b of M describes the strength of signal in each of the four channels for a unit emission of the FLN b. In principle the cross-talk is determined by the physics of system, and so it is assumed to be the same for all cycles.
where A is the interaction matrix, a (channel × cycle) × (base × position) matrix describing the effect that specific bases at each cycle have on the intensities for all cycles. Note that A allows for the cross-talk to vary between cycles and for the rate of phasing to depend on previous bases.
The statistical model described by equation 3 could be fitted to the raw intensity data using a variety of criteria (maximum likelihood, Bayesian techniques, etc.) but we chose a least squares criterion using an iterative approach. The major reasons for the use of least squares are that analytic solutions exist for many of the steps of the iteration, making it computationally efficient, and that the simple Iteratively reWeighted Least Squares (IWLS) technique can be used to fit the model in a manner robust to contamination . The IWLS approach is similar to Ordinary Least Squares (OLS), seeking to minimise the sum of squared errors over all the clusters, but the squared error for each cluster is weighted and the algorithm proceeds iteratively with the weights being updated between iterations; each iteration is equivalent to the Weighted Least Squares (WLS) criterion, which has an analytic solution. The weights are defined by a function of how well each cluster fits the model relative to the other clusters, so badly fitting clusters (high residual error) get progressively down-weighted. AYB uses the Cauchy function for weighting but many alternatives have been described and are summarised in the subsection on 'Robust Estimation' in Numerical Recipes .
Initialise estimates; set all weights to one.
Estimate interaction A and systematic noiseN.
Estimate cluster-specific luminescence λ i .
Call bases for each cluster, giving sequence S i .
Update weights for all clusters.
Iterate steps 2-5 to refine estimates.
Assess quality of calls.
Initialising to good values greatly helps the speed of the AYB algorithm, reducing the number of iterations needed until a good solution is found. A good starting value may be available from previous analyses using the same machine and protocol but, by default, AYB uses the more crude approach of ignoring phasing and dimming and assuming that the cross-talk is the same at all cycles: if M is a cross-talk matrix then the initial estimate of the interaction matrix is A 0 = ID ⊗ M where I D is the identity matrix of dimension cycle × position. An initial cross-talk matrix M can be found from the intensities of an early cycle of the run , making the implicit assumption that phasing does not contribute a significant amount to these observed intensities, but, since cross-talk is primarily determined by physics and has a similar form on different runs and machines, AYB instead initialises M to a fixed good value. Systematic noise is initially assumed to be absent.
Setting A to A 0 and solving equation 3 for λ i S i , assuming that the systematic and random noise (N and ∈ i ) are zero, gives a set of corrected intensities from which bases and luminescence can be estimated. The initial estimate of the base at each position is that which has the greatest intensity, and the luminescence of the cluster is the mean of the intensities of the called bases.
Estimation of the interaction matrix and noise
Equation 3, relating the sequence to the expected intensities, is linear but standard linear regression techniques produce unstable estimates for the interaction matrix; to see why, notice that any permutation of the columns of A and rows of S i leaves the intensities unchanged and so, when estimating the interaction is iterated with base-calling, the solution can jump between permutations. We use generalized Tikhonov regularisation  in a weighted least squares solution for A and N to favour there being no permutation of A and S i .
where ID is an identity matrix of the appropriate dimension, p is a constant specifying the strength of regularisation, and with 0 being a vector consisting of zeros. For convenience, the solution is regularised towards the value used to initialise the algorithm, although this is not a requirement and other choices may be more desirable.
There is an arbitrary scaling factor implicit in equation 3, corresponding to the scale that the luminescence is measured on. If all elements of interaction matrix are doubled and every λ i is halved, then the expected intensities are unchanged. This ambiguity is resolved by scaling the interaction matrix so that its determinant is one.
Estimation of luminescence
which is an ordinary least squares estimate since the weights used for the estimation of the interaction matrix are cluster-specific and so cancel.
As well as being linear in the interaction matrix, the observed intensities in equation 3 are also a linear function of the sequence. As described, equation 3 assumes that each element of the random error is independent and identically distributed (IID) but this is not found to be the case in real data (results not shown). The violation of the IID assumption is not a problem when estimating the interaction matrix, since these estimates are produced from a large number of independent clusters and so the random error is small, but is much more significant when trying to estimate the sequence since there are many fewer, dependent, observations. Forcing the IID assumption onto the random noise produces poor base calls (results not shown) and so correlation between the elements of ε i must be taken into account and the sequence that minimises the generalised least squared error must be found.
and the sequence that minimises the generalised least square error of equation 3 also minimises the generalised least square error of the same relation written in terms of the corrected intensities C i . This latter formulation is more convenient to work with. Finding the minimum generalised least square is a type of constrained binary quadratic programming problem and so difficult to solve exactly. Instead of solving directly, we make the additional assumption that each read position only depends on its immediate neighbours and so most positions are conditionally independent of each other. This dependence structure requires that the inverse of the covariance matrix is block tridiagonal; that is, it consists of a grid of base × base matrices and this grid is tridiagonal. The maximum likelihood estimate of the required covariance matrix is found by numerical optimisation (conjugate gradient algorithm) of the log-likelihood function parametrised in terms of the Cholesky factorisation of the matrix; full details are contained in Additional file 1.
The structure of the inverse covariance matrix means that the log-likelihood for a cluster i having the sequence s 0, s 1, . . . , s n can be written as for suitably chosen tensor a f,g,h and a constant k, and so is a one-dimensional Gibbs field. The classic Viterbi and Forward/Backward dynamic programming algorithms can be used to find the most likely sequence or the posterior distribution of bases at each position.
The weighting of the clusters plays an important part in making the AYB method robust to contamination and other misleading observations, reducing their influence on the parameter estimates. The weight for each cluster is calculated, after all model parameters have been fitted, using the Cauchy function so w i = 1/(1 + L i /2μ) where L i is the least square error for cluster i and μ is a measure of the central trend of the L i (their mean, for example). In contrast to OLS, where every cluster would receive a weight of one, this weighting function means that only perfect observations, those with a least square error of zero, receive full weight whereas worse-fitting observations receive progressively lower weights.
Iteration and termination
As the luminescence and individual bases are estimated (steps 3-4 above) using a different criterion to that used to estimate the interaction matrix (step 2), the least squares error is not guaranteed to decrease when the parameter estimation and base calling steps are iterated. Theoretically this could lead to problems with convergence but this was not found to be the case, a small number of cycles sufficing to get good estimates. Numerical experiments suggest that three to five iterations are sufficient, with little change in accuracy for addition iterations (results not shown).
Assessment of quality
where is the (scaled) probability density of the observed intensities, L i ; xj is the least square error for cluster i given that the base at position j is x, and π x is the prior probability of base x. The required least squared errors can be calculated for all bases and positions simultaneously using a Forwards/Backwards modification of the Viterbi algorithm used to find the best sequence of bases (step 4 above). The particular form of f i ; bj comes from assuming that the random error has an elliptical distribution defined by the Cauchy function, in keeping with our choice of weighting function for the IWLS estimation. The parameter n* should normally be equal to the dimension of the elliptical distribution but, instead, we use the median of the observed least squared errors as this helps to correct for skews in the distribution.
Reads can also be wrong for reasons other than base calling error, a good example of this being polymerase errors during sample preparation. The net effect of these 'generalised' errors is to bound the maximum possible quality of a call and they are incorporated into AYB's quality scores as a constant probability of error independent of the probability that the base was called incorrectly. The final probability that a base is correct, incorporating the notion of generalised error, is where ε is the constant probability of a generalised error. The corresponding quality score is . Despite the methods incorporated into AYB being robust and attempts being made to compensate for effects of unusual clusters, the quality scores produced may not be accurate because of differences between AYB's assumptions and how the sequencing machines actually operate: not every source of error can be incorporated into the model and the various distributional assumptions made can only be approximate. To improve concordance, quality scores may be calibrated to real data  using some form of table look-up or calibration function. There are many good methods to calibrate quality scores [21, 25, 31] but, for the purposes of the comparisons in this paper, we use a simple linear calibration function for AYB and note that it could be improved upon.
where the base × base × base table α and the constant β are chosen to agree with representative data from a real sequencing run. We expect these parameters to vary with difference machines, chemistries and experimental protocols, and typical values based on the sets of data analysed within this paper are distributed with the AYB software. We provide a tool to produce bespoke constants given a set of mapped data; however, we reiterate our recommendation that quality scores should be recalibrated using a reference whenever possible.
The development of AYB has benefited from the help of many people, for both advice and feedback about performance 'in the wild'. In particular we would like to thank Ewan Birney (EMBL-EBI); Jonathon Blake (EMBL); Gordon Brown (CRI); Kevin Howe (formerly CRI); Tony Cox, Klaus Maisinger, Lisa Murray (Illumina Inc.); Christophe Dessimoz and Christian Ledergerber (ETH Zürich); Richard Durbin, Tom Skelly (Sanger Institute); Nava Whiteford (formerly Sanger Institute); and David van Heel (Blizard Institute).
The 76 cycle ϕX174 and 76 cycle B. pertussis data sets were provided by the Sanger Institute. Other data as indicated in the text.
The development and maintenance of AYB was supported by a Wellcome Trust Technology Development grant WT088151MA.
- Varela I, Klijn C, Stephens PJ, Mudie LJ, Stebbings L, Galappaththige D, van der Gulden H, Schut E, Klarenbeek S, Campbell PJ, Wessels LFA, Stratton MR, Jonkers J, Futreal PA, Adams DJ: Somatic structural rearrangements in genetically engineered mouse mammary tumors. Genome Biology. 2010, 11: R100-10.1186/gb-2010-11-10-r100.PubMed CentralView ArticlePubMedGoogle Scholar
- Fuller CW, Middendorf LR, Benner SA, Church GM, Harris T, Huang X, Jovanovich SB, Nelson JR, Schloss JA, Schwartz DC, Vezenov DV: The challenges of sequencing by synthesis. Nature Biotechnology. 2009, 27 (11): 1013-1023. 10.1038/nbt.1585.View ArticlePubMedGoogle Scholar
- Metzker ML: Sequencing technologies -- the next generation. Nature Reviews Genetics. 2010, 11: 31-46. 10.1038/nrg2626.View ArticlePubMedGoogle Scholar
- Ledergerber C, Dessimoz C: Base-calling for next-generation sequencing platforms. Briefings in Bioinformatics. 2011, 12: 489-497. 10.1093/bib/bbq077.PubMed CentralView ArticlePubMedGoogle Scholar
- Li L, Speed T: An estimate of the crosstalk matrix in four-dye fluorescence-based DNA sequencing. Electrophoresis. 1998, 20: 1433-1442.View ArticleGoogle Scholar
- Kao WC, Stevens K, Song YS: BayesCall: a model-based basecalling algorithm for high-throughput short-read sequencing. Genome Research. 2009, 19: 1884-1895. 10.1101/gr.095299.109.PubMed CentralView ArticlePubMedGoogle Scholar
- Erlich Y, Mitra PP, de la Bastide M, McCombie WR, Hannon GJ: Alta-Cyclic: a self-optimizing base caller for next-generation sequencing. Nature Methods. 2008, 5: 679-682. 10.1038/nmeth.1230.PubMed CentralView ArticlePubMedGoogle Scholar
- Rougemont J, Amzallag A, Iseli C, Farinelli L, Xenarios I, Naef F: Probabilistic base calling of Solexa sequencing data. BMC Bioinformatics. 2008, 9: 431-10.1186/1471-2105-9-431.PubMed CentralView ArticlePubMedGoogle Scholar
- Kircher M, Stenzel U, Kelso J: Improved base calling for the Illumina Genome Analyzer using machine learning strategies. Genome Biology. 2009, 10: R83-10.1186/gb-2009-10-8-r83.PubMed CentralView ArticlePubMedGoogle Scholar
- Bravo HC, Irizarry RA: Model-based quality assessment and base-calling for second-generation sequencing data. Biometrics. 2010, 66: 665-674. 10.1111/j.1541-0420.2009.01353.x.PubMed CentralView ArticlePubMedGoogle Scholar
- Dohm JC, Lottaz C, Borodina T, Himmelbauer H: Substantial biases in ultra-short read sets from high-throughput DNA sequencing. Nucleic Acids Research. 2008, 36: e105-10.1093/nar/gkn425.PubMed CentralView ArticlePubMedGoogle Scholar
- Nakamura K, Oshima T, Morimoto T, Ikeda S, Yoshikawa H, Shiwa Y, Ishikawa S, Linak M, Hirai A, Takahashi H, Altaf-Ul-Amin M, Ogasawara N, Kanaya S: Sequence-specific error profile of Illumina sequencers. Nucleic Acids Research. 2011, 39: e90-10.1093/nar/gkr344.PubMed CentralView ArticlePubMedGoogle Scholar
- Kao WC, Song YS: naiveBayesCall: an efficient model-based base-calling algorithm for high-throughput sequencing. Proc 14th Annual Intl Conf on Research in Computational Molecular Biology. 2010, 233-247.View ArticleGoogle Scholar
- Li R, Li Y, Kristiansen K, Wang J: SOAP: short oligonucleotide alignment program. Bioinformatics. 2008, 24: 713-714. 10.1093/bioinformatics/btn025.View ArticlePubMedGoogle Scholar
- AYB website. [Accessed: 14 Jan. 2012], [http://www.ebi.ac.uk/goldman-srv/ayb/]
- Li H, Durbin R: Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009, 25: 1754-1760. 10.1093/bioinformatics/btp324.PubMed CentralView ArticlePubMedGoogle Scholar
- Ibis website. [Accessed: 12 Jan. 2012], [http://bioinf.eva.mpg.de/Ibis/]
- Zerbino DR, Birney E: Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Research. 2008, 18: 821-829. 10.1101/gr.074492.107.PubMed CentralView ArticlePubMedGoogle Scholar
- The 1000 Genomes Project Consortium: A map of human genome variation from population-scale sequencing. Nature. 2010, 467: 1061-1073. 10.1038/nature09534.PubMed CentralView ArticleGoogle Scholar
- Best Practice Variant Detection with the GATK v3. [http://www.broadinstitute.org/gsa/wiki/index.php/Best Practice Variant Detection with the GATK v3]. [Accessed: 12 Jan. 2012]
- DePristo M, Banks E, Poplin R, Garimella K, Maguire J, Hartl C, Philippakis A, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell T, Kernytsky A, Sivachenko A, Cibulskis K, Gabriel S, Altshuler D, Daly M: A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics. 2011, 43: 491-498. 10.1038/ng.806.PubMed CentralView ArticlePubMedGoogle Scholar
- Drmanac R, Sparks AB, Callow MJ, Halpern AL, Burns NL, Kermani BG, Carnevali P, Nazarenko I, Nilsen GB, Yeung G, Dahl F, Fernandez A, Staker B, Pant KP, Baccash J, Borcherding AP, Brownley A, Cedeno R, Chen L, Chernikoff D, Cheung A, Chirita R, Curson B, Ebert JC, Hacker CR, Hartlage R, Hauser B, Huang S, Jiang Y, Karpinchyk V, Koenig M, Kong C, Landers T, Le C, Liu J, McBride CE, Morenzoni M, Morey RE, Mutch K, Perazich H, Perry K, Peters BA, Peterson J, Pethiyagoda CL, Pothuraju K, Richter C, Rosenbaum AM, Roy S, Shafto J, Sharanhovich U, Shannon KW, Sheppy CG, Sun M, Thakuria JV, Tran A, Vu D, Zaranek AW, Wu X, Drmanac S, Oliphant AR, Banyai WC, Martin B, Ballinger DG, Church GM, Reid CA: Human genome sequencing using unchained base reads on self-assembling DNA nanoarrays. Science. 2010, 327: 78-81. 10.1126/science.1181498.View ArticlePubMedGoogle Scholar
- Complete Genomics pubic data repository. http://www.completegenomics.com/sequence-data/download-data/]. [Accessed: 12 Jan. 2012]
- Simon-Sanchez J, Scholz S, Fung HC, Matarin M, Hernandez D, Gibbs JR, Britton A, de Vrieze FW, Peckham E, Gwinn-Hardy K, Crawley A, Keen JC, Nash J, Borgaonkar D, Hardy J, Singleton A: Genome-wide SNP assay reveals structural genomic variation, extended homozygosity and cell-line induced alterations in normal individuals. Human Molecular Genetics. 2007, 16: 1-14. 10.1093/hmg/ddm004.View ArticlePubMedGoogle Scholar
- Abnizova I, Skelly T, Naumenko F, Whiteford N, Brown C, Cox T: Statistical comparison of methods to estimate the error probability in short-read Illumina sequencing. J Bioinform Comput Biol. 2010, 8: 579-591. 10.1142/S021972001000463X.View ArticlePubMedGoogle Scholar
- Langridge GC, Phan MD, Turner DJ, Perkins TT, Parts L, Haase J, Charles I, Maskell DJ, Peters SE, Dougan G, Wain J, Parkhill J, Turner AK: Simultaneous assay of every Salmonella typhi gene using one million transposon mutants. Genome Research. 2009, 19: 2308-2316. 10.1101/gr.097097.109.PubMed CentralView ArticlePubMedGoogle Scholar
- MiSeq product brochure. [http://www.illumina.com/documents//products/brochures/MiSeq_Brochure.pdf]. [Accessed: 12 Jan. 2012]
- Agresti A: Categorical Data Analysis. 2002, John Wiley & Sons, secondView ArticleGoogle Scholar
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C. 1992, Cambridge University Press, secondGoogle Scholar
- Tikhonov AN: Solution of incorrectly formulated problems and the regularization method. Soviet Math Dokl. 1963, 4: 1035-1038.Google Scholar
- Ewing B, Green P: Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Research. 1998, 8: 186-194.View ArticlePubMedGoogle Scholar
- Wilson EB: Probable inference, the law of succession, and statistical inference. J Amer Stat Assoc. 1927, 22: 209-212.View ArticleGoogle Scholar
This article is published under license to 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.