Merging RLBWTs Adaptively

The Burrows-Wheeler Transform (BWT) [8] is an important tool in modern genomics. As DNA sequencing technologies have advanced and researchers have started working with pangenomic references, they have begun using run-length compressed BWTs (RLBWTs). There are good algorithms for building RLBWTs for massive pangenomic references in which all the genomes are similar [4, 14] but the case when the reference contains many genomes from each of many species is still challenging. One possible solution is to build a separate RLBWT for each species and then merge them into a single RLBWT or run-length compressed extended BWT (eBWT) [19].

Merging BWTs is an established research topic [11, 12, 24] and merging RLBWTs [23, 17] is a natural extension of that. As far as we know, however, until recently all algorithms for merging RLBWTs were based on dynamic RLBWTs [22, 2] and used time at least roughly linear in the size of the uncompressed inputs, even if they used compressed space. Díaz-Domínguez et al. [9] gave an algorithm for merging BWTs that can be extended to RLBWTs and works faster when the inputs are individually repetitive but not similar to each other. Because it is based on prefix-free parsing, however, it does not have good worst-case bounds. In this paper we show how to merge two RLBWTs into a run-length compressed eBWT in $O(r)$ space and $O((r+L)\log (m+n))$ time, where $m$ and $n$ are the lengths of the uncompressed strings, $r$ is the number of runs in the final eBWT and $L$ is the sum of the longest common prefix (LCP) values at the beginnings of those runs (known as its irreducible LCP values). It is known that $L\in O((m+n)\log \delta )$ [15], where $\delta \le r$ is a powerful measure of the eBWT’s compressibility [16].

We describe here only how to merge the RLBWTs of two strings. Our algorithm can easily be extended to merge many RLBWTs or run-length compressed eBWTs, but we leave that to the full version of the paper. In Section 2 we review some preliminary concepts. As warm-ups, in Section 3 we present simple algorithms for merging positional BWTs (PBWTs) and RLBWTs. In Section 4 we present our main algorithm and in Section 5 we present a first analysis that yields a time bound of $O((r\log r+L)\log (m+n))$. In Section 6 we give the proof of a technical lemma we use, whose details we refer to in Section 7 when we optimize our algorithm slightly to reduce the $r\log r$ in our time bound to $r$, so the whole bound becomes $O((r+L)\log (m+n))$.

For the sake of brevity, we assume readers are familiar with run-length compression, the Burrows-Wheeler Transform (BWT), run-length compressed BWTs (RLBWTs), positional BWTs (PBWTs) [10], suffix arrays (SAs), the $\mathrm{\Psi }$ function and longest common prefix (LCP) arrays; see Mäkinen et al.’s [18] and Navarro’s [20] texts for an introduction. The extended BWT (eBWT) [19] of two strings $S[1..m]$ and $T[1..n]$ is an interleaving ${\mathrm{BWT}}_{S,T}[1..m+n]$ of the characters ${\mathrm{BWT}}_S[1..m]$ and ${\mathrm{BWT}}_T[1..n]$ of $S$ and $T$. Since ${\mathrm{BWT}}_S$ and ${\mathrm{BWT}}_T$ are subsequences of ${\mathrm{BWT}}_{S,T}$, if it has $r$ runs then they each have at most $r$ runs. For convenience, we write ${\mathrm{BWT}}_S$, ${\mathrm{BWT}}_T$ and ${\mathrm{BWT}}_{S,T}$ to denote both the BWTs and the RLBWTs, stating when they are run-length compressed.

For $1\le h\le m+n$, ${\mathrm{BWT}}_{S,T}[h]$ has the lexicographically $h$th smallest context. The contexts of ${\mathrm{BWT}}_S[i]$ and ${\mathrm{BWT}}_T[j]$ are

