Search Schemes for Approximate Pattern Matching: An Overview

Quickly locating a given string (pattern) in a long sequence is a canonical algorithmic problem with numerous practical applications. In bioinformatics, for example, an important task is to locate thousands of reads generated by DNA sequencing in a genomic sequence. In this example, as well as in most other practical settings, it is necessary to allow for a certain number of errors, i.e., differences between the pattern and its occurrence in the sequence. This leads to the problem of approximate pattern matching (APM) that we address in this paper. A common algorithmic formalization, which we follow in this paper, assumes that we are given an error threshold $k$ and that we are looking for all occurrences of a given pattern up to $k$ errors in a text. The error model is defined either by the Hamming distance, which considers only character substitutions, or by the edit (Levenshtein) distance, which also allows character insertions and deletions.

Pattern matching is a problem that has been extensively studied. One class of solutions, relevant for most practical applications, is indexed pattern matching, where the text is preprocessed into an index data structure that supports efficient pattern search. Among such data structures, the FM-index, also known as the BWT-index, is particularly important. It was proposed by Ferragina and Manzini in 2000 [7]. Unlike other popular data structures such as suffix trees [43, 28, 40] and suffix arrays [26], the FM-index is a succinct data structure whose size (in bits) is asymptotically the same as that of the underlying sequence. Interestingly, the FM-index is built upon the Burrows-Wheeler transform [5], which is rooted in word combinatorics [8] and was initially proposed for text compression [27]. However, the idea proved extremely useful for indexing: the FM-index and its variants became integral components of many bioinformatics tools, such as Bowtie [19, 18], BWA and BWA-SW [20, 21], SOAP2 [22], Masai [37], Centrifuge [12], FMAlign [23, 44], Centrifuger [38], or ProPhyle [3].

In contrast to exact pattern matching, which can be solved very efficiently, approximate pattern matching with $k$ errors is a more challenging problem, as known algorithms are, in the worst case, exponential in $k$ with respect to either time or space requirements [39]. Considerable effort has been made to develop practically efficient APM algorithms. One proposed heuristic is seeding, based on the observation that approximately matching fragments typically share exact patterns of certain size. The idea is then to use small exactly matching patterns, or their combinations, as seeds for potential enclosing approximate matches. Another perspective on this approach is filtering where seeds are viewed as indicators of potential match locations in the sequence, filtering out irrelevant regions and narrowing the search. To efficiently detect seeds, all patterns in the text that could serve as seeds are indexed in a data structure, typically a hash table. The most prominent example of such algorithms is the BLAST alignment method, which was a predominant tool in bioinformatics for many years [13]. While seeding is generally not used for exhaustive search, it can also be applied to APM in a non-trivial way through the concept of spaced seeds [4, 24, 6]. This approach allows for defining one or several seeds such that any pattern occurrence contains at least one of them. However, computing spaced seeds that ensure a lossless search is a challenging problem, even for the Hamming distance [14]. For the edit distance, this approach becomes even more difficult to implement effectively.

In this paper, we focus on a different approach that leverages efficient indexing of the text via an FM-index. Although the standard method of pattern search using an FM-index scans the pattern backward (right-to-left), extensions of the FM-index, known as the bidirectional FM-index, support searching in both directions [17, 35, 36, 2, 1, 30]. This feature can be advantageous for approximate pattern matching, as first demonstrated by Lam et al. [17] for restricted cases of one or two errors. The idea of combining forward and backward searches for APM was also explored in [19] and [20], although those works employed a separate index for each direction.

In general, the idea involves partitioning the pattern into $k+1$ or more parts and breaking the search into independent searches, which can potentially run in parallel. Each search begins with one of the parts, allowing either no errors or a small number of errors. This approach drastically reduces the search space compared to the naive backtracking method where the pattern is searched in one direction while exploring all possible errors at each step. This idea is formalized through the concept of search schemes which we study in this paper. Search schemes constitute a flexible and versatile formalism for specifying search strategies. They can also be used to simulate some existing filtering methods, such as the one based on the pigeonhole principle, the suffix filter [10] or the 01*0 filter [42].

