Burrows-Wheeler Transform, in depth: rotations, suffix arrays, LF-mapping inversion, bzip2 and FM-index search

The transform: sorted rotations

Append a sentinel character, written $, that occurs nowhere else and sorts before every other symbol. Write down all n rotations of the string, sort them lexicographically and read off the last column. For banana$ the sorted rotations are $banana, a$banan, ana$ban, anana$b, banana$, na$bana and nana$ba, so the BWT is annb$aa.

Why do similar characters cluster? Sorting the rotations groups them by what follows each position. The last column holds the character that precedes each of those sorted contexts. In English, the contexts beginning with he are preceded by t far more often than anything else, so all those t characters land in one stretch of L. A longer repetitive example makes it visible: the BWT of abracadabra_abracadabra$ is aarrdd_$rrccaaaaaaaabbbb, with a run of eight a characters and four b characters.

Computing it from a suffix array

Sorting n full rotations is far too slow for real inputs, but the sentinel makes rotations and suffixes sort identically: no two suffixes compare past the $. So the BWT is a direct read from the suffix array SA: L[i] is the character just before suffix SA[i], wrapping to $ when SA[i] is 0.

def bwt_naive(s):
    assert s.endswith("$") and s.count("$") == 1
    rotations = sorted(s[i:] + s[:i] for i in range(len(s)))
    return "".join(r[-1] for r in rotations)

def suffix_array(s):
    # O(n^2 log n) teaching version; use SA-IS or libsais in production
    return sorted(range(len(s)), key=lambda i: s[i:])

def bwt_from_sa(s, sa):
    return "".join(s[i - 1] for i in sa)   # s[-1] is '$' when i == 0

s = "banana$"
print(suffix_array(s))              # [6, 5, 3, 1, 0, 4, 2]
print(bwt_from_sa(s, suffix_array(s)))   # annb$aa

Production code builds the suffix array in linear time with SA-IS or a tuned library such as libsais or divsufsort, using roughly 4 to 8 bytes per input character during construction. The suffix array construction article walks through prefix doubling and SA-IS; everything below assumes you can obtain SA efficiently.

Inverting with the LF mapping

The surprising property is that L alone is enough to recover the text. Sorting L gives the first column F, because F is just every character in sorted order. The key lemma is rank preservation: the k-th occurrence of a character c in L and the k-th occurrence of c in F are the same position in the original text. Rows that begin with c are sorted by what follows c; rows that end with c are sorted by what follows that same c after rotating it to the front, which is the same order.

That gives the last-to-first mapping LF(i) = C[L[i]] + occ(L[i], i), where C[c] counts characters smaller than c and occ(c, i) counts occurrences of c in L before row i. Starting at row 0, whose row is $ followed by the text, L[0] is the last real character; following LF walks the text backwards one character per step.

def inverse_bwt(L):
    counts = {}
    for ch in L:
        counts[ch] = counts.get(ch, 0) + 1
    C, total = {}, 0
    for ch in sorted(counts):              # C[c]: characters smaller than c
        C[ch] = total
        total += counts[ch]
    seen, LF = {}, []
    for ch in L:                           # LF[i] = C[L[i]] + occ(L[i], i)
        LF.append(C[ch] + seen.get(ch, 0))
        seen[ch] = seen.get(ch, 0) + 1
    row, out = 0, []
    for _ in range(len(L) - 1):
        out.append(L[row])                 # L[row] precedes row's first char
        row = LF[row]
    return "".join(reversed(out)) + "$"

print(inverse_bwt("annb$aa"))              # banana$

Worked trace for annb$aa: C is $=0, a=1, b=4, n=5, and LF for rows 0 to 6 is 1, 5, 6, 4, 0, 2, 3. From row 0 we read a and jump to row 1, read n and jump to 5, read a and jump to 2, read n and jump to 6, read a and jump to 3, read b. The collected characters a, n, a, n, a, b reversed spell banana. Inversion is linear time, but each step is a random memory access, so on large blocks it is cache-miss bound, not compute bound.

From BWT to compression

A BWT compressor is a pipeline in which each stage converts the previous stage's structure into something the next stage likes. bzip2, the best-known example, splits input into blocks of 100 to 900 kB (the -1 to -9 flags), applies the BWT to each block, then move-to-front coding, run-length coding of the resulting zeros, and Huffman coding with several tables chosen per group of symbols.

Move-to-front (MTF) keeps a list of symbols and outputs each symbol's current index, then moves it to the front. A run of identical characters becomes one index followed by zeros. On the abracadabra example above, MTF yields 15 zeros out of 24 outputs, a skewed distribution that Huffman coding encodes in far fewer bits than the uniform-looking input.

def move_to_front(L):
    alphabet = sorted(set(L))
    out = []
    for ch in L:
        k = alphabet.index(ch)
        out.append(k)
        alphabet.insert(0, alphabet.pop(k))
    return out

