Efficient construction of the extended BWT from grammar-compressed DNA sequencing reads

Consider a string $T[1..n-1]$ over alphabet $\Sigma[2..\sigma]$, and the sentinel symbol $\Sigma[1]=\texttt{\$}$, which we append at the end of $T$. The suffix array (SA) [21] of $T$ is a permutation of $[n]$ that enumerates the suffixes $T[i..n]$ of $T$ in increasing lexicographic order, $T[SA[i]..n]<T[SA[i+1]..n]$. The BWT [4] is a permutation of the symbols of $T$ obtained by extracting the symbol that precedes each suffix in $SA$, that is, $BWT[i]=T[SA[i]-1]$ (assuming $T[0]=T[n]=\texttt{\$}$). A run-length compressed representation of the BWT [20] adds sublinear-size structures that compute, in logarithmic time, the so-called $\mathsf{LF}$ step and its inverse: if $BWT[j]$ corresponds to $T[i]$ and $BWT[j^{\prime}]$ to $T[i-1]$ (or to $T[n]=\textsf{\$}$ if $i=1$), then $\mathsf{LF}(j)=j^{\prime}$ and $\mathsf{LF}^{-1}(j^{\prime})=j$. Note that $\mathsf{LF}$ regards $T$ as a circular string.

Let $\mathcal{T}=\{T_{1},T_{2},...T_{m}\}$ be a collection of $m$ strings of average size $k$. We then define the string $T[1..n]=T_{1}\texttt{\$}T_{2}\texttt{\$}..T_{n}\texttt{\$}$. The extended BWT (eBWT) of $\mathcal{T}$ [6] regards it as a set of independent circular strings: the BWT of $T$ is slightly modified so that, if $eBWT[j]$ corresponds to $T_{i}[1]$ inside $T$, then $\mathsf{LF}(j)=j^{\prime}$, so that $eBWT[j^{\prime}]$ corresponds to the sentinel $ at the end of $T_{i}$, not of $T_{i-1}$.

LOUDS [15] is a succinct representation that encodes an ordinal tree $T^{\prime}$ with $t$ nodes into a bitmap $B[1..2t+1]$, by traversing its nodes in levelwise order and writing down its arities in unary. The nodes are identified by the position where their description start in $B$. Adding $o(t)$ bits on top of $B$ enables constant-time operations like $\textsf{parent}(u)$ (the parent of node $u$), $\textsf{child}(u,i)$ (the $i$-th child of $u$), $\textsf{psibling}(u)$ (the sibling preceding $u$), $\textsf{nodemap}(u)$ (the level-wise rank of node $u$), $\mathsf{leafrank}(u)$ (the number of leaves in level-order up to leaf $u$), $\mathsf{internalrank}(u)$ (the rank of the internal node $u$ in level-order), and $\mathsf{internalselect}(r)$ (the identifier of the r-th internal node in level order).

Grammar compression consists of encoding $T$ as a small context-free grammar $\mathcal{G}$ that only produces $T$. Formally, a grammar is a tuple $(V,\Sigma,\mathcal{R},\mathqhv{S})$, where $V$ is the set of nonterminals, $\Sigma$ is the set of terminals, $\mathcal{R}$ is the set of replacement rules and $\mathqhv{S}\in V$ is the start symbol. The right-hand side of $\mathqhv{S}\rightarrow C\in\mathcal{R}$ is referred to as the compressed form of $T$. The size of $\mathcal{G}$ is usually measured in terms of $r=|\mathcal{R}|$; the number of rules, $g$; the sum of the length of the right-hand sides of $\mathcal{R}$, and $c=|C|$; the size of the compressed string. The more repetitive $T$ is, the smaller are these values.

