Efficient Construction of a Complete Index for Pan-Genomics Read Alignment

by   Alan Kuhnle, et al.

While short read aligners, which predominantly use the FM-index, are able to easily index one or a few human genomes, they do not scale well to indexing databases containing thousands of genomes. To understand why, it helps to examine the main components of the FM-index in more detail, which is a rank data structure over the Burrows-Wheeler Transform (BWT) of the string that will allow us to find the interval in the string's suffix array (SA) containing pointers to starting positions of occurrences of a given pattern; second, a sample of the SA that --- when used with the rank data structure --- allows us access the SA. The rank data structure can be kept small even for large genomic databases, by run-length compressing the BWT, but until recently there was no means known to keep the SA sample small without greatly slowing down access to the SA. Now that Gagie et al. (SODA 2018) have defined an SA sample that takes about the same space as the run-length compressed BWT --- we have the design for efficient FM-indexes of genomic databases but are faced with the problem of building them. In 2018 we showed how to build the BWT of large genomic databases efficiently (WABI 2018) but the problem of building Gagie et al.'s SA sample efficiently was left open. We compare our approach to state-of-the-art methods for constructing the SA sample, and demonstrate that it is the fastest and most space-efficient method on highly repetitive genomic databases. Lastly, we apply our method for indexing partial and whole human genomes, and show that it improves over Bowtie with respect to both memory and time.


page 1

page 2

page 3

page 4


Strong link between BWT and XBW via Aho-Corasick automaton and applications to Run-Length Encoding

The boom of genomic sequencing makes compression of set of sequences ine...

A Memory-Efficient FM-Index Constructor for Next-Generation Sequencing Applications on FPGAs

FM-index is an efficient data structure for string search and is widely ...

MONI can find k-MEMs

Maximal exact matches (MEMs) have been widely used in bioinformatics at ...

Tailoring r-index for metagenomics

A basic problem in metagenomics is to assign a sequenced read to the cor...

The Read-Optimized Burrows-Wheeler Transform

The advent of high-throughput sequencing has resulted in massive genomic...

PFP Data Structures

Prefix-free parsing (PFP) was introduced by Boucher et al. (2019) as a p...

Computing Maximal Unique Matches with the r-index

In recent years, pangenomes received increasing attention from the scien...

1 Introduction

The FM-index, which is a compressed subsequence index based on Burrows Wheeler transform (BWT), is the primary data structure majority of short read aligners — including Bowtie [19], BWA [13] and SOAP2 [21]. These aligners build a FM-index based data structure of sequences from a given genomic database, and then use the index to perform queries that find approximate matches of sequences to the database. And while these methods can easily index one or a few human genomes, they do not scale well to indexing the databases of thousands of genomes. This is problematic in analysis of the data produced by consortium projects, which routinely have several thousand genomes.

In this paper, we address this need by introducing and implementing an algorithm for efficiently constructing the FM-index. This will allow for the FM-index construction to scale to larger sets of genomes. To understand the challenge and solution behind our method, it helps to examine the two principal components of the FM-index: first a rank data structure over the BWT of the string that enables us to find the interval in the string’s suffix array (SA) containing pointers to starting positions of occurrences of a given pattern (and, thus, compute how many such occurrences there are); second, a sample of the SA that, when used with the rank data structure, allows us access the SA (so we can list those starting positions). Searching with an FM-index can be summarized as follows: starting with the empty suffix, for each proper suffix of the given pattern we use rank queries at the ends of the BWT interval containing the characters immediately preceding occurrences of that suffix in the string, to compute the interval containing the characters immediately preceding occurrences of the suffix of length 1 greater; when we have the interval containing the characters immediately preceding occurrences of the whole pattern, we use a SA sample to list the contexts of the corresponding interval in the SA, which are the locations of those occurrences.

Although it is possible to use a compressed implementation of the rank data structure that does not become much slower or larger even for thousands of genomes, the same cannot be said for the SA sample. The product of the size and the access time must be at least linear in the length of the string for the standard SA sample. This implies that the FM-index will become much slower and/or much larger as the number of genomes in the databases grows significantly. This bottleneck has forced researchers to consider variations of FM-indexes adapted for massive genomic datasets, such as Valenzuela et al.’s pan-genomic index [33] or Garrison et al.’s variation graphs [7]. Some of these proposals use elements of the FM-index, but all deviate in substantial ways from the description above. Not only does this mean they lack the FM-index’s long and successful track record, it also means they usually do not give us the BWT intervals for all the suffixes as we search (whose lengths are the suffixes’ frequencies, and thus a tightening sequence of upper bounds on the whole pattern’s frequency), nor even the final interval in the suffix array (which is an important input in other string processing tasks).

