U-Index: A Universal Indexing Framework for Matching Long Patterns

Recall that we deal with a string $T[0..n)$ of length $n$ over an alphabet $\mathrm{\Sigma }=[0,\sigma )$. We assume that $\mathrm{\Sigma }$ is an integer alphabet of polynomial size in $n$, i.e., $\sigma =n^{𝒪(1)}$. A substring of $T$ of length $k>0$ is called a $k$-mer. Given a query pattern $P[0..m)$ with $m\ge \mathrm{\ell }$ for some lower bound $\mathrm{\ell }>0$, the goal is to support the following operation:

We refer to the number of occurrences of $P$ in $T$, that is, $|\text{Locate}(P,T)|$, with $\text{Count}(P,T)$.

We assume that we have random read-only access to $T$ and count the space (in number of words) occupied on top of the space occupied by $T$. We assume the standard word RAM model of computation with machine words of $\mathrm{\Omega }(\log n)$ bits.

We use a specific class of randomized methods to sketch a string, called minimizers [34, 35]. Minimizers are defined as the triple $(k,w,h)$: from a window of $w$ consecutive $k$-mers of $T$, the leftmost smallest $k$-mer according to the order $h$ is chosen as the minimizer of the window. Since at least one $k$-mer must be chosen every $w$ positions, the fraction of sampled $k$-mers – defined as the density of the sampling algorithm – is always at least $1/w$. Several minimizer sampling algorithms have been proposed in the literature. See Section 3 of [19] for a recent overview of different sampling strategies and orders that lead to different densities. Here, however, we use the folklore random minimizer sampling, which is as defined above and uses a pseudo-random hash function for the order $h$. We have the following result.

We fix $\mathrm{\ell }$ to be the minimum pattern length and let $w=\mathrm{\ell }-k+1$. Each substring of length $\mathrm{\ell }$ of $T$ therefore contains one minimizer. (In practice, we expect to have $|P|\gg \mathrm{\ell }$ and that the sketch of $P$ is a sequence of several minimizers.) Further, we let ${ℳ}_{\mathrm{\ell },k}(T)$ indicate the sorted list of positions in $T$ of the minimizers of $T$. Let $z=|{ℳ}_{\mathrm{\ell },k}(T)|$ be the number of minimizers. By Theorem 1, we have that $z\approx 2(|T|-\mathrm{\ell }+1)/(\mathrm{\ell }-k+2)$ in expectation.

Given a set $𝒳$ of strings over the alphabet $\mathrm{\Sigma }$, a trie $\text{Trie}(𝒳)$ is a rooted tree whose nodes represent the prefixes of the strings in $𝒳$. The edges of $\text{Trie}(𝒳)$ are labeled by letters from $\mathrm{\Sigma }$; the prefix corresponding to node $u$ is denoted by $\text{str}(u)$ and is given by the concatenation of the letters labeling the path (sequence of edges) from the root of $\text{Trie}(𝒳)$ to $u$. The node $u$ is called the locus of $\text{str}(u)$. The parent-child relationship in $\text{Trie}(𝒳)$ is defined as follows: the root node is the locus of the empty string $\epsilon$; and the parent $u$ of another node $v$ is the locus of $\text{str}(v)$ without the last letter. This letter is the edge label of $(u,v)$. The order on $\mathrm{\Sigma }$ induces an order on the edges outgoing from any node of the trie. A node $u$ is branching if it has at least two children and terminal if $\text{str}(u)\in 𝒳$.

A compacted trie is obtained from $\text{Trie}(𝒳)$ by removing all nodes except the root, the branching nodes, and the terminal nodes. A compacted trie is thus a trie where all unary paths are collapsed into a single edge, labeled by the string obtained by concatenating all the letters of the edges in the unary path. The compacted trie takes $𝒪(|𝒳|)$ space provided that the edge labels are implicitly represented as pointers to fragments of strings in $𝒳$. Given the lexicographic order on $𝒳$ along with the lengths of the longest common prefixes between any two consecutive elements (in this order) of $𝒳$, one can compute $\text{Trie}(𝒳)$ in $𝒪(|𝒳|)$ time [26].

Let $p$ be a prime number and choose $r\in [0,p)$ uniformly at random. The rolling hash value of $T[i..j]$ – which we term fingerprint – is defined as [25]:

The adjective “rolling” refers to the way the hash value is updated incrementally as a fixed-length window slides through the string $T$. The function $\varphi$ allows one to compute the fingerprint of a window just knowing the fingerprint of the previous window and the character that is being removed/added, instead of recalculating the fingerprint from scratch.

