The TAG Array of a Multiple Sequence Alignment

Genomic analyses and diagnostics is often reference-based, that is, samples are compared with a reference genome of e.g. humans to speed up the analysis. For a long time, such a reference consisted of a single genome. However, a single reference sequence cannot capture the genetic diversity of a population. As a consequence, modern tools focus on representations that include (common) variations. Such a representation of the genomes of a population is commonly called a pangenome.¹¹1The concept of a pangenome was originally developed for bacterial studies [38], but now can refer to the entire genomic variation of any population. In this paper, we focus on the pangenome of a species. A common use case for pangenome representations is read mapping, where one seeks to determine the locus in the genome of a short substring of a DNA sample called a read, while accounting for sequencing errors and natural variation. Using pangenome representations can significantly reduce reference bias and thus mapping errors [8, 23, 26, 27, 35].

Pangenome representations are mostly based on sequence graphs (see e.g. [1, 2]). However, no widely agreed-upon metric for the quality of such graphs exists because several desirable qualities are at odds with each other [2]. For example, using the same subgraph to represent a variant that occurs in multiple sequences may lead to a smaller graph, but may also induce chimeric sequences (i.e., sequences that are valid in the graph but do not exist in the pangenome), e.g. by allowing the combination of variants that do not occur together in nature. Additionally, under the strong exponential-time hypothesis (SETH) it is impossible to index a sequence graph in polynomial time such that string matching queries can be answered in sub-quadratic time [16]. For this reason, methods based on sequence graphs either resort to heuristic matching or cannot guarantee a good worst-case time complexity. In contrast, there are indexes for strings or sets of strings which achieve optimal construction and query time complexities [19]. While these data structures would require space proportional to the combined size of all haplotypes present in a sequence graph and thus be far less practical, almost as good time complexities can be achieved with indexes based on the run-length compressed Burrows-Wheeler Transform (BWT) [9] such as variants of the $r$-index [22] (see e.g. [4, 5, 12, 30]).

In combination with algorithms to construct these indexes for such huge datasets [6, 14, 27, 28, 31, 33], this opens the possibility to represent a population just by the set of sequences and not worry about e.g. graph indexing.

One problem of such a simple representation is that now a match may occur at thousands of positions in the index even when those occurrences all correspond to the same locus. However, in practice, it is often desirable to locate a match relative to a linear reference sequence. Recent works that use the BWT for pangenome indexing are for instance ropebwt3 [27] and Moni-align [39]. The former aligns samples using the BWT-based index, while the latter merely uses the index to search for seeds (maximal exact matches between the pangenome and the pattern) and then aligns only to the most promising sequences in the pangenome.

Recently, we showed how to compute a multiple sequence alignment (MSA) of even very long sequences (such as human chromosomes), given that those sequences are very similar [32]. In an MSA, equal substrings that correspond to the same locus are aligned. Therefore, we can identify the aforementioned thousands of positions in the index by just the corresponding column in the MSA. Note that it is trivial to directly map such a column in the MSA to a single reference sequence if that reference is part of the alignment. Tools such as ropebwt3 “cannot project the alignment [of a pattern] to a designated reference genome” [27], so “it seems [..] interesting to know which column of the alignment a character in the BWT comes from” [21].

For $i,j\in {\mathbb{N}}_0$ we denote the set $\{k\in {\mathbb{N}}_0\mid i\le k\le j\}$ by the interval notations $[i..j]=[i..j+1)=(i-1..j]=(i-1..j+1)$. A string (or array) $S$ of length $n$ over an alphabet $\mathrm{\Sigma }$ is a sequence of $n$ characters from $\mathrm{\Sigma }$. We denote the length $n$ of $S$ by $|S|$ and the $i$th symbol of $S$ by $S[i-1]$, i.e., strings and arrays are zero-indexed. The substring (or subarray) of $S$ from $i$ to $j$ is denoted by $S[i..j]=S[i..j+1)=S(i-1..j]=S(i-1..j+1)=S[i]S[i+1]\mathrm{\dots }S[j]$. For $i>j$, $S[i..j]$ is the empty string $ϵ$ of length $0$. A substring of the form $S[i..n)$ is a suffix of $S$ and is denoted by ${\mathrm{suf}}_i(S)$. A bit vector is an array over the binary alphabet $\{0,1\}$.

