Optimal-Time Mapping in Run-Length Compressed PBWT

The Positional Burrows–Wheeler Transform (PBWT) [5] is a data structure designed for efficiently representing and querying large collections of sequences, such as haplotype panels in genomics. Originally proposed by Durbin [5], the PBWT stores a set of $0pt$ haplotypes across $0pt$ variant sites in a $0pt\times 0pt$ binary matrix, where the rows at each column $j$ are arranged in co-lexicographic order according to the prefixes of the haplotypes up to column $j-1$. This ordering facilitates efficient querying between a given haplotype and the panel, enabling the identification of set-maximal exact matches (SMEMs).

Two fundamental operations are defined on the PBWT: forward stepping and backward stepping. Forward stepping, denoted by $\text{fore}[i][j]$, maps the position of a haplotype in the permutation of the PBWT at site $j$ to its position in the permutation at site $j+1$, while backward stepping, denoted by $\text{back}[i][j]$, performs the inverse mapping, tracing a haplotype’s position at site $j+1$ back to its position at site $j$. These operations play an essential role in many PBWT-based algorithms for haplotype matching and analysis.

Although the PBWT enables efficient computations, its memory usage grows rapidly with large haplotype datasets, posing a challenge for population-scale cohorts like the UK Biobank [8]. To address this limitation, the $\mu$-PBWT [4], also discussed in [2], was introduced as a compressed variant of the PBWT that leverages run-length encoding (RLE) to reduce space usage. A run in a sequence is defined as a maximal contiguous block of identical symbols, and the $\mu$-PBWT^′s storage requirement is $O(\tilde{r})$ words, where $\tilde{r}$ denotes the total number of runs in the corresponding PBWT across all $0pt$ sites. This run-length compression reduces the $\mu$-PBWT’s memory usage to only a fraction of that of the original PBWT, making it feasible to index and query massive haplotype panels such as those in the UK Biobank [8].

Beyond space efficiency, the $\mu$-PBWT’s structure has proven particularly useful for downstream applications such as computing matching statistics [4, 3], Minimal Positional Substring Covers (MPSC) [3], and SMEMs [4, 3]. Each of these applications relies on repeated forward and backward stepping operations. However, within the $O(\tilde{r})$-space bound of the $\mu$-PBWT, no data structure supporting both operations in constant time was known prior to our work. Achieving faster stepping operations would directly improve the performance of all these applications and further enhance the efficiency of run-length compressed PBWT indexing for large-scale genomic analyzes.

Durbin [5] noted that, due to the co-lexicographic ordering of the PBWT, the forward (resp, backward) stepping operation is the natural analogue of the LF- (resp, FL-) mapping in the classical Burrows–Wheeler Transform (BWT). Gagie et al. [7, Lemma 2.1] presented a data structure requiring $O(r)$ words of space, built over a text of length $n$, that supports LF- and FL-mappings on the BWT in $O(\log {\log }_w(n/r))$ time, where $r$ is the number of runs in the BWT of the text and $w$ is the number of bits in a machine word. By applying this data structure to each column of the PBWT, we can achieve, for every site $j$, an $O(r_j)$-word data structure that supports both forward and backward stepping queries in $O(\log {\log }_w(0pt/r_j))$ time, where $r_j$ denotes the number of runs at site $j$.

Nishimoto and Tabei [12] recently proposed the move structure to accelerate LF-mapping operations on the BWT. Their approach first divides the $r$ runs in the BWT into at most $2r$ sub-runs, and then constructs an $O(r)$-word data structure over these sub-runs. Given a position $i$ in the BWT and the index of the sub-run containing $i$, the structure can compute LF$(i)$ and determine the index of the sub-run containing LF$(i)$ in constant time.

Typically, in applications of the PBWT, forward and backward stepping are performed iteratively. For example, to traverse a haplotype in the PBWT from site $j$ to site $j^{\prime }$, the forward stepping procedure is invoked sequentially $j^{\prime }-j$ times. However, efficient (e.g., constant-time) forward and backward stepping in the run-length encoded PBWT cannot be achieved simply by applying the move structure independently at each site. Indeed, the move structure only returns the index of the sub-run at site $j$ that contains $\text{fore}[i][j]$, whereas continuing the forward step from site $j+1$ to site $j+2$ requires the index of the sub-run at site $j+1$ that contains $\text{fore}[i][j]$. See Figure 1 for an example. Hence, new methods are needed to accelerate fore and back queries while maintaining the run-length encoded PBWT.