Recently, Gagie, Navarro and Prezza [11] proposed a different approach to SA sampling, that takes space proportional to that of the compressed rank data structure while still allowing reasonable access times. While their result yields a potentially practical FM-index on massive databases, it does not directly lead to a solution since the problem of how to efficiently construct the BWT and SA sample remained open. In a direction toward to fully realizing the theoretical result of Gagie et al. [11], Boucher et al. [2] showed how to build the BWT of large genomic databases efficiently. We refer to this construction as prefix-free parsing. It takes as input string , and in one-pass generates a dictionary and a parse of with the property that the BWT can be constructed from dictionary and parse using workspace proportional to their total size and time. Yet, the resulting index of Boucher et al. [2] has no SA sample, and therefore, only supports counting and not locating. This makes this index not directly applicable to many bioinformatic applications, such as sequence alignment.

Our contributions. In this paper, we present a solution for building the FM-index111With the SA sample of Gagie et al. [11], this index is termed the -index. for very large datasets by showing that we can build the BWT and Gagie et al.’s SA sample together in roughly the same time and memory needed to construct the BWT alone. We note that this algorithm is also based on prefix-free parsing. Thus, we begin by describing how to construct the BWT from the prefix-free parse, and then show that it can be modified to build the SA sample in addition to the BWT in roughly the same time and space. We implement this approach, and refer to the resulting implementation as bigbwt. We compare it to state-of-the-art methods for constructing the SA sample, and demonstrate that bigbwt currently the fastest and most space-efficient method for constructing the SA sample on large genomic databases.

Next, we demonstrate the applicability of our method to short read alignment. In particular, we compare the memory and time needed by our method to build an index for collections of chromosome 19 with that of Bowtie. Through these experiments, we show that Bowtie was unable to build indexes for our largest collections (500 or more) because it exhausted memory, whereas our method was able to build indexes up to 1,000 chromosome 19s (and likely beyond). At 250 chromosome 19 sequences, the our method required only about 2% of the time and 6% the peak memory of Bowtie’s. Lastly, we demonstrate that it is possible to index collections of whole human genome assemblies with sub-linear scaling as the size of the collection grows.

Related work. The development of methods for building and the FM-index on large datasets is closely related to the development short-read aligners for pan-genomics — an area where there is growing interest [27, 5, 12]. Here, we briefly describe some previous approaches to this problem and detail its connection to the work in this paper. We note that majority of pan-genomic aligners requiring building the FM-index for a population of genomes and thus, can increase proficiency by using the methods described in this paper.

GenomeMapper [27], the method of Danek et al. [5], and GCSA [29] represent the genomes in a population as a graph, and then reduce the alignment problem to finding a path within the graph. Hence, these methods require all possible paths to be identified, which is exponential in the worst case. Some of these methods — such as the GCSA — use the FM-index to store and query the graph and could capitalize on our approach by building the index in the manner described here. Another set of approaches [24, 8, 12, 32] consider the reference pan-genome as the concatenation of individual genomes and exploits redundancy by using a compressed index. The hybrid index [8] operates on a Lempel-Ziv compression of the reference pan-genome. An input parameter sets the maximum length of reads that can be aligned; the parameter

has a large impact on the final size of the index. For this reason, the hybrid index is suitable for short-read alignment only, although there have been recent heuristic modifications to allow longer alignments

[9]. In contrast, the -index, of which we provide an implementation in this work, has no such length limitation. The most recent implementation of the hybrid index is CHIC [33]. Although CHIC has support for counting multiple occurrences of a pattern within a genomic database, it is an expensive operation, namely , where is the number of occurrences in the databases and is the length of the database. However, the -index is capable of counting all occurrences of a pattern of length in time up to polylog factors. There are a number of other approaches building off the hybrid index or similar ideas [5, 34]; for an extended discussion, we refer the reader to the survey of Gagie and Puglisi [12].

Finally, a third set of approaches [14, 23] attempts to encode variants within a single reference genome. BWBBLE by Huang et al. [14] follows this by supplementing the alphabet to indicate if multiple variants occur at a single location. This approach does not support counting of the number of variants matching a specific alignment; also, it suffers from memory blow-up when larger structural variations occur.

2 Background

2.1 Bwt and FM indexes