A multiple sequence alignment (MSA) of a set of strings $𝒮$ is obtained by inserting a number of gap characters ‘-’ into each string in $𝒮$ such that the resulting strings all have the same length. An example MSA of the strings $\{\text{ACGACT\$},\text{AAACT\$},\text{ACGCAGT\$}\}$ is given in the top-left of Figure 1. A “good” MSA has few gap-symbols and columns mostly contain the same character. In this paper, we are not concerned with the precise optimization objective or methods for construction, so we refer the interested reader to [11] for an overview. However, it is noteworthy that any “good” MSA of sufficiently similar sequences will exhibit contextual locality, that is, the suffixes starting in a column in the MSA are likely to be similar (i.e., share a long prefix) [3, 21]. For instance, the suffixes starting in the fifth column of the MSA in Figure 1 are ACT$, ACT$ and AGT$.

We assume totally ordered alphabets. This induces a total order on strings. Specifically, we say a string $S$ of length $n$ is lexicographically smaller than a string $T$ of length $m$ if and only if there is some $\mathrm{\ell }\le \min \{n,m\}$ such that $S[0..\mathrm{\ell })=T[0..\mathrm{\ell })$ and either ($S$ is a prefix of $T$) or and .

Given a string $S$, ${\mathrm{𝑟𝑎𝑛𝑘}}_c(S,i)$ is the number of $c$’s occurring in $S$ up to (but excluding) index $i$, i.e., ${\mathrm{𝑟𝑎𝑛𝑘}}_c(S,i)=\left|\{j\in [0..i)\mid S[j]=c\right\}|$. The select function returns the index of the $i$th occurrence of $c$ (zero-based) in $S$, i.e., ${\mathrm{𝑟𝑎𝑛𝑘}}_c(S,{\mathrm{𝑠𝑒𝑙𝑒𝑐𝑡}}_c(S,i))=i$.

The generalized suffix array for a collection of strings $S_1,\mathrm{\dots },S_n$ is an array $\mathrm{𝖦𝖲𝖠}$ where $\mathrm{𝖦𝖲𝖠}[i]=(k,j)$ indicates that there are $i$ lexicographically smaller suffixes than ${\mathrm{suf}}_j(S_k)$ among the suffixes of $S_1,\mathrm{\dots },S_n$. We assume the strings to be dollar-terminated, that is, the last character of each string $S_i$ is “$”, which is smaller than all other characters. In the case of a tie between equal suffixes, we define the one occurring earlier in the input to be smaller.²²2This is equivalent to using $n$ terminal symbols and terminating $S_i$ with ${\text{\$}}_i$. Hence the name “multidollar-$\mathrm{𝖤𝖡𝖶𝖳}$.”

Throughout this paper, we use the multidollar-$\mathrm{𝖤𝖡𝖶𝖳}$ as the BWT for string collections (see [10] for an overview of such BWT variants). In the remainder of this paper, we refer to the multidollar-$\mathrm{𝖤𝖡𝖶𝖳}$ just by $\mathrm{𝖡𝖶𝖳}$. It is defined as follows. Let $\mathrm{𝖦𝖲𝖠}[i]=(k,j)$. Then $\mathrm{𝖡𝖶𝖳}[i]=S_k[j-1]$ if $j>0$ and $\mathrm{𝖡𝖶𝖳}[i]=S_k[|S_k|-1]=\text{‘}{\text{\$}}_k\text{’}$ otherwise. An example can be seen in Figure 1. The $\mathrm{𝖫𝖥}$-mapping is a function such that $\mathrm{𝖫𝖥}[i]={\mathrm{𝖦𝖲𝖠}}^{-1}[(k,j-1)]$ for $\mathrm{𝖦𝖲𝖠}[i]=(k,j)$ if $j>0$ and $\mathrm{𝖫𝖥}[i]={\mathrm{𝖦𝖲𝖠}}^{-1}[(k,|S_k|-1)]$ otherwise. $\mathrm{𝖫𝖥}$ can thus be used to iterate over an input string in reverse order, given the index $p$ of the last character of the input string in $\mathrm{𝖡𝖶𝖳}$. Specifically, $S_k=\mathrm{𝖡𝖶𝖳}[{\mathrm{𝖫𝖥}}^{|S_k|-1}[p]]\mathrm{\dots }\mathrm{𝖡𝖶𝖳}[{\mathrm{𝖫𝖥}}^1[p]]\mathrm{𝖡𝖶𝖳}[{\mathrm{𝖫𝖥}}^0[p]]$ where $\mathrm{𝖦𝖲𝖠}[p]=(k,0)$.