The goal of this paper is to provide an overview of search schemes, relate it to other methods, and present its current state-of-the-art. We will also survey some techniques for designing efficient search schemes. Finally, we provide some experimental results demonstrating and comparing the efficiency of different search schemes.

A string is denoted as $A=a_1{a}_2\mathrm{\dots }{a}_{|A|}$ where $|A|$ represents its length, and $a_i$ refers to the $i$-th character of the string. An infix (or substring) of $A$ is denoted as $A_{i,j}=a_i{a}_{i+1}\mathrm{\dots }{a}_j$. Given an error threshold $k$, the approximate pattern matching problem for a query pattern $P=p_1{p}_2\mathrm{\dots }{p}_{|P|}$ and a text $T=t_1{t}_2\mathrm{\dots }{t}_n$ is defined as finding the set of all substrings $\{T_{s_1,e_1},T_{s_2,e_2},\mathrm{\dots }\}$ such that $d(T_{s_i,e_i},P)\le k$ for some distance measure $d(\cdot ,\cdot )$. In this work, we consider either the Hamming distance $d_{\mathrm{ham}}$ or the edit distance $d_{\mathrm{edit}}$. The Hamming distance $d_{\mathrm{ham}}(X,Y)$ between two strings $X$ and $Y$ of equal length is the number of positions at which the corresponding symbols are different. The edit distance $d_{\mathrm{edit}}(X,Y)$ between two strings $X$ and $Y$ is the minimum number of operations required to transform $X$ into $Y$, where an operation is an insertion, deletion, or substitution of a single character.

A bidirectional index of a text is a data structure for pattern matching that supports incremental search of the pattern in both directions. This allows a pattern to be searched starting from any position, extending either to the left or right, and possibly alternating directions. More specifically, given an already matched substring $M$, the index supports extending that match to either $Mc$ or $cM$, where $c$ is a character.

Assume a query pattern $P$ is divided into $p$ non-overlapping parts $P=P_1\mathrm{\dots }{P}_p$. A search is a triplet $S=(\pi ,L,U)$. Here, $\pi =(\pi [1],\mathrm{\dots },\pi [p])$ is a permutation of $\{1,\mathrm{\dots },p\}$ defining the search order in which these parts will be matched. To make it consistent with bidirectional indices, it is required that $\pi$ fulfills the connectivity property: for each $i>1$, ${\pi }_i$ is either , or . $L=L[1]\mathrm{\dots }L[p]$ and $U=U[1]\mathrm{\dots }U[p]$ are non-negative integer sequences of length $p$ that specify respectively lower and upper bounds on the cumulative number of errors when searching for consecutive parts in the order specified by $\pi$. Formally, $L$ and $U$ must satisfy $L[1]\le \mathrm{\dots }\le L[p]\le k$, $U[1]\le \mathrm{\dots }\le U[p]\le k$, and $L[i]\le U[i]$ for all $i\le p$. A search scheme $𝒮$ is a collection of searches $𝒮=\{S_1,S_2,\mathrm{\dots },S_{|𝒮|}\}$. For ease of reading, sequences $\pi$, $L$ and $U$ are written as strings without separators, for example search $((2,3,1),(0,0,0),(0,1,2))$ is written as $(231,000,012)$.

An error configuration $e=(e_1,\mathrm{\dots },e_p)$ is a distribution of at most $k$ errors over the $p$ parts, i.e. ${\sum }_{i=1}^p{e}_i\le k$. A search $S=(\pi ,L,U)$ covers $e$ if $L[i]\le {\sum }_{j\le i}{e}_{\pi [j]}\le U[i]$ for all $i$. A search scheme is called lossless if each possible error configuration is covered by at least one of its searches. Furthermore, if each error configuration is covered by only one search, it is called non-redundant. In the case of Hamming distance, a lossless and non-redundant search scheme will lead to each match being found exactly once. In contrast, in the case of edit distance, different error configurations can actually correspond to the same match.

