Prefix-free parsing for building big BWTs

Abstract

High-throughput sequencing technologies have led to explosive growth of genomic databases; one of which will soon reach hundreds of terabytes. For many applications we want to build and store indexes of these databases but constructing such indexes is a challenge. Fortunately, many of these genomic databases are highly-repetitive-a characteristic that can be exploited to ease the computation of the Burrows-Wheeler Transform (BWT), which underlies many popular indexes. In this paper, we introduce a preprocessing algorithm, referred to as prefix-free parsing, that takes a text T as input, and in one-pass generates a dictionary D and a parse P of T with the property that the BWT of T can be constructed from D and P using workspace proportional to their total size and O(|T|)-time. Our experiments show that D and P are significantly smaller than T in practice, and thus, can fit in a reasonable internal memory even when T is very large. In particular, we show that with prefix-free parsing we can build an 131-MB run-length compressed FM-index (restricted to support only counting and not locating) for 1000 copies of human chromosome 19 in 2 h using 21 GB of memory, suggesting that we can build a 6.73 GB index for 1000 complete human-genome haplotypes in approximately 102 h using about 1 TB of memory.

Keywords: Burrows-Wheeler Transform; Compression-aware algorithms; Genomic databases; Prefix-free parsing.

Fig. 1

The suffix trie for our example with the three strings GATTACAT, GATACAT and GATTAGATA. The input is shown at the bottom, in red because we do not need to store it

Fig. 2

The suffix tree for our example. We now also need to store the input

Fig. 3

The suffix array for our example is the sequence of values stored in the leaves of the tree (which we need not store explicitly). The LF mapping is shown as the arrows between two copies of the suffix array; the arrows to values i such that $T[\mathrm{SA}[\mathit{i}]]=\mathtt{A}$ are in red, to illustrate that they point to consecutive positions in the suffix array and do not cross. Since $\mathrm{\Psi }$ is the inverse of the LF mapping, it can be obtained by simply reversing the direction of the arrows

Fig. 4

The BWT and the sorted list of characters for our example. Drawing arrows between corresponding occurrences of characters in the two strings gives us the diagram for the LF-mapping

Fig. 5

Results versus various choices of parameters (w, p) on a collection of 1000 Salmonella genomes (2.7 GB)

Fig. 6

RLFM indexing efficiency for successively larger collections of genetically distinct human chr19s. Results for the prefix-free parsing step (“pfbwt”) are shown alongside the overall RLFM index-building (“rlfm_total”) and Bowtie (“bowtie”) results

Fig. 7

RLFM average “count” query time for successively larger collections of genetically distinct human chr19s

Fig. 8

Computational efficiency of the three stages of index building when SA sampling is included. Results are shown for the prefix-free parsing (“pfbwt”), back-stepping (“bwtscan”) and run-length encoding (“rle”) steps. “total” is the sum of the three steps