We consider the problem of on-line exact string matching of a pattern in a set of highly similar sequences. This can be useful in cases where indexing the sequences is not feasible. We present a preliminary study by restricting the problem for a specific case where we adapt the classical Morris-Pratt algorithm to consider borders with errors. We give an original algorithm for computing borders at Hamming distance 1. We exhibit experimental results showing that our algorithm is much faster than searching for the pattern in each sequences with a very fast on-line exact string matching algorithm.
High-throughput sequencing or Next Generation Sequencing (NGS) technologies allow to produce a great amount of DNA sequences with a high rate of similarity. For instance, the 1000 genomes project1 aimed at sequencing a large amount of individual whole human genomes. This generates massive amounts of sequences (3 billion letters A, C, G, T) which are identical more than 99% to the reference human genome. The generated data form a collection of sequences where each di↵ers from another by a few number of di↵erences such as substitutions or single nucleotide variants (SNVs), indels, copy number variations (CNVs) or translocations to name a few.
With the large mass of available data, storing, indexing and support for fast pattern matching have become important research topics.
Pattern matching can be carried out in two ways: o↵-line by using an index or on-line when indexing is not possible. Although the first kind of solutions seems Copyright c by the paper’s authors. Copying permitted only for private and academic purposes. to be more suitable, the index issue might not seem significant in some cases even if it is compressed. The main problem we can face is to have insucient storage space to build the index. Thus one may have to scan the whole sequence rather than index it.
In this paper, we focus on o↵ering a preliminary study that allows to find the exact occurrences of a given pattern of length m in a set of highly similar sequences. We propose a solution that follows a tight analysis of the Morris-Pratt algorithm. We point out occurrences of the pattern by performing a left to right traversal over the reference sequence and at the same time we take into account variations contained in other sequences. Our approach makes a simplistic assumption that sequences include variations only of type substitutions and that there exists at most only one variation in a window of length m.
The rest of the paper is organized as follows. Section 2 presents related works. We set up notations and formalize the problem in Sect. 3. We give our new algorithm in Sect. 4 while experimental results are exhibited in Sect. 5. Finally we give our conclusions in Sect. 6. 2
Related work Storing genetic sequences of many individual of the same species is a major challenge in biological research. Basic structures usually store redundant information which lead to a memory requirement proportional to the total length of input data. Recently, several works focusing on indexing similar sequences were implemented to allow building data structures taking advantage from high similarity between the considered data. These works aim to reduce the memory requirement from the length of all input sequences to the length of a single sequence (reference sequence) plus the number of variations.
Huang et al. [
All these results are concerned with o↵-line string matching but to the best of our knowledge there exists no result for on-line string matching in a set of highly similar sequences. 3
In what follows, we consider a finite alphabet ⌃= {A, C, T, G} for DNA sequences. A string or a sequence is a succession of zero or more symbols of the alphabet. The empty string is denoted by ". The set of all non empty strings over ⌃ is denoted by ⌃ +. All strings over the alphabet ⌃ are element of ⌃ ⇤ = ⌃ + [ { ✏}. The string w of length m is represented by w[0 . . m 1] where w[i] 2 ⌃ and 0 i m 1. The length of w is denoted by |w|.
A string x is a factor (substring) of y if there exist u and v such y = uxv, where u, v, x, y 2 ⌃ ⇤ . Let 0 j m 1 be the starting position of x in y, thus x = y[j . . j + |x| 1]. A factor x is a prefix of y if y = xv, v 2 ⌃ ⇤ . Similarly a factor x is a sux of y if y = ux, for u 2 ⌃ ⇤ . A factor u is a border of x, if it is both a prefix and a sux of x, then there exist v, w 2 ⌃ ⇤ such x = vu = uw. The reverse of the string x is denoted by x⇠ . The longest common prefix between two strings u and v is denoted by lcp(u, v).
The exact pattern matching problem consists in finding all the occurrences of a
pattern x in a string y. That is, all possible j such that y[j + i] = x[i] holds for
all 0 i m 1. This problem can be extended in a very interesting way by
considering a set of sequences and find whether a given pattern occurs distributed
horizontally where di↵erent parts of the pattern can be located in consecutive
positions of di↵erent texts. More formally, given a set of sequences Y = {y0, . . . , yr 1}
of equal length n, point out all positions 0 j n m + 1, such that for
0 i m 1 we have x[i] = yg[j + i] for some g 2 [0; r 1]. This latter problem
is known as distributed pattern matching [
The problem we focus on in this paper is formally defined as follows. Let
y0, y1, . . . , yr 1 be r highly similar sequences with the same length n defined over
the alphabet ⌃. Let y0 be the reference sequence. The sequences y1, y2, . . . , yr 1
are represented by variations over y0. Thus, we consider the set Z = {(G, j, c)},
such that c = yg[j] 6= y0[j] for all 0 j n 1, g 2 G where 1 g r 1
and c 2 ⌃. Furthermore, for ( G, j, c), (G0, j0, c0) 2 Z we have |j j0| > M for some
integer M . We wish to find all occurrences of an arbitrary pattern x of length
m M in yg where 0 g r 1. This problem can be viewed as an hybrid
between distributed pattern matching and approximate string matching with k
mismatches [
We o↵er an algorithm to solve the problem described above in the same fashion as
the Morris-Pratt (MP) algorithm [
Given the pattern x of length m, we consider three cases when a prefix x[0 . . i] for 0 i m 1, is recognized when scanning the r sequences at position j: Case 1 x[0 . . i] = y0[j . . j + i] and 6 9 (G, k, c) 2 Z such that j k j + i.
This means that x[0 . . i] matches on y0 and there is no variation in all others sequences in the current window then x[0 . . i] matches equally in all sequences. Case 2 x[0 . . i] = y0[j . . j + i] and 9 (G, k, c) 2 Z such that j k j + i. This means that x[0 . . i 1] matches all sequences except yg.
Case 3 x[0 . . i] = yg[j . . j + i] and 9 (G, k, c) 2 Z such that j g 2 G . Then x[0 . . i 1] matches only sequence yg. k j + i and 4.1
The preprocessing phase consists in precomputing arrays storing positions of borders for each prefix of the pattern. Borders at Hamming distance 0 are computed as in the MP algorithm and stored in an array called mpNext. For borders at Hamming distance 1 we will use the two following arrays.
The classical array prefx is the array of prefixes of the pattern x: prefx[i] is the length of the longest prefix of x starting at position i, for 0 i m 1.
The array pref ⇠x stores for each position i, the length of the longest common prefix starting at position i when reading the pattern from right to left with each position i0 < i where lcp(x[0 . . i]⇠ , x[0 . . i0]⇠ ) 6= 0. It is defined for all 0 i m 1 and i0 < i by pref ⇠x [i] = {(i0, `) | i0 < i, x[i] = x[i0] and 0 < ` = |lcp(x[0 . . i]⇠ , x[0 . . i0]⇠ |)}. Then pref ⇠x [i + 1] can be easily computed from pref ⇠x [i].
For 0 i m, B[i] contains the sux starting positions of borders of x[0 . . i] with one change. More formally B[i] = {i0 | Ham(x[0 . . i i0], x[i0 . . i]) = 1}. Proposition 4.1 For 1 i m 1, if i0 2 B[i] then x[0 . . i x[0 . . i0 + prefx[i0] 1]x[prefx[i0]]x[i0 + prefx[i0] + 1 . . i]. i0] is a border of
Then B[i] = {i0 | 9 (i i0, `) 2 pref ⇠x [i] and i i0 = prefx[i0]+`}[{ i | prefx[i] = 0}. Proposition 4.2 The number of elements in the array B is O(m). 4.2
The searching phase consists in scanning the sequences from beginning to end. A general situation is the following: a prefix x[0 . . i] of the pattern matches at least one sequence of Y at position j. Then position j + 1 is scanned and a decision has to be made in accordance with the three cases.
Case 1 If x[i + 1] = y0[j + i + 1] then if there is no variation at position j + i + 1 in the other sequences we remain in Case 1 otherwise we move to case 2. If x[i + 1] 6= y0[j + i + 1] then if there is no variation in the other sequences at 2.8 2.6 2.4 2.2 1.6 1.4 1.2 )s 2 ( e iTm 1.8 position j + i + 1 then a shift is performed as in the MP algorithm otherwise if the variant symbol c is equal to x[i + 1] we move to Case 3 otherwise a shift has to be performed using mpNext and B.
Case 2 If x[i + 1] = y0[j + i + 1] then we remain in Case 2. If x[i + 1] 6= y0[j + i + 1] then a shift has to be performed using B.
Case 3 If x[i + 1] = y0[j + i + 1] then we remain in Case 3. If x[i + 1] 6= y0[j + i + 1] then a shift has to be performed using B.
When a shift is performed using B, the case can change according to the length of the shift and the position of the variation.
Theorem 4.3 The searching phase runs in O(n) time. 5
Experiments
We use simulated data of length 150 MB and mutation rate 0.25% among them 30%
to 50% are at the same position in di↵erent sequences. We first generate a reference
text, then mutate it at random positions with random nucleotides. We compare
our solution with FJS [9] which is one of the most ecient exact string matching
algorithm in this setting (see [
As mentioned above our algorithm takes into account the reference sequence with positions of variations. We consider variations every 500 positions. For FJS algorithm we launch the execution texts one by one. We ran the FJS algorithm on each sequence successively. Our goal is to demonstrate that from a certain number of sequences, it is more ecient to use our solution than a classical exact string matching solution. Results are shown in Figure 1 and show that our solution is faster (in our settings) when considering more that 3 highly similar sequences. 6
Discussion and conclusions
We have presented a new algorithm for searching in a set of similar sequences. The
design of the algorithm follows a tight analysis of the Morris and Pratt algorithm.
Recall that we use a suitable representation of data such that we take into account
a reference sequence with only positions of variations of the other sequences. The
searching time of our algorithm depends on the size of the reference text. However,
our solution works with a particular model, since we are limited to a certain gap
between consecutive variations. We will to improve our algorithm to overcome this
point. On the other hand, we aim to use variants of the Boyer-Moore algorithm [