Although we use the dollar-$\mathrm{𝖤𝖡𝖶𝖳}$, the algorithms presented in this paper are applicable to all BWT variants where for each input string $S$ there is an index $i$ such that the indices ${\mathrm{𝖫𝖥}}^0[i],{\mathrm{𝖫𝖥}}^1[i],\mathrm{\dots },{\mathrm{𝖫𝖥}}^{|S|-1}[i]$ are all distinct and $S=\mathrm{𝖡𝖶𝖳}[{\mathrm{𝖫𝖥}}^{|S|-1}[i]]\mathrm{\dots }\mathrm{𝖡𝖶𝖳}[{\mathrm{𝖫𝖥}}^0[i]]$.³³3Notably, this excludes the original EBWT for non-primitive input strings, because there each root of such a string corresponds to a distinct cycle in $\mathrm{𝖫𝖥}$. As far as we can tell, all other BWT variants for string collections satisfy this requirement.

Note that this definition gives an immediate linear-time algorithm for computing $\mathrm{𝖳𝖠𝖦}$: For each character in the alignment, use the inverse of $\mathrm{𝖦𝖲𝖠}$ to find the corresponding position in $\mathrm{𝖳𝖠𝖦}$ and write the character’s column. However, this naïve approach requires holding the inverse of $\mathrm{𝖦𝖲𝖠}$ in main memory. Definition 1 is equivalent to the following definition based on the $\mathrm{𝖡𝖶𝖳}$ instead of the suffix array, which follows immediately from the definition of the $\mathrm{𝖡𝖶𝖳}$ given above.

We will use this definition from now on as it does not depend on $\mathrm{𝖦𝖲𝖠}$. In simpler terms, consider a character in column $i$ in row $s$ which corresponds to $\mathrm{𝖡𝖶𝖳}[j]$. Then $\mathrm{𝖳𝖠𝖦}[j]$ is $\mathrm{𝑐𝑜𝑙}$, where $\mathrm{𝑐𝑜𝑙}$ is the first column after $i$ where there is a non-gap character in row $s$. An example of the $\mathrm{𝖳𝖠𝖦}$ array of a multiple string alignment (MSA) is shown in Figure 1. For instance, the base G in the first string corresponds to $\mathrm{𝖡𝖶𝖳}[9]=\text{G}$, and the character following this G in the alignment occurs in column $4$. Therefore, we have $\mathrm{𝖳𝖠𝖦}[9]=4$.

In the next section, we present an algorithm that requires only access to the $\mathrm{𝖡𝖶𝖳}$, $\mathrm{𝖫𝖥}$ and the positions of the gaps in the alignment to compute the run-length encoded TAG array.

Before describing our algorithm we briefly recall how we can reconstruct an input string from the $\mathrm{𝖡𝖶𝖳}$ using $\mathrm{𝖫𝖥}$. For this, we start at the position $p$ of the string’s last character (a “$” in our case) in the $\mathrm{𝖡𝖶𝖳}$.⁴⁴4Equivalently, the rank of the input string in the sorted list of all rotations. Many tools for computing the $\mathrm{𝖡𝖶𝖳}$ directly output these indices [6, 31, 33]. $\mathrm{𝖫𝖥}[p]$ now gives us the position in the $\mathrm{𝖡𝖶𝖳}$ of the character preceding the “$” in the string, $\mathrm{𝖫𝖥}[\mathrm{𝖫𝖥}[p]]$ gives us the character preceding that, and so on.

Essentially, we compute $\mathrm{𝖳𝖠𝖦}$ using the inverse of Fact 2: for each TAG value $\mathrm{𝑡𝑎𝑔}$, we find the indices where the TAG array $\mathrm{𝖳𝖠𝖦}$ has value $\mathrm{𝑡𝑎𝑔}$. Consider a column $\mathrm{𝑐𝑜𝑙}$. By Fact 2, $\mathrm{𝖳𝖠𝖦}$ should have value $\mathrm{𝑐𝑜𝑙}$ at the indices where the characters preceding the characters in column $\mathrm{𝑐𝑜𝑙}$ are in $\mathrm{𝖡𝖶𝖳}$.