Two key applications of the back and fore queries on the PBWT are haplotype retrieval and prefix search [6]. In haplotype retrieval, the goal is to extract any haplotype from the PBWT. In prefix search, the objective is to find the longest common prefix $P[1..m^{\prime }]$ between any haplotype and a query pattern $P[1..m]$, and to list the indices of all haplotypes prefixed by $P[1..m^{\prime }]$. Using the data structure from [7, Lemma 2.1], built over each column of the PBWT, haplotype retrieval can be performed in $O(0pt\log {\log }_{w}h)$ time by invoking forward stepping iteratively $0pt$ times, with an overall space usage of $O(\tilde{r})$ words. Gagie et al. [6] present a textbook solution based on the PBWT that uses $O(0pt\cdot 0pt)$ words of space and supports each prefix query in $O(m\log 0pt+\text{occ})$ time, where occ denotes the number of reported indices. They also provide various time-space tradeoffs in the same work. Nishimoto and Tabei [12] propose a data structure for prefix search that combines a BWT built over all haplotypes, a compact trie [10], the move data structure, and at most $20pt-1$ marked positions. The structure supports computing the longest prefix $P[1..m^{\prime }]$ in $O(m^{\prime })$ time, using $O(r^{\ast }+0pt)$ words of space, where $r^{\ast }$ denotes the number of runs in the BWT of the concatenated sequence.

We consider the multi-allelic PBWT as described in [11]. We propose a “move data structure”$-$style solution specifically tailored for the PBWT. Our approach introduces a simple yet effective algorithm that partitions the $\tilde{r}$ runs of the PBWT into at most $2\tilde{r}$ sub-runs across all $0pt$ columns of the PBWT (Section 4). We then design an $O(\tilde{r})$-word data structure over the sub-runs (Theorem 13) that supports $O(1)$-time computation of both forward and backward stepping operations, allowing to be invoked iteratively (Section 5).

As a first application, we present a PBWT-based solution for prefix searches. Specifically, assuming all haplotypes have the same length, we design an $O(\tilde{r}+0pt)$-word data structure that finds the longest common prefix $P[1..m^{\prime }]$ between any haplotype and a query pattern $P[1..m]$ in $O(m^{\prime }\log {\log }_w\sigma +\text{occ})$ time. When the haplotypes are sorted lexicographically, we can further reduce the space complexity from $O(\tilde{r}+0pt)$ to $O(\tilde{r})$. Note that the problem setting considered by Nishimoto and Tabei [12] allows haplotypes of varying lengths. We also show that our PBWT-based solution can be adapted to this more general case. While our query time is not as efficient as that of Nishimoto and Tabei [12], our approach has the advantage of not requiring a compact trie. As a second application, we design an $O(\tilde{r})$-word data structure to represent the haplotypes, supporting retrieval of any haplotype $S_i$ ($1\le i\le 0pt$) in $O(\log {\log }_{w}0pt+0pt)$ time (Theorem 17).

Before introducing our “move”-like $O(\tilde{r})$-word data structure, we complement our results with lower and upper bounds on the measure $\tilde{r}$. We establish a connection between $\tilde{r}$ and the number $0p{t}^{\prime \prime }$ of adjacent haplotype pairs $(S_i,S_{i+1})$ such that $S_i\ne S_{i+1}$. Specifically, we show that $\tilde{r}\ge 0p{t}^{\prime \prime }+1$ and $\tilde{r}\le 0pt(0p{t}^{\prime \prime }+1)$.

In Section 2, we introduce the notation and preliminary results used throughout the paper. In Section 3, we establish lower and upper bounds on the size $\tilde{r}$. Section 4 defines the three-overlap constraint over two lists of intervals and presents an algorithm that divides intervals in one list into sub-intervals satisfying this constraint. Applying this algorithm, we describe in Section 5 how to divide PBWT runs into sub-runs for back and fore queries. In Section 6, we design data structures over these sub-runs and develop the corresponding algorithms for back and fore queries. Section 7 discusses two applications of back and fore queries: haplotype retrieval and prefix searches. Finally, Section 8 concludes the paper and outlines directions for future work.

