How BWA generate index files?

Question

I need to know how BWA generate bw and sa in less memory usage. I didn't find details in the original paper. Do they save all rotations temporary ? I need to know do BWA stores all rotations so that it can sort the rotations alphabetically ?

Thanks.

Answer 1

Papers

The space-efficient BWT construction algorithm is in bwt_gen.c. One of the first things you see there is the copyright line:

Copyright (C) 2004, Wong Chi Kwong.

With that name, you can find the relevant papers:

• Lam et al: Compressed indexing and local alignment of DNA (Bioinformatics, 2008). This paper describes the BWT-SW alignment algorithm.
• Hon et al: Practical Aspects of Compressed Suffix Arrays and FM-index in Searching DNA Sequences (ALENEX 2004). This was the first paper to demonstrate that it's feasible to build the FM-index for a human genome.
• Hon et al: A Space and Time Efficient Algorithm for Constructing Compressed Suffix Arrays (Algorithmica, 2007). This paper contains the theoretical description of the construction algorithm.

The third paper only mentions BWT briefly in the last section. It deals with compressed suffix arrays, which are space-efficient representations of the inverse function of LF-mapping. But if you have a representation of (the inverse of) the LF-mapping, you have a representation of the BWT. Hence the algorithm can be adapted for BWT construction.

Algorithm

When building BWT, you usually append a unique character $ to the end of the text, either implicitly or explicitly. This way you can reduce sorting rotations to sorting suffixes. You sort the suffixes in lexicographic order, and then you output the preceding character for each suffix in that order to form the BWT.

Assume that you already have the BWT for string X, and you want to transform it into the BWT of string cX, where c is a single character. You can do that in three steps (this assumes 0-based array indexing):

1. Find the position i of the full suffix X, and replace the $ in the BWT with the inserted character c (set BWT[i] = c).
2. Determine the number of suffixes j that are lexicographically smaller than cX with j = LF(i, c) = C[c] + BWT.rank(i, c). Here C[c] is the number of suffixes that start with a character smaller than c, and BWT.rank(i, c) is the number of occurrences of c in the prefix BWT[0..i-1].
3. Insert $ between BWT[j-1] and BWT[j].

By iterating this, you can build the BWT for any string slowly but space-efficiently. The algorithm of Hon et al. is based on the same idea, but they make it faster by inserting multiple characters at the same time. Given the BWT of string X and another string Y, they produce the BWT of the concatenation YX. The details are quite complicated, but you can find them in the third paper and in the source code.

Answer 2

The following is only a comment to @JouniS' answer, but it is too long to fit in 512 characters.

The algorithm @JouniS was describing has been rediscovered multiple times. The C program below constructs BWT with this algorithm. It is proof-of-concept only, not intended for long strings.

#include <string.h>
#include <stdio.h>
int str2bwt(char *s)
{
    int i, l, c = 0, k = 0;
    l = strlen(s); // length of input string
    for (i = l - 1; i >= 0; --i) {
        int j, r = 0, a = s[i];
        memmove(&s[i], &s[i+1], k); // make room for insertion
        s[i+k] = a; // insert the char
        for (j = i; j < i + k; ++j) r += (s[j] <= a);
        for (; j < l; ++j) r += (s[j] < a);
        k = r + 1, c = a;
    }
    return k; // the position of sentinel '$'
}
int main(int argc, char *argv[])
{
    int k;
    if (argc == 1) return 1;
    k = str2bwt(argv[1]);
    fwrite(argv[1], 1, k, stdout); putchar('$'); puts(&argv[1][k]);
    return 0;
}

PS: this is one of the simplest ways to construct BWT and does not require any additional dependencies. It is slow for long strings, though.