Consider a string of length from a totally ordered alphabet , such that the last character of is lexicographically less than any other character in . Let be the list of ’s characters sorted lexicographically by the suffixes starting at those characters, and let be the list of ’s characters sorted lexicographically by the suffixes starting immediately after those characters. The list is termed the Burrows-Wheeler Transform [3] of and denoted BWT. If is in position in then is in position in . Moreover, if then and have the same relative order in both lists; otherwise, their relative order in is the same as their lexicographic order. This means that if is in position in then, assuming arrays are indexed from 0 and denotes lexicographic precedence, in it is in position . The mapping is termed the LF mapping. Finally, notice that the last character in always appears first in . By repeated application of the LF mapping, we can invert the BWT, that is, recover from . Formally, the suffix array SA of the string is an array such that entry is the starting position in of the th largest suffix in lexicographical order. The above definition of the BWT is equivalent to the following:


The BWT was introduced as an aid to data compression: it moves characters followed by similar contexts together and thus makes many strings encountered in practice locally homogeneous and easily compressible. Ferragina and Manzini [10] showed how the BWT may be used for indexing a string : given a pattern of length , find the number and location of all occurrences of within . If we know the range occupied by characters immediately preceding occurrences of a pattern in , then we can compute the range occupied by characters immediately preceding occurrences of in , for any character , since

Notice is the number of occurrences of in . The essential components of an FM-index for are, first, an array storing for each character and, second, a rank data structure for BWT that quickly tells us how often any given character occurs up to any given position222Given a sequence (string) over an alphabet , a character , and an integer , is the number of times that appears in .. To be able to locate the occurrences of patterns in (in addition to just counting them), the FM-index uses a sampled333Sampled means that only some fraction of entries of the suffix array are stored. suffix array of

and a bit vector indicating the positions in

BWT of the characters preceding the sampled suffixes.

2.2 Prefix-free parsing

Next, we give an overview of prefix-free parsing, which produces a dictionary and a parse by sliding a window of fixed width through the input string . We refer the reader to Boucher et al. [2] for the formal proofs and Section 3.1 for the algorithmic details. A rolling hash function identifies when substrings are parsed into elements of a dictionary, which is a set of substrings of . Intuitively, for a repetitive string, the same dictionary phrases will be encountered frequently.

We now formally define the dictionary and parse . Given a string444For technical reasons, the string must terminate with copies of lexicographically least symbol. of length , window size and modulus , we construct the dictionary of substrings of and the parse as follows. We let be a hash function on strings of length , and let be the sequence of substrings such that or or , ordered by initial position in ; let .By construction the strings

form a parsing of in which each pair of consecutive strings and overlaps by exactly characters. We define ; that is, consists of the set of the unique substrings of such that and the first and last characters of form consecutive elements in . If has many repetitions we expect that . With a little abuse of notation we define the parsing as the sequence of lexicographic ranks of substrings in : . The parse indicates how may be reconstructed using elements of . The dictionary and parse may be constructed in one pass over in time if the hash function can be computed in constant time.

2.3 -index locating

Policriti and Prezza [26] showed that if we have stored for each value such that is the beginning or end of a run (i.e., a maximal non-empty unary substring) in BWT, and we know both the range occupied by characters immediately preceding occurrences of a pattern in and the starting position of one of those occurrences of , then when we compute the range occupied by characters immediately preceding occurrences of in , we can also compute the starting position of one of those occurrences of . Bannai et al [1] then showed that even if we have stored only for each value such that is the beginning of a run, then as long as we know we can compute .

Gagie, Navarro and Prezza [11] showed that if we have stored in a predecessor data structure for each value such that is the beginning of a run in BWT, with stored as satellite data, then given we can compute in time as where is a query to the predecessor data structure. Combined with Bannai et al.’s result, this means that while finding the range occupied by characters immediately preceding occurrences of a pattern , we can also find and then report in -time, that is, -time per occurrence.

Gagie et al. gave the name -index to the index resulting from combining a rank data structure over the run-length compressed BWT with their SA sample, and Bannai et al. used the same name for their index. Since our index is an implementation of theirs, we keep this name; on the other hand, we do not apply it to indexes based on run-length compressed BWTs that have standard SA samples or no SA samples at all.

3 Methods

Here, we describe our algorithm for building the SA or the sampled SA from the prefix free parse of a input string , which is used to build the -index. We first review the algorithm from [2] for building the BWT of from the prefix free parse. Next, we show how to modify this construction to compute the SA or the sampled SA along with the BWT.

