IBB: Fast Burrows-Wheeler Transform Construction for Length-Diverse DNA Data

We construct $BWT(W)$ without concatenating $W$ and without using a single-string BWT construction algorithm. Instead, like BCR and ropeBWT2, we use iterations $t=0,\mathrm{\dots },M$. In each iteration, we insert at most one character of each word to compute a partial BWT string $bwt(t)$. Furthermore, like BCR and ropebwt2, we insert the characters of each word in reverse order. However, we do not start to insert the last character of each word $S_j$ in iteration $t=0$. Instead, we begin inserting characters of longer words in earlier iterations with the goal to terminate the insertion after the insertion of the first symbol of each word $S_j$ in the second last iteration $t=M-1$ followed by the insertion of a $\$$ character for each word $S_j$ in the last iteration $t=M$.

As we will process the strings $S_j$ from the end to the beginning, in each iteration $t$, we access the symbols $S_j[M-1-t]$ of strings $S_j$ with $|S_j|+t>M$. In the following, we write the shortcut $W_j[t]$ for $S_j[M-1-t]$. We also set $W_j[M]=\$$ for all words $W_j$. 1(a) shows the words $W_0,\mathrm{\dots },W_8$ already in reverse, right-aligned order during the iterations $t=0,\mathrm{\dots },10$.

Our insertion order has two advantages: First, we keep the partial $bwt(t)$ smaller compared to starting all words in iteration $t=0$, which allows faster computation of $bwt(t+1)$ from $bwt(t)$. Second, we only have to consider the 4 characters $A$, $C$, $G$, and $T$ instead of the 5 characters $A$, $C$, $G$, $T$, and $\$$ in $bwt(t)$ for all $t=0,\mathrm{\dots },M-1$, which allows us to use smaller data structures.

The words $W_0,\mathrm{\dots },W_m$ are not necessarily ordered by their length and swapping their order would change $BWT(W)$. We start inserting a word $W_j$ in iteration $t=M-|W_j|$ and call the word active then. Let the number of currently active words in iteration $t$ be $\alpha (t)$. When we start inserting a word $W_j$, we need to compute the insert position $P_j(M-|W_j|)$ for $W_j$ in iteration $t=M-|W_j|$, because the insertion of shorter words with an index in has not yet started. Our solution to this insert-position-problem uses a bit-vector $SB(t)$ of length $m$. Initially, $SB(0)[j]=0$ for all and we set $SB(t)[j]=1$ when $W_j$ starts in iteration $t=M-|W_j|$. After updating $SB(t)$ in an iteration $t=M-|W_j|$, the insert position $P_j(t)$ of $W_j[t]$ is $ran{k}_{SB(t)}(j,1)$, which counts the number of active words $W_z$ with index $z\in [0,j-1]$. For example, in iteration $t=6$, the word $W_7$ is started. Its insert position shown in 1(b) is $P_7(6)=ran{k}_{SB(6)}(7,1)=2$.

As in Ferragina et al. [9], the relative order of the symbols in $bwt(t)$, $t\ge 0$, is unchanged when inserting the new symbols to get $bwt(t+1)$. Thus, given the last insert position $P_j(t)$ of character $W_j[t]$, $P_j(t+1)$ can be obtained from an adjusted LF-Mapping. The adjustment regards the number of active words $\alpha (t+1)$ in iteration $t+1$ and is necessary, because $bwt(t)$, , is not a valid BWT due to the missing $\$$ symbols. We need to adjust by $\alpha (t+1)$ instead of $\alpha (t)$ because within a non-partial BWT, every active word in iteration $t+1$ would need a $\$$ symbol.

The next insert position is $P_j(t+1)=LF(P_j(t))+\alpha (t+1)$ for $W_j$ in iteration $t+1$ for . In the last iteration $t=M$, we do not calculate any next insert positions, because we terminate after iteration $t=M$. Inserting the next character of each active word $W_j$ into $bwt(t)$ – sorted according to the insert positions $P_j(t)$ – yields $bwt(t+1)$.

In the example of 1(c), the next insert position $P_3(7)$ of $W_3$ in iteration $t+1=7$ given $P_3(6)=8$ and $bwt(6)=CTCCGAACCGCCG$ is