Unlike the regular BWT, the position of each $C[j]$ in the eBWT does not depend on the whole suffix $C[j+1..c]$, but on the string $S=C[j+1..j+p^{\prime}]C[j-p..j]$. This sequence is a circular permutation of the compressed string $T_{x}\in\mathcal{T}$ encoded in the range $C[j-p..j+p^{\prime}]$. Computing $S$ from the grammar tree $\mathcal{P}$ is simple as $\mathcal{G}$ is string independent (see Definition 2.6). This feature implies that every $T_{x}\in\mathcal{T}$ maps a specific range $C[k..k^{\prime}]$, with $k\leq k^{\prime}$. Therefore, we do not have to deal with border cases. For instance, when the compressed suffix or prefix of $T_{x}$ lies in between two consecutive symbols of $C$.

For constructing the eBWT of $C$ we require $\mathcal{P}$ and $(L^{h},R^{h},f^{h})$, the last triplet produced by $\mathsf{GLex}$. Given the definition of $\mathcal{P}$, we can easily obtain the root child $v$ encoding $C[j]$ as $v=\mathsf{nodeselect}(j+1)$. Once we retrieve $v$, we obtain $C[j]$ with $f^{h}(\mathsf{label}(v))$. For accessing the circular string $S$ to the right of $C[j]$ we define the function $\mathsf{cright}$. This procedure receives as input a position $j\in[1..c]$ and returns another position $j^{\prime}\in[1..c]$ such that $C[j^{\prime}]$ is the circular right context of $C[j]$. We use $\mathsf{cright}$ as the underlying operator for another function, $\mathsf{ccomp}$. This method compares lexicographically two circular permutations located at different positions of $C$. Similarly, we define a function $\mathsf{cleft}$ that returns the circular left context of $C[j]$. We use it to get the eBWT symbols once we sort the circular permutations. To support these operations, we consider the border cases $C[k^{\prime}+1]=C[k]$ and $C[k-1]=C[k^{\prime}]$ for every $T_{x}$. These exceptions require us to include a bitmap $U[1..c]$ that marks as $U[j]=1$ every root child in $\mathcal{P}$ whose recursive expansion is suffixed by a dummy symbol. The functions $\mathsf{cleft},\mathsf{cright}$ and $\mathsf{ccomp}$ are described in Algorithm 2.

We start the computation of the eBWT of $C$ by creating a table $A[1..c]$ with $|R^{h}|$ lexicographical buckets. Then, we scan the children of the root of $\mathcal{P}$ from left to right, and for every node $v$, we store its child rank in the leftmost available cell of bucket $f^{h}(\mathsf{label}(v))$. This process yields a partial sorting of circular permutations of $C$; every bucket $b$ contains the strings that start with symbol $b$. To finish the sorting, we apply a local $\mathsf{quicksort}$ in every bucket, using $\mathsf{ccomp}$ as the comparison function. Notice these calls to $\mathsf{quicksort}$ should be fast as most of the buckets have few elements, and the amount of right contexts we have to access in every comparison is small. Finally, we produce $B^{h}$ by scanning $A$ from left to right and pushing every symbol $f^{h}(\mathsf{label}(\mathsf{nodeselect}(\mathsf{cleft}(A[j])+1)))$ with $j\in[1..|A|]$.

We define a method called $\mathsf{nextBWT}$ for inferring $B^{i-1}$ from $B^{i}$ (Line 8 of Algorithm 1). This procedure requires us to have a mechanism to map a symbol $B^{i}[j]$ to its phrase $F\in\mathcal{D}^{i-1}\cup I^{i-1}$. We support this feature by replacing $f^{i}$ with an inverted function $f^{i}_{inv}$ that receives a rank $b\in R^{i}$ and returns its associated $\mathcal{P}$ label $l\in L^{i}$ (Line 6 of Algorithm 1). Thus, we obtain the grammar tree node $v$ that encodes $F$ with $\textsf{internalselect}(l-\sigma)$. To spell the sequence of $F$ we proceed as follows; we use $\mathcal{P}$ to simulate a pre-order traversal over the subtree of $v$ in the parse tree of $\mathcal{G}$. When we visit a new node $v^{\prime}$, we check if its label $l^{\prime}=\mathsf{label}(v^{\prime})$ belongs to $L^{i-1}$. If that so, then we push $f^{i-1}(l^{\prime})$ to $F$ and skip the parse subtree under $v^{\prime}$. The process stop when we reach $v$ again. We call this process the level-decompression of $F$, or $\mathsf{ldc}(v)=F$.

