The Burrows–Wheeler transform

The phage λ DNA sold to laboratories is cut with the enzyme HindIII into eight pieces of known length, which serve as a ruler for the lengths of other DNA pieces separated on a gel. But if you look for the word AAGCTT, where HindIII cuts, in the λ reference genome from the NCBI database, there are six sites and seven pieces. Laboratory λ is a special strain, cI857 ind1 Sam7, and it differs from the reference in four letters out of 48,502, one of which creates the seventh site. How do you find these letters? The way the differences of any genome are found: cut the DNA into short pieces, read them, which gives the reads, and find each one's place in the reference. Below are twenty thousand simulated reads of this strain, and an index built on the Burrows–Wheeler transform finds them all.

Under the lens is position 37,589. The reference has C there, while the reads that landed on this place have T. The letters around it read AAGCTT: this is the seventh HindIII site. The “next difference” button goes through all places where the reads disagree with the reference. With a read error rate of 1 %, of twenty thousand reads 12,005 land exactly, 6165 with one substitution, 1830 are not found, and there are four differences, all four real substitutions of the strain. Set the error rate to 3 %. Every read that landed on a substitution has already spent its only allowed error on it, and a read with one more error does not land at all. Now 9047 reads are not found, only three land on position 37,589, and the scene sees two substitutions of the four: searching with one error pulls the picture towards the reference.

The idea

Sort all suffixes of the text, as in a suffix array, and write out the letter that stands before each: that is the Burrows–Wheeler transform. This one string, as long as the text, lets you find any word: read it from the end, and after each letter narrow the range of table rows that begin with the tail read so far. Each step is two counts of “how many times has this letter occurred in the string up to here”.

The transform itself, and why it gathers identical letters into runs, is worked through in the course chapter “The packing contest”, along with how bzip2 compresses it and how the text is restored. Here is the other side: the transform as an index.

Append to a text T of length n a character $, smaller than all letters and occurring nowhere else, and write out all cyclic rotations of T$ in alphabetical order. This gives a table of n + 1 rows. Its first column F is the letters of the text in alphabetical order, and its last column L is the transform. Since $ stands at the end and is the smallest, the order of rotations matches the order of suffixes, and L[i] is the letter before suffix number SA[i]: L[i] = T[SA[i] − 1], and $ for the suffix that starts at the beginning of the text. For phage λ the first letters of L are GCGCGGAGAAAGGAG…: the 48,503 characters of the transform have 35,329 runs of identical letters. The DNA of a single genome has short runs, as it has few repeats.

The last column leads to the first

Let C[c] be the number of characters of T$ smaller than c, $ included: the rows beginning with c occupy the stretch of the table from C[c] to C[c] + (number of c) − 1. And let rank_c(i) be how many times the letter c occurs in L[0..i), before row i. Everything rests on one observation.

An LF step takes the row where the text starts at position p to the row where it starts at position p − 1: one letter to the left. If you start from row 0, which holds the rotation “$T” whose last letter is the last letter of the text, and take LF steps, column L spells out all the letters of the text from right to left. This is the inverse transform in n steps, and finding it yourself is the problem of the course chapter.

Backward search

The rows of the table that begin with a pattern P stand together, since they are sorted. Find this stretch by reading P from the end.

Start with the range [0, n + 1), with which all rows begin, and add the letters of the pattern one at a time from right to left. After m steps the length of the range is the number of occurrences of the pattern. Each step is two calls of rank, and if rank answers in constant time, a pattern is counted in O(m), independent of the length of the text. In the listing the EcoRI site, GAATTC, is found in the λ genome with 12 calls of rank: two per letter.

It remains to make rank fast. Keeping a ready answer for every letter and every row is σ numbers per row, more than the text itself. So answers are kept every 64 rows, and the rest is counted from L: at most 63 letters. Real implementations keep a block of L as bits, and the count is one or two population-count instructions on a machine word; for large alphabets there are wavelet trees, where rank costs O(log σ).

Here backward search runs through “War and Peace” in the Russian original. The word «наташа» (“natasha”) is read from the end. After «а» the range is huge: 183,707 suffixes begin with this letter. After «ша» there are 3908, after «аша» 972, after «таша» 827, and the number stops changing: every «таша» in the novel comes from the name «Наташа». Try the single letter «и» and compare with a long phrase: the number of rank calls depends only on the length of what you are looking for.

The main thing in this scene is what it lacks. The browser did not load the text of the novel: only the Burrows–Wheeler transform and the positions of every 256th row. The count, the places on the strip and even the excerpts around what was found it got from this one string. The file with the permutation and the samples, compressed with gzip, weighs 1.0 MB, while the novel itself, compressed the same way, weighs 1.4 MB: runs of identical characters compress better than text. The novel's permutation has 1,293,114 runs over 2,904,145 characters, 2.25 on average: prose has longer runs than phage DNA.

Where the match stands

The range says how many times the pattern occurs, but not where. The position of a row is SA[i], and keeping the whole suffix array is n numbers of 4 bytes, four times the text. So only every s-th position is kept: in the listing, positions that are multiples of 32, together with the numbers of their rows. To find the position of row i, take LF steps until you reach a marked row: each step moves one letter to the left, and if after k steps you reach the row of position q, the position sought is q + k.