BCR [3] and ropeBWT2 [14] segment $BWT(W)$ into buckets $B_c$ based on the number of a character $c$ in $W$. Formally, $B_c[i]=BWT(W)[i+coun{t}_W(c)]$ and $BWT(W)=B_{\$}{B}_A{B}_C{B}_G{B}_T(B_N)$ in case of ropeBWT2s alphabet that includes $N$ additionally. We call the buckets $B_c$ level-1 buckets and refine these buckets by looking at the first $k$ characters instead of one character $c$. Hereby, $k$ with is a fixed natural number. As shown in 1(c) for $k=2$, we have the level-k buckets $B_{c_1\mathrm{\cdots }{c}_k}$ and we call $c_1\mathrm{\cdots }{c}_k$ the predecessor sequence of the symbol to be inserted into the bucket. We omit special buckets for any predecessor sequence containing $\$$-symbols from $W$ and replace the character $\$$ and all following characters by $A$s for the definition of the buckets, so we use only characters $c_i\in \mathrm{\Sigma }$ as predecessor symbols. The reason for this decision is that $\$$-symbols do not occur within a word $W_i$ and thus, some level-k buckets, for example $A\$C$ for $k=3$, would be empty by construction.

Like BCR [3] and ropeBWT2 [14], we can properly define the buckets of the final $BWT(W)$. For that purpose, we define $A{W}_j=(A^k\cdot W_j)$ as the word $W_j$ prepended by $k$ $A$ symbols before the word $W_j$. We transfer the index $i$ of each symbol in $W_j$ to $A{W}_j$, so $A{W}_j[i]=W_j[i]$ for every word $W_j$. In 1(a), the grey $A$s before the words $W_j$ belong to $A{W}_j$. For each predecessor sequence $D\in {\mathrm{\Sigma }}^k$, let be the number of predecessor sequences that are lexicographically smaller than $D$, obtained from the $k$-length substrings of $A{W}_j$. Then, we get that $BWT(W)=B_{A^k}{B}_{A^{k-1}C}{B}_{A^{k-1}G}{B}_{A^{k-1}T}{B}_{A^{k-2}CA}\mathrm{\cdots }{B}_{T^{k-1}G}{B}_{T^k}$ is a well-defined partition into $4^k$ segments where a bucket is defined as $B_D[i]=BWT(W)[i+C_W(D)]$.

To properly insert the character $W_j[t]$ into its corresponding bucket, we need to know the predecessor sequence $D_j(t)$ for $W_j[t]$ in iteration $t$. Because the BWT can be obtained by taking the last column of the sorted rotations and the buckets are ordered, the predecessor sequence $D_j(t)$ occurs at the start of the rotation of the character $W_j[t]$, if the $\$$ is replaced with $k$ $A$ symbols. Thus, $D_j(t)=A{W}_j[t-1]\mathrm{\cdots }A{W}_j[t-k]$ is the proper predecessor sequence for character $W_j[t]$.

Since $D_j(t)=A{W}_j[t-1]\mathrm{\cdots }A{W}_j[t-k]$, $D_j(t)$ are the last $k$ inserted characters of word $W_j$ with the last inserted character $W_j[t-1]$ first or, if , $D_j(t)$ is prepended by $k-t$ $A$s. Especially, you can get $D_j(t+1)=W_j[t]\cdot D_j(t)[0,k-2]$ by shifting the previous predecessor sequence $D_j(t)$ and placing the last inserted character $W_j[t]$ at index $0$. For example, in 1(b), $D_2(7)=CC$, because we insert $W_2[7-1]=C$ in the $CG=D_2(6)$ bucket, so $D_2(7)=W_2[6]\cdot D_2(6)[0,k-2]=C\cdot (CG)[0,0]=CC$.

The idea is that knowing the bucket $B_{W_j[t-1]}$ to insert into and the insert position $P_j^{W_j[t-1]}(t)$ within that bucket is sufficient for insertion and that the global insert position $P_j(t)$ is not needed for insertion. To calculate the next insert position $P_j(t+1)$, we need $ran{k}_{bwt(t)}(P_j(t),W_j[t])$ and $coun{t}_{bwt(t)}(W_j[t])$ for the LF-mapping. Using buckets and predecessor sequences, we can omit computing (or maintaining) and adding $coun{t}_{bwt(t)}(W_j[t])$: Let $P_j^{W_j[t-1]}(t)$ be the insert position of character $W_j[t]$ where we insert the character $W_j[t]$ into the bucket $B_{W_j[t-1]}$. Bauer et al. [3] prove $P_j^{W_j[t-1]}(t)=ran{k}_{bwt(t-1)}(P_j(t-1),W_j[t-1])$ using a different notation. In case of $W_j[t]=A$, we get the next insert position $P_j^A(t+1)=ran{k}_{bwt(t)}(P_j(t),A)+\alpha (t+1)$ due to our merge of the $\$$ bucket into the $A$ bucket. In other words, $coun{t}_{bwt(t)}(W_j[t])$ is equal to the total number of characters in the buckets of smaller characters . If $W_j[t]=A$, $coun{t}_{bwt(t)}(A)$ is the number of active words $\alpha (t+1)$, which would be the size of the bucket $B_{\$}$ in $bwt(t+1)$ without the merge of the bucket $B_{\$}$ into the bucket $B_A$.