If we level-decompress all the phrases encoded in $B^{i}$, then we obtain all the symbols of $T^{i-1}$, although not sorted according the eBWT’s definition. The position of each decompressed symbol $F[u]\in R^{i-1}$ in $B^{i-1}$ depends, in most of the cases, only on the suffix $F[u+1..|F|]$, except when this suffix has length less than two (Lemma 2.5). In addition, if two symbols $F[u]$ and $F^{\prime}[u^{\prime}]$, level-decompressed from different positions $B^{i}[j]$ and $B^{i}[j^{\prime}]$ (respectively), are followed by the same suffix in their respective phrases $F$ and $F^{\prime}$, then their relative orders in $B^{i-1}$ only depend on the values of $j$ and $j^{\prime}$. This property is formally defined as follows:

If we generalize Lemma 3.1 to $x\geq 1$ occurrences of $S$, then we can use the following Lemma for building $B^{i-1}$:

Finding the range of $O$ in $B^{i-1}$ requires to sort the distinct suffixes in $\mathcal{D}^{i-1}\cup I^{i-1}$ in lexicographical order. Thus, if $S$ has rank $b$, then $O$ is the b-th range of $B^{i-1}$. We refer to these suffixes as the context strings of the symbols in $B^{i-1}$. The problem, however, is that the suffixes in $\mathcal{D}^{i-1}$ of length one are ambiguous as we cannot assign them a lexicographical order (Lemma 2.5). We solve this problem by extending the occurrences of the ambiguous context strings with the phrases in $\mathcal{D}^{i-1}\cup I^{i-1}$ that occur to their (circular) right in $T^{i-1}$. By doing this extension, we induce a partition over the $O$ list of an ambiguous $S$; the symbols in first range $O_{X}=O[z..z^{\prime}]$ precede the occurrences of the extended phrase $SX$, the symbols in the second range $O_{Y}=O[z^{\prime}+1..z^{\prime\prime}]$ precede the occurrences of another phrase $SY$, and so on. In this way, every segment of $O$ can be placed in some contiguous range of $B^{i-1}$.

We build an FM-index from $B^{i}$ to compute the string extensions, but without including the SA. The idea is that when we extract the occurrence of an ambiguous $S$ from $B^{i}[j]$, we use $\mathsf{LF}^{-1}(j)=j^{\prime}$ to get the position in $B^{i}$ with the symbol that occurs to the right of $B^{i}[j]$ in $T^{i}$. In this way, we concatenate $S$ to the phrase $F$ obtained from level-decompressing the symbol $B^{i}[j^{\prime}]$. Equivalently, we use $\mathsf{LF}$ to find the left symbols of the context strings that are not proper suffixes of their enclosing phrases.

When we sort the distinct context strings, we reference them by using nodes of $\mathcal{P}$. If $S$ has length two an appears as suffix in different right-hand sides of $\mathcal{R}$, then there is a SuffPair internal node in $\mathcal{P}$ for it (nodes $y2,y3$ and $y4$ of Figure 3). Additionally, if $S$ appears as suffix only in one right-hand side, then there is a node (leaf or internal) that we can use as identifier (node $y1$ of Figure 3). Finally, if $S$ is not a proper suffix in the right-hand side of the rule, then we use the internal node in $\mathcal{P}$ for that rule (node $v$ of Figure 3).

