Approximate String Matching using a Bidirectional
Index
Gregory Kucherov∗

Kamil Salikhov∗†

Dekel Tsur‡

arXiv:1310.1440v3 [cs.DS] 6 Sep 2015

Abstract
We study strategies of approximate pattern matching that exploit bidirectional text indexes, extending and generalizing ideas of [9]. We introduce a
formalism, called search schemes, to specify search strategies of this type, then
develop a probabilistic measure for the efficiency of a search scheme, prove
several combinatorial results on efficient search schemes, and finally, provide
experimental computations supporting the superiority of our strategies.

1

Introduction

Approximate string matching has numerous practical applications and has long been
a subject of extensive studies by algorithmic researchers [18]. If errors are allowed
in a match between a pattern string and a text string, most of fundamental ideas
behind exact string search algorithms become inapplicable.
The problem of approximate string matching comes in different variants. In this
paper, we are concerned with the indexed variant, when a static text is available for
pre-processing and storing in a data structure (index), before any matching query is
made. The challenge of indexed approximate matching is to construct a small-size
index supporting quick search for approximate pattern occurrences, within a worstcase time weakly dependent on the text length. From the theoretical perspective,
even the case of one allowed error turned out to be highly nontrivial and gave rise
to a series of works (see [10] and references therein). In the case of k errors, existing
solutions generally have time or space complexity that is exponential in k, see [22]
for a survey.
The quest for efficient approximate string matching algorithms has been boosted
by a new generation of DNA sequencing technologies, capable to produce huge
quantities of short DNA sequences, called reads. Then, an important task is to
map those reads to a given reference genomic sequence, which requires very fast
and accurate approximate string matching algorithms. This motivation resulted in
a very large number of read mapping algorithms and associated software programs,
we refer to [13] for a survey.
∗

CNRS/LIGM, Université Paris-Est Marne-la-Vallée, France
Mechanics and Mathematics Department, Lomonosov Moscow State University, Russia
‡
Department of Computer Science, Ben-Gurion University of the Negev, Israel
†

1

Broadly speaking, read mapping algorithms follow one of two main approaches,
or sometimes a combination of those. The filtration approach proceeds in two steps:
it first identifies (with or without using a full-text index) locations of the text
where the pattern can potentially occur, and then verifies these locations for actual matches. Different filtration schemes have been proposed [5, 7, 8, 17]. Filtration
algorithms usually don’t offer interesting worst-case time and space bounds but are
often efficient on average and are widely used in practice. Another approach, usually
called backtracking, extends exact matching algorithms to the approximate case by
some enumeration of possible errors and by simulating exact search of all possible
variants of the pattern. It is this approach that we follow in the present work.
Backtracking and filtration techniques can be combined in a hybrid approach [15].
Some approximate matching algorithms use standard text indexes, such as suffix
tree or suffix arrays. However, for large datasets occurring in modern applications,
these indexes are known to take too much memory. Suffix arrays and suffix trees
typically require at least 4 or 10 bytes per character respectively. The last years saw
the development of succinct or compressed full-text indexes that occupy virtually
as much memory as the sequence itself and yet provide very powerful functionalities [16]. For example, the FM-index [6], based on the Burrows-Wheeler Transform [3], may occupy 2–4 bits of memory per character for DNA texts. FM-index has
now been used in many practical bioinformatics software programs, e.g. [11, 12, 21].
Even if succinct indexes are primarily designed for exact string search, using them for
approximate matching naturally became an attractive opportunity. This direction
has been taken in several papers, see [19], as well as in practical implementations [21].
Interestingly, succinct indexes can provide even more functionalities than classical ones. In particular, succinct indexes can be made bidirectional, i.e. can perform
pattern search in both directions [2, 9, 19, 20]. Lam et al. [9] showed how a bidirectional FM-index can be used to efficiently search for strings up to a small number
(one or two) errors. The idea is to partition the pattern into k + 1 equal parts,
where k is the number of errors, and then perform multiple searches on the FMindex, where each search assumes a different distribution of mismatches among the
pattern parts. It has been shown experimentally in [9] that this improvement leads
to a faster search compared to the best existing read alignment software. Related
algorithmic ideas appear also in [19].
In this paper, we extend the search strategy of [9] in two main directions. We
consider the case of arbitrary k and propose to partition the pattern into more than
k + 1 parts that can be of unequal size. To demonstrate the benefit of both ideas,
we first introduce a general formal framework for this kind of algorithm, called
search scheme, that allows us to easily specify them and to reason about them
(Section 2). Then, in Section 3 we perform a probabilistic analysis that provides
us with a quantitative measure of performance of a search scheme, and give an
efficient algorithm for obtaining the optimal pattern partition for a given scheme.
Furthermore, we prove several combinatorial results on the design of efficient search
schemes (Section 4). Finally, Section 5 contains comparative analytical estimations,
based on our probabilistic analysis, that demonstrate the superiority of our search
strategies for many practical parameter ranges. We further report on large-scale
experiments on genomic data supporting this analysis.
2

2

Bidirectional search

In the framework of text indexing, pattern search is usually done by scanning the
pattern online and recomputing index points referring to the occurrences of the
scanned part of the pattern. With classical text indexes, such as suffix trees or
suffix arrays, the pattern is scanned left-to-right (forward search). However, some
compact indexes such as FM-index provide a search algorithm that scans the pattern
right-to-left (backward search).
Consider now approximate string matching. For ease of presentation, we present
most of our ideas for the case of Hamming distance (recall that the Hamming distance between two strings A and B of equal lengths is the number of indices i for
which A[i] 6= B[i]), although our algorithms extend to the edit distance as well.
Section 3.1.2 below will specifically deal with the edit distance.
Assume that k letter mismatches are allowed between a pattern P and a substring of length |P | of a text T . Both forward and backward search can be extended
to approximate search in a straightforward way, by exploring all possible mismatches
along the search, as long as their number does not exceed k and the current pattern still occurs in the text. For the forward search, for example, the algorithm
enumerates all substrings of T with Hamming distance at most k to a prefix of P .
Starting with the empty string, the enumeration is done by extending the current
string with the corresponding letter of P , and with all other letters provided that
the number of accumulated mismatches has not yet reached k. For each extension,
its positions in T are computed using the index. Note that the set of enumerated
strings is closed under prefixes and therefore can be represented by the nodes of a
trie. Similar to forward search, backward search enumerates all substrings of T with
Hamming distance at most k to a suffix of P .
Clearly, backward and forward search are symmetric and, once we have an implementation of one, the other can be implemented similarly by constructing the
index for the reversed text. However, combining both forward and backward search
within one algorithm results in a more efficient search. To illustrate this, consider
the case k = 1. Partition P into two equal length parts P = P1 P2 . The idea is to
perform two complementary searches: forward search for occurrences of P with a
mismatch in P2 and backward search for occurrences with a mismatch in P1 . In both
searches, branching is performed only after |P |/2 characters are matched. Then, the
number of strings enumerated by the two searches is much less than the number of
strings enumerated by a single standard forward search, even though two searches
are performed instead of one.
A bidirectional index of a text allows one to extend the current string A both left
and right, that is, compute the positions of either cA or Ac from the positions of A.
Note that a bidirectional index allows forward and backward searches to alternate,
which will be crucial for our purposes. Lam et al. [9] showed how the FM-index can
be made bidirectional. Other succinct bidirectional indexes were given in [2, 19, 20].
Using a bidirectional index, such as FM-index, forward and backward searches can
be performed in time linear in the number of enumerated strings. Therefore, our
main goal is to organize the search so that the number of enumerated strings is
minimized.

3