In order to save computation time for rank operations on previous buckets, we introduce four balanced binary trees, A-tree, C-tree, G-tree, and T-tree, of depth $2(k-1)$, as shown in 1(c) for $k=2$. We save the balanced binary trees in an array $T[0,2^{2k}-1]$ of length $4^k=𝒪(|W|)$ using the Eytzinger Layout: The root of the balanced binary trees, A-tree to T-tree, are saved in $T[4]$ to $T[7]$; the left child of a node $T[i]$ is stored at index $2i$ and the right child position at index $2i+1$ in $T$. Each internal node has $|\mathrm{\Sigma }|=4$ counters, where the counter $T[i].c$ associated with $c\in \mathrm{\Sigma }$ represents the number of characters $c$ in the left subtree of node $i$. The leaf nodes are the buckets, which we discuss in Subsection 4.7.

The trees $T[4]$, $T[5]$, $T[6]$, and $T[7]$ represent the level-1 buckets $B_A$, $B_C$, $B_G$, and $B_T$, respectively. We continue to refer to each of these buckets as the concatenation of all level-k buckets of the corresponding tree.

We speed up the computation of $ran{k}_{bwt(t)}(P_j(t),W_j[t])$ by storing additional information in $T[0]$ to $T[3]$. In $T[0]$, we store the total number of characters in the leaves of the A-tree. In $T[1]$, we store the number of characters of the C-tree plus the A-tree, so we have an accumulate count of the number of characters of the whole part of $bwt(t)$ prior to the start of the G-tree. In the same way, $T[2]$ and $T[3]$ represent the total number of characters up to the end of the G-tree and T-tree. 1(c) visualizes the tree structure $T$.

We use $T[0]$ to $T[3]$ to determine the number of characters $W_j[t]$ up to the beginning of tree corresponding to $W_j[t-1]$. Thereby, $ran{k}_{bwt(t)}(P_j(t),W_j[t])$ can be determined by processing only the tree corresponding to $W_j[t-1]$. We define the level-1 rank as $rank1(j,t)=ran{k}_{B_{W_j[t-1]}}(P_j^{W_j[t-1]}(t),W_j[t])$, which is the rank of symbol $W_j[t]$ inserted at its insert position $P_j^{W_j[t-1]}(t)$ on the tree of symbol $W_j[t-1]$.

To get the local insertion position $P_j^{D_j(t)}(t)$ within level-k bucket $B_{D_j(t)}$, we navigate through the balanced binary trees in $T$ using the predecessor sequence $D_j(t)=c_1{c}_2\mathrm{\dots }{c}_k$: While the predecessor symbol $c_1$ of each word determines the selection of one of the 4 binary trees, the predecessor symbols $c_2,\mathrm{\dots },c_k$ of $W_j[t]$ determine the navigation direction inside the $c_1$-tree towards the bucket $B_{D_j(t)}$. Explicitly, for the navigation within the first two levels of the $c_1$-tree, the symbol $c_2$ is used, for the next two levels $c_3$, and so on. We navigate twice left if $c_i=A$; first left and second right if $c_i=C$; first right and second left if $c_i=G$; and twice right if $c_i=T$. If a left navigation step at node $i$ is issued for $c=W_j[t]$, we need to increase $T[i].c$ by one, because we will insert the character $c$ into a leaf in left subtree of node $i$.

If we want to insert symbols from different words at the same time and some symbols have to be inserted into a bucket left of an inner node $i$ while other symbols have to be inserted into a bucket right of the inner node $i$, the left steps must increase $T[i].c$ at a node $i$ before we can proceed with the right steps at node $i$. Otherwise, $T[i].c$ would not be the correct number of occurrences of the character $c$ in the left subtree of node $i$.