We consider each column of the alignment from right to left and maintain the BWT positions by iterating over the input strings from right to left using $\mathrm{𝖫𝖥}$ as is done when reconstructing the strings from $\mathrm{𝖡𝖶𝖳}$. By incorporating the positions of the gap symbols, this is synchronized such that we consider a column of the alignment at each time step.

Note that, given the indices of the characters in column $\mathrm{𝑐𝑜𝑙}$ in the $\mathrm{𝖡𝖶𝖳}$, the indices of the preceding characters can be found with a single $\mathrm{𝖫𝖥}$ step each. If a symbol in the alignment in column $\mathrm{𝑐𝑜𝑙}$ is a gap (‘-’), the first preceding non-gap character has a TAG value greater than $\mathrm{𝑐𝑜𝑙}$. Therefore, in this case we just ignore the current column in this string.

Algorithm 1 shows the procedure. Note that any two iterations of the second for-loop of Algorithm 1 concern two different $\mathrm{𝑐𝑜𝑙}$ values (columns). Therefore, no two iterations of this loop can concern the same TAG run. It is thus possible to immediately compute the TAG runs without having to store the entire $\mathrm{𝖳𝖠𝖦}$ array.

Modifying Algorithm 1 accordingly, we obtain an algorithm that outputs the TAG runs in order of decreasing TAG value. To obtain the run-length compressed TAG array, it is hence necessary to sort the TAG runs by (start or end) index afterwards.

Because both the MSA and BWT possess contextual locality, the TAG array of an MSA of similar strings is likely run-length compressible [21]. One can facilitate random access to the run-length compressed TAG array using a sparse bit vector indicating the boundaries of the TAG runs. Assuming $r_{\mathrm{𝖳𝖠𝖦}}$ TAG runs and an alignment with length $m$ and $N$ non-gap characters, this would require at least $r_{\mathrm{𝖳𝖠𝖦}}\cdot {\log }_{2}m+r_{\mathrm{𝖳𝖠𝖦}}\cdot {\log }_2\frac{N}{r_{\mathrm{𝖳𝖠𝖦}}}$ bits for the TAG run labels and the run boundaries.

However, $r_{\mathrm{𝖳𝖠𝖦}}$ is likely to be larger than the number $r$ of runs of the $\mathrm{𝖡𝖶𝖳}$ (and guaranteed to be at least the length $m$ of the alignment). Additionally, the alphabet used for $\mathrm{𝖳𝖠𝖦}$ is much larger than that for $\mathrm{𝖡𝖶𝖳}$. Therefore, a naïve run-length compressed TAG array would require many times as much memory as the run-length compressed $\mathrm{𝖡𝖶𝖳}$, with the majority of the memory used for the run labels. In this section we present a technique for sampling the TAG array where we store the labels of only a fraction $\frac{1}{s}$ of the TAG runs, for some user-defined parameter $s$. For this, $\mathrm{𝖫𝖥}$ is assumed to be available.⁵⁵5This is not a major restriction because efficient access to $\mathrm{𝖫𝖥}$ can be facilitated using space proportional to the number of runs $r$ of the $\mathrm{𝖡𝖶𝖳}$ [22], and $\mathrm{𝖫𝖥}$ is available anyway in many tools using the $\mathrm{𝖡𝖶𝖳}$.

Note that there is no obvious method to just sample every $s$-th TAG value and use $\mathrm{𝖫𝖥}$ to walk to the next sampled TAG run for a query. This is because, regardless of which cells are chosen as the starting points for $\mathrm{𝖫𝖥}$ steps in the TAG runs, the $\mathrm{𝖫𝖥}$ steps may “jump over” sampled TAG runs.⁶⁶6In any sensible MSA, there is no column which contains only gap symbols. Note however, that there may very well be TAG runs where each corresponding symbol in the alignment is preceded by a gap symbol. As a consequence, it would not be possible to guarantee that a sampled TAG run can be reached with a given number of $\mathrm{𝖫𝖥}$ steps. In the following, we present a more involved sampling strategy which guarantees that we can always reach a sampled TAG run with fewer than $s$ $\mathrm{𝖫𝖥}$ steps.

The following concepts and data structures are illustrated in Figure 2 for our running example.