We denote by $[i,j]$ the interval of integers $i,i+1,\mathrm{\dots },j$, and define $[i,j]=\mathrm{\varnothing }$ if . Given an interval $𝖨=[i,j]$, we denote its left and right endpoints by $𝖨.b$ and $𝖨.e$, respectively, so that $𝖨.b=i$ and $𝖨.e=j$. For any matrix $A$, we denote by ${\text{col}}_j(A)$ its $j$-th column and by ${\text{col}}_j(A)[i]$ the entry $A[i][j]$. Let $M$ be the matrix that stores $S_1,S_2,\mathrm{\dots },S_{0}pt$ in rows $1,2,\mathrm{\dots },0pt$. We note that the input haplotypes might not be pairwise distinct.

Consider an ordered alphabet $\{0,\mathrm{\dots },\sigma -1\}$ with . A string $\alpha$ over this alphabet is a finite sequence of symbols from $\{0,\mathrm{\dots },\sigma -1\}$, that is, $\alpha =\alpha [1]\alpha [2]\mathrm{\dots }\alpha [|\alpha |]$, where $|\alpha |$ denotes the length of $\alpha$, and $\alpha [i]\in \{0,\mathrm{\dots },\sigma -1\}$ for all $1\le i\le |\alpha |$. The empty string is denoted by $\epsilon$. For indices $1\le i\le j\le |\alpha |$, we denote by $\alpha [i..j]$ the substring of $\alpha$ spanning positions $i$ through $j$ (and define $\alpha [i..j]=\epsilon$ if $i>j$). The string $\alpha [1..i]$ is referred to as the $i$-th prefix of $\alpha$, and $\alpha [i..|\alpha |]$ as the $i$-th suffix, for $1\le i\le |\alpha |$. A proper prefix (respectively, proper suffix) of $\alpha$ is a prefix (respectively, suffix) $\beta$ such that $\beta \ne \alpha$.

Given two strings $\alpha$ and $\beta$ over the alphabet set $\{0,\mathrm{\dots },\sigma -1\}$, we say that $\alpha$ is lexicographically smaller than $\beta$, denoted by $\alpha \prec \beta$, if and only if one of the following holds: i) there exists an index $k$ such that $\alpha [i]=\beta [i]$ for all , and ; ii) or $\alpha$ is a proper prefix of $\beta$. We say that $\alpha$ is co-lexicographically smaller than $\beta$, denoted by $\alpha {\prec }_{colex}\beta$, if and only if one of the following holds: i) There exists an index $k$ such that $\alpha [|\alpha |-i+1]=\beta [|\beta |-i+1]$ for all , and ; ii) or $\alpha$ is a proper suffix of $\beta$.

Given a sorted list $S$ of integers, a predecessor query takes an integer $x$ as input and returns $\max \{y\in S\mid y\le x\}$ and its rightmost position in $S$.

Given a sequence $A[1..n^{\prime }]$ over an alphabet $\{0,\mathrm{\dots },\sigma -1\}$, the operation ${\text{rank}}_c(A,j)$ returns the number of occurrences of $c$ in $A[1..j]$, for $c\in \{0,\mathrm{\dots },\sigma -1\}$. The operation ${\text{select}}_c(A,j)$ returns the position of the $j$-th occurrence of $c$ in $A$, for $1\le j\le {\text{rank}}_c(A,n^{\prime })$, and returns $n^{\prime }+1$ if $j>{\text{rank}}_c(A,n^{\prime })$.

Closely related to the PBWT is the Prefix Array (PA), a matrix that records, for each column, the permutation of haplotype indices induced by the PBWT. Formally, the Prefix Array PA built for the matrix $M$ is an $0pt\times 0pt$ matrix, in which ${\text{col}}_1(\text{PA})$ is simply the list $1,2,\mathrm{\dots },0pt$, and ${\text{col}}_j(\text{PA})$, for $j>1$, stores the permutation of the set $\{1,\mathrm{\dots },0pt\}$ induced by the co-lexicographic ordering of prefixes of $\{S_1,\mathrm{\dots },S_{0}pt\}$ up to column $j-1$, that is, ${\text{col}}_j(\text{PA})[i]=k$ if and only if $S_k[1..j-1]$ is ranked $i$ in the co-lexicographic order of $S_1[1..j-1],\mathrm{\dots },S_{0pt}[1..j-1]$.