3.1 Construction of Bwt from prefix-free parse

We assume we are given a prefix-free parse of with window size consisting of a dictionary and a parse . We represent the dictionary as a string where ’s are the dictionary phrases in lexicographic order and is a unique separator. We assume we have computed the SA of , denoted by in the following, and the suffix array of , denoted , and the array such that stores the number of occurrences in the parse of the dictionary phrase . These preliminary computations take time.

By the properties of the prefix-free parsing, each suffix of is prefixed by exactly one suffix of a dictionary phrase with . We call the representative prefix of the suffix . From the uniqueness of the representative prefix we can partition ’s suffix array into ranges

with , for , and , such that for all suffixes

have the same representative prefix . By construction .

By construction, any suffix of the dictionary is also prefixed by the suffix of a dictionary phrase. For , let denote the longest prefix of which is the suffix of a phrase (i.e. ). By construction the strings ’s are lexicographically sorted . Clearly, if we compute and discard those such that , the remaining ’s will coincide with the representative prefixes ’s. Since both ’s and ’s are lexicographically sorted, this procedure will generate the representative prefixes in the order . We note that more than one can be equal to some since different dictionary phrases can have the same suffix.

We scan , compute and use these strings to find the representative prefixes. As soon as we generate an we compute and output the portion corresponding to the range associated to . To implement the above strategy, assume there are exactly entries in prefixed by . This means that there are distinct dictionary phrases that end with . Hence, the range contains elements. To compute we need to: 1) find the symbol immediately preceding each occurrence of in , and 2) find the lexicographic ordering of ’s suffixes prefixed by . We consider the latter problem first.

Computing the lexicographic ordering of suffixes.

For consider the -th occurrence of in and let denote the position in the parsing of of the phrase ending with the -th occurrence of . In other words, is a dictionary phrase ending with and . By the properties of the lexicographic ordering of ’s suffixes prefixed by coincides with the ordering of the symbols in . In other words, precedes in if and only if ’s suffix prefixed by the -th occurrence of is lexicographically smaller than ’s suffix prefixed by the -th occurrence of .

We could determine the desired lexicographic ordering by scanning and noticing which entries coincide with one of the dictionary phrases that end with but this would clearly be inefficient. Instead, for each dictionary phrase we maintain an array of length containing the indexes such that . These sorts of “inverted lists” are computed at the beginning of the algorithm and replace the which can be discarded.

Finding the symbol preceding .

Given a representative prefix from we retrieve the indexes of the dictionary phrases that end with . Then, we retrieve the inverted lists and we merge them obtaining the list of the positions such that is a dictionary phrase ending with . Such list implicitly provides the lexicographic order of ’s suffixes starting with .

To compute the BWT we need to retrieve the symbols preceding such occurrences of . If is not a dictionary phrase, then is a proper suffix of the phrases and the symbols preceding in are those preceding in that we can retrieve from and . If coincides with a dictionary phrase , then it cannot be a suffix of another phrase. Hence, the symbols preceding in are those preceding in that we store at the beginning of the algorithm in an auxiliary array along with the inverted list .

3.2 Construction of Sa and Sa sample along with the Bwt

We now show how to modify the above algorithm so that, along with BWT, it computes the full SA of or the sampled SA consisting of the values and , where is the number of maximal non-empty runs in BWT and and are the starting and ending positions in BWT of the -th such run, respectively. Note that if we compute the sampled SA the actual output will consist of start-run pairs and end-run pairs since the values alone are not enough for the construction of the -index.

We solve both problems using the following strategy. Simultaneously to each entry , we compute the corresponding entry . Then, if we need the sampled SA, we compare and and if they differ, we output the pair among the end-runs and the pair among the start-runs. To compute the SA entries, we only need additional arrays (one for each dictionary phrase), where , and contains the ending position in of the dictionary phrase which is in position of .

Recall that in the above algorithm for each occurrence of a representative prefix , we compute the indexes of the dictionary phrases that end with . Then, we use the lists to retrieve the positions of all the occurrences of in , thus establishing the relative lexicographic order of the occurrences of the dictionary phrases ending with . To compute the corresponding entries, we need the starting position in of each occurrence of . Since the ending position in of the phrase with relative lexicographic rank is , the corresponding SA entry is . Hence, along with each BWT entry we obtain the corresponding SA entry which is saved to the output file if the full SA is needed, or further processed as described above if we need the sampled SA.

4 Time and memory usage for Sa and Sa sample construction