Let $e[i]$ be the index of the TAG run containing $\mathrm{𝖫𝖥}[b_i]$, where $b_i$ is the start of TAG run $i$. Let $ℝ𝔹$ be a (sparse) bit vector indicating the TAG run boundaries, i.e., $ℝ𝔹[i]=1$ if and only if $i$ is the start of a TAG run. We have $b_i={\mathrm{𝑠𝑒𝑙𝑒𝑐𝑡}}_1(ℝ𝔹,i)$ and $e[i]={\mathrm{𝑟𝑎𝑛𝑘}}_1(ℝ𝔹,\mathrm{𝖫𝖥}[b_i]+1)-1$. Therefore, $e[i]$ can be computed with one $\mathrm{𝖫𝖥}$ computation and a ${\mathrm{𝑠𝑒𝑙𝑒𝑐𝑡}}_1$ and a ${\mathrm{𝑟𝑎𝑛𝑘}}_1$ query and we do not have to store $e$ explicitly.

Now let ${\mathrm{𝖳𝖠𝖦}}^{\prime }$ be such that ${\mathrm{𝖳𝖠𝖦}}^{\prime }[i]$ is the TAG value of the $i$th TAG run. Note that there are at most $n$ indices $i$ where ${\mathrm{𝖳𝖠𝖦}}^{\prime }[e[i]]\ge {\mathrm{𝖳𝖠𝖦}}^{\prime }[i]$, because the run head of such a TAG run must correspond to the first character in a string of the alignment. Let the set of these indices be $R$. Now consider the digraph $G=(V,E)$ with $V=\{0,\mathrm{\dots },r_{\mathrm{𝖳𝖠𝖦}}-1\}$ and $E=\{(i,e[i])\mid i\in [0..r_{\mathrm{𝖳𝖠𝖦}}-1]\setminus R\})$. That is, $G$ is a graph where each TAG run is a node and there is an edge from a node $u$ to $v$ if and only if an $\mathrm{𝖫𝖥}$ step from the run-head of TAG run $u$ results in an index in TAG run $v$. By definition, for each edge $(i,e[i])\in E$ we have . Therefore, $G$ is acyclic. Additionally, the out-degree of each node is at most $1$. Thus, $G$ is not only a DAG but a rooted forest with $R$ as the roots and each edge pointing “upwards” towards a root (cf. Figure 2).

We require a set of sampled TAG runs to satisfy one constraint: Given a sampling rate $s$, the distance (in $G$) from an unsampled node to a sampled node should be less than $s$. Note that this implies that the roots $R$ of the forest are sampled.

Sampled TAG runs are marked in a bit vector $𝕊$. Rank support on this bit vector is used to access an array $𝖫$ storing the TAG values of the sampled TAG runs. We can thus decide in constant time whether a TAG run is sampled and retrieve the TAG value if it is.

To retrieve the TAG value of an unsampled TAG run $i$, we compute $e[i]$ as described above, recursively find the TAG value of TAG run $e[i]$, and add the difference ${\mathrm{𝖳𝖠𝖦}}^{\prime }[i]-{\mathrm{𝖳𝖠𝖦}}^{\prime }[e[i]]$ to the result. The distance-constraint ensures that we can get any TAG value using fewer than $s$ recursion steps. For this to work, we encode the difference ${\mathrm{𝖳𝖠𝖦}}^{\prime }[i]-{\mathrm{𝖳𝖠𝖦}}^{\prime }[e[i]]$ for all unsampled $i$ with unary encoding in a bit vector $𝔻$. For the sampled TAG runs, we also encode the value $1$ for simplicity. Specifically, if the $i$th one-bit in $𝔻$ is at position $p$, the $(i+1)$th one-bit is at position $p+({\mathrm{𝖳𝖠𝖦}}^{\prime }[i]-{\mathrm{𝖳𝖠𝖦}}^{\prime }[e[i]])$ if TAG run $i$ is sampled and at position $p+1$ otherwise. Two select queries on this bit vector then suffice to extract this difference and thus to compute ${\mathrm{𝖳𝖠𝖦}}^{\prime }[i]$ given ${\mathrm{𝖳𝖠𝖦}}^{\prime }[e[i]]$. An additional $1$ is appended for an easier implementation. Algorithm 2 shows how the TAG access is performed.