(unless ${\mathrm{SA}}_S[i]=1$, in which case ${\mathrm{context}}_S(i)$ is just $S[1..m]$, or ${\mathrm{SA}}_T[j]=1$, in which case ${\mathrm{context}}_T(j)$ is just $T[1..n]$), where $\circ$ denotes concatenation. If we assume $S$ and $T$ are each terminated by an end-of-string symbol that occurs nowhere else in $S$ and $T$ then two characters in the same BWT cannot have the same context, and a character in ${\mathrm{BWT}}_S$ and a character in ${\mathrm{BWT}}_T$ have the same context if and only if $S=T$ and the characters are in the same positions (in which case we can break ties arbitrarily without affecting ${\mathrm{BWT}}_{S,T}$).

We can store a move structure [21] (see also [6]) for the $\mathrm{\Psi }$ function ${\mathrm{\Psi }}_S$ of $S$ in $O(r)$ space such that, given ${\mathrm{BWT}}_S$ and position $i$ in it, we can extract ${\mathrm{BWT}}_S[i]$’s context

character by character by iteratively evaluating ${\mathrm{\Psi }}_S$ in $O(\log r)$ time plus constant time per iteration or, equivalently, per character extracted. Brown [5] and Bertram et al. [3] (see also [7]) showed how we can build such a move structure in $O(r\log r)$ time. An analogous result holds for $T$.

It follows that we can compare ${\mathrm{context}}_S(i)$ and ${\mathrm{context}}_T(j)$ in $O(\log r)$ time plus time proportional to the length of their longest common prefix. We give a proof of Lemma 1 in Section 6 and then explain in Section 7 how we can pay the $O(\log r)$ overhead only three times during a binary search to find, for example, the smallest value $j$ with ${\mathrm{context}}_T(j)$ lexicographically larger than ${\mathrm{context}}_S(i)$, denoted ${\mathrm{context}}_T(j)\succ {\mathrm{context}}_S(i)$, rather than paying it at every step of the search.

Suppose we have the PBWTs for two sets of haplotypes from the same species, with the same number of columns (but not necessarily the same number of rows) and with the columns representing the same variation sites, and we want the PBWT for all the haplotypes. For example, consider the PBWTs shown in Figure 1, with the bottom PBWT being the merge of the top two (and the example from [13] with indexing from 1 for consistency with the rest of this paper). The entries of the prefix arrays (showing what haplotype each PBWT entry corresponds to) are shown in parentheses. Cells are red if their information is for one of the first 10 haplotypes and blue if it is for one of the second 10. This makes it easy to see that, by the definition of the PBWT, each column in the merged PBWT is an interleaving of the corresponding columns in the two input PBWTs.

The first column from the left of the merged PBWT is the concatenation of the the first columns of the two input PBWTs. Suppose we have already merged the $(j-1)$st columns of the input PBWTs into the $(j-1)$st column of the merged PBWT, and now we want to merge their $j$th columns into its $j$th column. In particular, suppose we know how the $(j-1)$st column of the merged PBWT is divided into maximal blocks of consecutive bits that all come from the same input PBWT. In our example, the blocks in the 9th column have lengths $2,2,2,1,5,4,1,3$ and alternate between blue and red with the first block being blue.

To compute the $j$th column of the merged PBWT after computing the $(j-1)$st column, we first consider all the 0s in the $(j-1)$st column and list the next bits in their haplotypes (in the 0s’ order in the $(j-1)$st column), then consider all the 1s in the $(j-1)$st column and list the next bits in their haplotypes (in the 1s’ order in the $(j-1)$st column); the $j$th column of the merged PBWT is the concatenation of the two lists. By induction, the bits in the $j$th column are in the co-lexicographic order of the preceding prefixes in their haplotypes.

Because the $j$th columns of the two input PBWTs already contain all those next bits and in the correct order, we can use the input PBWTs instead of the explicit haplotypes. To do this, we start with an empty draft of the $j$th column of the merged PBWT and scan the $(j-1)$st column of the merged PBWT twice. During the first scan, whenever we see a 0 we copy a bit (the next bit in that 0’s haplotype) from the $j$th column of the input PBWT that is the source of that 0. During the second scan, whenever we see a 1 we copy a bit (the next bit in that 1’s haplotype) from the $j$th column of the input PBWT that is the source of that 1. However, our algorithm may work more efficiently considering blocks rather than literally bit by bit.