Lam et al. [9] gave a new search algorithm, called bidirectional search, that
utilizes the bidirectional property of the index. Consider the case k = 2, studied
in [9]. In this case, the pattern is partitioned into three equal length parts, P =
P1 P2 P3 . There are now 6 cases to consider according to the placement of mismatches
within the parts: 011 (i.e. one mismatch in P2 and one mismatch in P3 ), 101,
110, 002, 020, and 200. The algorithm of Lam et al. [9] performs three searches
(illustrated in Figure 1):
1. A forward search that allows no mismatches when processing characters of P1 ,
and 0 to 2 accumulated mismatches when processing characters of P2 and P3 .
This search handles the cases 011, 002, and 020 above.
2. A backward search that allows no mismatches when processing characters of
P3 , 0 to 1 accumulated mismatches when processing characters of P2 , and 0
to 2 accumulated mismatches when processing characters of P1 . This search
handles the cases 110 and 200 above.
3. The remaining case is 101. This case is handled using a bidirectional search.
It starts with a forward search on string P 0 = P2 P3 that allows no mismatches
when processing characters of P2 , and 0 to 1 accumulated mismatches when
processing the characters of P3 . For each string A of length |P 0 | enumerated by
the forward search whose Hamming distance from P 0 is exactly 1, a backward
search for P1 is performed by extending A to the left, allowing one additional
mismatch. In other words, the search allows 1 to 2 accumulated mismatches
when processing the characters of P1 .
We now give a formal definition for the above. Suppose that the pattern P is
partitioned into p parts. A search is a triplet of strings S = (π, L, U ) where π
is a permutation string of length p over {1, . . . , p}, and L, U are strings of length
p over {0, . . . , k}. The string π indicates the order in which the parts of P are
processed, and thus it must satisfy the following connectivity property: For every
i > 1, π(i) is either (minj<i π(j)) − 1 or (maxj<i π(j)) + 1. The strings U and L
give upper and lower bounds on the number of mismatches: When the j-th part is
processed, the number of accumulated mismatches between the active strings and
the corresponding substring of P must be P
between L[j] and U [j]. Formally, for a
string A over integers, the weight of A is i A[i]. A search S = (π, L, U ) covers
Pi
a string A if L[i + 1] ≤
j=1 A[j] ≤ U [i] for all i (assuming L[p + 1] = 0). A
k-mismatch search scheme S is a collection of searches such that for every string
A of weight k, there is a search in S that covers A. For example, the 2-mismatch
scheme of Lam et al. consists of searches Sf = (123, 000, 022), Sb = (321, 000, 012),
and Sbd = (231, 001, 012). We denote this scheme by SLam .
In this work, we introduce two types of improvements over the search scheme of
Lam et al.
Uneven partition. In SLam , search Sf enumerates more strings than the other
two searches, as it allows 2 mismatches on the second processed part of P , while
the other two searches allow only one mismatch. If we increase the length of P1 in
the partition of P , the number of strings enumerated by Sf will decrease, while the
4

(a) Forward search

(b) Backward search

(c) Bidirectional search

Figure 1: The tries representing the searches of Lam et al. for binary alphabet {a, b},
search string P = abbaaa, and number of errors 2. Each trie represents one search
and assumes that all the enumerated substrings exist in the text T . In an actual
search on a specific T , each trie contains of a subset of the nodes, depending on
whether the strings of the nodes in the trie appear in T . A vertical edge represents
a match, and a diagonal edge represents a mismatch.

5

number of strings enumerated by the two other searches will increase. We show that
for some typical parameters of the problem, the decrease in the former number is
larger than the increase of the latter number, leading to a more efficient search.
More parts. Another improvement can be achieved using partitions with k + 2 or
more parts, rather than k + 1 parts. We explain in Section 3.2 why such partitions
can reduce the number of enumerated strings.

3

Analysis of search schemes

In this section we show how to estimate the performance of a given search scheme S.
Using this technique, we first explain why an uneven partition can lead to a better
performance, and then present a dynamic programming algorithm for designing an
optimal partition of a pattern.

3.1

Estimating the efficiency of a search scheme

To measure the efficiency of a search scheme, we estimate the number of strings
enumerated by all the searches of S. We assume that performing single steps of
forward, backward, or bidirectional searches takes the same amount of time. It is
fairly straightforward to extend the method of this section to the case when these
times are not equal. Note that the bidirectional index of Lam et al. [9] reportedly
spends slightly more time (order of 10%) on forward search than on backward search.
For the analysis, we assume that characters of T and P are randomly drawn
uniformly and independently from the alphabet. We note that it is possible to
extend the method of this section to a non-uniform distribution. For more complex
distributions, a Monte Carlo simulation can be applied which, however, requires
much more time than the method of this section.
3.1.1

Hamming distance

Our approach to the analysis is as follows. Consider a fixed search S, and the trie
representing this search (see Figure 1). The search enumerates the largest number
of strings when the text contains all strings of length m as substrings. In this
case, every string that occurs in the trie is enumerated. For other texts, the set of
enumerated strings is a subset of the set of strings that occurs in trie. The expected
number of strings enumerated by S on random T and P is equal to the sum over
all nodes v of the trie of the probability that the corresponding string appears in T .
We will first show that this probability depends only on the depth of v (Lemmas 1
and 2 below). Then, we will show how to count the number of nodes in each level
of the trie.
Let prob n,l,σ denote the probability that a random string of length l is a substring
of a random string of length n, where the characters of both strings are randomly
chosen uniformly and independently from an alphabet of size σ. The following
lemma gives an approximation for prob n,l,σ with a bound on the approximation
error.
6