Let PBWT be the matrix representing the positional BWT of $M$. Then PBWT is also an $0pt\times 0pt$ matrix in which ${\text{col}}_j(\text{PBWT})[i]={\text{col}}_j(M)[{\text{col}}_j(\text{PA})[i]]$ for all $i\in [\mathrm{1..0}pt]$ and $j\in [\mathrm{1..0}pt]$. We refer to a maximal substring of identical characters in ${\text{col}}_j(\text{PBWT})$ as a run. Throughout the paper, let $r_j$ denote the number of runs for ${\text{col}}_j(\text{PBWT})$. We define $\tilde{r}$ as ${\sum }_{1\le j\le 0pt}{r}_j$.

We define $\text{fore}[i][j]$, for $i\in [\mathrm{1..0}pt]$ and $j\in [\mathrm{1..0}pt)$, as the (row) index of ${\text{col}}_j(\text{PA})[i]$ in ${\text{col}}_{j+1}(\text{PA})$, and $\text{back}[i][j]$, for $i\in [\mathrm{1..0}pt]$ and $j\in (\mathrm{1..0}pt]$, as the (row) index of ${\text{col}}_j(\text{PA})[i]$ in ${\text{col}}_{j-1}(\text{PA})$.

Previously, the operation $\text{fore}[i][j]$ could be implemented as follows. Let $C_c$ denote the number of occurrences of symbols in column $j$ of the PBWT, that is, in ${\text{col}}_j(\text{PBWT})$. Then, $\text{fore}[i][j]=C_c+{\text{rank}}_c({\text{col}}_j(\text{PBWT}),i)$, where $c={\text{col}}_j(\text{PBWT})[i]$. The backward stepping operation can be implemented in a symmetric manner.

Proposition 3 states the key properties of fore.

When constructing the $j$-th column of PBWT, note that the prefixes of the haplotypes up to column $j-1$ are stably sorted. Therefore, statements (a) and (b) follow directly from this stable ordering. $◀$

Consider any pair $(S_i,S_{i+1})$ such that $S_i\ne S_{i+1}$. Let $j^{\prime \prime }$ denote the length of the longest common prefix between $S_i$ and $S_{i+1}$. Observe that the indices $i$ and $i+1$ remain consecutive in ${\text{col}}_j(\text{PA})$ for all $1\le j\le j^{\prime \prime }+1$. Let $\tau$ be the row index such that ${\text{col}}_{j^{\prime \prime }+1}(\text{PA})[\tau ]=i$. By the observation above, we have ${\text{col}}_{j^{\prime \prime }+1}(\text{PA})[\tau +1]=i+1$. Since ${\text{col}}_{j^{\prime \prime }+1}(\text{PBWT})[\tau ]\ne {\text{col}}_{j^{\prime \prime }+1}(\text{PBWT})[\tau +1]$ by the definition of $j^{\prime \prime }$, there must be a run boundary between rows $\tau$ and $\tau +1$ in the $(j^{\prime \prime }+1)$-st column of PBWT. Therefore, each pair $(S_i,S_{i+1})$ with $S_i\ne S_{i+1}$ corresponds to a run boundary in PBWT. The mapping from such a pair to a run boundary at $(\tau ,j^{\prime \prime }+1)$ is injective, since each haplotype index $i$ appears at a unique row $\tau$ in PA of any column. This establishes that $\tilde{r}\ge 0p{t}^{\prime \prime }+1$. $◀$

The statement holds because the number of distinct haplotypes in $\{S_1,S_2,\mathrm{\dots },S_{0pt}\}$ is at most $0p{t}^{\prime \prime }+1$. $◀$

We call an interval $[b,e]$ canonical with respect to column $j$ of the PBWT, for $1\le j\le 0pt$, if $[b,e]$ is the maximal interval such that $S_{{\text{col}}_j(\text{PA})[b]}=S_{{\text{col}}_j(\text{PA})[b+1]}=\mathrm{\cdots }=S_{{\text{col}}_j(\text{PA})[e]}$. Let ${\mathrm{\ell }}_j$ denote the number of canonical intervals with respect to column $j$ for $1\le j\le 0pt$. Lemma 6 describes the relationship between ${\mathrm{\ell }}_j$ and $r_j$.

Consider column $j$ for any $1\le j\le 0pt$. Let $[b,e]$ denote any canonical interval with respect to this column. Observe that ${\text{col}}_j(\text{PBWT})[b]={\text{col}}_j(\text{PBWT})[b+1]=\mathrm{\cdots }={\text{col}}_j(\text{PBWT})[e]$. Hence, $[b,e]$ corresponds to a contiguous block of identical symbols in ${\text{col}}_j(\text{PBWT})$, although the block might not be maximal. This implies that $r_j\le {\mathrm{\ell }}_j$. $◀$