Suppose the $(j-1)$st column of the merged PBWT consists of $b$ blocks. For $k$ from 1 to $b$, if there are $c$ copies of 0 in the $k$th block in the $(j-1)$st column then we append to our draft of the $j$th column of the merged PBWT the next $c$ bits of the $j$th column of the input PBWT that is the source of the bits in the $k$th block. After that, for $k$ from 1 to $b$, if there are $c$ copies of 1 in the $k$th block of the $(j-1)$st column then we append to our draft of the $j$th column of the merged PBWT the next $c$ bits of the $j$th column of the input PBWT that is the source of the bits in the $k$th block.

In our example, there are no copies of 0s in the first block of the 9th column, which is blue, so we append no bits from the 10th column of the blue PBWT; there is 1 copy of 0 in the second block of the 9th column, which is red, so we append 1 bit (0) from the 10th column of the red PBWT (from haplotype (1)); then 1 bit (0) from the blue PBWT (from haplotype (17)); then 1 bit (1) from the red PBWT (from haplotype (9)); then 3 bits (000) from the blue PBWT (from haplotypes (12), (19) and (18)); then 4 bits (0000) from the red PBWT (from haplotypes (5), (6), (7) and (8)); then 1 bit (0) from the blue PBWT (from haplotype (20)); then 3 bits (000) from the red PBWT (from haplotypes (2), (3) and (4)).

There are 2 copies of 1 in the first block of the 9th column, which is blue, so we append 2 bits (00) from the 10th column of the blue PBWT (from haplotypes (15) and (16)); then 1 bit (0) from the from the red PBWT (from haplotype (10)); then 1 bit (0) from the blue PBWT (from haplotype (11)); then no bits from the red PBWT; then 2 bits (00) from the blue PBWT (from haplotypes (13) and (14)). When we are finished, we know the 10th column of the merged PBWT is 00100000000000000000 and we know its blocks have lengths 1, 1, 3, 4, 1, 3, 2, 1, 3; its prefix array is (1) (17) (9) (12) (19) (18) (5) (6) (7) (8) (20) (2) (3) (4) (15) (16) (10) (11) (13) (14).

If a block in the $(j-1)$st column of the merged PBWT overlaps $t$ runs of bits in that column, then computing the numbers of 0s and 1s in that block takes $O(t)$ time. We perform one append operation for the 0s in that block and another for the 1s. If the input PBWTs are run-length compressed and we want only the $j$th column and its block structure – which is all we need to continue merging the two input PBWTs – but not its prefix array, then each append operation takes time proportional to the number of runs in the substring of bits we append. It follows that we can merge two PBWTs in time $O(r+B)$, where $r$ and $B$ are the total numbers of runs and blocks in all the columns of the merged PBWT, respectively.

Theorem 2 may be of independent interest but we have proven it mainly as a warm-up. As a second warm-up, we now give a simple algorithm for merging RLBWTs that uses $O(r)$ space and $O\left(\text{<span class="ltx\_rule" style="width:0.0pt;height:8.6pt;--ltx-bg-color:black;display:inline-block;"></span>}r\log r+(b\log r+B+\sigma )\log (m+n)\right)$ time, where $r$ is the number of runs in the resulting eBWT, $b$ is the number of maximal blocks of characters in that eBWT that all come from the same input RLBWT, $B$ is the sum of the LCP values at the beginnings of those blocks, $\sigma$ is the size of the alphabet, and $m$ and $n$ are the lengths of the uncompressed input BWTs.

Suppose we are to merge the RLBWTs ${\mathrm{BWT}}_S[1..m]$ and ${\mathrm{BWT}}_T[1..n]$ of strings $S[1..m]$ and $T[1..n]$ into the run-length compressed eBWT ${\mathrm{BWT}}_{S,T}[1..m+n]$ of $S$ and $T$ together, and we have already built the data structures for $\mathrm{\Psi }$ for $S$ and $T$ described in Lemma 1 so that we can compare characters’ contexts. For example, if $S$ is the concatenation of the $-terminated strings GATTACAT, GATACAT and GATTAGATA and $T$ is the concatenation of those three strings’ $-terminated reverse complements, then

