{'fraction_non_alphanumeric': 0.04050514742110313, 'fraction_numerical': 0.013925265150532894, 'mean_word_length': 4.220522539596589, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FAMOUS: Fast Approximate string Matching using OptimUm search Schemes

Finding approximate occurrences of a pattern in a text using a full-text index is a central problem in bioinformatics. Bidirectional indices have opened new possibilities for solving the problem as they allow the search to be started from anywhere within the pattern and extended in both directions. In particular, use of search schemes (partitioning the pattern into several pieces and searching the pieces in certain orders with bounds on the number of errors in each piece) has shown significant potential in speeding up approximate matching. However, finding the optimal search scheme to maximize the search speed is a difficult combinatorial optimization problem. In this paper, we propose, for the first time, a method to solve the optimal search scheme problem for Hamming distance with given number of pieces. Our method is based on formulating the problem as a mixed integer program (MIP). We show that the optimal solutions found by our MIP significantly improve upon previously published ad-hoc solutions. Our MIP can solve problems of considerable size to optimality in reasonable time and has the attractive property of finding near-optimal solutions for much larger problems in a very short amount of time. In addition, we present FAMOUS (Fast Approximate string Matching using OptimUm search Schemes), a bidirectional search (for Hamming and edit distance) implemented in SeqAn that performs the search based on the optimal search schemes from our MIP. We show that FAMOUS is up to 35 times faster than standard backtracking and anticipate that it will improve many tools as a new core component for approximate matching and NGS data analysis. We exemplify this by searching Illumina reads completely in our index at a speed comparable to or faster than current read mapping tools. Finally, we pose several open problems regarding our MIP formulation and use of its solutions in bidirectional search.


Introduction
Finding approximate occurrences of a string in a large text is a fundamental problem in computer science with numerous applications. The approximate string matching (ASM) problem for Hamming distance considered in this paper is defined as follows: Given a number of mismatches K, a string (here referred to as a read ) of length R, and a text of length T , composed of characters from an alphabet of size σ, find a substring of the text whose Hamming distance to the read is at most K. A similar definition can be provided for ASM for edit distance where in addition to mismatches, insertions and deletions, are also considered.
Solving the ASM problem has become especially important in bioinformatics due to the advances in sequencing technology during the last years. The mainstream second generation sequencing techniques like Illumina produce reads of length 150-250 with an error rate of about 1%, mostly substitutions caused by the sequencing technology. Other sequencing technologies, e.g., Pacific Bioscience or Oxford Nanopore, produce much longer reads but with a higher error rate (in the range of 15%) containing both substitutions and insertions/deletions. A standard problem is to map the reads back to a reference genome while taking into account the errors introduced by the sequencing technology as well as those caused by biological variation, such as SNPs or small structural variations. Such a problem is almost always modeled as the ASM problem for Hamming or edit distance.
There are two main algorithmic strategies to address the ASM problem for large input sizes (in number of reads and size of the text): filtering and indexing. In this work, we focus on the indexing approach. Here, the main idea is to preprocess the reference sequence, the set of reads, or both, in a more intricate way. Such preprocessing into full-text string indices has the benefit that we usually do not have to scan the whole reference, but can conduct queries much faster at the expense of larger memory consumption. String indices that are currently used are the suffix array [10] and enhanced suffix array [1] as well as the FM-index [3], a data structure based on the Burrows-Wheeler Transform [2] and some auxiliary tables. For an in-depth discussion see [14]. Such indices are used to compute exact matches between a query and a text as a subroutine in backtracking approaches. For the ASM problem for Hamming or edit distance, the existing algorithms all have exponential complexity in K (e.g. [5,16]), and are hence only suited for small K.
Lam et al. [7] introduced bidirectional indices to speed up ASM for Hamming distance. For the cases K = 1 and 2, they partitioned the read into K + 1 equal pieces, and argued that performing approximate matching on a certain combination of these pieces in a bidirectional index amounts to faster approximate matching of the whole read. This combination is such that all possible mismatch patterns, i.e., all possible distributions of K mismatches among the pieces, are covered. The main idea behind improved speed is that a bidirectional index not only can start the search from the beginning (or end) of the read, but also from the beginning (or end) of any of the pieces. Therefore, we can start the search from a middle piece and then expand it to the left or right into adjacent pieces in any order we like. By choosing multiple appropriate orderings of pieces for this purpose, we can perform a much faster ASM compared to a unidirectional search because we can enforce exact or near-exact searches on the first pieces in the partition, significantly reducing the number of backtrackings, while using different orderings of pieces to ensure all possible mismatch patterns are still covered.
Kucherov et al. [6] formalized and generalized this idea by defining the concept of search schemes. Assume a read can be partitioned into a given number of pieces, denoted by P (not necessarily equal to K + 1). The pieces are indexed from left to right. A search scheme S = {(π s , L s , U s ), s = 1, . . . , S} is a collection of S searches, where each search s is designated by a triplet (π s , L s , U s ). π s is a permutation of 1, . . . , P and denotes the order in which the pieces of the partition are searched in search s. If π s,i = j, then piece j is searched at position i in the order (shortly referred to as iteration i in this paper). Due to the way a bidirectional index works, the permutation π s must satisfy the so-called connectivity condition, i.e, a piece j can appear at iteration i > 1 in the permutation only if at least one of pieces j − 1 or j + 1 have appeared at an iteration before i. L s and U s each are strings of P numbers. L s,i , is the lower bound on the cumulative number of mismatches allowed at iteration i of search s, and U s,i is the upper bound on this value.
Having this formal framework, the answer to the Optimal Search Scheme problem, defined as follows, can potentially have a great impact on improving the running time of ASM using a Bidirectional index (ASM-B).
Optimal Search Scheme Problem: What is the search scheme that minimizes the number of steps in ASM-B while ensuring all possible mismatch patterns are covered?
It turns out that this is a very difficult combinatorial optimization problem due to several reasons: There are a large number of attributes that define a solution (including S, P , size of each piece, and (π s , L s , U s ) for each search) with a large number of possibilities for each attribute; the solution must satisfy complex combinatorial constraints; and, calculating the objective function, i.e., number of steps in the ASM-B algorithm, for a given solution is complicated.
Kucherov et al. [6] presented some interesting results contributing initial insight into this key problem. More specifically, they assumed the number of steps in the ASM-B algorithm with a given search scheme, is a constant factor of the (weighted) total number of substrings enumerated by the algorithm in all searches. Assuming that a randomly generated read is to be matched to a randomly generated text, they presented a method to calculate this objective function for a given search scheme. They then showed that unequal pieces in the partition can potentially improve the objective function compared to equal pieces, and presented a dynamic programming (DP) algorithm that for a single pre-specified search, with given P and (π, L, U ), finds the optimal sizes of pieces assuming that we only calculate the objective function as the total number of substrings up to a limited length (justified by total randomness of the read and the text); see [6] for more details. In fact, the superiority of this DP over explicit enumeration is only due to this assumption. Nevertheless, this DP is very inefficient, and most importantly, it only finds the optimal piece sizes for a pre-specified search. In other words, it does not address the problem of finding an optimal search scheme which calls for determining S and all attributes of each search in the search scheme, and ensuring that they cover all mismatch patterns.
Kucherov et al. [6] also presented solutions for another limited problem, i.e., lexicographically minimizing the lexicographically maximal U string (critical U string) in a search scheme, only for P = K + 1 or K + 2 and assuming that the L strings for all searches contain only zeros. The usefulness of these solutions is justified by the high probability that the search with the critical U string has the largest share in the objective function; see [6] for details. Again from the perspective of finding an optimal search scheme, this result has similar limitations. Only one of the attributes (U ) of one of the searches for two specific values of P are optimized by fixing all L strings, which is far from designing a globally optimal search scheme as defined above. Consequently, in their computational experiments, Kucherov et al. [6] use a greedy algorithm based on this limited result to construct search schemes with unknown quality and only optimize the piece sizes for these schemes using their DP.
In this paper, for the first time, we propose a method to solve the optimal search scheme problem for ASM-B with Hamming distance, for any given P and equal-size pieces. Our method is based on a novel and powerful mixed integer linear program (MIP) that gets K, R, P , and an upper on S, denoted by S, as input, and provides, as its solution, all the attributes of the exact optimal search scheme (MIP methodology for optimization has been addressed in many references such as [11,17]). We present our MIP and provide the results of our computational study on the characteristics of its optimal solution and its running time, for different values of its input parameters.
Next we present FAMOUS (Fast Approximate string Matching using OptimUm search Schemes), a bidirectional search (for Hamming and edit distance) implemented in SeqAn [13] that performs the search based on optimal search schemes we have obtained from our MIP. The bidirectional search in FAMOUS uses a recent fast implementation of bidirectional indices [12] based on EPR dictionaries. We show that, for practical ranges of various input parameters, the number of substrings for the optimal search schemes found by our MIP can reduce to as small as half the number of substrings in the unidirectional complete backtracking. More importantly, we show that with our optimal search schemes, FAMOUS improves the search time for approximate string matching dramatically. Searching for all occurrences of Illumina reads on human genome sequence using FAMOUS is up to 3, 14, and 35 times faster than unidirectional search for K = 1, 2, and 3 Hamming distance errors, respectively. Although our MIP finds the optimal search schemes for Hamming distance, when we used its optimal schemes with the edit distance, we got similar improvements.
Additionally, we perform two other experiments using FAMOUS: In the first one, we search for all approximate occurrences of Illumina reads in the genome (all-mapping), and in the second one, instead of all-mapping, we do the search according to the typical strata-based strategy, i.e., we search for all best occurrences plus a number of errors up to a global error upper bound. We show that FAMOUS in strata mode is already competitive with Yara [15]. The all-mapping results suggest that a combination of our optimal search schemes and a verification strategy might lead to a tool outperforming all current ASM approaches for genomic reads. We note that the optimal searches found by our MIP are implemented in the SeqAn library for efficient sequence analysis, and hence, can in the future speed up any tool part of which includes searching for approximate occurrences of a read.
We will introduce our MIP for solving the optimal search scheme problem in Section 2. In Section 3, we first discuss our computational studies on solving our MIP, and then present the impact of using its optimal search schemes on solving ASM-B with FAMOUS using various realistic scenarios. We will conclude in Section 4 by raising several open problems and possible extensions, which provide exciting directions for future research stemming from this work.

Solving Optimal Search Scheme Problem using MIP
In this section we provide our MIP-based methodology for finding optimal search schemes after presenting some preliminaries. Our MIP is for Hamming distance, but as mentioned before, based on our computational experiments (Section 3.2), its optimal schemes for Hamming distance are also very good (but not necessarily optimal) search schemes for the edit distance as well.

Preliminaries
Our MIP presented in Section 2.2 will solve the optimal search scheme problem assuming P is given as an input (is not a decision variable in optimization) and all P pieces of the partition are equal in length, i.e., R = mP , where m denotes the length of any piece. Note that these assumptions pose no practical restrictions. Solving problems which include P as a decision variable and allow unequal pieces is part of our future research plan. Given the upper bound on Hamming distance (maximum number of mismatches) K as an input, a mismatch pattern is a particular distribution of h mismatches among the P pieces, for any h ≤ K. Specifically, the mismatch pattern q is a string of P integers a q,1 . . . a q,P such that a q,j ∈ {0, . . . , min{m, K}} for j = 1, . . . , P , and P j=1 a q,j = h. For given K and P , we denote the set of all possible mismatch patterns by M. Note that if K ≤ m then |M| = K h=0 h+P −1 h . Given a search s = (π s , L s , U s ), a mismatch pattern q is said to be covered by s if at every iteration i = 1, . . . , P of s, L s,i ≤ i t=1 a q,πs,t ≤ U s,i , i.e., the cumulative number of mismatches up to iteration i is between the allowed lower and upper bounds of search s. A search scheme S is feasible if and only if every mismatch pattern in M is covered by at least one search in S.
A search scheme can be visualized by representing each of its searches as a trie that captures all substrings enumerated by the search. Each edge at a level of the trie corresponds to a character of the alphabet at that level of search. A vertical edge represents a match, and a diagonal edge represents a mismatch. Fig. 1(a) shows the tries associated with the search scheme presented by Lam et al. [7] for K = 2 and P = 3, S Lam , applied on the six-character read "abbaaa" from alphabet {a, b} (note that the tries are slightly different from the ones given in [6], which contains a small error). Fig. 1(b) shows a search scheme with a single unidirectional search (complete backtracking), S Uni , for the same problem, and Fig. 1(c) shows the optimal search scheme, S Opt , found by our MIP, for the same problem. Each one of the three schemes in Fig  The total number of edges in SOpt (optimal number of edges) is 59, which is less than that of SUni, and significantly less than that of SLam. As shown in Section 3.2, for bigger problems, the reduction in the total number of edges of the optimal search scheme found by our MIP compared to the unidirectional search is much more significant (up to 50%).
of length R, the search enumerates all substrings in the tries; hence, the performance of the search scheme can be measured by the total number of edges in the search scheme. Otherwise, only a subset of the substrings in the tries will be enumerated depending on whether they occur in the text or not.
To address the performance measure in this latter case, Kucherov et al. [6] assumed the read and the text are randomly and independently drawn from the alphabet according to a uniform distribution, and hence, calculated the expected number of substrings enumerated by the scheme as the sum, over all non-root nodes of the tries, of the probability that the corresponding substring appears in the text. As a result, they presented a weighted sum of number of edges as the measure of performance. Due to the assumption of complete randomness and independence of the read and the text, they show that the weights of the edges at levels lower than log σ T + c σ of the tries, where c σ is that ((σ − 1)/σ) cσ is sufficiently small, are almost zero meaning that they can be dropped from the weighted summation.
For the main application of our interest, i.e. ASM of DNA sequence reads to reference genomes, the assumption of randomness and independence of the read and the text is far from reality. Calculating the expected number of substrings enumerated by a scheme calls for significant more study on determining probabilities that DNA sequence reads of particular length from a sample occur in the reference genomes. As currently there is no trivial answer to this problem, in this paper, we use the same performance measure of total number of edges in the tries of the search scheme even for the case where not all substrings occur in the text. Of course, our MIP can be easily modified to incorporate any other weighting scenario which might be proposed in the future.
Adapting the method from [6], the total number of edges in the search scheme is calculated by where n s,l,d is defined as the number of edges at level l of the trie of search s that end at nodes corresponding to substrings with d cumulative mismatches up to that level. The value of n s,l,d can be calculated using the following recursive equation, which is an adaptation of the formula in [6]: where, by definition, n s,0,0 = 1, n s,0,−1 = 0 and n s,0,d = 0, for d ≥ 1, s = 1, . . . , S, and L s,l and U s,l denote the smallest and largest cumulative number of mismatches that can occur at level l of the trie of search s, respectively, calculated as L s,l = max{L s, l/m −1 , L s, l/m − m l/m + l} and U s,l = min{U s, l/m , U s,l−1 + 1}. Here l/m , the smallest integer greater than or equal to l/m, would be the index of the iteration in which level l falls, and by definition, L s,0 = U s,0 = 0, for s = 1, . . . , S. For example, for search s bi of S Opt , we have L s bi = (0, 0, 0, 1, 1, 1) and U s bi = (0, 0, 1, 1, 2, 2).

MIP Formulation of Optimal Search Scheme Problem
Our MIP formulation, presented below, solves the optimal search scheme problem assuming P is given as an input and pieces are all equal in length. More specifically, for given K, R, P , and S, this MIP finds the search scheme with minimum total number of edges among all feasible search schemes that have at most S searches. The optimal solution to the MIP provides the (π, L, U ) of all searches in the optimal search scheme. The objective value of this optimal solution provides the minimum total number of edges (substrings) achievable among all feasible search schemes.
subject to The objective function (3) minimizes the total number of edges as calculated by (1) with n s,l,d as defined before. The binary variables x s,i,j capture the assignment of pieces to iterations, i.e., x s,i,j = 1 if piece j is searched at iteration i of search s, and x s,i,j = 0 otherwise. We define x s,i,0 = x s,i,P +1 = 0 to simplify presentation of constraints. At optimality, these variables determine the π s values for the optimal search scheme. Constraints (4a) and (4b) make sure that for any search s, only one piece is assigned to an iteration and only one iteration is assigned to a piece.
Constraints (5a)-(5b) ensure the connectivity of the pieces and are in fact linearization of the following constraint using auxiliary binary variables t + s,i,j and t − s,i,j : x s,h,j−1 = 2 for all s and i = 2, . . . , P − 1, which is one way to enforce connectivity of pieces. The term i h=1 x s,h,j will have a binary value which denotes whether or not piece j has been searched at any of iterations 1 to i of search s. The term i h=1 x s,h,j−1 captures the same notion for piece j − 1. If at any iteration all searched pieces form a connected block on the read, the value of i h=1 x s,h,j − i h=1 x s,h,j−1 will be equal to 1 only for one j, −1 for another j, and 0 for all other j's, which is ensured by (10), and hence its linearization.
Constraints (6a)-(6c) enforce calculation of n s,l,d based on the recursive equation (2) with the help of binary variables z s,l,d and z s,l,d . Due to (6a), if d ≥ L s, l/m − m l/m + l, then z s,l,d = 1, and due to (6b), if d ≤ U s, l/m , then z s,l,d = 1. Calculation of equation (2) is then enforced by (6c). When z s,l,d = z s,l,d = 1, (6c) reduces to n sld − n s,l−1,d − (σ − 1)n s,l−1,d−1 ≥ 0, which implies n sld − n s,l−1,d − (σ − 1)n s,l−1,d−1 = 0 since the objective function is to be minimized. If any of z s,l,d or z s,l,d is equal to 0, (6c) does not enforce anything as − l d (σ − 1) d is a lower bound on the right-hand side of (6c). Constraints (7a)-(7b) ensure L s,i and U s,i are non-decreasing as they are cumulative values. Constraints (8a)-(8b) ensure feasibility of the search scheme. λ q,s is a binary variable designating whether or not mismatch pattern q is covered by search s. Constraint (8a) forces λ q,s = 0 if search s does not cover mismatch pattern q and constraint (8b) ensures every mismatch pattern q is covered by at least one search, for q = 1, . . . , |M|.
Constraints (4a)-(9c) are enough to formulate the MIP; however, we have noticed that imposing the additional constraints x t,1,k ≤ (S − s + 1)(1 − x s,1,j ) for all s and j = 2, . . . , P (11b) x sij + P j=i x sij = 1 for all s and i ≥ P/2 + 1 (12) strengthens the formulation while preserving at least one optimal solution, resulting in faster solution time for the MIP. Constraints (11a) and (11b) eliminate some symmetry in the solution space. For every search scheme, there is an equivalent search scheme obtained by reversing all π s , s = 1, . . . , S. Constraint (11a) eliminates one of these two equivalent solutions in each pair by forcing piece P to be assigned to iteration P in the first search, eliminating the solutions in which piece 1 is assigned to iteration P . For any search scheme, another equivalent search scheme can be obtained by permuting the indexes of searches within the scheme. Existence of only one of the search schemes obtained by this index permutation in the feasible solution set is enough. This can be achieved by sorting (in ascending order) the searches based on the piece assigned to their first iteration. This is done by constraint (11b), which does not allow pieces 1, . . . , j − 1 to be assigned to the first iteration of searches s, . . . , S if piece j is assigned to the first iteration of search s. In addition to symmetry elimination, notice that the connectivity condition of pieces implies that the piece assigned to iteration P is either piece 1 or piece P , and in general, the piece assigned to iteration i ≥ P/2 + 1 is one of pieces 1, . . . , P − i + 1, i, . . . , P . Constraint (12) enforces this property, which strengthens the formulation, and according to our computational tests, reduces the running time of the MIP.

Remark 1.
A powerful feature of our MIP is that S is an upper bound on the number of searches, i.e., our MIP finds the optimal search scheme among all schemes with at most S searches. In our MIP, variables are defined for S searches, and if the optimal search scheme has S * < S searches, our MIP generates S − S * empty searches, i.e. searches in which L s,i > U s,i for some i.

Computational Evaluation
In this section, we present results of our computational studies on solving our MIP and using FAMOUS for solving ASM. Section 3.1 provides results on the MIP optimal objective value as a function of its input parameters, the running time of solving MIP, and speed of its convergence to the optimal solution. Section 3.2 presents the computational advantages gained in ASM by using FAMOUS.

Solving MIP
We used CPLEX 12.7.1 solver [4] to solve our MIP by implementing the code in C++ using CPLEX Callable Library. All instances were run over four 28-core nodes (2.4 GHz Intel Broadwell) with 64GB of memory per node. We ran our MIP solver for instances generated for a broad range of parameters K, R, S, and P and gave each instance a 3-hour time limit. Fig. 2(a) is a small representative of our results. It shows the optimal objective value (total number of edges) for R = 100, K = 1, . . . , 4, P = 5, 6, and S = 1, . . . , 5. If the problem is not solved to optimality in 3 hours, the best solution found within this time limit is shown. The optimal objective value does not show a consistent change pattern in terms of change in P ; however, as expected, it increases as K increases, as R increases (not shown), and as S decreases. In all instances, the optimal objective value shows a sharp drop from S = 1 to S = 2, then a modest drop to S = 3, and negligible change beyond S = 3, generating empty searches in many cases. Therefore, as long as S = 5, it is advisable to use S = 3 instead if we would like to reduce the MIP running time and still find an optimal or near-optimal solution for S = 5. We also noticed that the optimal search scheme obtained by our MIP is not sensitive to the value of R (see open problems in Section 4). Therefore, when R is large, it is advisable to solve the MIP for a much smaller reasonable value of R, e.g., R = KP , in order to get a solution that is most probably optimal for the large R in a much shorter amount of time.
Using the MIP formulation, we were able to solve considerable size problems to optimality. For instance, we were able to solve a problem with K = 4, R = 100, P = 3, and S = 3 to optimality in 5802 seconds. However, more complicated cases reached the time limit of 3 hours without proving solution (b) Fig. 2. (a). Sensitivity of optimal objective value to parameters R, K, S, and P . For some cases due to memory overflow, there is no data point. (b) Rapid convergence of feasible solutions to the optimal solution.
optimality. Consequently, it is important to investigate the rate of convergence of the solutions found during execution of MIP to the the optimal solution. Fig. 2(b) illustrates the ratio of the best solutions found by MIP during its execution to the final optimal objective value plotted against running time for some instances which reached optimality. We observe that in all cases, within 0.1% to 1% of the total running time, the MIP finds a solution which is finally proved to be optimal or very close to the optimal after MIP execution is complete. In other words, the MIP solver finds optimal or near optimal solutions very early on and spends the rest of its time ensuring that no better solution exists. This can be partly due to the remaining symmetry in the solution space. Nevertheless, from practical perspective, this is an attractive property because, when the input parameters are much larger, we can run the MIP for a short time and find solutions which are most probably optimal or near-optimal.

Computational Power of FAMOUS
In this section, we discuss the computational advantages achieved by FAMOUS. While FAMOUS is completely independent of alphabet size and read length, we chose to concentrate on parameter values relevant to standard sequencing reads, e.g., Illumina reads.
In Table 1, for a number of relevant parameter values, we have shown how the total number of edges using the optimal search schemes found by our MIP is reduced compared to the unidirectional backtracking scheme. It can be seen that the reduction is between 42% and 49%. Also, for K = 2 and K = 3, the optimal search scheme with P = K + 2 has fewer edges than the one with P = K + 1. Table 1. Total number of edges in the optimal search schemes found by our MIP for K = 1, 2, 3 and P = K + 1 compared to full backtracking. The factor column shows the ratio of total number of edges in each scheme to that in backtracking. The optimal search schemes are as follows:

Distance
Search Scheme Although the reduction factors in total number of edges obtained by our optimal search schemes in Table 1 are very significant in themselves, due to the stochastic nature of occurrence of errors in sequencing reads and occurrence of approximate matches in the reference genome, the real-case ASM speed-up factors achieved by these optimal search schemes compared to backtracking can be yet much more significant. To gain insight into this speed-up and also compare FAMOUS with other read mapping tools, we performed two experiments searching real Illumina reads in the human genome hg38 and compared the running time of FAMOUS with unidirectional backtracking for Hamming and edit distance, and with the read mappers Bowtie2 [8], BWA [9], and Yara [15], as well as an available implementation of the 01 * 0-filter scheme (named Bwolo) [16], for edit distance only as all these tools work for edit distance. All of our tests were conducted on Debian GNU/Linux 7.1 with Intel Xeon E5-2667V2 CPUs at fixed frequency of 3.3 GHz to prevent dynamic overclocking effects. All data was stored on tmpfs, a virtual file system in main memory to prevent loading data just on demand during the search and thus effecting the speed of the search by I/O operations. All tools were run with a single thread to make the results comparable.
Before addressing the results of our two experiments, we would like to note that while Bowtie2, BWA, and Yara are standard read mapping tools for Illumina reads, Bwolo is interesting in the context of our approach. Vroland et al. [16] presented this fast method for searching in an index by partitioning the read into K+2 pieces for K errors and then exploiting the fact that the K+2 pieces must contain the mismatch pattern 01 * 0. Their method searches for an occurrence of the 01 * 0 mismatch pattern in the read and then verifies the remaining pieces of the read partially in the index (they use a unidirectional index) and partially in the text. Interestingly, their 01 * 0 method can be formulated as a search scheme. For example, for K = 2, and hence P = 4, their method can be expressed as the following search scheme: S 01 * 0 = {(4321, 0000, 0122), (3214, 0000, 0122), (2134, 0000, 0022)}. The optimal search scheme found by our MIP for these parameters is S Opt = {(4321, 0002, 0122), (3214, 0000, 0112), (2134, 0011, 0022)}, which is quite similar, but has yet fewer edges. Indeed, we verified the practical superiority of S Opt to S 01 * 0 by searching 100, 000 real Illumina reads in the human genome with K = 2 Hamming and edit distance and noticed that S Opt took less time than S 01 * 0 (18.8s vs. 22.0s for Hamming, and 112s vs. 160s for edit distance, respectively).
In our first experiment, we searched for all approximate matches of 100, 000 Illumina reads of length R = 101 (SRA accession number ERX1959065) for K = 1, 2, 3. The results are shown in Table  2. For Hamming distance, we compared FAMOUS with optimal search scheme obtained by our MIP for P = K + 1 and P = K + 2 with the unidirectional backtracking scheme. For the edit distance, we did a similar comparison as well as comparison with the aforementioned read mapping tools. BWA was run with the options -N -n <K>, and Yara was run with the options -e <K> -s <K> -y full -t 1. Bowtie2 is not designed with all-mapping in mind (for our dataset, it did not terminate in 3 hours with default configuration and -a option). Moreover, imposing an all-mapping with maximum K errors in Bowtie2 in a way that its results are comparable to FAMOUS and other tools is difficult. Bowtie2 settings used in [16] do not enforce this, and nonetheless, led to a very long running time for our dataset (did not terminate in 13 hours). Consequently, we removed Bowtie2 from further study. Table 2. Running time comparison of searching all approximate matches of 100, 000 Illumina reads (R = 101) using FAMOUS with P = K + 1 and P = K + 2 versus backtracking for Hamming distance, and versus backtracking scheme, BWA, Bwolo, Yara for edit distance. The factor column is the speed-up ratio versus backtracking in each category.

Distance
Search Tool We can see in Table 2 that for both Hamming and edit distance, FAMOUS, which uses our optimal search schemes, is much faster than backtracking, verifying our expectation. The respective speed-ups for K = 1, 2, and 3 are 3.1, 14.3, and 35.2 for Hamming distance, and 4.1, 11.1, and 21.3 for edit distance, much more significant than reduction in the total number of edges reported in Table 1.
Comparing FAMOUS with BWA, Bwolo, and Yara, we see that FAMOUS is faster for K = 1, comparable for K = 2, and much worse for K = 3. This shows that for larger K and this read length the verification in the index is too costly, especially if the number of successful verifications is low (in our case almost all reads occur only once in the genome). Although Vroland et al. [16] raised the option of using a bidirectional index for verification, they only used standard verification for the last pieces as they had only a unidirectional index at hand. Interestingly, for their dataset (40bp, exactly 3 errors), we outperform Bwolo by a factor of almost 1.5 (data not shown), so the read length matters.
In our second experiment, we performed our search using the strata-based strategy as it is done in many read mappers. The 0-strata search means we first search the reads with 0 errors, then search all the reads with no exact match, with 1 error, and so on, until K is reached. This strategy can be generalized to s-strata, where s ≤ K. This means that, for b = 0 to K − s, for all reads with a b-error best match, we compute all occurrences with up to b + s errors.
In Table 3, we compare (for edit distance) the running time of FAMOUS for P = K + 2 in 1-strata mode with Yara in the same mode using -e <K> -s 0 -y full -t 1. We can see that FAMOUS is about as fast as Yara for K = 1 and 2, and slightly worse for K = 3, which is consistent with the results of our first experiment. The difference in the 1-strata mode is that we only have to map about 7% of the 100, 000 reads with 3 errors, which does not affect the overall running time as heavily as in the first experiment.

Conclusions
In this paper, we proposed for the first time a method to solve the optimal search scheme problem for ASM-B for Hamming distance (using a MIP formulation). In addition, we presented FAMOUS (Fast Approximate string Matching using OptimUm search Schemes), a bidirectional search (for Hamming and edit distance) implemented in SeqAn [13] that performs the search based on optimal search schemes obtained from our MIP and showed that FAMOUS is up to 35 times faster than standard backtracking. This speed-up makes it, for example, possible to search Illumina reads directly in a bidirectional index with up to 2 errors with times comparable with standard read mappers. Although FAMOUS includes an optimal search scheme that is superior to the 01 * 0 scheme of Vroland et al. [16] for search in the index, for a larger number of errors, one has to take into account the number of remaining verifications versus the number of edges in the trie for the remaining pieces. If that ratio is low, it does not pay to verify in the index as our comparison with Bwolo showed. We see the implementation of such a realistic read mapper based on FAMOUS as future work. Our approach in this research raises some interesting open problems: 1. Our computational experiments in Section 3.1 showed that our current MIP has two attractive properties: the early solutions it finds are optimal or near-optimal, and its optimal search scheme is insensitive to the value of R (we ask: "is this insensitivity to R a theoretically provable fact?"). This makes our current MIP quite powerful in practice because, even if all input parameters K, R, P , S are quite large, we can run the MIP for a short time with a much smaller R to get a solution that is most probably optimal or near-optimal for the original problem. Nevertheless, solving the MIP completely to ascertain optimality is of great interest and currently consumes considerable computational resources for large instances, especially when S > 5, K > 4, P > 6, R > 100. We ask "can the solution time be improved by introducing other MIP formulations, or strengthening the current formulation using strong cutting planes or further elimination of symmetric solutions?" 2. Our current MIP is restricted to optimizing over equal-size pieces, with P and S given as part of the input. We ask "are there (MIP) approaches to optimize over the number of pieces in the partition, and/or unequal piece lengths, and/or with no upper bound on the number of searches?" 3. Our optimal search scheme, S Opt , for the example of Lam et al. [7] covers every possible mismatch pattern only once, ensuring that no duplicate computational effort is spent on the same mismatch pattern. We, however, believe that this is not always the case for every optimal search scheme. We ask "how would enforcing this requirement on the solution of MIP affect its solution time and the performance of optimal search schemes it obtains?" 4. We demonstrated that solutions found by our MIP, although for Hamming distance, perform very well for the edit distance as well. We ask "are there (MIP) approaches to find the actual optimal search scheme for the edit distance?" 5. We demonstrated that the verification of few occurrences with high errors in the index is worse than a verification in the text. We ask "what is the best point to stop verification in the index and start verifying in the text instead?," which could be individually decided for each pattern.