The proof proceeds by induction on $j$. Let $H$ be the array consisting of in increasing order. It follows that $|H|=1+0p{t}^{\prime \prime }$.

For the base case $j=1$, observe that each interval $[H[t-1]+1,H[t]]$ for $1\le t\le |H|$ (under the convention $H[0]=0$) forms a distinct canonical interval with respect to the first column. Since ${\bigcup }_t[H[t-1]+1,H[t]]=[1,0pt]$, we have ${\mathrm{\ell }}_1=|H|$.

Assume inductively that ${\mathrm{\ell }}_{j-1}\le 0p{t}^{\prime \prime }+1$. Let $[b_{j-1},e_{j-1}]$ denote any canonical interval with respect to column $j-1$, so we have ${\text{col}}_{j-1}(\text{PBWT})[b_{j-1}]=\mathrm{\cdots }={\text{col}}_{j-1}(\text{PBWT})[e_{j-1}]$. By Proposition 3(b), the integers $\text{fore}[b_{j-1}][j-1],\text{fore}[b_{j-1}+1][j-1],\mathrm{\dots },\text{fore}[e_{j-1}][j-1]$ are consecutive, forming an interval $[\text{fore}[b_{j-1}][j-1],\text{fore}[e_{j-1}][j-1]]$. Furthermore, by the definitions of fore and $[b_{j-1},e_{j-1}]$, the haplotypes $S_{{\text{col}}_j(\text{PA})[t]}$ for all $t\in [\text{fore}[b_{j-1}][j-1],\text{fore}[e_{j-1}][j-1]]$ are identical. Thus, this interval is contained within some canonical interval with respect to column $j$. This containment may be strict whenever multiple canonical intervals with respect to column $j-1$, corresponding to identical haplotype sequences, are mapped to adjacent positions in column $j$, thereby merging into a single canonical interval with respect to column $j$. Therefore, each canonical interval with respect to column $j-1$ is mapped by fore to a subinterval of a canonical interval with respect to column $j$. Since the union of these mapped intervals covers all rows $[1,0pt]$, the number of canonical intervals in column $j$ must be less than or equal to the number in column $j-1$ (${\mathrm{\ell }}_j\le {\mathrm{\ell }}_{j-1}$). By the induction hypothesis, ${\mathrm{\ell }}_j\le {\mathrm{\ell }}_{j-1}\le 0p{t}^{\prime \prime }+1$, completing the proof. $◀$

By Lemmas 6 and 7, we have $r_j\le 0p{t}^{\prime \prime }+1$ for all $j$. Therefore, $\tilde{r}={\sum }_{1\le j\le 0pt}{r}_j\le {\sum }_{1\le j\le 0pt}(0p{t}^{\prime \prime }+1)=0pt(0p{t}^{\prime \prime }+1).$ $◀$

Let $x=|I_p|$ and $I_p=[p_1,p_1^{\prime }],[p_2,p_2^{\prime }],\mathrm{\dots },[p_x,p_x^{\prime }]$ be $x$ pairwise-disjoint intervals into which $[1,n]$ is partitioned, where $p_1=1,p_x^{\prime }=n$, and $p_i^{\prime }+1=p_{i+1}$ for . Let $y=|I_q|$. Without loss of generality, assume that $y\ge 4$; otherwise, $I_p$ immediately satisfies the three-overlap constraint. Let $I_q=[q_1,q_1^{\prime }],[q_2,q_2^{\prime }],\mathrm{\dots },[q_y,q_y^{\prime }]$ be $y$ pairwise-disjoint intervals into which $[1,n]$ is partitioned, where $q_1=1,q_y^{\prime }=n$, and $q_j^{\prime }+1=q_{j+1}$ for .

The algorithm $\text{normalization}(I_p,I_q)$ is described as follows. Create a variable $k$ and initiate $k$ to $1$. Iterate and process each interval $[p_i,p_i^{\prime }]$ in $I_p$ from left to right as follows: If $[p_i,p_i^{\prime }]$ overlaps with at most three intervals of $I_q$, then skip $[p_i,p_i^{\prime }]$ and move on to the next interval; otherwise, divide $[p_i,p_i^{\prime }]$ into two sub-intervals, $[p_i,d_k]$ and $[d_k+1,p_i^{\prime }]$, where $1\le d_k\le n$ is the largest integer such that $[p_i,d_k]$ overlaps three intervals of $I_q$. Then, substitute $[p_i,p_i^{\prime }]$ for $[p_i,d_k]$ and $[d_k+1,p_i^{\prime }]$, increment $k$ by one, and move on to the next interval, that is, $[d_k+1,p_i^{\prime }]$.