with exponents indicating run lengths, and ${\mathrm{BWT}}_S$ and ${\mathrm{BWT}}_T$ are interleaved in ${\mathrm{BWT}}_{S,T}$ as shown in Figure 2. In this toy example, $r=27$, $b=18$, $B=18$, $\sigma =5$ and $m+n=54$.

We first set $i$ and $j$ to 1, then compare ${\mathrm{context}}_S(i)$ and ${\mathrm{context}}_T(j)$ to check whether the first block of characters in ${\mathrm{BWT}}_{S,T}$ is from ${\mathrm{BWT}}_S$ or ${\mathrm{BWT}}_T$, and set a Boolean flag to true in the former case and false in the latter. As long as $i\le m$ and $j\le n$, if the flag is true then we use a doubling search in ${\mathrm{BWT}}_S[i..m]$ to find the largest value $i^{\prime }$ with ${\mathrm{context}}_S(i^{\prime })\prec {\mathrm{context}}_T(j)$, output ${\mathrm{BWT}}_S[i..i^{\prime }]$, reset $i$ to $i^{\prime }+1$, and reset the flag to false. If the flag is false then we use a doubling search in ${\mathrm{BWT}}_T[j..n]$ to find the largest value $j^{\prime }$ with ${\mathrm{context}}_T(j^{\prime })\prec {\mathrm{context}}_S(i)$, output ${\mathrm{BWT}}_T[j..j^{\prime }]$, reset $j$ to $j^{\prime }+1$, and reset the flag to true. Eventually, when $i=m+1$ or $j=n+1$, we output either ${\mathrm{BWT}}_T[j..n]$ or ${\mathrm{BWT}}_S[i..m]$, respectively. The pseudocode for our algorithm is shown in Algorithm 1.

By induction, during each pass through the while loop, we output the next block of ${\mathrm{BWT}}_{S,T}$, from ${\mathrm{BWT}}_S$ if the flag is true and from ${\mathrm{BWT}}_T$ if the flag is false. If the flag is true then all of the $O(\log m)$ comparisons we perform during the current pass are against ${\mathrm{context}}_T(j)$, which is in the block from ${\mathrm{BWT}}_T$ we will output during the next pass, so the number of characters we extract for each of those comparisons is at most 1 plus the maximum of the LCP values at the boundaries of that block from ${\mathrm{BWT}}_T$. If the flag is false then all of the $O(\log n)$ comparisons we perform during the current pass are against ${\mathrm{context}}_S(i)$, which is in the block from ${\mathrm{BWT}}_S$ we will output during the next pass, so the number of characters we extract for each of those comparisons is at most 1 plus the maximum of the LCP values at the boundaries of the block from ${\mathrm{BWT}}_S$.

It follows that we can charge all the comparisons to blocks of ${\mathrm{BWT}}_{S,T}$ such that each block is charged $O(\log (m+n))$ comparisons and for each comparison charged to a block, the number of characters we extract is at most 1 plus the maximum of the LCP values at the boundaries of that block. Since at most $\sigma$ of those LCP values are 0, this means we use $O\left(\text{<span class="ltx\_rule" style="width:0.0pt;height:8.6pt;--ltx-bg-color:black;display:inline-block;"></span>}r\log r+(b\log r+B+\sigma )\log (m+n)\right)$ total time, including the $O(r\log r)$ time to build the data structures for comparing contexts and the $O(\log r)$ overhead for each of the $O(b\log (m+n))$ comparisons. Our algorithm is adaptive in the sense that it takes advantage of cases when $b$ and $B$ are small, which is likely when $S$ and $T$ are individually repetitive but dissimilar.

Again, suppose we are to merge the RLBWTs ${\mathrm{BWT}}_S[1..m]$ and ${\mathrm{BWT}}_T[1..n]$ of strings $S[1..m]$ and $T[1..n]$ into the run-length compressed eBWT ${\mathrm{BWT}}_{S,T}[1..m+n]$ of $S$ and $T$ together, and we have already built the data structures for $\mathrm{\Psi }$ for $S$ and $T$ described in Lemma 1 so that we can compare characters’ contexts. The pseudocode for our algorithm is shown in Algorithm 2.