We compare the running time and memory usage of bigbwt555Our implementation of the algorithm in Section 3, available here: https://gitlab.com/manzai/Big-BWT. with the following methods, which represent the current state-of-the-art.

  • Once the BWT has been computed, the SA or SA sample may be computed by applying the LF mapping to invert the BWT and the application of Eq. 1. Therefore, as a baseline, we use bigbwt to construct the BWT only, as in Boucher et al. [2]; next, we load the BWT as a Huffman-compressed string with access, rank, and select support to compute the LF mapping. We step backwards through the BWT and compute the entries of the SA in non-consecutive order. Finally, these entries are sorted in external memory to produce the SA or SA sample. This method may be parallelized when the input consists of multiple strings by stepping backwards from the end of each string in parallel.

  • A second baseline is to compute the SA directly from the input; for this computation, we use the external-memory algorithm pSAscan [17], with available memory set to the memory required by bigbwt on the specific input; with the ratio of memory to input size obtained from bigbwt, pSAscan is the current state-of-the-art method to compute the SA. Once pSAscan has computed the full SA, the SA sample may be constructed by loading the input text into memory, streaming the SA from the disk, and the application of Eq. 1 to detect run boundaries. We denote this method of computing the SA sample by pSAscan+.

We compared the performance of all the methods on two datasets: (1) Salmonella genomes obtained from GenomeTrakr [31]; and (2) chromosome 19 haplotypes derived from the 1000 Genomes Project phase 3 data [4]. The Salmonella strains were downloaded from NCBI (NCBI BioProject PRJNA183844) and preprocessed by assembling each individual sample with IDBA-UD [25] and counting -mers (=32) using KMC [6]. We modified IDBA by setting kMaxShortSequence to per public advice from the author to accommodate the longer paired end reads that modern sequencers produce. We sorted the full set of samples by the size of their -mer counts and selected 1,000 samples about the median. This avoids exceptionally short assemblies, which may be due to low read coverage, and exceptionally long assemblies which may be due to contamination.

Next, we downloaded and preprocessed a collection of chromosome 19 haplotypes from 1000 Genomes Project. Chromosome 19 is 58 million base pairs in length and makes up around 1.9% of the total human genome sequence. Each sequence was derived by using the bcftools consensus tool to combine the haplotype-specific (maternal or paternal) variant calls for an individual in the 1KG project with the chr19 sequence in the GRCH37 human reference, producing a FASTA record per sequence. All DNA characters besides A, C, G, T and N were removed from the sequences before construction.

(a) Salmonella, 1 thread
(b) chr19, 1 thread
(c) chr19, 16 threads
(d) Peak memory usage
Figure 1: Runtime and peak memory usage for construction of full SA.
(a) Salmonella, 1 thread
(b) chr19, 1 thread
(c) chr19, 16 threads
(d) Peak memory usage
Figure 2: Runtime and peak memory usage for construction of SA sample.

We performed all experiments in this section on a machine with Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHz and 324 GB RAM. We measured running time and peak memory footprint using /usr/bin/time -v, with peak memory footprint captured by the Maximum resident set size (kbytes) field and running time by the User Time and System Time field.

We witnessed that the running time of each method to construct the full SA is shown in Figs. 1(a)1(c). On both the Salmonella and chr19 datasets, bigbwt ran the fastest, often by more than an order of magnitude. In Fig. 1(d), we show the peak memory usage of bigbwt as a function of input size. Empirically, the peak memory usage was sublinear in input size, especially on the chr19 data, which exhibited a high degree of repetition. Despite the higher diversity of the Salmonella genomes, bigbwt remained space-efficient and the fastest method for construction of the full SA. Furthermore, we found qualitatively similar results for construction of the SA sample, shown in Fig. 2. Similar to the results on full SA construction, bigbwt outperformed both baseline methods and exhibited sublinear memory scaling on both types of databases.

5 Application to many human genome sequences

We studied how the -index scales to repetitive texts consisting of many similar genomic sequences. Since an ultimate goal is to improve read alignment, we benchmark against Bowtie (version 1.2.2) [19] . We ran Bowtie with the -v 0 and --norc options; -v 0 disables approximate matching, while --norc causes Bowtie (like -index) to perform the locate query with respect to the query sequence only and not its reverse complement.

5.1 Indexing chromosome 19s