For each right step, our goals are to transform the insert position $P_j^{W_j[t-1]}(t)$ into the local insert position $P_j^D(t)$, and to determine $rank1(j,t)$. The number of characters $c$ in the local buckets before the bucket $B_D$ is equal for all words of which characters are inserted into the bucket $B_D$. This equality is true both in the case of every single $c\in \mathrm{\Sigma }$ and for the sum of all characters $c\in \mathrm{\Sigma }$. We use accumulators $R_D.c$ to reduce the number of necessary additions in the right step: For each sequence $D$, $R_D.c$ is the sum of all counters $T[i].c$ where $i$ is a node with a right step issued by $D$. Hereby, $D$ must not be of length $k$ during computation and can contain a $\leftarrow$ or $\to$ symbol at the end to indicate the step of the last half-processed predecessor symbol. We only compute $R_D$ for predecessor sequences $D$ where at least one character $W_j[t]$ has $D$ as predecessor sequence in iteration $t$. Thus, we have at most $4\cdot \min (4^k,m)$ different $R_D.c$ accumulators instead of exactly $4m$ if we would use four variables for each individual word.

Using the accumulator $R_D$ for the character $W_j[t]$, so $D=D_j(t)$, the local insert position of $W_j[t]$ is $P_j^{D_j(t)}(t)=P_j^{W_j[t-1]}(t)-{\sum }_{c\in \mathrm{\Sigma }}{R}_D.c$. Because ${\sum }_{c\in \mathrm{\Sigma }}{R}_D.c$ is again constant for all words that insert into the bucket $B_D$, the sum ${\sum }_{c\in \mathrm{\Sigma }}{R}_D.c$ is only computed once per predecessor sequence $D$. Furthermore, we do not subtract the sum ${\sum }_{c\in \mathrm{\Sigma }}{R}_D.c$ from each $P_j(t)$. Instead, we initialize the position counter with ${\sum }_{c\in \mathrm{\Sigma }}{R}_D.c$. An example for the insertion using $R_D.c$ is shown in Figure 2c.

To compute $rank1(j,t)$ for the insert position $P_j(t+1)$, we use $R_D.c$ with $c=W_j[t]$, which is the number of symbols $c$ up to the bucket $B_D$. Then, we compute the level-k rank $rank\text{k}(j,t)=ran{k}_{B_{D_j(t)}}(P_j^{D_j(t)}(t),W_j[t])$ during insertion into the bucket $B_{D_j(t)}$. Thus, we get the overall expression of the insert position $P_j^{W_j[t]}(t+1)$ for iteration $t+1$ from the iteration $t$ as:

We parallelize the processing of the children each time both children need to be processed, as the further steps of the left and right child do not interfere with each other. We limit the number of threads close to the number of available processors on the used system, as this limit is most efficient.

If symbols from multiple words are inserted into the same bucket $B_D$, we must insert them in correct relative order. We insert the symbols from left to right; for that purpose, we sort the words $W_j$ from which the characters $W_j[t]$ are taken according to their predecessor sequence $D_j(t)$ and then according to their local insert position $P_j^{D_j(t)}(t)$ in bucket $B_{D_j(t)}$. The resulting sort order is equal if we sort according to the first predecessor symbol $D_j(t)[0]=W_j[t-1]$ first and then to the insert position $P_j^{W_j(t-1)}(t)$. We will use the latter ones.

In iteration $0$, all new started words are sorted by their insert position $P_j(0)$ by construction; thus they are also sorted according to their local insert position $P_j^{A^k}(0)$ because $P_j(0)=P_j^A(0)=P_j^{A^k}(0)$. Next, at the start of iteration $t$, we prepend all new active words $W_z$ to the sorted list of active words, because their predecessor sequence $D$ are $A^k$ and their insert positions $P_z^A(t)$ are lower than the insert positions $P_j^A(t)$ of the older active words $W_j$ due to the addition of $\alpha (t+1)$.

For iteration $t+1>0$, we sort the active words from iteration $t$ at the end of iteration $t$ by a single stable radix-sort step on $W_j[t]$. In 1(b), the radix sort sorts the words according to the characters $W_j[6]$ resulting in the same order of active words in iteration $t=7$. The next insert positions $P_j(t+1)$ for one character $W_j[t+1]$ depends only on $ran{k}_{bwt(t)}(P_j(t),W_j[t])$ for a group all words $W_j$ with the same character $c$ at position $t$, or on the constant term $\alpha (t+1)$ in case of $W_j[t]=A$. As the $rank$ is a monotonic non-decreasing function of the position, the sort-order of the group of all $W_j$ with $W_j[t]=c$ is stable. Thus, a single radix-sort step on $W_j[t]$ is sufficient to achieve the correct order of the active words $W_j$ of the iteration $t$.