print(move_to_front("aarrdd_$rrccaaaaaaaabbbb"))
# [2, 0, 6, 0, 6, 0, 4, 4, 3, 0, 6, 0, 5, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0]

Block size is the main tuning knob. Larger blocks expose longer-range repetition and compress better, but cost more memory and time per block, and a block must be complete before any of it can be output, so BWT compressors cannot stream at byte granularity. Blocks are independent, which makes parallel compression straightforward, as pbzip2 does.

The FM-index and backward search

The FM-index, introduced by Ferragina and Manzini in 2000, stores L, the C array and an occ structure, and answers pattern queries by backward search: process the pattern from its last character to its first, maintaining the half-open range of sorted rows that begin with the suffix of the pattern seen so far. Each step is two LF computations.

class FMIndex:
    def __init__(self, s):
        self.sa = suffix_array(s)
        self.L = bwt_from_sa(s, self.sa)
        chars = sorted(set(self.L))
        self.C, total = {}, 0
        for ch in chars:
            self.C[ch] = total
            total += self.L.count(ch)
        # occ[ch][i] = count of ch in L[:i]; real indexes sample this
        self.occ = {ch: [0] * (len(self.L) + 1) for ch in chars}
        for i, x in enumerate(self.L):
            for ch in chars:
                self.occ[ch][i + 1] = self.occ[ch][i] + (x == ch)

    def count(self, pattern):
        lo, hi = 0, len(self.L)
        for ch in reversed(pattern):
            if ch not in self.C:
                return 0, (0, 0)
            lo = self.C[ch] + self.occ[ch][lo]
            hi = self.C[ch] + self.occ[ch][hi]
            if lo >= hi:
                return 0, (lo, hi)
        return hi - lo, (lo, hi)

    def locate(self, pattern):
        n, (lo, hi) = self.count(pattern)
        return sorted(self.sa[lo:hi])

fm = FMIndex("banana$")
print(fm.count("ana"), fm.locate("ana"))   # (2, (2, 4)) [1, 3]
print(fm.count("nab"))                     # (0, (2, 2))

Trace for ana: start with rows [0, 7). Processing a gives [1 + 0, 1 + 3) = [1, 4), the rows starting with a. Processing n gives [5 + occ(n, 1), 5 + occ(n, 4)) = [5, 7), the rows starting with na. Processing a again gives [1 + occ(a, 5), 1 + occ(a, 7)) = [2, 4). Two rows, so two occurrences; SA[2] = 3 and SA[3] = 1 give positions 1 and 3. Counting cost depends on the pattern length, not the text length.

Making the index small

The toy index above stores a full occ table and a full suffix array, which is larger than the text. Real FM-indexes make three substitutions. Occ is stored as checkpoints every 64 to 256 rows plus a popcount over the packed L bits between checkpoint and query row; for large alphabets a wavelet tree answers occ in time logarithmic in the alphabet size. The suffix array is sampled, keeping SA only for text positions that are multiples of k; locate walks LF from an unsampled row until it hits a sampled one, adding the step count. That trades up to k LF steps per reported hit for a k-fold smaller SA. Finally L itself can be entropy-compressed, which is where the compressed in FM-index comes from.

For a DNA reference, with a four-letter alphabet packed at two bits per base, this brings a genome-scale index into a few gigabytes. The choice of k is a direct latency versus memory dial: count queries do not touch SA at all, so applications that only need counts or presence can sample very sparsely. For a comparison with the other text indexes, see suffix arrays in depth.

Failure modes

• Sentinel collisions. The input already contains the byte you chose as $. Use an out-of-band sentinel, a virtual end-of-string position, or escape the input; never assume text will not contain a zero byte.
• Wrong sentinel order. If $ does not sort before every symbol, rotation order and suffix order diverge and the BWT from SA is wrong. Unit-test against the naive rotation version on random strings, as the code here was.
• Naive construction in production. Sorting suffixes with string comparison is quadratic on repetitive input and falls over on logs and genomes.
• Pathological repetition. Highly repetitive inputs can be slow for some construction algorithms; bzip2 runs a first run-length pass partly to blunt this.
• Locate blow-up. A pattern with millions of hits, such as a repeat in DNA, makes locate cost hits times k. Cap reported hits or report counts first.
• Off-by-one in occ. Mixing inclusive and exclusive occ definitions shifts every range by one. Pick half-open ranges and keep them everywhere.

Trade-offs

What to do next

1. Type in the code above and reproduce annb$aa, the inversion trace and the ana search.
2. Write a property test comparing bwt_naive with bwt_from_sa and inverse_bwt on random strings over a small alphabet, where repetition is common.
3. Replace the toy suffix array with a library call such as pydivsufsort and time a 10 MB file.
4. Compress the same file with bzip2 -9 and with gzip -9 and compare ratio and time.
5. Convert the occ table to checkpoints every 64 rows plus counting, and sample SA every 32 positions; measure memory and locate latency.
6. Read the BWA or Bowtie documentation to see how a production FM-index handles reverse complements and mismatches.