We performed our experiments on collections of one or more versions of chromosome 19. These versions were obtained from 1000 Genomes Project haplotypes in the manner described in the previous section. We used 10 collections of chromosome 19 haplotypes, containing 1, 2, 10, 30, 50, 100, 250, 500, and 1000 sequences, respectively. Each collection is a superset of the previous. Again, all DNA characters besides A, C, G, T and N were removed from the sequences before construction. All experiments in this section were ran on a Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50GHz machine with 512GB memory. We measured running time and peak memory footprint as described in the previous section.

First we constructed -index and Bowtie indexes on successively larger chromosome 19 collections (Figure 3, 3). The -index’s peak memory is substantially smaller than Bowtie’s for larger collections, and the gap grows with the collection size. At 250 chr19s, the -index procedure takes about 2% of the time and 6% the peak memory of Bowtie’s procedure. Bowtie fails to construct collections of more than 250 sequences due to memory exhaustion.

Next, we compared the disk footprint of the index files produced by Bowtie and -index (Figure 3). The -index currently stores only the forward strand of the sequence, while the Bowtie index stores both the forward sequence and its reverse as needed by its double-indexing heuristic [19]. Since the heuristic is relevant only for approximate matching, we omit the reverse sequence in these size comparisons. We also omit the 2-bit encoding of the original text (in the *.3.ebwt and *.4.ebwt files) as these too are used only for approximate matching. Specifically, the Bowtie index size was calculated by adding the sizes of the forward *.1.ebwt and *.2.ebwt files, which contain the BWT, SA sample, and auxiliary data structures for the forward sequence. The size of the -index increased more slowly than Bowtie’s, though the -index was larger for the smallest collections. This is because, unlike Bowtie which samples a constant fraction of the SA elements (every 32nd by default), the density of the -index SA sample depends on the ratio . When the collection is small, is small and more SA samples must be stored per base. At 250 sequences, the -index index takes 6% the space of the Bowtie index.

Figure 3: Scalability of r-index and bowtie indexes against chr19 haplotype collection size and total sequence length (megabases) with respect to index construction time (seconds) (a), index construction peak memory (megabytes) (b), index disk space (megabytes) (c), and locate time (seconds) of 100,000 100bp queries (d).

We then compared the speed of the locate query for -index and Bowtie. We extracted 100,000 100-character substrings from the chr19 collection of size 1, which is also contained in all larger collections. We queried these against both the Bowtie and -indexes. We used the --max-hits option for -index and the -k option for Bowtie to set the maximum number of hits reported to be equal to the collection size. The actual number of hits reported will often equal this number, but could be smaller (if the substring differs between individuals due to genetic variation) or larger (if the substring is from a repetitive portion of the genome). Since the source of the substrings is present in all the collections, every query is guaranteed to match at least once. As seen in Figure 3, the -index locate query was faster for the collection of 250 chr19s. No comparison was possible for larger collections because Bowtie could not build the indexes.

5.2 Indexing whole human genomes

Lastly, we used -index to index many human genomes at once. We repeated our measurements for successively larger collections of (concatenated) genomes. Thus, we first evaluated a series of haplotypes extracted from the 1000 Genomes Project [4] phase 3 callset (1KG). These collections ranged from 1 up to 10 genomes. As the first genome, we selected the GRCh37 reference itself. For the remaining 9, we used bcftools consensus to insert SNVs and other variants called by the 1000 Genomes Project for a single haplotype into the GRCh37 reference.

Second, we evaluated a series of whole-human genome assemblies from 6 different long-read assembly projects (“LRA”). We selected GRCh37 reference as the first genome, so that the first data point would coincide with that of the previous series. We then added long-read assemblies from a Chinese genome assembly project [28], a Korean genome assembly project [16] a project to assemble the well-studied NA12878 individual [15], a hydatidiform mole (known as CHM1) assembly project [30] and the Celera human genome project [20]. Compared to the series with only 1000 Genomes Project individuals, this series allowed us to measure scaling while capturing a wider range of genetic variation between humans. This is important since de novo human assembly projects regularly produce assemblies that differ from the human genome reference by megabases of sequence (12 megabases in the case of the Chinese assembly [28]), likely due to prevalent but hard-to-profile large-scale structural variation. Such variation was not comprehensively profiled in the 1000 Genomes Project, which relied on short reads.

The 1KG and LRA series were evaluated twice, once on the forward genome sequences and once on both the forward and reverse-complement sequences. This accounts for the fact that different de novo assemblies make different decisions about how to orient contigs. The -index method achieves compression only with respect to the forward-oriented versions of the sequences indexed. That is, if two contigs are reverse complements of each other but otherwise identical, -index achieves less compression than if their orientations matched. A more practical approach would be to index both forward and reverse-complement sequences, as Bowtie 2 [18] and BWA [22] do.