An example of a search scheme is given in Fig. 1. The scheme contains three searches for a pattern of length $8$ over a binary alphabet, with up to $k=2$ errors under the Hamming distance. The pattern is partitioned into three parts of length $3,3,2$ respectively. Note that scheme $S_2$ requires bidirectional search: after initially searching part $P_2$, part $P_1$ is searched backward and then $P_3$ is searched forward. The scheme simulates the well-known pigeonhole principle for APM when the pattern is partitioned into $k+1$ parts and an exact occurrence of each part is then extended into a potential match (see Section 3). Observe that the scheme is lossless but redundant.

Fig. 2 depicts another lossless search scheme for the search setting of Fig. 1. This scheme follows the algorithm of Lam et al. [17]. The third search $S_{\mathrm{bi}}$ covers error configuration $(1,0,1)$, the only one not covered by $S_{\mathrm{fwd}}$ and $S_{\mathrm{bwd}}$. One can observe that ${𝒮}_{\mathrm{Lam}}$ has a smaller search space than the search scheme from Fig. 1.

To gain intuition about what constitutes an efficient search scheme, we first consider naive backtracking as a baseline algorithm. Naive backtracking can be expressed as a search scheme with a single search ${𝒮}_{\mathrm{bt}}=\{(1,0,k)\}$, where pattern $P$ is not partitioned $(p=1)$ and is matched in its entirety, with a lower bound of $0$ and an upper bound of $k$ errors. Using an index such as the FM-index, candidate occurrences $O$ of the pattern $P$ in the text $T$ are incrementally enumerated character by character in a depth-first manner. For simplicity, we assume that $O$ is spelled from left to right. The number of errors between a (partial) candidate occurrence $O$ and (a prefix of) $P$ is tracked using a dynamic programming matrix $D$, where each matrix element $D(i,j)$ represents $d(O_{1,i},P_{1,j})$, i.e., the distance between the first $i$ characters of $O$ and the first $j$ characters of $P$. As soon as all values in row $|O|$ of $D$ exceed the threshold $k$, it is impossible to further extend $O$ to match $P$ within $k$ errors. In this case, the algorithm backtracks to the most recent branching point with unexplored characters and continues the search from there. If $P$ can be completely matched (i.e., $d(O,P)\le k$), then $O$ constitutes an approximate occurrence of $P$ in $T$ and its position(s) in $T$ can be reported. Note that in the case of the edit distance, maintaining a banded matrix $D$ with $2k+1$ entries per row suffices. In the case of the Hamming distance, it is sufficient to track only of the diagonal elements of $D$.

The naive backtracking algorithm induces a search tree where each node corresponds to a single character extension and, therefore, to a unique substring $O$ of the text that is enumerated. Since the FM-index enables character extensions in constant time, the runtime of naive backtracking is proportional to the number of nodes in this search tree. The size of the search tree, i.e., the total number of strings $O$ enumerated, is governed by two factors: (i) the string $O$ should exist in $T$, and (ii) either $d(O,P)\le k$ (and then $O$ should be reported as an occurrence), or it must still be possible to further extend $O$ to some $O^{\prime }$, such that $d(O^{\prime },P)\le k$. Criterion (i) is controlled by the FM-index which spells out only candidate occurrences that exist in $T$. For (ii), the (banded) dynamic programming (DP) matrix $D$ is used to track the distance between a (partial) candidate occurrence $O$ and pattern $P$. Using naive backtracking, the search tree grows rapidly with the maximum number of allowed errors $k$ and becomes impractically large even for modest values of $k$. This is because the region near the root of the search tree is densely branched. This is illustrated in Fig. 3 showing the trie of all substrings of the human reference genome. Each search enumerates a part of this trie, and a large number of short strings $O$ are enumerated only to find out that the vast majority of them cannot be extended to an actual (approximate) occurrence of $P$ in $T$. The reason that many short strings are enumerated is because short strings $O$ are likely to be present in $T$ and not (yet) exceed the threshold of $k$ errors.

