Burrows–Wheeler Transform

The Burrows–Wheeler transform (BWT) is a reversible block-sorting string permutation within computation and algorithms that rearranges a character sequence into runs of identical symbols without losing the ability to reconstruct the original text. Invented by Michael Burrows and David J. Wheeler at the Digital Equipment Corporation (DEC) Systems Research Center in 1994 (building on an unpublished 1983 permutation discovered by Wheeler), the transform was originally designed for lossless block data compression (as in the bzip2 utility), where it outperforms sliding-window Lempel–Ziv complexity coders by approaching higher-order empirical Shannon entropy at modest computational cost.

In 2000, Paolo Ferragina and Giovanni Manzini discovered that combining the Burrows–Wheeler transform with succinct rank data structures yields a compressed full-text substring index—the FM-index—capable of counting and locating occurrences of a query string of length m within a text of length N in O(m) operations, independent of N. Beginning in 2009 with the development of Bowtie (Ben Langmead et al.) and BWA (Heng Li and Richard Durbin), the BWT and FM-index revolutionized bioinformatics by reducing the memory footprint of a indexed 3.1-gigabase human reference genome from > 15 gigabytes (for uncompressed suffix trees) to ~1.5–2.5 gigabytes of RAM while accelerating short-read sequence alignment by orders of magnitude.

The Burrows–Wheeler Matrix (M)

Let T be a string of length N − 1 over an ordered finite alphabet Σ (for example, the nucleotide alphabet Σ = {A, C, G, T} with lexicographical ordering A < C < G < T). To ensure that every cyclic shift is distinct and that the end of the string is unambiguously marked, a unique sentinel character $ that is lexicographically smaller than every symbol in Σ ($ < A < C < G < T) is appended to form T$ of length N.

The conceptual construction of BWT(T$) proceeds in three steps:

1. Form Cyclic Rotations: Generate all N cyclic shifts of T$ by successively moving the terminal character to the front of the string, placing the N rotations as rows of an N × N conceptual matrix M_raw.
2. Lexicographical Sort: Sort the rows of M_raw in ascending lexicographical order to obtain the Burrows–Wheeler matrix M. The first column of M is denoted F (First), and the last (N-th) column of M is denoted L (Last).
3. Extract Last Column (L): The Burrows–Wheeler transform BWT(T$) is the string of length N formed by reading down the last column L[i] for i = 0, 1, ..., N − 1.

Because every row in M is a cyclic rotation of T$, the character in the last column L[i] directly precedes the character in the first column F[i] in the original circular string T$. Moreover, because the rows of M are sorted by their right-hand suffixes (starting at F[i]), positions in T that are followed by identical right-hand contexts (for instance, the suffix "ANA" or a conserved genomic repeat) are grouped into adjacent rows of M. Consequently, the characters directly preceding those contexts cluster into long homogeneous runs in L, making L readily compressible by run-length encoding and move-to-front coding.

Concrete Worked Example (T = "BANANA")

For the 7-character string T = "BANANA" (N = 7, alphabet ordering $ < A < B < N), sorting all 7 cyclic rotations produces the Burrows–Wheeler matrix M:

Reading down the final column L yields BWT("BANANA") = "ANNBAA", in which the two 'N' characters and two of the 'A' characters have been brought together into contiguous runs.

Linear-Time Construction via Suffix Arrays (SA)

Explicitly materializing and sorting an N × N character matrix M would require O(N²) memory and O(N² log N) comparison time—an impossible ~10¹⁹ bytes for a 3.1 × 10⁹ bp eukaryotic genome. However, because sentinel $ is smaller than all other characters, sorting the N cyclic rotations of T$ produces the identical row permutation as sorting the N suffixes of T$.

Let the suffix array *SA*[0..N − 1] (introduced by Udi Manber and Gene Myers in 1990) be the permutation of starting indices {0, 1, ..., N − 1} such that the suffix T[*SA*[i]..N − 1] is the i-th lexicographically smallest suffix of T$. Because the character in the last column L[i] of row i is simply the character directly preceding the start of suffix *SA*[i], BWT[i] is determined directly from *SA*[i] in O(1) time per entry:

• BWT[i] = T[*SA*[i] − 1] if *SA*[i] > 0
• BWT[i] = $ if *SA*[i] = 0

Using linear-time induced-sorting suffix array algorithms—such as the Karkkainen–Sanders DC3 skew algorithm (2003) or Nong, Zhang, and Chan's SA-IS (Suffix Array Induced Sorting, 2009)—both SA and the Burrows–Wheeler transform are constructed in optimal O(N) linear time and O(N log |Σ|) bits of working memory.

Why Relative Rank Order Is Preserved Between L and F

A remarkable algebraic property of the Burrows–Wheeler matrix M is that for any character c ∈ Σ, the k-th occurrence of character c in the last column L corresponds to the identical text character as the k-th occurrence of character c in the first column F.

To see why this rank invariance holds, consider any two rows i < j in M that both end with the same character L[i] = L[j] = c. In the sorted matrix M, row i precedes row j because row i's prefix string α_i is lexicographically smaller than row j's prefix string α_j (α_i < α_j). Cyclic right-rotation of both rows shifts their terminal character c to the front, yielding rows cα_i and cα_j in the set of rows that begin with c in column F. Because prepending the identical symbol c to both strings preserves their lexicographical inequality (cα_i < cα_j), the two instances of c appear in the identical relative top-to-bottom order in column F as they do in column L.

LF Mapping Formula and O(N) String Reconstruction

This order-preserving bijection defines the Last-to-First (LF) mapping, which maps any row index i in column L (where c = L[i]) to the unique row index LF(i) in column F that starts with that same character instance:

LF(i) = C[c] + Occ(c, i)

where:

1. C[c] (Cumulative Character Count Array): A tiny lookup table of length |Σ| storing the total number of characters in T$ that are lexicographically smaller than c. In sorted column F, the contiguous block of rows starting with character c begins at row index C[c].
2. Occ(c, i) (Occurrence / Rank Function): The number of times character c appears in the prefix L[0..i − 1] prior to index i (0-based rank).

Using LF(i), the original string T is reconstructed from right to left in N steps starting at row r₀ = 0 (since row 0 of M starts with $, its last character L[0] is T[N − 2], the final character before $):

For step k = 1 to N − 1: output L[r_k−1] and set r_k = LF(r_k−1).

In our "BANANA" example (C = { : 0, A : 1, B : 4, N : 5}), starting at r₀ = 0 outputs L[0] = 'A', jumps to r₁ = C['A'] + Occ('A', 0) = 1 + 0 = 1, outputs L[1] = 'N', jumps to r₂ = C['N'] + Occ('N', 1) = 5 + 0 = 5, outputs L[5] = 'A', and continues through 'N', 'A', 'B', reconstructing "BANANA$" without storing the full matrix M.

Ferragina–Manzini Backward Search Algorithm

In 2000, Paolo Ferragina and Giovanni Manzini realized that the LF mapping turns the Burrows–Wheeler transform into a self-contained substring index (FM-index). In a suffix array SA, all occurrences of a query pattern P[0..m − 1] share the common prefix P and therefore occupy a single contiguous interval of rows [sp, ep] in matrix M.

Rather than binary-searching SA forward in O(m log N) time, the FM-index narrows the suffix array interval [sp, ep] backward from the rightmost character P[m − 1] down to P[0], updating the interval bounds at each character c = P[k] via two LF evaluations:

sp_new = C[c] + Occ(c, sp_old)

ep_new = C[c] + Occ(c, ep_old + 1) − 1

If sp_new ≤ ep_new after processing all m characters of P, then pattern P occurs ep − sp + 1 times in the text T, and every row i ∈ [sp, ep] corresponds to a genomic coordinate *SA*[i].

Wavelet Trees, Occurrence Checkpoints, and Sampled Suffix Arrays