We store each level-k bucket $B_D$ on external memory and insert into the bucket $B_D$ by alternating two files. We denote the two files by $F_0(B_D)$ and $F_1(B_D)$. Let $L_D(t)\in \{0,1\}$ be the index of the file $F_{L_D(t)}(B_D)$ that we suppose to write to. Initially, we set $L_D(0)=0$ for all $D$. If no word $W_j$ has $D$ as predecessor sequence, we skip the bucket $B_D$. Otherwise, we input $F_{1-L_D(t)}(B_D)$ as a stream. We output the symbols of the input stream to $F_{L_D(t)}(B_D)$ until we reach the next lowest insert position $P_j^D(t)$ of a word $W_j$. We then output $W_j[t]$ before continuing with the input stream and the next lowest insert position $P_z^D(t)$ of some word $W_z$. If no such word $W_z$ with $D$ as predecessor sequence exists, we output the rest of the input stream. Both, the input and output stream, are buffered for performance reasons. At the end, we flip $L_D(t)$ for the next iteration $t+1$ if and only if bucket $B_D$ was changed.

To get $rank\text{k}(j,t)$ for the next insert position $P_j^{W_j[t]}(t+1)$, we count the number of characters $W_j[t]$ that are written to file $F_{L_D(t)}(B_D)$ until the current local insert position $P_j^D(t)$. As $rank\text{k}(j,t)$ is added to the value of $R_D.c$ with $c=W_j[t]$, we can directly increment the counter $R_D.c$ during writing $c$ to the output file $F_{L_D(t)}(B_D)$. Thereby, we get $rank1(j,t)$ from $R_D.c$ with $c=W_j[t]$ just before outputting the symbol $W_j[t]$ at position $P_j^D(t)$.

In Figure 2, the characters $W_2[6]$ and $W_3[6]$ are inserted into $bwt(5)$ in iteration $6$. First, in Figure 2a, we update $T[0]$ to $T[3]$ by increasing $T[0].C$ by $2$ and $T[1].C$, $T[2].C$, and $T[3].C$ by $4$ by using an accumulator $h$ with four counters $h.c$. Second, in Figure 2b, we process the C-tree, because the first predecessor symbol is $D_2(6)[0]=D_3(6)[0]=C$. We do a right and then a left child step because the second predecessor symbol is $D_2(6)[1]=D_3(6)[1]=G$. In Figure 2c, the right step is marked with (1) and the left step is marked with (2). We initialize the accumulators $R_C.c$ with $0$ for all $c$. Because there is no word that requires a left child step at $T[5]$, we omit cloning $R_C$ into $R_{C\leftarrow }$ and just keep $R_C$ as $R_{C\to }$. Additionally, we add the counters of the root of the C tree at $T[5]$ to $R_{C\to }$ to $R_C$, so we properly have the number of characters in the left subtree of $T[5]$ in $R_C$. Then, we process the right child $T[2\ast 5+1]=T[11]$. First, we increase for $W_2$ and $W_3$ $T[11].C=0$ twice by 1 to $2$, because both are left steps. As no word requires a right step at $T[11]$, we can omit that branch. By navigating to the child with index $2\ast 11=22$, we reach a leaf of the tree, thus we insert the characters $W_2[6]$ and $W_3[6]$ to that bucket $B_{CG}$.

In Figure 2c, we perform the insertion into the bucket $CG$. The position counter $i$ is initialized with value $1$, which is the sum of the accumulators $R_C.c$. Because $P_j^C(6)=1=i$, we insert the character $W_2[6]=C$ at the beginning of bucket $B_{CG}$. We increase $R_{CG}.C$ by $1$ for the inserted $C$ of word $W_2$. In Figure 2d, we compute the next insert position $P_2^C(7)$ for $W_2$ before incrementing $R_{CG}.C$. The next insert position $P_2^C(7)=3$ for $W_2$ is the sum of the number of $C$ before the C-tree $T[5]$, which is $T[0].C=3$, and the number of $C$ within the $C$ tree up to that position, which is $R_{CG}.C=0$. In the same way we insert $W_3[6]$ and get the next insert position $P_3^C(7)=4$ for $W_3$ in iteration $t=7$.