One way to avoid exploring the densely branched region near the root of the substring trie is to use the pigeonhole scheme ${𝒮}_{\mathrm{ph}}$ illustrated in Fig. 1. If $k$ errors are allowed, the pattern $P$ is partitioned into $p=k+1$ parts $P=P_1\mathrm{\dots }{P}_p$. Regardless of how the $k$ errors are distributed across the parts, at least one part must be error-free. Consequently, each approximate occurrence $O$ of $P$ with at most $k$ errors must contain at least one such error-free part of $P$. The key idea is to first match this error-free part of $P$, and then to extend this match with the remaining parts of $P$ using backtracking, allowing up to $k$ errors in total. To identify all occurrences, this procedure must be performed $p$ times, each time assuming a different part of $P$ to be error-free. In each such search, the initial exact matching bypasses the densely branched upper levels of the substring trie, and the backtracking procedure is performed only in deeper levels which are more sparsely branched (nodes have fewer possible character extensions, see Fig. 3). This pigeonhole-based approach can be formally expressed as a search scheme with $p=k+1$ searches: ${𝒮}_{\mathrm{ph}}=\{S_1,\mathrm{\dots },S_{k+1}\}$, where

Formula (1) expresses that during search $S_i$, part $P_i$ is matched first with no errors allowed. Next, this seed is extended with parts $P_{i+1},\mathrm{\dots },P_p$ to the right, followed with parts $P_{i-1},\mathrm{\dots }{P}_1$ to the left, allowing up to $k$ errors in total. Collectively, searches $S_i$ of ${𝒮}_{\mathrm{ph}}$ will identify all approximate occurrences $O$ of $P$. However, the same occurrence may be reported by multiple searches. Indeed, any occurrence with more than one error-free part will be reported by multiple searches. The pigeonhole-based search has been applied in many earlier works (see e.g. [29]). Note that searches $S_i$ with $i=2\mathrm{\dots }p-1$ require bidirectional search that supports extending a partial occurrence $O$ with a single character $c$ either to the left ($cO$) or to the right ($Oc$).

Although conceptually simple, the pigeonhole search scheme can already outperform naive backtracking by a significant margin in practical applications. Yet, each search in the pigeonhole search scheme already admits the maximum number of $k$ errors starting from the second part that is searched. Intuitively, performance could benefit from only gradually allowing more errors when additional parts are matched. Indeed, deeper levels in the substring trie are associated with fewer branches (see Fig. 3), so admitting more errors is best postponed as much as possible.

Kärkkäinen and Na [10] proposed a suffix filter that enhances the pigeonhole method. The pattern is still partitioned into $p=k+1$ parts $P=P_1\mathrm{\dots }{P}_p$, but the search is now performed for strongly matching suffixes $P_i\mathrm{\dots }{P}_p$. A suffix $P_i\mathrm{\dots }{P}_p$ is said to strongly match a string $O=O_i\mathrm{\dots }{O}_p$ (suffix of a potential approximate match $O=O_1\mathrm{\dots }{O}_p$) if for any $j$ with $i\le j\le p$, the condition $d(P_i\mathrm{\dots }{P}_j,O_i\mathrm{\dots }{O}_j)\le j-i$ holds. It was proven in [10] that this filter is lossless, meaning that for any error configuration, there always exists a strongly matching suffix. Note that the distance function $dist$ can be either $d_{\mathrm{ham}}$ or $d_{\mathrm{edit}}$.

It is easily seen that the suffix filter can be simulated by a search scheme that first searches forward for a strongly matching suffix and then extends it backward to the beginning of $P$. The idea of combining suffix filter with the search on the FM-index was exploited in [25, 41, 16] in application to DNA read alignment problems. Formally, the search scheme for suffix filter is defined by ${𝒮}_{\mathrm{suf}}=\{S_1,\mathrm{\dots },S_{k+1}\}$, where