To make Occ(c, i) execute in O(1) constant time without storing an uncompressed |Σ| × N integer table, the FM-index divides the BWT string L into fixed-size blocks (typically every B = 32, 64, or 128 characters) and stores checkpoint counts at block boundaries. Evaluating Occ(c, i) adds the precomputed checkpoint count at block ⌊i / B⌋ to a hardware POPCNT bit-parallel population count over the remaining i mod B bit-packed 2-bit DNA characters—completing in a single CPU cache-line access. For larger alphabets (such as 20 amino acids or ASCII text), Roberto Grossi, Ankur Gupta, and Jeffrey Scott Vitter (2003) introduced the wavelet tree, a balanced binary tree over Σ that answers Occ(c, i) in O(log |Σ|) bit-vector rank operations.

Similarly, to convert matched BWT row indices i ∈ [sp, ep] back into reference genome coordinates *SA*[i] without storing all N 32-bit suffix array integers (which would cost 12 GB for the human genome), the FM-index stores a sampled suffix array at every s-th position (typically s = 32). If row i is not sampled directly, the aligner steps backward t < s times via i ← LF(i) until hitting a sampled row i′, returning coordinate *SA*[i′] + t.

Seed-and-Extend Alignment with Mismatches and Gaps

Biological sequencing reads differ from a reference genome due to single-nucleotide polymorphisms (SNPs), insertions/deletions (indels), and sequencing errors. While pure backward search finds zero-mismatch matches in O(m) time, Bowtie (Langmead et al., 2009) and BWA-backtrack (Li and Durbin, 2009) accommodate mismatches by bounded backtracking across the FM-index trie, pruning branches using Phred quality ceilings and bidirectional BWT bounds.

For longer Illumina, PacBio, and Oxford Nanopore reads, Bowtie 2 (2012) and BWA-MEM (Li, 2013) employ a seed-and-extend architecture: BWA-MEM uses a bidirectional FMD-index (indexing both forward and reverse-complement DNA strands simultaneously) to identify Super-Maximal Exact Matches (SMEMs) in O(m) time, chains collinear seeds across the genome, and executes SIMD-vectorized banded Smith–Waterman / Gotoh dynamic programming only within the small genomic windows flanking high-scoring seed chains. More recently, Gagie, Navarro, and Prezza's r-index (2020) indexes run-length encoded BWTs in O(r) space (where r is the number of BWT runs), allowing thousands of human pangenome assemblies to be indexed simultaneously in a few gigabytes.

Why does backward search in an FM-index proceed from the end of the query pattern to the beginning?

In the Burrows–Wheeler matrix M, column L contains the character directly preceding each sorted suffix in column F. Given the suffix-array interval [sp, ep] of all text positions matching suffix P[k + 1..m − 1], the LF mapping directly identifies the sub-interval of those suffixes that are preceded by character P[k], extending the match leftward by one character to P[k..m − 1] in O(1) time.

How does the Burrows–Wheeler transform compress data if it only permutes characters without shortening the string?

The BWT itself has the identical length N as the input string T$, so it does not compress data directly. Instead, by sorting characters according to their right-hand contexts, the BWT transforms local context dependencies into long contiguous runs of identical symbols. Passing BWT(T$) through a Move-to-Front (MTF) transform and Run-Length + Huffman/Arithmetic coding compresses those homogeneous runs to near the k-th order empirical entropy H_k(T).

1. Burrows, M., & Wheeler, D. J. (1994). A block-sorting lossless data compression algorithm. DEC Systems Research Center Research Report 124. https://doi.org/10.1007/978-0-387-30162-4_57
2. Ferragina, P., & Manzini, G. (2000). Opportunistic data structures with applications. Proceedings 41st Annual Symposium on Foundations of Computer Science (FOCS), 390–398. https://doi.org/10.1109/SFCS.2000.892127
3. Langmead, B., Trapnell, C., Pop, M., & Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 10(3), R25. https://doi.org/10.1186/gb-2009-10-3-r25
4. Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 25(14), 1754–1760. https://doi.org/10.1093/bioinformatics/btp324
5. Gagie, T., Navarro, G., & Prezza, N. (2020). Fully functional suffix trees and optimal text searching in BWT-runs bounded space. Journal of the ACM, 67(1), 1–54. https://doi.org/10.1145/3375890