We begin $\mathsf{nextBWT}$ by initializing two empty semi-external lists, $Q$ and $Q^{\prime}$. Then, we scan $B^{i}$ from left to right. For every $B^{i}[j]$, we first perform $\mathsf{LF}(j)=j^{\prime}$ to find the position in $B^{i}$ with the symbol that occur to the left of $B^{i}[j]$ in $T^{i}$. Then, we use the functions $f^{i}_{inv}$ to obtain the internal nodes $u$ and $v$ in $\mathcal{P}$ that encode the phrases of $B^{i}[j^{\prime}]$ and $B^{i}[j]$, respectively. Subsequently, we level-decompress the phrase $W\in\mathcal{D}^{i-1}\cup I^{i-1}$ of $u$ and push the pair $(b,v)$ to $Q$, where $b\in R^{i-1}$ is the last symbol of $W$. The next step is to level-decompress the phrase of $v$. During the process, when we reach a node whose label is in $L^{i-1}$, we get its associated symbol $b\in R^{i-1}$ and its right sibling $y$. The next action depends on the value of $l_{y}=\mathsf{label}(y)$. If $l_{y}$ belongs to $L^{i-1}$ and $y$ is not the rightmost child of its parent, then we push the pair $(b,y)$ to $Q$. If, on the other hand, $y$ is the rightmost child of its parent, but $l_{y}$ does not belong to $L^{i-1}$, then we update its value to $y=\mathsf{internalselect}(l_{y}-\sigma)$ and then push $(b,y)$ into $Q$. The purpose of the $\mathsf{internalselect}$ operation is to get the internal node with the first occurrence of $l_{y}$ in $\mathcal{P}$. The last option is that $l_{y}$ belongs to $L^{i-1}$ and $y$ is the rightmost child of its parent. In such situation, we extend the context of $y$. We use $\mathsf{LF}^{-1}(j)=j^{\prime\prime}$ to get the position with the symbol to the right of $B^{i}[j]$. As before, we get the internal node $z$ associated to $B^{i}[j^{\prime\prime}]$ and finally push the triplet $(b,y,z)$ to $Q^{\prime}$.

Once we complete the scan of $B^{i}$, we group the pairs in $Q$ according to $y$, without changing their relative orders. In $Q^{\prime}$ we do the same, but we sort the elements according to $(y,z)$. Subsequently, we form a sorted set $U$ with the distinct $y$ symbols of $Q$ plus the distinct pairs $(y,z)$ of $Q^{\prime}$. To get the relative order of two elements in $U$, we level-decompress their associated phrases and compare them in lexicographical order. For the pairs $(y,z)$, their phrases are obtained by concatenating the level-decompression of $y$ and $z$. Finally, if a given value $y\in Q$ has rank $b$ among the other elements of $U$, then we place the left symbols of its range in $Q$ at the b-th range of $B^{i-1}$. The same idea applies for the pairs $(y,z)$ of $Q^{\prime}$. Algorithm 3 implements $\mathsf{nextBWT}$ and Figure 3 depicts our method for processing $B^{i}[j]$.

The $\mathsf{nextBWT}$ algorithm works provided all the occurrences of a context string $S$ appear as suffixes in the phrases of $\mathcal{D}^{i-1}\cup I^{i-1}$. Still, this is not always the case as sometimes they occur in between nonterminals.

An implicit occurrence appears when $F[1]$ is classified as LMS-type during the parsing. We use the following lemma to detect this situation:

Before running $\mathsf{nextBWT}$, we scan $L^{i}$ from left to right to find the internal nodes of $\mathcal{P}$ that meet Lemma 3.5. For every node $v^{\prime}$ that meets the lemma, we create a pair $p=(l_{l},l_{r})$, with the labels in $\mathcal{P}$ of its left and right children, respectively. Then, we record $p$ in a hash table $\mathcal{H}$, with $v^{\prime}$ as its value. During the execution of $\mathsf{nextBWT}$, when we level-decompress the symbol to the left of $B^{i}[j]$ (Line 8 of Algorithm 3), we check if the pair formed by the grammar tree label of $W[|W|]$ and the label $f^{i}_{inv}(B^{i}[j])$ has an associated value $v^{\prime}$ in $\mathcal{H}$. If that happens, then we insert $(W[|W|-1],v^{\prime})$ to $Q$.