Let ${\widehat{I}}_p$ denote the list of intervals outputted by the above algorithm. Figure 2 illustrates an example of the algorithm. Clearly, we have $|{\widehat{I}}_p|=x+k$ and it follows that every interval in ${\widehat{I}}_p$ overlaps at most three intervals in $I_q$. It remains to prove that $k\le \lfloor \frac{y}{2}\rfloor$. To this end, we define $\mathrm{\Lambda }({\widehat{I}}_p[i])$ for any $1\le i\le x+k$ to be the set of intervals in $I_q$, overlapping ${\widehat{I}}_p[i]$.

Observe that we have $d_{\tau }\in \{q_1^{\prime },\mathrm{\dots },q_y^{\prime }\}$ for any , and thus the integers $d_{\tau }$ and $d_{\tau }+1$ are always in different intervals in $I_q$. Moreover, $d_{t_1}+1\notin \widetilde{𝖨}$ for any $\widetilde{𝖨}\in \mathrm{\Lambda }({𝖨}_1)$, in view of the algorithm.

Assume that there exists an interval $𝖨\in \mathrm{\Lambda }({𝖨}_1)\cap \mathrm{\Lambda }({𝖨}_2)$. Let $z$ denote the left endpoint of the interval ${𝖨}_2$. Then, we have $d_{t_1}\in 𝖨$ and $z\in 𝖨$. Since , we have $d_{t_1}+1\in 𝖨\in \mathrm{\Lambda }({𝖨}_1)$ as well, a contradiction. Hence, the assumption is false. $\mathrm{⊲}$

In view of the algorithm, the interval in ${\widehat{I}}_p$ ending at $d_{\tau }$ for any overlaps exactly three intervals in $I_q$. In view of Claim 11, we have $k-1\le \lfloor \frac{y}{3}\rfloor$; therefore, $x+k\le x+1+\lfloor \frac{y}{3}\rfloor \le x+\lfloor \frac{y}{2}\rfloor$, for $y\ge 4$.

Clearly, the algorithm runs in $O(|{\widehat{I}}_p|+|I_p|+|I_q|)\subseteq O(|I_p|+|I_q|)$ time. $◀$