# Genomes Length (MB)
1 6,072 6,072 1.86 1.86
2 12,144 12,484 3.70 3.58
3 18,217 17,006 5.38 4.83
4 24,408 22,739 7.13 6.25
5 30,480 28,732 8.87 7.80
6 36,671 34,420 10.63 9.28
Figure 4: Sequence length and statistic with respect to number of whole genomes for the first 6 collections in the 1000 Genomes (1KG) and long-read assembly (LRA) series.
Figure 5: Peak index-building memory for -index when indexing successively large collections of 1000-Genomes individuals (1KG) and long-read whole-genome assemblies (LRA).

We measured the peak memory footprint when indexing these collections (Figure 5). We ran these experiments on an Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz system with 256GB memory. Memory footprints for LRA grew more quickly than those for 1KG. This was expected due to the greater genetic diversity captured in the assemblies. This may also be due in part to the presence of sequencing errors in the long-read assembles; long-read technologies are more prone to indel errors than short-read technologies, for examples, and some may survive in the assemblies. Also as expected, memory footprints for the LRA series that included both forward and reverse complement sequences grew more slowly than when just the forward sequence was included. This is due to sequences that differ only (or primarily) in their orientation between assemblies. All series exhibit sublinear trends, highlighting the efficacy of -index compression even when indexing genetically diverse whole-genome assemblies. Indexing the forward and reverse complement strands of 10 1KG individuals took about 6 hours and 20 minutes and the final index size was 36GB.

We also measured lengths and ratios for each collection of whole genomes (Table 5). Consistent with the memory-scaling results, we see that the ratios are somewhat lower for the LRA series than for the 1KG series, likely due to greater genetic diversity in the assemblies.

6 Conclusions and Future Work

We give an algorithm for building the SA and SA sample from the prefix-free parse of an input string , which fully completes the practical challenge of building the index proposed by Gagie et al. [11]. This leads to a mechanism for building a complete index of large databases — which is the linchpin in developing practical means for pan-genomics short read alignment. In fact, we apply our method for indexing partial and whole human genomes, and show that it scales better than Bowtie with respect to both memory and time. This allows for an index to be constructed for large collections of chromosome 19s (500 or more); a task that is out of reach of Bowtie — as exceeded our limit of 512 GB of memory.

Even though this work opens up doors to indexing large collections of genomes, it also highlights problems that warrant further investigation. For example, there still remains a significant amount of work in adapting the index to work well on large sets of sequence reads. This problem not only requires the construction of the