By definition, if $T[i..i+\mathrm{\ell }]=T[j..j+\mathrm{\ell }]$, then $\varphi (T[i..i+\mathrm{\ell }])=\varphi (T[j..j+\mathrm{\ell }])$. On the other hand, if $T[i..i+\mathrm{\ell }]\ne T[j..j+\mathrm{\ell }]$, then $\varphi (T[i..i+\mathrm{\ell }])\ne \varphi (T[j..j+\mathrm{\ell }])$ with probability at least $1-\mathrm{\ell }/p$ [10]. Since we are comparing only substrings of equal length, the number of different possible substring comparisons is less than $n^3$. Thus, for any constant $c\ge 1$, we can set $p$ to be a prime larger than $\max (|\mathrm{\Sigma }|,n^{c+3})$ to make the function $\varphi$ perfect (i.e., no collisions) with probability at least $1-n^{-c}$ (this means with high probability). Any fingerprint of $T$ or $P$ fits in one machine word, so that comparing any two fingerprints takes $𝒪(1)$ time. In particular, we will use the following well-known fact.

All experiments were run on an Intel Core i7-10750H running at a fixed frequency of 2.6 GHz with hyperthreading disabled and cache sizes of 32 KiB (L1), 256 KiB (L2), and 12 MiB (shared L3). Code was compiled using rustc 1.85.0-nightly and GCC 14.2.1 on Arch Linux.

We use three textual datasets of different nature and alphabet size: (1) chromosome 1 of CHM13v2.0³³3https://www.ncbi.nlm.nih.gov/nuccore/NC_060925.1, which contains repetitive regions and consists of 248 million characters over the DNA alphabet ($\sigma =4$), thus 59 MiB when each character is coded using 2 bits; (2) the 200 MiB protein sequences available from the Pizza & Chili site⁴⁴4https://pizzachili.dcc.uchile.cl/texts/protein ($\sigma =27$), or 125 MiB when each character is coded using 5 bits; and (3) the 200 MiB English collection, also available from the Pizza & Chili site⁵⁵5https://pizzachili.dcc.uchile.cl/texts/nlang ($\sigma =239$). In the following, we discuss experimental results referring to Figure 4 (on page 4) for the DNA alphabet. The results for proteins and English have very similar shapes and are deferred to Appendix B due to space constraints.

For each dataset, we test ${10}^5$ positive queries that are uniformly sampled from the text. For DNA, queries consist of $512$ characters, while for the protein and English datasets, we use queries of $128$ characters.

We compare the suffix array and the FM-index, indicated with SA and FM-index in our results, against their U-index variants, with parameters $(k,\mathrm{\ell })\in \{(4,32),(8,64),(16,128),(28,256)\}$. For the suffix array construction, we use libsais⁶⁶6https://github.com/IlyaGrebnov/libsais [33, 23]. For the FM-index, we use the implementations in SDSL-lite (SDSL-v2) [16] and in AWRY [1].

For each index, we use the smallest $\tau \ge \lceil {\log }_{2}c\rceil$ that is supported by the index. In practice, this means that for the suffix array and SDSL FM-index we use $\lceil {\log }_{2}c/8\rceil$ bytes to represent each ID, so that each ID is a single character. This way, these indexes get built on exactly $z=|{ℳ}_{\mathrm{\ell },k}(T)|$ characters. In practice, this is strictly better than using a smaller $\tau$ and building an index on $2z$ or more characters. The AWRY FM-index does not support generic alphabets, and thus we had to consistently use $\tau =2$ to encode the IDs as DNA bases. We note that the AWRY FM-index uses a multi-threaded parallel construction algorithm, while all other methods are single-threaded. Further details can be found in Appendix A.

Finally, we also compare against our own implementation of the sparse suffix array at minimizer positions [17], that we call the sparse SA index in our results.

In all cases, the U-index variants take less space and are faster to construct than the classic indexes. Both space and construction time tend to be at least $8\times$ less. The time spent searching the suffix array (bottom, shaded) goes down as the suffix array becomes sparser, but the total query time goes up. This is due to the increasing number of false positive matches in minimizer space, starting at $10$ false positives per query for $(k,\mathrm{\ell })=(4,32)$ and going up to $100$ for $(k,\mathrm{\ell })=(28,256)$. These are caused by highly repetitive regions in the centromere [18]. Nevertheless, the query time is usually less than $2\times$ slower than the classic suffix array, while the space and construction time are greatly reduced.

When $k\in \{16,28\}$, most minimizers are unique, and the map $H$ mapping the minimizers to IDs has size linear in their number. This can be seen in that the shaded area in the rightmost two columns of the top-left group is accountable for most of the space used by the index at over $2^6=64$ MiB, while the input only takes 59 MiB when encoded using $2$ bits per character. Looking at SA -H, which omits $H$, we see that the index is much smaller for larger $k$, while the query time is unaffected. The construction time does increase though, since the alphabet is larger and hence requires more memory. The high spike for $k=16$ seems to be caused by SA (libsais) being particularly slow for 32-bit integers.