The construction of ${\text{SubIB}}_j$ proceeds by induction. In the base case, we set ${\text{SubIB}}_1={\text{intervals}}_1$. For each $j$ from $2$ to $0pt$, we construct ${\text{SubIB}}_j$ as follows: we apply the algorithm normalization$(I_p,I_q)$ (see Lemma 10) with $I_p:={\text{intervals}}_j$ and $I_q:={\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$. Recall that the function ${\text{foreL}}_{j-1}$ maps sub-intervals in column $j-1$ to sub-intervals in column $j$. The list of sub-intervals output by this algorithm is assigned to ${\text{SubIB}}_j$. See Figure 3 for an illustration of the algorithm and an example.

In view of the algorithm in Lemma 10, each interval in ${\text{SubIB}}_j$ for $1\le j\le 0pt$ is a sub-interval of some interval in ${\text{intervals}}_j$ and overlaps at most three intervals in the list ${\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$. Recall that every interval in ${\text{intervals}}_j$ corresponds to a run in ${\text{col}}_j(\text{PBWT})$. Hence, every interval in ${\text{SubIB}}_j$ corresponds to a sub-run in ${\text{col}}_j(\text{PBWT})$.

The construction of ${\text{SubIF}}_j$ is also performed by induction. In the base case, we set ${\text{SubIF}}_{0pt}={\text{intervals}}_{0pt}$. Then, for $j$ from $0pt-1$ down to $1$, we construct ${\text{SubIF}}_j$ as follows: we apply the algorithm normalization$(I_p,I_q)$ (see Lemma 10), where $I_p={\text{foreL}}_j({\text{intervals}}_j)$ and $I_q={\text{SubIF}}_{j+1}$. Let ${\text{SubIF}}_{j+1}^{\prime }$ be the list of intervals output by the algorithm; we then assign ${\text{SubIF}}_j={\text{backL}}_{j+1}({\text{SubIF}}_{j+1}^{\prime })$. Figure 6 (in Appendix A) illustrates the construction scheme along with an example.

By the normalization algorithm shown in Lemma 10, each interval in ${\text{SubIF}}_{j+1}^{\prime }$ for $1\le j\le 0pt$ is a sub-interval of some interval in ${\text{foreL}}_j({\text{intervals}}_j)$ and overlaps at most three intervals in ${\text{SubIF}}_{j+1}$. Recall that every interval in ${\text{foreL}}_j({\text{intervals}}_j)$ corresponds to a run in ${\text{col}}_j(\text{PBWT})$. Likewise, every interval in ${\text{SubIF}}_{j+1}^{\prime }$ corresponds to a sub-run in ${\text{col}}_j(\text{PBWT})$. Since ${\text{SubIF}}_j={\text{backL}}_{j+1}({\text{SubIF}}_{j+1}^{\prime })$, every interval in ${\text{SubIF}}_j$ is a sub-run in ${\text{col}}_j(\text{PBWT})$.

Observe that $|{\text{foreL}}_j({\text{SubIB}}_j)|=|{\text{SubIB}}_j|$. Recall that $|{\text{intervals}}_j|=r_j$ is always true. In the base case, we have $|{\text{SubIB}}_1|=r_1$. For , we have the recursion $|{\text{SubIB}}_j|\le |{\text{intervals}}_j|+\lfloor \frac{|{\text{SubIB}}_{j-1}|}{2}\rfloor \le r_j+\frac{|{\text{SubIB}}_{j-1}|}{2}$, in view of Lemma 10 and the observation above. By solving the recursion, it follows that $|{\text{SubIB}}_j|=\frac{r_1}{2^{j-1}}+\frac{r_2}{2^{j-2}}+\mathrm{\cdots }+\frac{r_j}{2^{j-j}}={\sum }_{1\le \tau \le j}\frac{r_{\tau }}{2^{j-\tau }}$. Therefore, we have . The bound on the total number of $|{\text{SubIF}}_j|$ for all $j$ can be computed similarly. This concludes the proof. $◀$

We first describe the data structure for back queries. For each column $j$ of the PBWT, we construct the list ${\text{SubIB}}_j$ as described in Section 5 (see Constructing sub-runs for back queries). Observe that each interval in ${\text{SubIB}}_j$ overlaps at most three intervals in ${\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$. Moreover, each interval $\widetilde{𝖨}\in {\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$ corresponds to the interval $[\text{back}[\widetilde{𝖨}.b][j],\text{back}[\widetilde{𝖨}.e][j]]$ belonging to ${\text{SubIB}}_{j-1}$. For every and each interval ${\text{SubIB}}_j[\tau ]$ with $1\le \tau \le |{\text{SubIB}}_j|$, we store a list $B_j^{\tau }$ of quadruples $(\tilde{s},\tilde{t},s,\lambda )$, where:

• $\mathrm{■}$

  $[\tilde{s},\tilde{t}]$ is a distinct interval in ${\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$ that overlaps ${\text{SubIB}}_j[\tau ]$;

• $\mathrm{■}$

  $s=\text{back}[\tilde{s}][j]$; and

• $\mathrm{■}$

  $\lambda$ is the index of the interval in ${\text{SubIB}}_{j-1}$ whose left endpoint is $s$.

Since each ${\text{SubIB}}_j[\tau ]$ overlaps at most three intervals in ${\text{foreL}}_{j-1}({\text{SubIB}}_{j-1})$, the list $B_j^{\tau }$ contains at most three tuples. By Lemma 12, , so the total space usage is $O(\tilde{r})$ words.

Let $x$, given in a query, denote the index of the sub-run in ${\text{SubIB}}_{j^{\prime }}$ containing $i^{\prime }$. To answer a query $\text{back}[i^{\prime }][j^{\prime }]$, we search in the list $B_{j^{\prime }}^x$ for the quadruple $(\tilde{s},\tilde{t},s,\lambda )$ satisfying $i^{\prime }\in [\tilde{s},\tilde{t}]$. Then $\text{back}[i^{\prime }][j^{\prime }]=i^{\prime }-\tilde{s}+s$, and $\lambda$ is the index of the interval in ${\text{SubIB}}_{j^{\prime }-1}$ containing $\text{back}[i^{\prime }][j^{\prime }]$. We return $(i^{\prime }-\tilde{s}+s,\lambda )$ as the result. Since $|B_{j^{\prime }}^x|\le 3$, each query takes $O(1)$ time. See Figure 4 for an illustration of the data structure and algorithm.

We store, for each $1\le j\le 0pt$, an array ${\text{ValBack}}_j$ of length $|{\text{SubIB}}_j|$, where each entry ${\text{ValBack}}_j[\tau ]$ is set to ${\text{col}}_j(\text{PBWT})[{\text{SubIB}}_j[\tau ].b]$. This requires $O(r_j)$ words per column, and $O(\tilde{r})$ words overall. Given the index $x$ of the sub-run containing $i^{\prime }$ in ${\text{SubIB}}_{j^{\prime }}$, we have ${\text{col}}_{j^{\prime }}(\text{PBWT})[i^{\prime }]={\text{ValBack}}_{j^{\prime }}[x]$, which can be retrieved in $O(1)$ time.

We construct the lists ${\text{SubIF}}_j$ for each column $j$ as described in Section 5 (see Constructing sub-runs for fore queries). Each interval $𝖨\in {\text{SubIF}}_j$ corresponds to $[\text{fore}[𝖨.b][j],\text{fore}[𝖨.e][j]]\in {\text{foreL}}_j({\text{SubIF}}_j)$, and each interval in ${\text{foreL}}_j({\text{SubIF}}_j)$ overlaps at most three intervals in ${\text{SubIF}}_{j+1}$. For every and each interval ${\text{SubIF}}_j[\tau ]$, we store a list $F_j^{\tau }$ of quintuplets $(s^{\prime },\tilde{s},s,t,\lambda )$, where:

• $\mathrm{■}$

  $s^{\prime }={\text{SubIF}}_j[\tau ].b$, $\tilde{s}=\text{fore}[s^{\prime }][j]$,

• $\mathrm{■}$

  $[s,t]$ is an interval in ${\text{SubIF}}_{j+1}$ that overlaps $[\text{fore}[{\text{SubIF}}_j[\tau ].b][j],\text{fore}[{\text{SubIF}}_j[\tau ].e][j]]$,

• $\mathrm{■}$

  $\lambda$ is the index of $[s,t]$ in ${\text{SubIF}}_{j+1}$.

Since each such interval, i.e., $[\text{fore}[{\text{SubIF}}_j[\tau ].b][j],\text{fore}[{\text{SubIF}}_j[\tau ].e][j]]$, overlaps at most three intervals in ${\text{SubIF}}_{j+1}$, each $F_j^{\tau }$ has at most three tuples. By Lemma 12, the total space is again $O(\tilde{r})$ words.

Let $x$ denote the index of the sub-run containing $i^{\prime }$ in ${\text{SubIF}}_{j^{\prime }}$. To answer $\text{fore}[i^{\prime }][j^{\prime }]$, we search $F_{j^{\prime }}^x$ for the quintuple $(s^{\prime },\tilde{s},s,t,\lambda )$ satisfying $i^{\prime }-s^{\prime }+\tilde{s}\in [s,t]$. Then $\text{fore}[i^{\prime }][j^{\prime }]=i^{\prime }-s^{\prime }+\tilde{s}$, and $\lambda$ is the index of the interval in ${\text{SubIF}}_{j^{\prime }+1}$ containing $\text{fore}[i^{\prime }][j^{\prime }]$. We return $(i^{\prime }-s^{\prime }+\tilde{s},\lambda )$ as the result. Since $|F_{j^{\prime }}^x|\le 3$, each query takes $O(1)$ time. An illustration of the data structure and the algorithm is depicted in Figure 7 (in Appendix B).

Analogously, for each $1\le j\le 0pt$, we store an array ${\text{ValFore}}_j$ of length $|{\text{SubIF}}_j|$, where ${\text{ValFore}}_j[\tau ]={\text{col}}_j(\text{PBWT})[{\text{SubIF}}_j[\tau ].b]$. This requires $O(r_j)$ words per column, and $O(\tilde{r})$ words in total. Given the index $x$ of the sub-run containing $i$ in ${\text{SubIF}}_j$, we have ${\text{col}}_j(\text{PBWT})[i]={\text{ValFore}}_j[x]$, retrievable in $O(1)$ time. $◀$