Another improvement of the pigeonhole-based approach was proposed by Vroland et al. [42]. Here, pattern $P$ is partitioned into $p=k+2$ parts instead of $k+1$. This guarantees that any approximate occurrence $O$ of $P$ contains at least two error-free parts. Among these pairs, there exists one separated by a sequence of zero or more parts each containing exactly one error [42, Lemma 2]. This (variable-length) sequence of parts is referred to as a 01*0 seed. Compared to the suffix filter, the seeding of the 01*0 filter requires searching for parts with no more than one error, which contributes to its efficiency. Unfortunately, a direct simulation of the 01*0 filter with a search scheme requires $(k+2)(k+1)/2$ searches. However, as observed by Pockrandt [31], one can simulate a weaker version by searching for an exactly matching part $P_i$ followed by part $P_{i+1}$ matching with at most one error. This requires only $k+1$ searches, where the last $(k+1)$-th search can be restricted to be followed by an exactly matching part $P_{k+2}$. In summary, this search scheme is defined by ${𝒮}_{01}=\{S_1,\mathrm{\dots },S_{k+1}\}$, where

Compared to the pigeonhole scheme (2) which uses $k+1$ parts, ${𝒮}_{01}$ uses $k+2$ parts which are thus slightly shorter, making the starting exact match of $P_i$ a slightly weaker filter. However, the requirement that the following part is matched with at most one error outweighs this weakness in many practical scenarios, as will be shown below in Section 4.

In their seminal paper, Kucherov et al. [15] provided search schemes with $k+1$ and $k+2$ parts for up to $k=4$ errors. They differ from ${𝒮}_{\mathrm{ph}}$, ${𝒮}_{\mathrm{suf}}$ and ${𝒮}_{01}$ in several aspects. First, they often contain a larger number of searches. For example, for $k=4$ errors and $p=k+1=5$ parts, ${𝒮}_{\mathrm{Kuch}}=\{S_1,\mathrm{\dots },S_8\}$ contains eight searches¹¹1Note that [15] has a different definition of lower bound, therefore those were adjusted to the definition we use in this paper. Furthermore, lower bounds in $S_6$ and $S_7$ were slightly corrected compared to [15].:

Increasing the number of searches reduces the search space per search and enables a more gradual increase in the allowed number of errors as additional parts are added. Indeed, in search scheme (12), most searches allow the maximum number of $k=4$ errors only from the fourth part that is matched. Kucherov et al. [15] define the critical string of a search scheme $𝒮$ as the lexicographically maximal $U$-string of a search in $𝒮$. The critical $U$-string in ${𝒮}_{\mathrm{Kuch}}$ (for $k=4$ errors and $p=5$ parts) is 02244 from search $S_1$. This is more favorable than 04444, the critical $U$-string for the pigeonhole-based search scheme ${𝒮}_{\mathrm{ph}}$. Second, the search schemes ${𝒮}_{\mathrm{Kuch}}$ were the first to exploit the lower bound $L$. By carefully introducing a minimum number of errors on certain parts of each search, overlap between searches is minimized, reducing the search space and reducing the number of redundantly reported occurrences. Nevertheless, the searches in ${𝒮}_{\mathrm{Kuch}}$ collectively still cover all error configurations and therefore maintain their lossless character. Third, the authors of [15] observed that the search pattern $P$ does not necessarily need to be divided into equal-length parts. In fact, the part lengths can be chosen arbitrarily, provided that the same part lengths are used across all searches. Since search $S_1$ is associated with the critical $U$-string, it will likely correspond to the largest search space (and thus the longest runtime) if $P$ is partitioned into equally sized parts. Therefore, it may be beneficial to slightly increase the size of part $P_1$ while slightly decreasing the sizes of the other parts. This adjustment reduces the search space for searches that match $P_1$ as their first, error-free part (i.e., searches $S_1$, $S_4$, and $S_8$), while increasing the search space for the remaining searches. Although some searches will become slower while others speed up, the overall effect is a net performance gain. A more formal justification of the benefit of uneven partitioning is given in [15], Section 3.2. Overall, the search schemes proposed in [15] proved to be very efficient for practical bioinformatics applications [34].