We implemented the algorithms and data structures described in this paper in C++. In particular, our implementation computes the sampled TAG array as described in Sections 3 and 4 and then computes the sampled TAG index as described in Section 4 together with the data structure for listing the distinct TAGs described in Section 4.2. For the bit vector and $\mathrm{𝑟𝑎𝑛𝑘}$ and $\mathrm{𝑠𝑒𝑙𝑒𝑐𝑡}$ implementations, we used the Succinct Data Structure Library 2.0 [24]. The source code of our implementation is publicly available.⁷⁷7https://gitlab.com/qwerzuiop/msatag

All experiments were conducted on a Linux-6.8.0 machine with an Intel Xeon Gold 6338 CPU and 512 GB of RAM. As the compiler we used GCC 13.3.0. Currently, our construction algorithm is single-threaded. We only test our query algorithm with one thread. However, note that concurrent access to our index is trivially possible.

As test data, we use $1000$ human Chromosome 19 haplotypes from [7] and the corresponding alignment from [32]. The MSA uses $1.67\cdot {10}^8$ gaps and has $5.91\cdot {10}^{10}$ non-gap characters.

To test mapping performance, we extracted ${10}^7$ “reads” with 100bp each, chosen uniformly at random from the sequences used for the alignment. We mutated each base with a probability of $1\%$ and reverse-complemented each read with a probability of $50\%$. We then used ropebwt3 to find the corresponding MEMs and queried our index with the resulting BWT-interval. Finally, we projected each TAGs to the first sequence in the data set as a designated linear reference. ropebwt3 finds $\mathrm{17\hspace{0.17em}953\hspace{0.17em}756}$ MEMs, which cumulatively correspond to $\mathrm{29\hspace{0.17em}185\hspace{0.17em}037}$ TAGs (columns in the MSA).

The dollar-$\mathrm{𝖤𝖡𝖶𝖳}$ of the sequences was computed using lg [31]. Ropebwt3 requires that both forward and reverse-complemented strands are in the index, and we of course need to compute the TAG runs with the same index to ensure that the TAG runs match the BWT. However, we only compute the TAGs for the forward strands. This means that our TAG runs do not cover all positions in the BWT. Space between two TAG runs is treated like any other (unsampled) TAG run, except marked as “no TAG available”.

The dollar-$\mathrm{𝖤𝖡𝖶𝖳}$ has $\mathrm{91\hspace{0.17em}081\hspace{0.17em}437}$ runs, and there are $\mathrm{137\hspace{0.17em}981\hspace{0.17em}814}$ runs in our index, $\mathrm{97\hspace{0.17em}979\hspace{0.17em}057}$ of which have a TAG value.

Table 1 shows construction and query performance for varying sampling rates. The increased construction time with increasing sample rate stems entirely from the construction of the index for reporting the distinct TAG values, constructing the TAG index always took about $175\mathrm{s}$. This is because, for constructing this index, we need to access all TAG values, which we do using our TAG index. Since the average query time of the TAG index is proportional to the sampling rate, a larger sampling rate necessarily slows this down. By using e.g. the non-sampled TAG array or using a more clever (semi-external) algorithm, this could be remedied in the future. The memory used for constructing the index is roughly the same for all tested sampling rates.

Note that with a sampling rate of $4$, our program takes less time for mapping the MEMs to columns in the MSA with a single thread ($76.0\mathrm{s}$) than ropebwt3 uses for finding the MEMs with $8$ threads ($81.7\mathrm{s}$).

We described an algorithm that can compute the TAG array of a multiple sequence alignment (MSA) using the $\mathrm{𝖫𝖥}$ function of the BWT whose working memory is proportional to the number of sequences. Additionally, we described a TAG index that is able to report the unique tags corresponding to a BWT interval (i.e., the columns where matches corresponding to the BWT interval occur in the MSA) in time proportional to the size of the output using standard document listing techniques. We also gave a non-trivial sampling strategy for the TAG index and showed that our TAG construction algorithm can be adapted to output just the sampled TAG values. Finally, we demonstrated experimentally that our construction algorithm and index perform well on real-world data.

Our work enables e.g. the efficient mapping of matches in BWT-based indices to a designated reference genome, which programs such as ropebwt3 currently lack [27]. Our techniques could also be used to obtain more effective chaining heuristics in programs using a BWT-based index to find seeds for a seed-and-extend approach (e.g. Moni-align [39]).