Abstract

1 INTRODUCTION

FM-index plays an important role in DNA sequence alignment, de novo assembly (Simpson and Durbin, 2012) and compression (Cox et al., 2012). Fast and lightweight construction of FM-index for a large dataset is the key to these applications. In this context, a few algorithms (Bauer et al., 2013; Liu et al., 2014) have been developed that substantially outperform earlier algorithms. However, they are only efficient for short reads. A fast and practical algorithm for long sequence reads is still lacking. This work aims to fill this gap.

2 METHODS

Let Σ = {A, C, G, T, N} be the alphabet of DNA with a lexicographical order A < C < G < T < N. Each element in Σ is called a symbol and a sequence of symbols called a string over Σ. Given a string P, |P| is its length and P[i] the symbol at position i. A sentinel $ is smaller than all the other symbols. For simplicity, we let P[−1] = P[|P|] = $. We also introduce $\stackrel{\sim }{P}$ as the reverse of P and $\overline{P}$ as the reverse complement of P.

Given a list of strings over Σ, (P_i)_0≤i<m, let T = P₀$₀ … P_{m − 1}$_{m − 1} with $₀ <  ⋯  <  $_m−1 < A < C < G < T < N. The suffix array of T is an integer array S such that S(i), 0 ≤ i < |T|, is the starting position of the i-th smallest suffix in the collection T. The Burrows-Wheeler Transform, or BWT, of T can be computed as B[i] = T[S(i)−1]. For the description of the algorithm, we segment B into B = B_$B_AB_CB_GB_TB_N, where B_a[i] = B[i + C(a)] with C(a) = |{j:T[j] < a}| being the array of accumulative counts. By the definition of suffix array and BWT, B_a consists of all the symbols with their next symbol in T being a.

The above defines BWT for an ordered list of strings. We next seek to define BWT for an unordered set of strings 𝒞 by imposing an arbitrary sorting order on 𝒞. We say list (P_i)_i is in the reverse lexicographical order or RLO, if ${\stackrel{\sim }{P}}_i\le {\stackrel{\sim }{P}}_j$ for any i < j; say it is in the reverse-complement lexicographical order or RCLO, if ${\overline{P}}_i\le {\overline{P}}_j$ for any i < j. The RLO-BWT of 𝒞, denoted by B^RLO(𝒞), is constructed by sorting strings in 𝒞 in RLO and then applying the procedure in the previous paragraph on the sorted list. RCLO-BWT B^RCLO(𝒞) can be constructed in a similar way. In $B^{\text{RCLO}}({\left\{P_i\right\}}_i\cup {\left\{{\overline{P}}_j\right\}}_j)$, the k-th smallest sequence is the reverse complement of the k-th sequence in the FM-index. This property removes the necessity of keeping an extra array to link the rank and the position of a sequence in the FM-index, and thus helps to reduce the memory of some FM-index–based algorithms (Simpson and Durbin, 2012). For short reads, RLO/RCLO-BWT is also more compressible (Cox et al., 2012).

As a preparation, we further define two string operations: rank(c, k; B) and insert(c, k; B), where rank(c, k; B) = |{i < k:B[i] = c}| gives the number of symbols c before the position k in B, and insert(c, k; B) inserts symbol c after k symbols in B with all the symbols after position k shifted to make room for c. We implemented the two operations by representing each B_c in a B+-tree in memory, where a leaf keeps a run-length encoded string and an internal node keeps the count of each symbol in the leaves descended from the node.

Algorithm 1 appends a string to an existing index by inserting each of its symbol from the end of P. It was first described by Chan et al. (2004). Algorithm 2 constructs RLO/RCLO-BWT in a similar manner to Algorithm 1 except that it inserts P[i] to [l, u), the suffix array interval of P’s suffix starting at i + 1, and that BWT symbols in this interval are already sorted. This process implicitly applies a radix sort from the end of P, sorting it into the existing strings in the BWT in RLO/RCLO. Note that if we change line 1 to “l ← u ← |{i:B[i] = $}|”, Algorithm 2 will be turned into Algorithm 1. Recall that the BCR algorithm (Bauer et al., 2013) is, to some extent, the multi-string version of Algorithm 1. Following similar reasoning, we can extend Algorithm 2 so as to insert multiple strings at the same time, which gives Algorithm 3. We use an array A(j) to keep the state of the j-th sequence after inserting its d-long suffix. At line 2, A(j).c is the previously inserted symbol and [A(j).l, A(j).u) is the interval to which the new symbol is inserted. In implementation, we may speed up the sorting mode by inserting multiple symbols at line 3.

When B is represented by a balanced tree structure, the time complexity of all three algorithms is O(n log⁡n), where n is the total number of symbols in the input. However, we will see later that for short strings, Algorithm 3 is substantially faster than the first two algorithms, due to the locality of memory accesses, the possibility of cached B+-tree update and the parallelization of the ‘for’ loop at line 1. These techniques are more effective for a larger batch of shorter strings.

Disregarding RLO/RCLO, Algorithm 3 is similar to BCR except that BCR keeps B in monolithic arrays. As a result, the time complexity of BCR is O(nl), where l is the maximum length of reads, not scaling well to l.

3 RESULTS AND DISCUSSION

We implemented the algorithm in ropeBWT2 and evaluated its performance together with BEETL (http://bit.ly/beetlGH), the original on-disk implementation of BCR and BCRext, ropeBWT-BCR (https://github.com/lh3/ropebwt), an in-memory reimplementation of BCR by us, and NVBio (http://bit.ly/nvbioio), a GPU-based algorithm inspired by CX1 (Liu et al., 2014). Table 1 shows that for ~100 bp reads, ropeBWT2 has comparable performance to others. For the ~875 bp Venter dataset, NVBio aborted due to insufficient memory under various settings. We did not apply BCR because it is not designed for long reads of unequal lengths. Only ropeBWT2 works with this data set and the even longer moleculo reads.

Funding: NHGRI U54HG003037; NIH GM100233.

Conflict of Interest: none declared.

REFERENCES