Fast and accurate short read alignment with Burrows-Wheeler transform

Abstract

Motivation: The enormous amount of short reads generated by the new DNA sequencing technologies call for the development of fast and accurate read alignment programs. A first generation of hash table-based methods has been developed, including MAQ, which is accurate, feature rich and fast enough to align short reads from a single individual. However, MAQ does not support gapped alignment for single-end reads, which makes it unsuitable for alignment of longer reads where indels may occur frequently. The speed of MAQ is also a concern when the alignment is scaled up to the resequencing of hundreds of individuals.

Results: We implemented Burrows-Wheeler Alignment tool (BWA), a new read alignment package that is based on backward search with Burrows-Wheeler Transform (BWT), to efficiently align short sequencing reads against a large reference sequence such as the human genome, allowing mismatches and gaps. BWA supports both base space reads, e.g. from Illumina sequencing machines, and color space reads from AB SOLiD machines. Evaluations on both simulated and real data suggest that BWA is approximately 10-20x faster than MAQ, while achieving similar accuracy. In addition, BWA outputs alignment in the new standard SAM (Sequence Alignment/Map) format. Variant calling and other downstream analyses after the alignment can be achieved with the open source SAMtools software package.

Availability: http://maq.sourceforge.net.

Figures

Fig. 1.

Prefix trie of string ‘GOOGOL’. Symbol ∧ marks the start of the string. The two numbers in a node give the SA interval of the string represented by the node (see Section 2.3). The dashed line shows the route of the brute-force search for a query string ‘LOL’, allowing at most one mismatch. Edge labels in squares mark the mismatches to the query in searching. The only hit is the bold node [1, 1] which represents string ‘GOL’.

Fig. 2.

Constructing suffix array and BWT string for X=googol$. String X is circulated to generate seven strings, which are then lexicographically sorted. After sorting, the positions of the first symbols form the suffix array (6, 3, 0, 5, 2, 4, 1) and the concatenation of the last symbols of the circulated strings gives the BWT string lo$oogg.

Fig. 3.

Algorithm for inexact search of SA intervals of substrings that match W. Reference X is $ terminated, while W is A/C/G/T terminated. Procedure InexactSearch(W, z) returns the SA intervals of substrings that match W with no more than z differences (mismatches or gaps); InexRecur(W, i, z, k, l) recursively calculates the SA intervals of substrings that match W[0, i] with no more than z differences on the condition that suffix Wi+1 matches interval [k, l]. Lines started with asterisk are for insertions to and deletions from X, respectively. D(i) is the lower bound of the number of differences in string W[0, i].

Fig. 4.

Equivalent algorithm to calculate D(i).