This is how memory and time trade places: with a step of 32 the numbers take 32 times less room than the whole array, and each position costs up to 31 steps. The novel scene marks every 256th row, to keep the file smaller, and there the steps are about 256 on average, but without a guarantee; that is why the strip shows at most three thousand places. With the same sampling in reverse, “a row for every 256th position”, the scene reads excerpts: it starts from the nearest marked position to the right of the piece needed and goes left by LF, writing out letters from L.

Reads: searching with one error

Reads almost never match the reference letter for letter: the machine makes errors, and the organism itself has its own substitutions, the ones we are looking for. Backward search is easily taught to forgive substitutions. At the step where letter P[k] is read, continue not only with it but with each of the three other letters as well, noting that the only allowed substitution has been spent. This gives backtracking over a tree of variants, but the tree is small: a branch ends as soon as its range is empty, and in a genome of 48 thousand letters almost any ten random letters in a row do not occur, so branches with a substitution near the end of a read die within a few steps.

The listing does this with two thousand 35-letter reads from both strands of the genome, each letter wrong with probability 1 %. 1409 reads landed exactly, 499 with one substitution, and 92 were not found: they have two errors or more. That is what was expected: 0.99³⁵ ≈ 70 % of reads stay error-free, and 35 · 0.01 · 0.99³⁴ ≈ 25 % have one error. All 1908 reads that landed did so at the place they were taken from. The search took 287 calls of rank per read on average, counting both strands and reporting positions; a plain pass over the genome would compare a read with 97,004 places.

The HindIII story is there too. In the reference the word AAGCTT occurs six times, and the pieces are 23,130, 2027, 2322, 9416, 564, 6682 and 4361 letters long. Change one letter, C to T at position 37,589, and there are seven sites, with the 6682 piece splitting into 125 and 6557. That is exactly the eight pieces of the laboratory marker. With EcoRI and BamHI there is no discrepancy: six EcoRI pieces, 21,226, 4878, 5643, 7421, 5804 and 3530, the same ladder laboratories get.

Code

The listing builds the permutation through a suffix array by doubling, as in the article on it: the ranks of pairs are packed into one number so as to sort integers. The index is the string L, the table C, counters every 64 rows, and the positions that are multiples of 32. The function interval is backward search, locate the LF steps to a mark, one_mismatch the search with one substitution.

Note that after building, the program never touches the text itself again: search, positions and the read check all go through L, C, the counters and the sample. In Python counting a letter within a block is done by str.count, a fast loop inside the interpreter; in Go by an ordinary loop over 64 bytes.

What it costs

For DNA σ = 4, and L takes two bits per letter. The human genome has about three billion letters: a suffix array of 4-byte numbers would take about 12 GB, while in 2009 the Bowtie program aligned more than 25 million reads an hour to the human genome on one processor core and used about 1.3 GB of memory.

Building is the heaviest part, but it is done once per genome. The doubling in the listing takes a fraction of a second for λ; the human genome needs linear algorithms and building in pieces, without holding the whole suffix array in memory.

Variants

Compression. Runs of identical letters in L compress well: bzip2 permutes blocks of up to 900 KB, recodes the letters as “how many different ones there have been since last time”, and compresses the result with a Huffman code; this is in the course chapter “The packing contest”. A compressed index keeps L already compressed and can compute rank without decompressing it whole.

The r-index. If the text is thousands of almost identical genomes of one species, the number of runs r in L grows slowly: a new genome adds almost no new runs. In 2018 Travis Gagie, Gonzalo Navarro and Nicola Prezza showed how to search and locate while keeping O(r) numbers rather than O(n). A single phage λ genome, as we saw, has almost as many runs as letters; the gain appears only on collections.

The bidirectional index. Backward search extends a pattern only to the left. If you also build the transform of the reversed text and keep the ranges in both in step, a pattern can be extended in both directions. This is how searches with several errors are done: a read is cut into parts, one of which is surely error-free, that part is found exactly and extended to the sides, spending errors only there.

Alignment with gaps. Reads can have inserted and missing letters, and one substitution will not find them. Programs like BWA-MEM find exact pieces, seeds, by backward search, and then align the read around them by dynamic programming, as in the article on sequence alignment.

History

David Wheeler came up with the permutation in 1983 but did not publish it. In 1994 it appeared in the report by Michael Burrows and Wheeler, “A block-sorting lossless data compression algorithm”, from DEC's research centre, as the basis of compression, and in 1996 Julian Seward released bzip2 built on it. Paolo Ferragina and Giovanni Manzini turned the permutation into an index: in 2000 they showed that a compressed permutation can be searched in time that depends only on the length of the pattern. They called their structure the FM-index, “Full-text index in Minute space”.

The FM-index came to biology when machines began producing millions of short reads per run. In 2009 Bowtie appeared, by Ben Langmead, Cole Trapnell, Mihai Pop and Steven Salzberg, with a backtracking search that takes the quality of each read letter into account, and so did BWA, the Burrows–Wheeler Aligner by Heng Li and Richard Durbin. Their new versions, Bowtie 2 and BWA-MEM, are still in use, and the r-index carries the same idea to collections of thousands of genomes.

Problems

In both problems only the permutation is given, with no text. In the first you have to count many words faster than you could restore the text and read it. In the second, find the places of reads with one error, as in the listing.