Authors: Nawar Malhis, Yaron S. N. Butterfield, Martin Ester, Steven J. M. Jones
Categories: Original Papers
Source: Bioinformatics
Motivation: A plethora of alignment tools have been created that are designed to best fit different types of alignment conditions. While some of these are made for aligning Illumina Sequence Analyzer reads, none of these are fully utilizing its probability (prb) output. In this article, we will introduce a new alignment approach (Slider) that reduces the alignment problem space by utilizing each read base's probabilities given in the prb files.
Results: Compared with other aligners, Slider has higher alignment accuracy and efficiency. In addition, given that Slider matches bases with probabilities other than the most probable, it significantly reduces the percentage of base mismatches. The result is that its SNP predictions are more accurate than other SNP prediction approaches used today that start from the most probable sequence, including those using base quality.
Contact: nmalhis@bcgsc.ca
**Supplementary information and ** http://www.bcgsc.ca/platform/bioinfo/software/slider
Novel parallel sequencing technologies including sequencing by synthesis (SBS) such as Illumina Sequence Analyzer (for a review, refer to Holt and Jones, 2008), 454 Life Sciences sequencing and Applied Biosystems SOLiD sequencing, coupled with ever decreasing costs, has provided new opportunities to study genomes. This includes the ability to identify sequence aberrations and chromosomal abnormalities, such as single nucleotide polymorphisms (SNPs), insertions, deletions and inversions, on a genome-wide scale. The amount of sequence produced has also provided the computational challenge of developing accurate and reliable software approaches for sequence alignment that can complete withina useful time frame using affordable computer hardware resources.
Here, we focus on the short sequence reads generated by the Illumina Sequence Analyzer G1 and G2 platforms. Each experiment run typically produces in excess of 1 Gb of nucleotide sequence as short reads (ranging from 27 to 42 nt in length). The Illumina sequencing by synthesis technology works by first randomly fragmenting DNA of interest, ligating adapters to each end, and then annealing the DNA fragments to a flow cell surface coated with complementary oligonucleotides. The DNA fragments are then amplified at millions of unique locations across the flow cell surface resulting in a random array of clonal clusters of DNA molecules. Sequencing of the amplified DNA occurs in cycles. At each cycle, all four fluorescently labelled reversibly terminated nucleotides are added to the reaction allowing the single correct nucleotide to be incorporated by the polymerase at each cluster. A subsequent laser excitation allows for four images to be generated each of which capture a specific excited wavelength from each of the four fluorescent labels corresponding to a single base type. Comparing the colour signal intensity of each base for each cluster, probabilities corresponding to the most likely base incorporation are generated and reported in a prb file. The most probable identity of that base is called according to these probabilities resulting in a final read sequence which is reported in a seq file. It is important to note that with each subsequent cycle, the sequencing error rate increases primarily due to the fact that the chemical reactions at each cycle do not complete to 100%, resulting in molecules within in each cluster becoming out of synchronization. This phenomenon, known as ‘dephasing’ increases the background noise and reduces the accuracy of each base called as the sequencing proceeds.
An important task is to align sequence reads to the appropriate reference genome for the detection of genetic variants. However, the alignment of these short reads presents a number of challenges due to their short lengths and error rate. Traditional methods of sequence alignment have worked well with sequences that are long enough to provide high location specificity whilst tolerating a number of mismatching bases (sequencing errors and/or biological SNPs). For short reads a correct interpretation of the quality of every base call is essential in the correct alignment of reads to the reference and the accurate identification of SNPs. For example, one base-calling error can result in a no alignment match to the reference genome, a match reporting a false SNP, or a misaligned match to a different region of the genome. The quality of the alignment and SNP prediction decreases as the number of inaccurate base calls increase. In general, if a base in the reference sequence has a high enough coverage with a consistent difference from its reference, this base likely represents a polymorphism. This redundant coverage must be high enough to confidently discriminate true variants from sequencing errors, and to account for the fact that different copies of the same read location might be generated from biological aliases (e.g. diploid or mitochondrial chromosomes).
A plethora of alignment tools have been created employing various methods (some of which are described here). A few make use of the quality information from the Illumina prb files and these typically use the base probability values to evaluate the probability of the called base only.
In this study, we will introduce a new alignment tool, Slider. Slider is an application for the Illumina Sequence Analyzer output that uses the probability files instead of the sequence files as an input for alignment to a reference sequence or a set of reference sequences. We will show that this approach is more accurate and more efficient than the currently available software packages.
First, we introduce existing tools and their approach and then describe Slider's algorithm.
Different nucleotide alignment programs are made to best fit different alignment properties which sequence length of references and reads, quality of reads with respect to sequencing errors, sequence similarity to a reference sequence (SNPs, micro-indels), and the number of reads to be aligned with the reference. Some of the more general use aligners
Given the two input sets, the set of reads and the set of reference sequences, most of the available aligners are based on the following two
However, when the aligned sequences have specific properties, building aligners to best fit these properties will reduce the problem space and improve both the efficiency and accuracy of the alignment process. Some aligners are designed especially for aligning the output of Illumina's Sequence Analyzer. These aligners take advantage of the short read lengths for speeding up the alignment (by using hash tables for speeding up seed identification) which allows the processing of a larger number of reads in a reasonable time
As of today, these aligners still have some common disadvantages with respect to the analysis of our target
To overcome these problems, it is important to first understand the properties of the Illumina Sequence Analyzer data, and then to build an approach that maximizes the utilization of these properties in order to optimize the alignment (efficiency and accuracy) and SNP prediction accuracy.
In contrast to other aligners, the Slider algorithm uses not only the most probable base, but also all possible bases with a probability above a certain base probability threshold (baseMinPrb) provided by the Illumina probability files in order to generate all possible reads with probability above a certain read probability threshold (read_0_MinPrb). This reduces thelevel of alignment approximation needed and improves accuracy (reducing the number misaligned reads and the number of base mismatches). For its core alignment, Slider sorts all these generated reads in lexicographical order and then crosses it sequentially with a presorted table of windows of reference sequence(s) and their reverse complement(s). This approach eliminates the need for an indexed structure by replacing random I/O with sequential I/O which allows Slider to scale to large input datasets. Our test results show that Slider efficiently provides more accurate alignments and SNP calls.
Today, the use of Paired End Tag (PET) data in alignments is an important aspect of extracting more information from these reads. Slider currently processes PET data and we are now expanding on this capability. In this article, we demonstrate Slider's concept of aligning a read which is applicable for both single-end and paired-end tag data.
There are three main properties of the Illumina Sequence Analyzer output reads that are important to
The base probability values are provided in the prb files in the form of Q-values, where probabilities can then be calculated
where,
Qbs (i): the probability of the base bs, in the form of a Q-value, at location i in a read in the range of [−Qmax, Qmax]. Qmax is usually set to 40.
Pbs (i): the probability of base bs at location i in a read, is in the range of [Pmin, Pmax].
At a given sequencing cycle, if a cluster is visible in only one of the four images, then the base for that cluster is as well resolved as it can be. We call this base a ‘crisp base’. In other words, we define a base as a crisp base if three of its four prb values are equal to −Qmax, and one is equal to Qmax. A crisp read is a read where all the constituent bases are crisp.
For example, when a base is crisp, the base will have such Q and corresponding P-values:
But, if the base is not crisp, the Q-values for a base might look like
Utilizing these properties, Slider provides major improvements in efficiency and accuracy over other traditional alignment algorithms.
We present some more definitions that will be used in later sections
Definitions and
Rsrc: the actual source read of nucleotide sequence.
SZr: the size (length) of a read.
SZd: the size (length) of a single database sequence (which is a subsequence of the reference).
Rps: a probable sequence of a read given its prb line; crisp reads have only one Rps.
Rmps: the most probable sequence (MPS) of a read; this sequence is the one given in the seq files which are generally used by many aligners as the only input.
Ral: a substring of the reference sequence where an Rps is aligned by an aligner. If Ral = Rsrc, we say that the read is accurately aligned. Otherwise, the read is misaligned.
Dmps: the edit distance (or the number of different bases) between Rmps and Rsrc.
Dal: the edit distance between an aligned Rps and Ral.
Pps: the probability, ‘weight’, of an Rps.
Pmps: equals Pps when Rps is Rmps. Crisp reads have only one Rps, which is the Rmps and this value is about one.
For example, assume a read of 8 bases, with two non-crisp bases (an example of a non-crisp read):
We first convert the Q-values to
Then, we generate all possible sequences considering every possible nucleotide with a probability higher than a baseMinPrb (a 0.2% value is used by default), and assign a probability for each sequence Pps. Those reads with a cumulative probability or ‘weight’ less than the read_0_MinPrb threshold (a 0.1% value is used by default) are ignored.
We will loosely define a read specificity as the ability of accurately aligning the read to a unique location on the reference sequence after allowing a limited number of changes on the read from a set of possible changes (Sch). In general, as the size of Sch increases, the read specificity decreases.
Unlike the Roche (454) Genome Sequencer 20 where errors are dominated by missing bases (Brockman et al., 2008), errors in Illumina Sequence Analyzer, according to our experience, are mostly base miscalls caused by PCR cluster generation (especially in the early stages), cross-talk between clusters (a considerable amount of such errors being generated when the cluster density is high) and optical noise. Most of these base miscalls canbe treated according to its probabilities provided in the prb files. Each base mismatch in an accurately aligned read can be either a sequencing error or a source mutation (SNP). By using the probability files as input instead of the sequence files, Slider reduces base mismatches for non-crisp bases which eliminates most of the sequencing errors and results in more accurate SNP predictions.
Slider works in five
Given that our datasets are large (hundreds of millions of records), in order for sorting to not impose main memory limitations on Slider, we used External Merge sorting, which has a runtime complexity of O(n logn) for n records (Knuth, 1998). As each of P0_Reads and P0_Reads tables account for sequencing errors, the size of these tables vary as the sequencing quality changes. The number of records in P0_Reads is about 5–10 times the number of prb lines, and each record is 29 bytes, so, if we have an input set of 100 million prb lines, P0_reads size is about 15–30 GB. P1_Reads size is mainly determined by the number of reads that did not match the reference and are not LQ. The P1_Reads table has about 10–30 times the number of prb lines, so for a 100 million reads input set P1_Reads table is about 30–90 GB. Since mismatches in low probability sequences are less likely to account for real SNPs, by increasing the read_1_MinPrb threshold value, we can reduce the size of the P1_Reads table without significantly affecting SNPs prediction.
While steps 1, 2 and 4 are straightforward, we will now cover in more detail the alignment (3) and the SNPs prediction (5) steps.
The alignment step is a key feature of Slider. It is fast, accurate and does not have a large memory requirement (a 4 GB machine is enough for the human genome). The slider alignment is based on two steps of exact match
While some aligners match sequences with up to a small number of mis-matches (typically, this is limited to two as increasing the number of possible mismatches reduces the read specificity), if Dmps, the edit distance between Rsrc and Rmps, is greater than this value, the read will either fail to align or will misalign. Eland, for example, uses only the most probable sequence to align and others, such as RMAP (Smith et al., 2008) make use of quality information (of the most probably base) in a discreet way in determining a degenerate or ‘wild card’ position. However, because the threshold value for non-crisp bases in a read can be large with respect to the number of possible mismatches used by other aligners, this only partially solves the problem.
Fig. 1. Slider scans both the lexographically sorted reference database and the lexicographically sorted Px_Reads (s.ol0) input table once to generate all exact matches. Exact matches are stored in the sorted s.m0 table. In this example, the set of input sequences is 6 bp (SZ
r= 6), which is aligned to a reference database of 10 bp (SZd= 10) oligos created with a sliding window across the reference. Reads that match are indicated in bold and underlined with an example of a unique match indicated by a solid line and that of a multiple match with a dashed line.
The analysis of the run time for each aligner must take into consideration a large number of factors read size, read quality, reference accuracy and contamination. In general, the computational complexity for both sorting and indexing is of an order of O(n * log n). Let us have an input set of M reads, and a reference sequence of size N. Given that our interest is when both M and N are too large to fit in main memory, the input set of reads are divided into k subsets of size m = M/k, such that each of these subset of sequences and their associated parameters can fit in main memory (for most aligners, these subsets are up to 15 million reads). In comparing the computational complexity of sorting with indexing, we can see
This gives the sorting approach an advantage when the input set of reads and the reference genomesare large such as for human.
SNP prediction is another key feature of Slider. The goal is to locate places in the reference sequence that are covered with read bases that are not equal to the reference sequence base.
Given that Slider attempts to align most of the possible reads for a prb line as a zero-off (U0) match, the final number of Slider's U0 matches is higher than other aligners, and its number of one-off (U1) matches is much lower. Given that most probability values are considered in the U0 and U1 phases, no U2 alignment is needed for Slider (Table 1). As a result,there are a lower number of base mismatches, and consequently more accurate SNP predictions. This parsimonious approach by Slider means that before a base mismatch is called, the probabilities of all the bases are utilized to ensure that the base is not one of the predicted nucleotides provided by the read prb values. Given the large number of sequence reads that are generated by a machine run, it will be expected that many of the actual bases in sequence reads will not be thosereported as the most probable.
A reference base can be covered by an aligned read in one of three different
Mismatch bases are predicted to be SNPs based on two main threshold
Other factors that affect the accuracy of SNP prediction
Fig. 2. Probability that a given base mismatch is a true SNP as a function of the read sequence weight.
Three counters for each base in the reference sequence are cmCount, mCount and nmCountwhich hold the number of crisp matches, non-crisp matches and non-matches, respectively. Based on these three counts, aligned reads' weight, and reads' base probability, three scores that reflect the likelihood of a SNP at each location in the reference are
SNPs are called when the wSNP% and the number of mismatch bases that cover a specific location in the reference sequence are higher than some user provided thresholds.
We tested Slider's alignments for both efficiency and accuracy and compared it with competing aligners. In addition to that, SNP prediction of Slider was tested for accuracy using Illumina sequence derived from high-coverage bacterial artificial chromosomes (BACs).
Three different datasets were used to evaluate Slider's performance and compare with other
Two aligners were chosen for comparing
We aligned a set of 12 flowcell lanes, 82.6 million reads in total, to the human reference genome using Eland, RMAP and Slider. Slider aligned the reads for all lanes in 28 h and 30 min with <2 GB of memory (using the –Xmx2G option of the Java VM). This is about two-thirds of the time that Eland took for the same job. RMAP by far took a lot longer to align this dataset (Table 2).
We should mention here that due to their memory requirements, there is a maximum number of reads that can be aligned simultaneously with Eland and RMAP. Therefore, results of aligning for more than one lane are accomplished by aligning to each lane separately and joining their output. Slider does not have such a limitation because it keeps both the input set of reads and the reference sequence stored in files and reads these files sequentially. This reduces Slider memory requirements to a small constant regardless of the data size. Given today's large datasets—such asthe Yoruba dataset that has been sequenced by Illumina and consists of 606 lanes (303 PET lanes), with a total of more than 4 billion reads (ftp.ncbi.nih.gov/pub/TraceDB/ShortRead/SRA000271) and those datasets in the 1000 Genomes project (www.1000genomes.org)—Slider's scalability makes it an attractive choice.
Starting with the most probable sequence of an edit distance Dmps from its actual source Rsrc, most aligners use approximation by introducing all possible x mutations trying to align this sequence to the reference. If the sequence did align with x<Dmps, then it will be aligned inaccurately, and x false mismatches are generated. Also, if the location where it was misaligned did have a SNP, this SNP will be covered with a match.
For measuring the alignment accuracy, we aligned a set of reads against both its reference BAC, RefBAC and an extra sequence from a different source, RefEX. Since the RefBAC is the accurate reference for the reads, most reads will be mapped to it; however, if RefEX is large enough, it will have two effects on the
This will make these reads less useful by moving them from the category of unique matches to the category of multiple matches.
As measuring the effect of some of these factors is a complicated issue and beyond the scope of this work, alignments of variable read and reference lengths were used in comparing reads misalignments of Eland, RMAP and Slider.
Using the same set of reads, and the same set of references (RefBAC and RefEX), we can compare the accuracy of different aligners by
In a hypothetical ideal case where the target sequence matches exactly the reference sequence, and read sequencing is error free, we will not have any misaligned reads; Pmis=0, and Puq is maximized (multiple matches are only resulted from reference similarities). In reality, the target sequence does not always match the reference sequence and there are sequencing errors. In order to perform the alignment, approximation is required. However, the less approximation we use, the less number of reads will be misaligned but the total number of accurately aligned reads will be reduced also. Different aligners use different approaches for approximation which result in a different number of reads that align and different values for Pmis, and Puq.
With the use of prb files as input and considering all base probability information, Slider is able to significantly reduce approximation compared with other aligners. The result is that less reads are misaligned.
To evaluate Slider's alignment accuracy and compare it with other aligners, alignments were done to their BAC reference sequence and three extra reference sequences RefEX:
These RefEX sequences are used to measure the level of false alignments for Eland, RMAP and Slider as a function to the RefEX length. Full results are available in the Supplementary Material on the Slider website. In Table 3, values for Pmis and Puq for the CT302 BAC with SZr in {27, 32, 36} aligned against its RefBAC and human genome excluding chromosome 6 (RefEX) shows that the percentage of reads misaligned Pmis for each of Eland and RMAP is more than double that of Slider.
As a result of its low mismatch bases (Section 3.2 and Table 1), Slider is able to give reasonably accurate SNP prediction at low coverage. We measured Slider's SNP prediction accuracy by first taking advantage of the high coverage of CT302 and T02. We generated a list of high-confidence SNPs for each dataset which are easily identified given their high coverage. These lists are then used to evaluate the accuracy of Slider's SNPs prediction by computing the probability of a single base mismatch to reflect a true SNP as a function of its read weight, Figure 2. We can see that mismatches that are in reads with higher weight are more likely to represent actual SNPs. Many factors can affect the probability of a base mismatch to reflect a true SNP including the reference sequence complexity, the length of the read and reference, and the accuracy of the prb values. A more complete statistical analysis of SNP prediction accuracy that model these factors is beyond the scope of this article, however, the higher the coverage, the more confident the SNPs call is.
There are three major advantages of using Slider for the alignment of sequences generated by the Illumina Genome Analyzer and the detection of SNPs based on these
S.J.M.J. is a senior scholar of the Michael Smith Foundation for Health Research. We thank Anthony Fejes for introducing the initial alignment problem.
Funding: IBM Canada Ltd. (in part).
Conflict of Interest: none declared.