Assume we have already merged ${\mathrm{BWT}}_S[1..i-1]$ and ${\mathrm{BWT}}_T[1..j-1]$ into ${\mathrm{BWT}}_{S,T}[1..i+j-2]$ correctly. When we first reach the while loop (lines 3–43 in Algorithm 2) $i=j=1$ so this assumption is trivially true. If $i>m$ then ${\mathrm{BWT}}_{S,T}[i+j-1..m+n]={\mathrm{BWT}}_T[j..n]$ (line 45) and if $j>n$ then ${\mathrm{BWT}}_{S,T}[i+j-1..m+n]={\mathrm{BWT}}_S[i..m]$ (line 47), so assume $i\le m$ and $j\le n$ and consider only the while loop. Let $k$ and $\mathrm{\ell }$ be the lengths of the leading runs of ${\mathrm{BWT}}_S[i..m]$ and ${\mathrm{BWT}}_T[j..n]$ respectively (lines 4 and 5).

First suppose ${\mathrm{BWT}}_S[i]={\mathrm{BWT}}_T[j]$ (lines 7–21), meaning the leading runs ${\mathrm{BWT}}_S[i..i+k-1]$ and ${\mathrm{BWT}}_T[j..j+\mathrm{\ell }-1]$ of ${\mathrm{BWT}}_S[i..m]$ and ${\mathrm{BWT}}_T[j..n]$ consist of copies of the same character. If there exists a character ${\mathrm{BWT}}_S[i+k]$ after the run ${\mathrm{BWT}}_S[i..i+k-1]$ and ${\mathrm{context}}_S(i+k)\prec {\mathrm{context}}_T(j+\mathrm{\ell }-1)$ (lines 8–11) then ${\mathrm{BWT}}_S[i+k]$ precedes some non-empty suffix ${\mathrm{BWT}}_T[j^{\prime }..j+\mathrm{\ell }-1]$ of the run ${\mathrm{BWT}}_T[j..j+\mathrm{\ell }-1]$ in ${\mathrm{BWT}}_{S,T}$, with $j^{\prime }$ being the smallest value with ${\mathrm{context}}_T(j^{\prime })\succ {\mathrm{context}}_S(i+k)$. We use binary search to find $j^{\prime }$ in the range $[j,j+\mathrm{\ell })$ (line 8), output

(line 9) and reset $i$ to $i+k$ (line 10) and $j$ to $j^{\prime }$ (line 11). Since the characters we output are all equal, we need not consider their order to merge the RLBWTs. (We note, however, that we would need to consider the characters’ order if we were also merging the RLBWTs’ SAs.) The case when there exists a character ${\mathrm{BWT}}_T[j+\mathrm{\ell }]$ after the run ${\mathrm{BWT}}_T[j..j+\mathrm{\ell }-1]$ and ${\mathrm{context}}_T(j+\mathrm{\ell })\prec {\mathrm{context}}_S(i+k-1)$ (lines 13–16) is symmetric.

By the definition of the eBWT, it is impossible for ${\mathrm{BWT}}_S[i+k]$ to precede a non-empty suffix of the run ${\mathrm{BWT}}_T[j..j+\mathrm{\ell }-1]$ in ${\mathrm{BWT}}_{S,T}$ and for ${\mathrm{BWT}}_T[j+\mathrm{\ell }]$ to precede a non-empty suffix of the run ${\mathrm{BWT}}_S[i..i+k-1]$ in ${\mathrm{BWT}}_{S,T}$ in ${\mathrm{BWT}}_{S,T}$, since then ${\mathrm{BWT}}_S[i+k]$ would have to both precede and follow ${\mathrm{BWT}}_T[j+\mathrm{\ell }]$. Therefore, the only remaining case with ${\mathrm{BWT}}_S[i]={\mathrm{BWT}}_T[j]$ is when either $i+k=m+1$ or ${\mathrm{BWT}}_S[i+k]$ follows the entire run ${\mathrm{BWT}}_T[j..j+\mathrm{\ell }-1]$ in ${\mathrm{BWT}}_{S,T}$ and either $j+\mathrm{\ell }=n+1$ or ${\mathrm{BWT}}_T[j+\mathrm{\ell }]$ follows the entire run ${\mathrm{BWT}}_S[i..i+k-1]$ in ${\mathrm{BWT}}_{S,T}$. In this case, we output