-index but also an efficient means to update the index as new datasets become available. Moreover, there is interest in supporting more sophisticated queries than just pattern matching, which would allow for more complex searches of large databases.


  • [1] H. Bannai, T. Gagie, and T. I. Online LZ77 parsing and matching statistics with RLBWTs. In Proceedings of tjhe 29th Annual Symposium on Combinatorial Pattern Matching (CPM), volume 105, pages 7:1–7:12, 2018.
  • [2] C. Boucher, T. Gagie, A. Kuhnle, and G. Manzini. Prefix-free parsing for building big BWTs. In Proceedings of 18th International Workshop on Algorithms in Bioinformatics (WABI), volume 113, pages 2:1–2:16, 2018.
  • [3] M. Burrows and D.J. Wheeler. A block sorting lossless data compression algorithm. Technical Report 124, Digital Equipment Corporation, 1994.
  • [4] The 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature, 526(7571):68–74, October 2015.
  • [5] A. Danek, S. Deorowicz, and S. Grabowski. Indexes of large genome collections on a PC. PLoS ONE, 9(10), 2014.
  • [6] S. Deorowicz, M. Kokot, S. Grabowski, and A. Debudaj-Grabysz. KMC 2: Fast and resource-frugal k-mer counting. Bioinformatics, 31(10):1569–1576, 2015.
  • [7] Garrison E. et al. Variation graph toolkit improves read mapping by representing genetic variation in the reference. Nature Biotechnology, 36(9):875–879, 2018.
  • [8] H. Ferrada, T. Gagie, T. Hirvola, and S. J. Puglisi. Hybrid indexes for repetitive datasets. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372(2016):1–9, 2014.
  • [9] H. Ferrada, D. Kempa, and S.J. Puglisi. Hybrid Indexing Revisited. In Proceedings of the 21st Algorithm Engineering and Experiments (ALENEX), pages 1–8, 2018.
  • [10] P. Ferragina and G. Manzini. Opportunistic data structures with applications. In Proceedings of the 41st Annual Symposium on Foundations of Computer Science (FOCS), pages 390–398, 2000.
  • [11] T. Gagie, G. Navarro, and N. Prezza. Optimal-time text indexing in bwt-runs bounded space. In Proceedings of the 29th Annual Symposium on Discrete Algorithms (SODA), pages 1459–1477, 2018.
  • [12] T. Gagie and S.J. Puglisi. Searching and Indexing Genomic Databases via Kernelization. Frontiers in Bioengineering and Biotechnology, 3:10–13, 2015.
  • [13] Li H. and Durbin R. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 25:1754–60, 2009.
  • [14] L. Huang, V. Popic, and S. Batzoglou. Short read alignment with populations of genomes. 29(13), 2013.
  • [15] M. Jain et al. Nanopore sequencing and assembly of a human genome with ultra-long reads. Nature Biotechnology, 36(4):338–345, April 2018.
  • [16] S. Jeong-Sun et al. De novo assembly and phasing of a Korean human genome. Nature, 538(7624):243–247, October 2016.
  • [17] J. Kärkkäinen, D. Kempa, and S. J. Puglisi. Parallel external memory suffix sorting. In Proceedings of the 26th Annual Symposium on Combinatorial Pattern Matching (CPM), pages 329–342, 2015.
  • [18] B. Langmead and S.L. Salzberg. Fast gapped-read alignment with Bowtie 2. Nature Methods, 9(4):357, March 2012.
  • [19] B. Langmead, C. Trapnell, M. Pop, and S. L. Salzberg. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 10, 2008.
  • [20] S. Levy et al. The Diploid Genome Sequence of an Individual Human. PLoS Biology, 5(10):e254, September 2007.
  • [21] C. Li, R.and Yu, Y. Li, T.-W. Lam, S.-M. Yiu, K. Kristiansen, and J. Wang. Soap2: an improved tool for short read alignment. Bioinformatics, 25(15):1966–1967, 2009.
  • [22] Heng Li. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997 [q-bio], March 2013. arXiv: 1303.3997.
  • [23] S. Maciuca, C. del Ojo Elias, G. McVean, and Z. Iqbal. A natural encoding of genetic variation in a Burrows-Wheeler transform to enable mapping and genome inference. In Proceedings of the 16th Annual Workshop on Algorithms in Bioinformatics (WABI), pages 222–233, 2016.
  • [24] V. Mäkinen, G. Navarro, J. Sirén, and N. Välimäki. Storage and retrieval of highly repetitive sequence collections. Journal of Computational Biology, 17(3):281–308, 2010.
  • [25] Yu Peng, Henry C M Leung, S. M. Yiu, and Francis Y L Chin. IDBA-UD: A de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics, 28(11):1420–1428, 2012.
  • [26] A. Policriti and N. Prezza. LZ77 computation based on the run-length encoded BWT. Algorithmica, 80(7):1986–2011, 2018.
  • [27] K. Schneeberger et al. Simultaneous alignment of short reads against multiple genomes. Genome Biology, 10(9), 2009.
  • [28] L. Shi et al. Long-read sequencing and de novo assembly of a Chinese genome. Nature Communications, 7:12065, June 2016.
  • [29] J. Sirén, N. Välimäki, and V. Mäkinen. Indexing graphs for path queries with applications in genome research. 11(2):375–388, 2014.
  • [30] K.M. Steinberg et al. Single haplotype assembly of the human genome from a hydatidiform mole. Genome Research, page gr.180893.114, November 2014.
  • [31] E.L. Stevens, R. Timme, E.W. Brown, M.W. Allard, E. Strain, K. Bunning, and S. Musser. The public health impact of a publically available, environmental database of microbial genomes. Frontiers in Microbiology, 8:808, 2017.
  • [32] D. Valenzuela and V. Mäkinen. CHIC: a short read aligner for pan-genomic references. BMC Bioinformatics, 19(Suppl 2):87, 2018.
  • [33] D. Valenzuela, T. Norri, N. Välimäki, E. Pitkänen, and V. Mäkinen. Towards pan-genome read alignment to improve variation calling. BMC Genomics, 19(2):87, 2018.
  • [34] S. Wandelt, J. Starlinger, M. Bux, and U. Leser. RCSI: Scalable similarity search in thousand(s) of genomes. Proceedings of the VLDB Endowment, 6(13):1534–1545, 2013.