Run-length decompression. As we descend over the iterations of $\mathsf{infBWT}$, the size of $B^{i}$ increases, but its alphabet decreases. This means that in the last steps of the algorithm we end up decompressing a small number of phrases several times. To reduce the time overhead produced by this monotonous decompression pattern, we exploit the runs of equal symbols in $B^{i}$. More specifically, if a given range in $B^{i}[j..j^{\prime}]$ of length $x=j^{\prime}-j+1$ has the same symbol $\mathqhv{X}$, then we level-decompress the associated phrase once and multiply the result by $x$ instead of performing the same operations $x$ times. By Lemma 3.1, we already now that the repeated symbols resulted from this run-length decompression will be placed in consecutive ranges of $B^{i-1}$. Therefore, when we process a given run $(x,\mathqhv{X})$ of $B^{i}$, we do not push every pair $(b,y)$ $x$ times to $Q$, but insert $(x,b,y)$ only once. We do the same for the triplets $(b,y,z)$ of $Q^{\prime}$. Figures 3B and 3C show an example of how the runs of $B^{i}$ are processed.

Unique context strings. Another way of exploiting the repetitions in $B^{i}$ is with the unique context phrases. Consider a nonterminal $\mathqhv{X}$ encoded in the grammar tree node $v$. If a given child $y$ of $v$ has a label in $L^{i-1}$, then the context phrase $S$ decompressed from $y$ only appears under the subtree of $v$. This means that if $\mathqhv{X}$ has $x$ occurrences in $B^{i}$, then $S$ has $x$ occurrences in $T^{i-1}$, all with the same left context symbol. Computing the elements in $B^{i-1}$ with unique context strings is simple; for each label $l\in L^{i}$, we get the internal node $v=\mathsf{internalselect}(l-\sigma)$, and its number of children $a=\mathsf{nchildren}(v)$. If $a>2$, then we count the $x$ occurrences of $f^{i}(l)\in R^{i}$ in $B^{i}$ by calling $\mathsf{rank}$ over the FM-index. Then, for every child $y$ of $v$ whose child rank is in the range $[2..a-1]$, we push the pair $(x,b,y)$ to $Q$, where $b\in R^{i-1}$ is the symbol obtained from the left sibling of $y$. Then, during the scan of $B^{i}$, we just skip these $y$ nodes. In $\mathsf{infBWT}$, we carry out this process before inverting $f^{i}$ (Line 6 of Algorithm 1). In Figure 3A, node $y_{1}$ in the subtree of $v$ encodes a unique context string.

Inferring the eBWT in parallel. Our algorithm for building $Q$ and $Q^{\prime}$ can run in parallel because of Lemma 3.3. For doing so, we divide $B^{i}$ in $p$ chunks and run the for loop between lines 4-33 of Algorithm 3 independently for every one of them. This process yields $p$ $Q$ and $Q^{\prime}$ lists, which we concatenate as $Q=Q_{1}\cup Q_{2},..\cup Q_{p}$ afterward (we do the same with the $Q^{\prime}$ lists). Then we continue with the inference of $B^{i-1}$ in serial.

Sorting the context strings. There is a trade-off between time and space in sorting $U$ (Line 35 of Algorithm 3). On hand side, decompressing the strings once and maintain them in plain form to compare them can produce a considerable space overhead. On the other, decompressing them every time we access them might be slow in practice. We opted for something in between; we only access the strings when we compare them, but we amortize the decompression time overhead by doing the sorting in parallel. Our approach is simple; we use a serial $\mathsf{countingsort}$ to reorder the strings of $U$ by their first symbols. We logically divide the resulting $U$ in buckets; all the strings beginning with the same symbol are grouped in a bucket. To finish the sorting, we just perform $\mathsf{quicksort}$ on parallel in every bucket, decompressing the strings on demand.