Lemma 1. |prob n,l,σ − (1 − e

−n/σ l

(
4nl/σ 2l
)| ≤
4l/σ l

if l ≥ logσ n
.
otherwise

Proof. Let A and B be random strings of length l and n, respectively. Let Ei be
the event that A appears in B at position i. The event Ei is independent of the
events {Ej : j ∈ {1, 2, . . . , n−l +1}\Fi }, where Fi = {i−l +1, i−l +2, . . . , i+l −1}.
By the Chen-Stein method [1, 4],
prob n,l,σ − (1 − e

−n/σ l

n−l+1
1 − e−λ X X
(Pr[Ei ] Pr[Ej ] + Pr[Ei ∩ Ej ]),
) ≤
λ
i=1 j∈F
i

where λ = n/σ l . Clearly, Pr[Ei ] = Pr[Ej ] = 1/σ l . It is also easy to verify that
l
Pr[Ei ∩ Ej ] = 1/σ 2l . Therefore, |prob n,l,σ − (1 − e−n/σ )| ≤ ((1 − e−λ )/λ) · 4nl/σ 2l .
The lemma follows since (1 − e−λ )/λ ≤ min(1, 1/λ) for all λ.
The bound in Lemma 1 on the error of the approximation of prob n,l,σ is large if l
is small, say l < 12 logσ n. In this case, we can get a better bound by observing that
l
prob n,l,σ ≥ prob n,l0 ,σ , where l0 = 43 logσ n. Since prob n,l0 ,σ ≥ 1 − e−n/σ 0 − 4l0 /σ l0 , we
l
l
l
obtain that |prob n,l,σ − (1 − e−n/σ )| ≤ max(e−n/σ , e−n/σ 0 + 4l0 /σ l0 ).
Let strings(S, X, σ, n) denote the expected number of strings enumerated when
performing a search S = (π, L, U ) on a random text of length n and random pattern
of length m, where X is a partition of the pattern and σ is the alphabet size (note
that m is not a parameter for strings since
P the value of m is implied from X). For
a search scheme S, strings(S, X, σ, n) = S∈S strings(S, X, σ, n).
Fix S, X, σ, and n. Let Al be the set of enumerated strings of length l when
performing search S on a random pattern of length m, partitioned by X, and a
text T̂ containing all strings of length at most m as substrings. Let Al,i be the i-th
element of Al (an order on Al will be defined in the proof of the next lemma). Let
nodes l = |Al |, namely, the number of nodes at depth l in the trie that represents
the search S. Let P ∗ be the string containing the characters of P according to the
order they are read by the search. In other words, P ∗ [l] is the character such that
every node at depth l − 1 of the trie has an edge to a child with label P ∗ [l].
Lemma 2. For every l and i, the string Al,i is a random string with uniform distribution.
Proof. Assume that the alphabet is Σ = {0, . . . , σ − 1}. Consider the trie that
represents the search S. We define an order on the children of each node of the
trie as follows: Let v be a node in the trie with depth l − 1. The label on the
edge between v and its leftmost child is P ∗ [l]. If v has more than one child, the
labels on the edges to the rest of the children of v, from left to right, are (P ∗ [l] +
1) mod σ, . . . , (P ∗ [l] + σ − 1) mod σ. We now order the set Al according to the
nodes of depth l in the trie. Namely, let v1 , . . . , vnodes l be the nodes of depth l in
the trie, from left to right. Then, Al,i is the string that corresponds to vi . We have
that Al,i [j] = (P ∗ [j] + ci,j − 1) mod σ for j = 1, . . . , l, where ci,j is the rank of the
node of depth j on the path from the root to vi among its siblings. Now, since each
letter of P is randomly chosen uniformly and independently from the alphabet, it
7

follows that each letter of Al,i has uniform distribution and the letters of Al,i are
independent.
By the linearity of the expectation,
strings(S, X, σ, n) =

X nodes
Xl
l≥1

i=1

Pr [Al,i is a substring of T ].

T ∈Σn

By Lemma 2 and Lemma 1,
strings(S, X, σ, n) =

m
X

nodes l · prob n,l,σ ≈

l=1

m
X

l

nodes l (1 − e−n/σ ).

(1)

l=1

We note that the bounds on the approximation errors of prob n,l,σ are small, therefore
even when these bounds are multiplied by nodes l and summed over all l, the resulting
bound on the error is small.
In order to compute the values of nodes l , we give some definitions. Let nodes l,d
be the number of strings in Al of length l with Hamming distance d to the prefix
of P ∗ of length l. For example, consider search Sbd = (231, 001, 012) and partition
of a pattern of length 6 into 3 parts of length 2, as shown in Figure 1(c). Then,
P ∗ = baaaba, nodes 5,0 = 0, nodes 5,1 = 2 (strings baabb and babab), and nodes 5,2 =
2 (strings baaba and babaa).
Let πX be a string obtained from π by replacing each character π(i) of π by a run
of π(i) of length X[π(i)], where X[j] is the length of the j-th part in the partition
X. Similarly, LX is a string obtained from L by replacing each character L[i] by a
run of L[i] of length X[π(i)], and UX is defined analogously. In other words, values
LX [i], UX [i] give lower and upper bounds on the number of allowed mismatches for
an enumerated string of length i. For example, for Sbd and the partition X defined
above, πX = 223311, LX = 000011, and UX = 001122.
Values nodes l are given by the following recurrence.
UX [l]

nodes l =

X

nodes l,d

(2)

d=LX [l]

nodes l,d



nodes l−1,d + (σ − 1) · nodes l−1,d−1
= 1


0

if l ≥ 1 and LX [l] ≤ d ≤ UX [l]
if l = 0 and d = 0
otherwise
(3)

For a specific search, a closed formula can be given for nodes l . If a search scheme
S contains two or more searches with the same π-strings, these searches can be
merged in order to eliminate the enumeration of the same string twice or more. It is
straightforward to modify the computation of strings(S, X, σ, n) to account for this
optimization.
l
Consider equation (1). The value of the term 1 − e−n/σ is very close to 1 for l ≤
logσ n − O(1). When l ≥ logσ n, the value of this term decreases exponentially. Note
that nodes l increases exponentially, but the base of the exponent of nodes l is σ − 1
8

l

whereas the base of 1 − e−n/σ is 1/σ. We can then approximate strings(S, X, σ, n)
with function strings 0 (S, X, σ, n) defined by
dlogσ ne+cσ
0

strings (S, X, σ, n) =

X

l

nodes l · (1 − e−n/σ ),

(4)

l=1

where cσ is a constant chosen so that ((σ − 1)/σ)cσ is sufficiently small.
From the above formulas we have that the time complexities for computing
strings(S, X, σ, n) and strings 0 (S, X, σ, n) are O(|S|km) and O(|S|k logσ n), respectively.
3.1.2

Edit distance

We now show how to estimate the efficiency of a search scheme for the edit distance.
We define strings edit analogously to strings in the previous section, except that
edit distance errors are allowed. Fix a search S = (π, L, U ) and a partition X. We
assume without loss of generality that π is the identity permutation. Similarly to
the Hamming distance case, define Al to be the set of enumerated strings of length
l when performing the search S on a random pattern of length m, partitioned by X,
and a text T̂ containing all the strings of length at most m + k as substrings. Unlike
the case of Hamming distance, here the strings of Al are not P
distributed uniformly.
Thus, we doPnot have the equality strings edit (S, X, σ, n) = m
l=1 nodes l · prob n,l,σ .
m
We will use l=1 nodes l · prob n,l,σ as an approximation for strings edit (S, X, σ, n), but
we do not have an estimation on the error of this approximation. Note that in the
Hamming distance case, the sizes of the sets Al are the same for every choice of the
pattern, whereas this is not true for edit distance. We therefore define nodes l (P ) to
be the number of enumerated strings of length l when performing the search S on
a pattern P of length m, partitioned by X, and a text T̂ . We also define nodes l to
be the expectation of nodes l (P ), where P is chosen randomly.
We next show how to compute values nodes l . We begin by giving an algorithm
for computing nodes l (P ) for some fixed P . Build a non-deterministic automaton AP
that recognizes the set of strings that are within edit distance at most k to P , and
the locations of the errors satisfy the requirements of the search [7, 14] (see Figure 2
for an example). For a state q and a string B, denote by δ̂P (q, B) the set of all
states q 0 for which there is a path in AP from q to q 0 such that the concatenation
of the labels on the path is equal to B. For a set of states Q and a string B,
δ̂P (Q, B) = ∪q∈Q δ̂P (q, B). Clearly, nodes l (P ) is equal to the number of strings B of
length l for which δ̂P (q0 , B) 6= ∅, where q0 is the initial state. Let nodes l,Q (P ) be the
number of strings B of length l for which δ̂P (q0 , B) = Q. The values of nodes l,Q (P )
can be computed using dynamic programming and the following recurrence.
X
X
nodes l,Q (P ) =
nodes l−1,Q0 (P ).
c∈Σ Q0 :δ̂P (Q0 ,c)=Q

The values nodes l,Q (P ) gives the values of nodes l (P ), since by definition,
X
nodes l (P ) =
nodes l,Q (P ),
Q

9

Figure 2: Non-deterministic automaton corresponding to the search S = (12, 00, 02)
and pattern P = bbabab over the alphabet Σ = {a, b}. A path from the initial
state q0 to the state in the i-th row and j-column of the automaton correspond to
a string with edit distance i − 1 to P [1..j − 1]. The nodes of the set Q4 are marked
by gray.
where the summation is done over all non-empty sets of states Q.
Note that for a string B of length l, set δ̂P (q0 , B) is a subset of a set of (k + 1)2
states that depends on l. This set, denoted Ql , includes the l + 1-th state in the first
row of the automaton, states l, l+1, l+2 on the second row, states l−1, l, . . . , l+3 on
the third row, and so on (see Figure 2). The size of Ql is 1 + 3 + 5 + · · · + (2k + 1) =
(k + 1)2 . Therefore, the number of sets Q for which nodes l,Q (P ) > 0 is at most
2
2(k+1) . If (k + 1)2 is small enough, a state can be encoded in one machine word,
and the computation of δ̂P (Q0 , c) can be done in constant time using precomputed
2
tables. Thus, the time for computing all values of nodes l,Q (P ) is O(2k σm).
Now consider the problem of computing the values of nodes l . Observe that for
Q ⊆ Ql , the value of δ̂P (Q, c) depends on the characters of P [l − k + 1..l + k + 1],
and does not depend on the rest of the characters of P . Our algorithm is based on
this observation. For an integer l, a set Q ⊆ Ql , and a string P 0 of length 2k + 1,
define
X
nodes l,Q,P 0 =
nodes l,Q (P ).
P :P [l−k+1..l+k+1]=P 0

Then,
nodes l,Q,P 0 =

XX
c0 ∈Σ

X

nodes l−1,Q0 ,Pc0 ,

c∈Σ Q0 :δ̂P (Q0 ,c)=Q
c

where Pc0 = c0 P 0 [1..2k], and Pc is a string satisfying Pc [(l−1)−k+1..(l−1)+k+1] = Pc0
(the rest of the characters of Pc can be chosen arbitrarily).
From the above, the time complexity for computing strings edit (S, X, σ, n) is
2
O(|S|2k σ 2k+3 m). Therefore, our approach is practical only for small values of k.

3.2

Uneven partitions

In Section 2, we provided an informal explanation why partitioning the pattern into
unequal parts may be beneficial. We now provide a formal justification for this. To

10

this end, we replace (4) by an even simpler estimator of strings(S, X, σ, n):
dlogσ ne
00

strings (S, X, σ, n) =

X

nodes l .

(5)

l=1

As an example, consider scheme SLam . Denote by x1 , x2 , x3 the lengths of the
parts in a partition X of P into 3 parts. It is straightforward to give closed formulas
for strings 00 (S, X, σ, n) for each search of SLam . For example,
(
N
if N ≤ x1
strings 00 (Sf , X, σ, n) =
3
2
c1 (N − x1 ) + c2 (N − x1 ) + c3 (N − x1 ) + N otherwise
where N = dlogσ ne, c1 = (σ−1)2 /6, c2 = (σ−1)/2, and c3 = −(σ−1)2 /6+(σ−1)/2.
Similar formulas can be given for Sb and Sbd . If x1 , x2 , and x3 are close to m/3
and N < m/3 then strings 00 (SLam , X, σ, n, m) = 3N and an equal sized partition is
optimal in this case. However, if m/3 < N < 2m/3, then
strings 00 (SLam ,X, σ, n) = c1 (N − x1 )3 + c2 (N − x1 )2 + c3 (N − x1 )
+ c01 (N − x3 )2 + c02 (N − x3 ) + c001 (N − x2 )2 + c002 (N − x2 ) + 3N.
It is now clear why the equal sized partition is not optimal in this case. The degree
of N − x1 in the above polynomial is 3, while the degrees of N − x2 and N − x3 are
2. Thus, if x1 = x2 = x3 = m/3, decreasing x2 and x3 by, say 1, while increasing x1
by 2 reduces the value of the polynomial.

3.3

Computing an optimal partition

In this Section, we show how to find an optimal partition for a given search scheme
S and a given number
 of parts p. An optimal partition can be naively found by
m−1
enumerating all p−1 possible partitions, and for each partition X, computing
strings 0 (S, X, σ, n). We now describe a more efficient dynamic programming algorithm.
We define an optimal partition to be a partition that maximizes strings(S, X, σ, n).
Let N = dlogσ ne + cσ . If m ≥ pN , then any partition in which all parts are of size
at least N is an optimal partition. Therefore, assume for the rest of this section
that m < pN . We say that a partition X is bounded if the sizes of the parts of X
are at most N . If X is not bounded, we can transform it into a bounded partition
by decreasing the sizes of parts which are larger than N and increasing the sizes of
parts which are smaller that N . This transformation can only decrease the value of
strings(S, X, σ, n). Therefore, there exists an optimal partition which is bounded.
Throughout this section we will consider only bounded partitions. For brevity, we
will use the term partition instead of bounded partition.
Our algorithm takes advantage of the fact that the value of strings 0 (S, X, σ, n)
does not depend on the entire partition X, but only on the partition of a substring
of P of length N induced by X. More precisely, consider a fixed S = (π, L, U ) ∈ S.
By definition, strings 0 (S, X, σ, n) depends on the values nodes 1 , . . . , nodes N (the
number of nodes in levels 1, . . . , N in the trie that correspond to the search S). From
11

Section 3.1, these values depend on the strings L and U which are fixed, and on the
string πX [1..N ]. The latter string depends on π[1..iX,π ], where iX,π is the minimum
PiX,π
index such that j=1
X[π(j)] ≥ N and on the values X[π(1)], . . . , X[π(iX,π )].
The algorithm works by going over the prefixes of P in increasing length order.
For each prefix P 0 , it computes a set of partitions of P 0 such that at least one partition
in this set can be extended to an optimal partition of P . In order to reduce the
time complexity, the algorithm needs to identify partitions of P 0 that cannot be
extended into an optimal partition of P . Consider the following example. Suppose
that m = 13, p = 5, N = 4 and S = {S1 , S2 , S3 }, where the π-strings of S1 , S2 , S3
are π 1 = 12345, π 2 = 32451, and π 3 = 43215, respectively. Consider a prefix
P 0 = P [1..8] of P , and let Y1 , Y2 be two partitions of P 0 , where the parts in Y1 are of
sizes 3,3,2, and the parts in Y1 are of sizes 4,2,2. Note that Y1 and Y2 have the same
number of parts, and they induce the same partition on P [8 − N + 1..8] = P [5..8].
We claim that one of these two partitions is always at least as good as the other for
every extension of both partitions to a partition of P . To see this, let Z denote a
partition of P [9..13] into two parts, and consider the three searches of S.
1. For search S1 we have that πY11 ∪Z [1..N ] = 1112 for every Z, and πY12 ∪Z [1..N ] =
1111 for every Z. It follows that the value of strings 0 (S1 , Y1 ∪ Z, σ, n) is the
same for every Z, and the value of strings 0 (S1 , Y2 ∪ Z, σ, n) is the same for
every Z. These two values can be equal or different.
2. For the search S2 we have that πY21 ∪Z [1..N ] = πY22 ∪Z [1..N ] = 3322. It follows
that strings 0 (S2 , Y1 ∪ Z, σ, n) = strings 0 (S2 , Y2 ∪ Z, σ, n) for all Z and this
common value does not depend on Z.
3. For the search S3 we have that πY31 ∪Z [1..N ] = πY32 ∪Z [1..N ] for every Z. For
example, if Z is a partition of P [9..13] into parts of sizes 2,2 then πY31 ∪Z [1..N ] =
πY32 ∪Z [1..N ] = 4433. It follows that strings 0 (S3 , Y1 ∪ Z, σ, n) = strings 0 (S3 , Y2 ∪
Z, σ, n) for every Z. This common value depends on Z.
We conclude that either strings 0 (S, Y1 ∪ Z, σ, n) < strings 0 (S, Y2 ∪ Z, σ, n) for every
Z, or strings 0 (S, Y1 ∪ Z, σ, n) ≥ strings 0 (S, Y2 ∪ Z, σ, n) for every Z.
We now give a formal description of the algorithm. We start with some definitions. For a partition Y of a substring P 0 = P [m00 ..m0 ] of pattern P , we define the
following quantities: mY is the length of P 0 , lY is the length of the last part of Y ,
pY is the number of parts in Y , and rY is the left-to-right rank of the part of Y
containing P 0 [m0 − N + 1]. Let prefix (Y ) be the partition of P [m00 ..m0 − lY ] of P 0
that is composed from the first pY − 1 parts of Y . For the example above, mY1 = 8,
lY1 = 2, pY1 = 3, rY1 = 2, and prefix (Y1 ) is a partition of P [1..6] with parts sizes 3, 3.
For a partition Y of a prefix P 0 of P , S(Y ) is a set containing every search S ∈ S
such that rY appears before pY + 1 in the π-string of S. If the length of P 0 is less
than N we define S(Y ) = ∅, and if P 0 = P we define S(Y ) = S. For the example
above, S(Y1 ) = {S1 , S2 }.
Let Y1 be a partition of a substring P1 = P [i1 ..j1 ] of P , and Y2 be a partition of
a substring P2 = P [i2 ..j2 ]. We say that Y1 and Y2 are compatible if these partitions
induce the same partition on the common substring P 0 = P [max(i1 , i2 ).. min(j1 , j2 )].
12

For example, the partition of P [4..6] into parts of sizes 1, 2 is compatible with the
partition of P [1..6] into parts of sizes 2, 2, 2.
Lemma 3. Let Y be a partition of a prefix of P of length at least N . Let S ∈ S(Y )
be a search. The value strings 0 (S, X, σ, n) is the same for every partition X of P
whose first pY parts match Y .
Proof. Let i0 be the index such that π(i0 ) = pY + 1. Since rY appears before
pY + 1 in string π, from the connectivity property of π we have that (1) Every
value in π that appears before pY + 1 is at most pY . In other words, π(i) ≤ pY
0
for
pY + 1 in π. By the definition of rY ,
PpYevery i < i . (2) rY , . . . , pY appear before
0
X[j]
≥
N
.
Therefore,
i
<
i
and
π(1),
. . . , π(iX,π ) ≤ pY . Thus, string
X,π
j=rY
π[1..iX,π ] and values X[π(1)], . . . , X[π(iX,π )] are the same for every partition X that
satisfies the requirement of the lemma.
Y of a prefix of P of length at least N , define v(Y ) to be
P For a partition
0
S∈S(Y ) strings (S, X, σ, n), where X is an arbitrary partition of P whose first pY
parts match Y (the choice of X does not matter due to Lemma 3). For a partition
Y of a prefix of P of length less than N , v(Y ) = 0. Define
X
∆(Y ) = v(Y ) − v(prefix (Y )) =
strings 0 (S, X, σ, n).
S∈S(Y )\S(prefix (Y ))

Lemma 4. Let Z be a partition of a substring P [m00 ..m0 ] such that pZ ≥ 2 and
mprefix (Z) = min(N, m0 − lY ). Let p0 ≥ pZ be an integer. The value of ∆(Y ) is the
same for every partition Y of P [1..m0 ] with p0 parts that is compatible with Z.
Proof. We assume N < m0 −lY (the case N ≥ m0 −lY is similar). Since mprefix (Z) =
min(N, m0 − lY ), the set S(Y ) \ S(prefix (Y )) is the same for every partition Y of
P [1..m0 ] with p0 parts that is compatible with Z. For a search S = (π, L, U ) in this
set, rY appears before pY + 1 in π, and pY appears before rprefix (Y ) . Let i = iX,π ,
where X is an arbitrary partition of P whose first pY parts are the parts of Y . We
obtain that rprefix (Y ) ≤ π(1), . . . , π(i) ≤ pY , and the lemma follows.
For Z, p0 that satisfy the requirements of Lemma 4, let ∆(Z, p0 ) denote the value
of ∆(Y ), where Y is an arbitrary partition of P [1..m0 ] with p0 parts that is compatible
with Z.
For m0 ≤ m, p0 ≤ p, and a partition Z of P [max(m0 − N + 1, 1)..m0 ] with at
most p0 parts, let v(m0 , p0 , Z) be the minimum value of v(Y ), where Y is a partition
of P [1..m0 ] into p0 parts that is compatible with Z.
Lemma 5. For m0 ≤ m, 2 ≤ p0 ≤ p, and a partition Z of P [max(m0 − N + 1, 1)..m0 ]
with at most p0 parts,
v(m0 , p0 , Z) = min
(v(m0 − lZ 0 , p0 − 1, prefix (Z 0 )) + ∆(Z 0 , p0 ))
0
Z

where the minimum is taken over all partitions Z 0 of a substring P [m00 ..m0 ] of P that
satisfy the following: (1) Z 0 is compatible with Z, (2) 2 ≤ pZ 0 ≤ p0 , (3) mprefix (Z 0 ) =
min(N, m0 − lZ 0 ), (4) pZ = p0 if m00 = 1.
13

An algorithm for computing the optimal partition
follows from Lemma
 5. The
Pmin(p−1,N
)
−1
time complexity of the algorithm is O (|S|kN + m) j=1
(p−j) Nj−1
, where

Pmin(p−1,N )
P
min(p−1,N )
−1
|S|kN j=1
(p−j) Nj−1
is time for computing ∆ values, and O m j=1
(p−

N −1
j) j−1 is time for computing v values.

4

Properties of optimal search schemes

Designing an efficient search scheme for a given set of parameters consists of (1) choosing a number of parts, (2) choosing searches, (3) choosing a partition of the pattern.
While it is possible to enumerate all possible choices, and evaluate the efficiency of
the resulting scheme using Section 3.1, this is generally infeasible due to a large number of possibilities. It is therefore desirable to have a combinatorial characterization
of optimal search schemes.
The critical string of a search scheme S is the lexicographically maximal U -string
of a search in S. A search of S is critical if its U -string is equal to the critical string
of S. For example, the critical string of SLam is 022, and Sf is the critical search. For
typical parameters, critical searches of a search scheme constitute the bottleneck.
Consider a search scheme S, and assume that the L-strings of all searches contain
only zeros. Assume further that the pattern is partitioned into equal-size parts. Let
` be the maximum index such that for every search S ∈ S and every i ≤ `, U [i] of S
is no larger than the number in position i in the critical string of S. From Section 3,
the number of strings enumerated by a search S ∈ S depends mostly on the prefix
of the U -string of S of length ddlogσ ne/(m/p)e. Thus, if ddlogσ ne/(m/p)e ≤ `, a
critical search enumerates an equal or greater number of strings than a non-critical
search.
We now consider the problem of designing a search scheme whose critical string
is minimal. Let α(k, p) denote the lexicographically minimal critical string of a kmismatch search scheme that partitions the pattern into p parts. The next theorems
give the values of α(k, k + 2) and α(k, k + 1). We need the following definition. A
string over the alphabet of integers is called simple if it contains a substring of the
form 01j 0 for j ≥ 0.
Lemma 6.

(i) Every string A of weight k and length at least k + 2 is simple.

(ii) If A is a non-simple string of weight k and length k+1 then A[1] ≤ 1, A[k+1] ≤
1, and A[i] ≤ 2 for all 2 ≤ i ≤ k. Moreover, there are no two consecutive 2’s
in A.
Proof. (i) The proof is by induction on k. It is easy to verify that the lemma holds
for k = 0. Suppose we proved the lemma for k 0 < k. Let A be a string of weight k
and length p ≥ k + 2. If A[1] ≥ 1 then by the induction hypothesis A[2..p] is simple,
and therefore A is simple. Suppose that A[1] = 0. Let i > 1 be the minimum index
such that A[i] 6= 1 (i must exist due to the assumption that p ≥ k + 2). If A[i] = 0
then we are done. Otherwise, we can use the induction hypothesis on A[i + 1..p] and
obtain that A is simple.
(ii) Let A be a non-simple string of weight k and length k + 1. If A[1] ≥ 2 then
0
A = A[2..k + 1] has weight k − A[1] ≤ k − 2 and length k, and thus by (i) we obtain
14

that A0 is simple, contradicting the assumption that A is non-simple. Similarly,
A[k + 1] cannot be greater than 1. For 2 ≤ i ≤ k, if A[i] ≥ 3 then either A[1..i − 1]
or A[i + 1..k + 1] satisfies the condition of (i). Similarly, if A[i] = A[i + 1] = 2 then
either A[1..i − 1] or A[i + 2..k + 1] satisfies the condition of (i).
We use the following notation. For two integers i and j, [i, j] denotes the string
i(i + 1)(i + 2) · · · j if i ≤ j, and the empty string if i > j. Moreover, [i, j] denotes
the string i(i − 1)(i − 2) · · · j if i ≥ j, and the empty string if i < j.
Theorem 7. α(k, k + 1) = 013355 · · · kk for every odd k, and α(k, k + 1) =
02244 · · · kk for every even k.
Proof. We first give an upper bound on α(k, k + 1) for odd k. We build a search
scheme as follows. The scheme contains searches Sk,i,j = ([i, k+2][i − 1, 1], 0 · · · 0, [0, j]jk · · · k)
for all i and j, which cover all simple strings of weight k and length k + 1. In order
to cover the non-simple strings, the scheme contains the following searches.
1
= ([i, k + 1][i − 1, 1], 0 · · · 0, 013355 · · · jj(j + 1)k · · · k) for every odd 3 ≤
1. Sk,i,j
j ≤ k (for j = k, the U -string is 013355 · · · kk).
2
2. Sk,i,j
= ([i, 1][i + 1, k + 1], 0 · · · 0, 013355 · · · jj(j + 1)k · · · k) for every odd 3 ≤
j ≤ k (for j = k, the U -string is 013355 · · · kk).

Let A be a non-simple string of weight k and length k + 1. By Lemma 6,
A = X0A1 0A2 0 · · · 0Ad 0Y where each of X and Y is either string 1 or empty string,
and each Ai is either 2, 12, 21, or 121. A string Ai is called a block of type 1, 2, or 3
if Ai is equal to 12, 21, or 121, respectively. Let B1 , . . . , Bd0 be the blocks of type 1
and type 2, from left to right.
We consider several cases. The first case is when X and Y are empty strings,
and B1 is of type 1. Since the weight of A is odd, it follows that d0 is odd. If A has
1
no other blocks, A is covered by search Sk,i,k
, where i + 1 is the index in A in which
1
B1 starts. Otherwise, if B2 is of type 1, then A is covered by search Sk,i,j
, where
i + 1 is the index in A in which B1 starts, and i + j + 1 is the index in which the first
block to the right of B1 starts (this block is either B2 , or a block of type 3). Now
2
suppose that B2 is of type 2. If B3 is of type 2, then A is covered by search Sk,i,j
,
where i − 1 is the index in A in which B3 ends, and i − j − 1 is the index in which
the first block to the left of B3 ends. By repeating these arguments, we obtain that
A is covered unless the types of B1 , . . . , Bd0 alternate between type 1 and type 2.
1
However, since d0 is odd, Bd0 is of type 1, and in this case A is covered by Sk,i,j
,
where i + 1 is the index in A in which B1 starts, and k − j is the index in which the
first block to the left of B1 ends.
Now, if X is empty string and
Y = 1, define a string A0 = A20. By
the above,
j0
j0
j0
0
.
A is covered by some search Sk+2,i,j . Then, A is covered by either Sk,i,j
or Sk,i,j−2
The same argument holds for the case when X = 1. The proof for the case when
B1 is of type 2 is analogous and thus omitted.
The lower bound on α(k, k+1) for odd k is obtained by considering the string A =
012020 · · · 20. The U -string of a search that covers A must be at least 013355 · · · kk.
We next give an upper bound on α(k, k + 1) for even k. We define k-mismatch
search schemes Sk recursively. For k = 0, S0 consists of a single search S0,1 = (1, 0, 0).
For k ≥ 2, Sk consists of the following searches.
15

1. For every search Sk−2,i = (π, 0 · · · 0, U ) in Sk−2 , Sk contains a search Sk,i =
(π · k(k + 1), 0 · · · 0, U · kk).
2. A search Sk,k = ([k + 1, 1], 0 · · · 0, 01kk · · · k).
3. A search Sk,k+1 = (k(k + 1)[k − 1, 1], 0 · · · 0, 01kk · · · k).
Note that the critical string of Sk is 02244 · · · kk corresponding to item 1 above. We
now claim that all number strings of length k + 1 and weight at most k are covered
by the searches of Sk . The proof is by induction on k. The base k = 0 is trivial.
Suppose the claim holds for k − 2. Let A be a number string of length k + 1 and
weight k 0 ≤ k. If A[k] + A[k + 1] ≤ 1, then A is covered by either Sk,k or Sk,k+1 .
Otherwise, the weight of A0 = A[1..k − 1] is at most k 0 − 2. By induction, A0 is
covered by some search Sk−2,i . Then search Sk,i covers A.
To prove that α(k, k + 1) ≥ 02244 · · · kk for even k, consider the string A =
0202 · · · 020. It is easy to verify that the U -string of a search that covers A must be
at least 02244 · · · kk.
Theorem 8. α(k, k + 2) = 0123 · · · (k − 1)kk for every k ≥ 1.
Proof. We first give an upper bound on α(k, k + 1). We build a k-mismatch search
scheme S that contains searches Sk,i,j = ([i, k + 2][i − 1, 1], 0 · · · 0, [0, j]jk · · · k) for
all i and j. Let A be a string of weight k and length k + 2. By Lemma 6 there are
indices i and j such that A[i..i + j + 1] = 01j 0, and therefore A is covered by Sk,i,j .
The lower bound is obtained from the string A = 011 · · · 110. It is easy to verify
that the U -string of a search that covers A must be at least 0123 · · · (k − 1)kk.
An important consequence of Theorems 7 and 8 is that for some typical cases,
partitioning the pattern into k + 2 parts brings an advantage over k + 1 parts. For
k = 2, for example, we have α(2, 3) = 022 while α(2, 4) = 0122. Since the second
element of 0122 is smaller than that of 022, a 4-part search scheme potentially
enumerates less strings than a 3-part scheme. On the other hand, the average
length of a part is smaller when using 4 parts, and therefore the branching occurs
earlier in the searches of a 4-part scheme. The next section shows that for some
parameters, (k + 2)-part schemes outperform (k + 1)-part schemes, while for other
parameters the inverse occurs.

5

Case studies

In this Section, we provide results of several computational experiments we have
performed to analyse practical applicability of our techniques.
We designed search schemes for 2, 3 and 4 errors (given in Appendix) using a
greedy algorithm. The algorithm iteratively adds searches to a search scheme. At
each step, the algorithm considers the uncovered string A of weight k such that the
lexicographically minimal U -string that covers A is maximal. Among the searches
that cover A with minimal U -string, a search that covers the maximum number of
uncovered strings of weight k is chosen. The L-string of the search is chosen to
be lexicographically maximal among all possible L-string that do not decrease the
16

number of uncovered strings. For each search scheme and each choice of parameters,
we computed an optimal partition.

5.1

Numerical comparison of search schemes

We first performed a comparative estimation of the efficiency of search schemes using
the method of Section 3.1.1 (case of Hamming distance). More precisely, for a given
search scheme S, we estimated the number of strings strings(S, X, σ, n) enumerated
during the search.
Results for 2 mismatches are given in Table 1 and Table 2 for 4-letter and 30letter alphabets respectively. Table 3 contains estimations for nonuniform letter
distribution. Table 4 contains estimations for 3 mismatches for 4-letter alphabet.
We first observe that our method provides an advantage only on a limited range
of pattern lengths. This conforms to our analysis (see Section 3.2) that implies
that our schemes can bring an improvement when m/(k + 1) is smaller than logσ n
approximately. When m/(k + 1) is small, Tables 1–4 suggest that using more parts
of unequal size can bring a significant improvement. For big alphabets (Table 2), we
observe a larger gain in efficiency, due to the fact that values nodes l (see equation (2))
grow faster when the alphabet is large, and thus a change in the size of parts can
have a bigger influence on these values. Moreover, if the probability distribution
of letters in both the text and the pattern is nonuniform, then we obtain an even
larger gain (Table 3), since in this case, the strings enumerated during the search
have a larger probability to appear in the text than for the uniform distribution.
For 3 mismatches and 4 letters (Table 4), we observe a smaller gain, and even a
loss for pattern lengths 36 and 48 when shifting from 4 to 5 parts. This is explained
by Theorem 7 showing the difference of critical strings between odd and even numbers of errors. Thus, for 3 mismatches and 4 parts, the critical string is 0133 while
for 5 parts it is 01233. When patterns are not too small, the latter does not lead
to an improvement strong enough to compensate for the decrease of part length.
Note that the situation is different for even number of errors, where incrementing
the number of parts from k + 1 to k + 2 leads to transforming the critical strings
from 0224 · · · to 0123 · · · .
Another interesting observation is that with 4 parts, obtained optimal partitions
have equal-size parts, as the U -strings of all searches of the 4-part scheme are all
the same (see Appendix).
These estimations suggest that our techniques can bring a significant gain in
efficiency for some parameter ranges, however the design of a search scheme should
be done carefully for each specific set of parameters.

5.2

Experiments on genomic data

To perform large-scale experiments on genomic sequences, we implemented our
method using the 2BWT library provided by [9] (http://i.cs.hku.hk/2bwt-tools/).
We then experimentally compared different search schemes, both in terms of running
time and average number of enumerated substrings. Below we only report running
time, as in all cases, the number of enumerated substrings produced very similar

17

Table 1: Values of strings(S, X, 4, 416 ) for 2-mismatch search schemes, for different
pattern lengths m. Second column corresponds to search scheme SLam with three
equal-size parts, the other columns show results for unequal partitions and/or more
parts. The partition used is shown in the second sub-column.
m 3 equal
3 unequal
4 unequal
5 unequal
24
1197
1077
9,7,8
959 7,4,4,9 939
7,1,6,1,9
36
241
165 15,10,11 140 12,5,7,12 165 11,1,9,1,14
48
53
53 16,16,16 51 16,7,9,16 53 16,1,15,1,15
Table 2: Values of strings(S, X, 30, 307 ) for 2-mismatch search schemes.
m 3 equal 3 unequal
4 unequal
5 unequal
15
846
286 6,4,5 231 5,2,3,5 286 5,1,3,1,5
18
112
111 7,6,5 81 6,2,4,6 111 6,1,4,1,6
21
24
24 7,7,7 23 7,3,4,7 24 7,1,6,1,6
Table 3: Values of strings(S, X, 4, 416 ) for 2-mismatch search schemes, using a nonuniform letter distribution (one letter with probability 0.01 and the rest with probability 0.33 each).
m 3 equal
3 unequal
4 unequal
5 unequal
24
3997
3541 10,8,6 3592 6,7,1,10 3541
6,1,7,1,9
36
946
481 16,10,10 450 11,6,6,13 481 10,1,9,1,15
48
203
157 18,15,15 137 16,7,9,16 157 15,1,14,1,17
Table 4: Values of strings(S, X, 4, 416 ) for 3-mismatch search schemes. Best partitions obtained for 4 parts are equal.
m
4 equal/unequal
5 unequal
24 11222
6,6,6,6
8039
4,6,5,1,8
36 416
9,9,9,9
549
6,11,5,1,13
48 185 12,12,12,12 213 11,11,11,1,14
results.
The experiments were done on the sequence of human chromosome 14 (hr14 ).
The sequence is 88 · 106 long, with nucleotide distribution 29%, 21%, 21%, 29%.
Searched patterns were generated as i.i.d. sequences. For every search scheme and
pattern length, we ran 105 pattern searches for Hamming distance and 104 searches
for the edit distance.
5.2.1

Hamming distance

For the case of 2 mismatches, we implemented the 3-part and 4-part schemes (see
Appendix), as well as their equal-size-part versions for comparison. For each pattern length, we computed an optimal partition, taking into account a non-uniform
distribution of nucleotides. Results are presented in Table 5.
Using unequal parts for 3-part schemes yields a notable time decrease for patterns
18

Table 5: Total time (in sec) of search for 105 patterns in hr14, up to 2 mismatches.
2nd column contains time obtained on partition into three equal-size parts. The 3rd
(respectively 4th and 5th) column shows the running time respectively for the 3unequal-parts, 4-equal-parts and 4-unequal-parts searches, together with their ratio
(%) to the corresponding 3-equal-parts value.
m 3 equal
3 unequal
4 equal
4 unequal
15
24.8
25.4 (102%)
6,6,3
25.3 (102%) 25.3 (102%) 3,5,1,6
24
5.5
4.2 (76%)
10,7,7
5.2 (95%)
4.0 (73%)
7,4,4,9
33
1.73
1.45 (84%) 13,10,10 2.07 (120%) 1.25 (72%) 11,5,6,11
42
0.71
0.71 (100%) 14,14,14 1.24 (175%) 0.82 (115%) 14,6,8,14
Table 6: Total time (in sec) of search for 105 patterns in hr14, up to 3 mismatches.
m 4 equal
5 equal
5 unequal
15
241
211 (86%)
206 (85%)
2,3,5,1,4
24
19.7
26.7 (136%) 19.6 (99%)
2,9,3,1,9
33
4.3
6.9 (160%) 4.7 (109%)
6,9,6,1,11
42
1.85
2.52 (136%) 2.05 (111%) 10,10,9,1,12
51
1.07
1.57 (147%) 1.06 (99%) 12,13,12,1,13
Table 7: Total time (in sec) of search for
distance).
m 3 equal
3 unequal
15
11.5
11.4 (99%)
6,6,3
24
2.1
1.3 (62%)
11,5,8
33
0.34
0.22 (65%) 13,10,10
42
0.08
0.08 (100%) 14,14,14

104 patterns in hr14, up to 2 errors (edit
4 equal
10.9 (95%)
1.5 (71%)
0.35 (103%)
0.18 (225%)

4 unequal
11.1 (97%)
3,5,1,6
1.0 (48%)
7,4,4,9
0.19 (56%) 11,5,6,11
0.08 (100%) 14,6,8,14

of length 24 and 33 (respectively, by 24% and 16%). Furthermore, we observe that
using unequal part lengths for 4-part schemes is beneficial as well. For pattern
lengths 24 and 33, we obtain a speed-up by 27% and 28% respectively. Overall, the
experimental results are consistent with numerical estimations of Section 5.1.
For the case of 3 mismatches, we implemented 4-part and 5-part schemes from
Appendix, as well as their equal part versions for comparison. Results (running
time) are presented in Table 6. In accordance with estimations of Section 5.1,
here we observe a clear improvement only for pattern length 15 and not for longer
patterns.
5.2.2

Edit distance

In the case of edit distance, along with the search schemes for 2 and 3 errors from the
previous section, we also implemented search schemes for 4 errors (see Appendix).
Results are shown in Table 7 (2 errors), Table 8 (3 errors) and Table 9 (4 errors).
For 2 errors, we observe up to two-fold speed-up for pattern lengths 15, 24 and
33. For the case of 3 errors, the improvement is achieved for pattern lengths 15

19

Table 8: Total time (in sec) of search for
distance).
m 4 equal
5 equal
15
233
174 (75%)
24
13.5
13.2 (98%)
33
0.74
1.81 (245%)
42
0.28
0.45 (161%)
51
0.13
0.24 (185%)
Table 9: Total
distance).
m 5 equal
15
4212
24
145
33
6.5
42
1.66
51
0.60
60
0.28

104 patterns in hr14, up to 3 errors (edit
5 unequal
168 (72%)
2,2,6,1,4
10.8 (80%)
3,8,3,1,9
1.07 (145%) 5,10,5,1,12
0.37 (132%) 9,10,9,1,13
0.14 (108%) 12,12,12,1,14

time (in sec) of search for 104 patterns in hr14, up to 4 errors (edit
5 unequal
6 equal
3222 (76%)
3,1,8,1,2
4028 (96%)
133 (92%)
7,3,5,1,8
131 (90%)
5.8 (89%)
8,7,5,8,5
6.6 (102%)
1.16 (70%)
12,8,7,8,7
1.51 (91%)
0.49 (82%)
13,11,9,9,9
0.74 (123%)
0.24 (86%) 14,13,11,11,11 0.44 (157%)

6 unequal
3401 (81%)
2,2,1,7,1,2
113 (78%)
2,7,3,4,5,3
5.1 (78%)
4,8,6,3,5,7
1.17 (70%)
7,8,8,5,2,12
0.54 (90%)
9,10,9,9,1,13
0.28 (117%) 11,12,11,11,1,14

Table 10: Total time (in sec) of search for 105 reads in hr14, up to 4 errors. First
row corresponds to read set with constant error rate 0.03. Second row corresponds
to read set with error rate increasing from 0.0 to 0.03.
m 5 equal
6 equal
6 unequal
100
247
250 (101%) 283 (115%) 20,20,20,19,1,20
100
415
367 (88%) 350 (84%) 20,20,20,19,1,20
and 24 (respectively 28% and 20%). Finally, for 4 errors, we obtain a significant
speed-up (18% to 30%) for pattern lengths between 15 and 51.
5.2.3

Experiments on simulated genomic reads

Experiments of Section 5.2 have been made with random patterns. In order to
make experiments closer to the practical bioinformatic setting occurring in mapping genomic reads to their reference sequence, we also experimented with patterns
simulating reads issued from genomic sequencers. For that, we generated realistic
single-end reads of length 100 (typical length of Illumina reads) from hr14 using
dwgsim read simulator (https://github.com/nh13/DWGSIM). Two sets of reads
were generated using two different error rate values (parameter -e of dwgsim):
0.03 for the first dataset and 0.0-0.03 for the second one. This means that in
the first set, error probability is uniform over the read length, while in the second
set, this probability gradually increases from 0 to 0.03 towards the right end of the
read. The latter simulates the real-life situation occurring with current sequencing
technologies including Illumina.
The results are shown in Table 10. As expected, due to a large pattern length, our

20

schemes did not produce a speed-up for the case of constant error rate. Interestingly
however, for the case of non-uniform distribution of errors, our schemes showed a
clear advantage. This illustrates another possible benefit of our techniques: they
are better adapted to a search for patterns with non-uniform distribution of errors,
which often occurs in practical situations such as mapping genomic reads.

6

Conclusions

This paper can be seen as the first step towards an automated design of efficient
search schemes for approximate string matching, based on bidirectional indexes.
More research has to be done in order to allow an automated design of optimal
search schemes. It would be very interesting to study an approach when a search
scheme is designed simultaneously with the partition, rather than independently as
it was done in our work.
We expect that search schemes similar to those studied in this paper can be
applied to hybrid approaches to approximate matching (see Introduction), as well
as possibly to other search strategies.
Acknowledgements. GK has been supported by the ABS2NGS grant of the
French government (program Investissement d’Avenir ) as well as by a EU MarieCurie Intra-European Fellowship for Carrier Development. KS has been supported
by the co-tutelle PhD fellowship of the French government. DT has been supported
by ISF grant 981/11.

References
[1] A. D. Barbour, L. Holst, and S. Janson. Poisson approximation. Clarendon
Press Oxford, 1992.
[2] D. Belazzougui, F. Cunial, J. Kärkkäinen, and V. Mäkinen. Versatile succinct
representations of the bidirectional burrows-wheeler transform. In Proc. 21st
European Symposium on Algorithms (ESA), pages 133–144, 2013.
[3] M. Burrow and D. Wheeler. A block-sorting lossless data compression algorithm. Technical report 124, Digital Equipment Corporation, California, 1994.
[4] L. H. Y. Chen. Poisson approximation for dependent trials. The Annals of
Probability, 3(3):534–545, 1975.
[5] M. Farach-Colton, G. M. Landau, S. C. Sahinalp, and D. Tsur. Optimal spaced
seeds for faster approximate string matching. In Proc. 32nd International Colloquium on Automata, Languages and Programming (ICALP), LNCS 3580, pages
1251–1262, 2005.
[6] P. Ferragina and G. Manzini. Opportunistic data structures with applications.
In Proc. 41st Symposium on Foundation of Computer Science (FOCS), pages
390–398, 2000.
21

[7] J. Kärkkäinen and J. C. Na. Faster filters for approximate string matching. In
Proc. 9th Workshop on Algorithm Engineering and Experiments (ALENEX),
pages 84–90, 2007.
[8] G. Kucherov, L. Noé, and M. Roytberg. Multi-seed lossless filtration.
IEEE/ACM Transactions on Computational Biology and Bioinformatics,
2(1):51–61, January-March 2005.
[9] T. W. Lam, R. Li, A. Tam, S. C. K. Wong, E. Wu, and S.-M. Yiu. High
throughput short read alignment via bi-directional BWT. In Proc. IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pages 31–36,
2009.
[10] T. W. Lam, W. K. Sung, and S. S. Wong. Improved approximate string matching using compressed suffix data structures. In Proc. 16th International Symposium on Algorithms and Computation (ISAAC), pages 339–348, 2005.
[11] B. Langmead, C. Trapnell, M. Pop, and S. Salzberg. Ultrafast and memoryefficient alignment of short DNA sequences to the human genome. Genome
Biology, 10(3):R25, 2009.
[12] H. Li and R. Durbin. Fast and accurate short read alignment with BurrowsWheeler transform. Bioinformatics, 25(14):1754–1760, 2009.
[13] H. Li and N. Homer. A survey of sequence alignment algorithms for nextgeneration sequencing. Briefings in Bioinformatics, 11(5):473–483, 2010.
[14] S. Mihov and K. U. Schulz. Fast approximate search in large dictionaries.
Computational Linguistic, 30(4):451–477, 2004.
[15] G. Navarro and R. Baeza-Yates. A hybrid indexing method for approximate
string matching. J. of Discrete Algorithms, 1(1):205–239, 2000.
[16] G. Navarro and V. Mäkinen. Compressed full-text indexes. ACM Computing
Surveys, 39(1), 2007.
[17] G. Navarro and M. Raffinot. Flexible Pattern Matching in Strings. Cambridge
University Press, 2002.
[18] Gonzalo Navarro. A guided tour to approximate string matching. ACM Comput. Surv., 33(1):31–88, March 2001.
[19] L.M.S. Russo, G. Navarro, A.L. Oliveira, and P. Morales. Approximate string
matching with compressed indexes. Algorithms, 2(3):1105–1136, 2009.
[20] T. Schnattinger, E. Ohlebusch, and S. Gog. Bidirectional search in a string with
wavelet trees and bidirectional matching statistics. Information and Computation, 213:13–22, 2012.
[21] J.T. Simpson and R. Durbin. Efficient de novo assembly of large genomes using
compressed data structures. Genome Research, 22(3):549–556, 2012.
22

[22] W.-K. Sung. Indexed approximate string matching. In Ming-Yang Kao, editor,
Encyclopedia of Algorithms, pages 1–99. Springer US, 2008.

23

Appendix
The following search schemes were used in experiments described in Section 5.
For 2 mismatches or errors:
1. Slightly modified scheme SLam . The searches are: Sf = (123, 000, 022), Sb =
0
0
is 213
= (213, 001, 012). Note that the π-string of Sbd
(321, 000, 012), and Sbd
0
and not 231 as in Sbd . While Sbd and Sbd have the same efficiency for equal-size
partitions, this in not the case for unequally sized parts.
2. 4-part scheme with searches (1234, 0000, 0112),
(2341, 0001, 0012), and (1234, 0002, 0022).

(4321, 0000, 0122),

For 3 mismatches or errors:
1. 4-part scheme with searches (1234, 0000, 0133),
(3421, 0000, 0133), and (4321, 0011, 0133).

(2134, 0011, 0133),

2. 5-part scheme with searches (12345, 00000, 01233), (23451, 00000, 01223),
(34521, 00001, 01133), and (45321, 00012, 00333).
For 4 mismatches or errors:
1. 5-part scheme with searches (12345, 00000, 02244), (54321, 00000, 01344),
(21345, 00133, 01334),
(12345, 00133, 01334),
(43521, 00011, 01244),
(32145, 00013, 01244), (21345, 00124, 01244) and (12345, 00034, 00444).
2. 6-part scheme with searches (123456, 00000, 012344), (234561, 00000, 012344),
(654321, 000001, 012244), (456321, 000012, 011344), (345621, 000023, 011244),
(564321, 000133, 003344), (123456, 000333, 003344), (123456, 000044, 002444),
(342156, 000124, 002244) and (564321, 000044, 001444).

24