We can also omit the sketched text $S$ and instead reconstruct it on-the-fly, as explained in Section 4. This saves significant space when $k$ is large (where also $H$ can be omitted altogether). The construction time is not affected since $S$ is discarded after the construction. The time spent in searching the suffix array does increase significantly though, up to $4\times$ (shaded), but as most of the time is spent verifying potential matches, the total query time only goes up slightly.

The SDSL implementation takes around 90 MiB for the plain input text, and goes significantly below this when using the U-index⁷⁷7For some technical reason, SDSL does not accept NULL bytes, hence minimizers must necessarily be mapped by $H$ in the range $[1,c+1)$. Thus, SDSL takes a significant amount of space for $k\ge 16$.. The construction time improves as well, almost $2\times$ each time we double $\mathrm{\ell }$ because the number of sampled positions nearly halves.

A main drawback of the SDSL implementation is its significantly larger query times compared to all other methods, starting at $32\times$ slower for the plain index and increasing up to $1000\times$ slower for $\mathrm{\ell }=256$. We suspect this is due to the inherent complexity in the wavelet tree data structures used to represent the Burrows-Wheeler transform of the text [7], that has increasingly more levels as $k$ increases and hence the number of bits $\tau$ in each $k$-mer ID grows. This makes SDSL an impractical choice in this scenario or, more in general, for applications where the final size is not a bottleneck but fast queries are the primary concern.

AWRY uses around twice as much space as SDSL, and has a similar construction time, likely because both are limited by the internal suffix array construction. On the other hand, query times for AWRY are significantly faster than SDSL because AWRY is optimized for the DNA alphabet, and the U-index version with $k=4$ is slightly faster than AWRY on the plain text. However, as $k$ and $\mathrm{\ell }$ increase, both SDSL and AWRY slow down by roughly a factor of 2 overall. For AWRY we can also omit $H$ and this almost halves the size for large $k$, but negates the previously seen speedup in construction time since $S$ is significantly larger.

Lastly, the sparse SA index takes strictly less space than the U-index variant of the suffix array, since it stores strictly less information: only one permutation of the minimizer positions in the original text. It is also significantly faster to query, since no sketching is needed. Additionally, it has much fewer false positives, since comparisons are made in the plain input space rather than in sketch space. The construction time is worse than SA, however, since there is no available linear-time implementation for sparse suffix array construction⁸⁸8Indeed, our implementation is very simple and relies on sorting the minimizer positions using the Rust sort_by function that returns the suffix associated with each position as needed.. Nevertheless, this is still faster than constructing an FM-index.

We run the following experiment. As input, we take the full human genome (CHM13v2.0) and a prefix of 450 PacBio HiFi long reads from the HG002-rep1 dataset⁹⁹9Downloaded from https://downloads.pacbcloud.com/public/revio/2022Q4/HG002-rep1/.. These reads are approximately 99.9% accurate and have an average length of 16kbp.

We partition (chunk) each read into patterns of length 256, and search each of these patterns in our index. Short leftover suffixes are ignored. Ideally, each read has then at least one pattern that exactly matches the text, which can then be used to anchor an alignment.

We build the U-index on top of the libsais suffix array in the configuration where the sketched text $S$ and map $H$ are both stored. We use $k=8$ and $\mathrm{\ell }=128$, so that at least 2 and on average 4 minimizers are sampled from every pattern. This results in 53M minimizers, and the entire U-index is built in 12 seconds.

Out of the 450 reads, 445 have at least one matching pattern, and in total, 14 824 of the 28 243 patterns match (52%). As observed before, an issue with DNA is that it contains many long repetitive regions. In particular, there are 160 patterns that match around 3820 times each for a total of 611k matches, while the 28k remaining patterns only match 27k in total. Worse, there are 721M mismatches, i.e., candidate matches in sketch space that turn out not to be matches in the genome. Verifying these candidate matches takes over 98% of the time.

Only very few (<200) patterns have more than 10 matches, and thus we stop searching once we hit 10 matches. Further, most patterns have relatively few mismatches, while a few patterns have a lot of mismatches. To also avoid those negative effects, we generally only consider the first 100 matches in sketch space. We still match 445 of 450 reads, while the number of matched patterns goes down to 13 717. On the other hand, the number of mismatches is now 663k (23 per query) and the number of matches is 25k (0.9 per query).

The result is a query time of 8.7$\mu$s per pattern or 550$\mu$s per read, of which 33% is sketching the input, 18% is locating the sketch, and the remaining 48% is verification.