(line 18) and reset $i$ to $i+k$ (line 19) and $j$ to $j+\mathrm{\ell }$ (line 20).

Now suppose ${\mathrm{BWT}}_S[i]\ne {\mathrm{BWT}}_T[j]$ (lines 23–42). Because $S[{\mathrm{SA}}_S[i]-1]={\mathrm{BWT}}_S[i]$ (or $S[m]={\mathrm{BWT}}_S[i]$ if ${\mathrm{SA}}_S[i]=1$ so $S[1..{\mathrm{SA}}_S[i]-1]$ is empty) and $T[{\mathrm{SA}}_T[j]-1]={\mathrm{BWT}}_T[j]$ (or $T[n]={\mathrm{BWT}}_T[j]$ if ${\mathrm{SA}}_T[j]=1$ so $T[1..{\mathrm{SA}}_T[j]-1]$ is empty), ${\mathrm{context}}_S(i)$ and ${\mathrm{context}}_T(j)$ differ on their last characters. Therefore, we consider only two cases: when ${\mathrm{context}}_S(i)\prec {\mathrm{context}}_T(j)$ (lines 24–31) and when ${\mathrm{context}}_S(i)\succ {\mathrm{context}}_T(j)$ (lines 33–40).

Suppose ${\mathrm{context}}_S(i)\prec {\mathrm{context}}_T(j)$, meaning ${\mathrm{BWT}}_S[i]$ precedes ${\mathrm{BWT}}_T[j]$ in ${\mathrm{BWT}}_{S,T}$. If ${\mathrm{context}}_S(i+k-1)\prec {\mathrm{context}}_T(j)$ then the entire run ${\mathrm{BWT}}_S[i..i+k-1]$ precedes ${\mathrm{BWT}}_T[j]$ in ${\mathrm{BWT}}_{S,T}$, so we output that run (line 25) and reset $i$ to $i+k$ (line 26). Otherwise, some non-empty prefix ${\mathrm{BWT}}_S[i..i^{\prime }-1]$ of the run ${\mathrm{BWT}}_S[i..i+k-1]$ precedes ${\mathrm{BWT}}_T[j]$ in ${\mathrm{BWT}}_{S,T}$ and some non-empty suffix ${\mathrm{BWT}}_S[i^{\prime }..i+k-1]$ follows it, with $i^{\prime }$ being the smallest value with ${\mathrm{context}}_S(i^{\prime })\succ {\mathrm{context}}_T(j)$. We use binary search (line 28) to find $i^{\prime }$ in the range $(i,i+k-1]$, comparing the context of a character in ${\mathrm{BWT}}_S[i+1..i+k-1]$ to ${\mathrm{context}}_T(j)$ at each step. We output ${\mathrm{BWT}}_S[i..i^{\prime }-1]$ (line 29) and reset $i$ to $i^{\prime }$ (line 30). The case when ${\mathrm{context}}_S(i)\succ {\mathrm{context}}_T(j)$ is symmetric.

We exit the while loop only after having consumed all of at least one of ${\mathrm{BWT}}_S$ and ${\mathrm{BWT}}_T$. If we have not consumed all of the other, we output the rest of it (lines 44 to 48).

Building the structures in Lemma 1 takes $O(r\log r)\subseteq O(r\log (m+n))$ time, where $r$ is the number of runs in ${\mathrm{BWT}}_{S,T}$. Since ${\mathrm{BWT}}_S$ and ${\mathrm{BWT}}_T$ both consist of at most $r$ runs, outputting ${\mathrm{BWT}}_S[i..m]$ (line 45) or ${\mathrm{BWT}}_T[j..n]$ (line 47) in run-length compressed form takes $O(r)$ time. Therefore, we can focus on the complexity of the while loop.