# an efficient algorithm for mismatch profile

We have M=10000 binary sequences of length N=1000.

given length L=15, for each pair of sequences, $S_1$ and $S_2$, we define the mismatch profile, mp($S_1$,$S_2$,$L$)[$m$], for m=0,1,...,L as following:

mp($S_1$,$S_2$,$L$)[$m$] = number of all substrings of length $L$ in $S_1$ and $S_2$ that exactly differ in m positions.

for example:

$\begin{array}{ccc} L=5 \\ N=7 \\ S1='1010101' \\ S2='1011000' \end{array}$

set of all L-mers in $S_1$ is { 10101, 01010, 10101} and set of all L-mers in $S_2$ is { 10110, 01100, 11000}

All the pairs of L-mers (one from S1 and the other from S2) are (with the number of mismatches for each pair in the parenthesis):

$\begin{array}{ccc} 10101 & 10110 & (2) \\ 10101 & 01100 & (3) \\ 10101 & 11000 & (3) \\ 01010 & 10110 & (3) \\ 01010 & 01100 & (2) \\ 01010 & 11000 & (2) \\ 10101 & 10110 & (2) \\ 10101 & 01100 & (3) \\ 10101 & 11000 & (3) \end{array}$

hence the mismatch profile for $S_1$ and $S_2$ for this example is: mp($S_1$,$S_2$,$L$)[0]= 0
mp($S_1$,$S_2$,$L$)[1]= 0 mp($S_1$,$S_2$,$L$)[2]= 4 mp($S_1$,$S_2$,$L$)[3]= 5 mp($S_1$,$S_2$,$L$)[4]= 0 mp($S_1$,$S_2$,$L$)[5]= 0

what is a fast algorithm to calculate the mismatch profile for every pair of sequences. The trivial algorithm would be $O(N^2 M^2)$ . can we make this better ?

• I think you mean "substring" not "subsequence". 11111 is a subsequence of 0010010001100000100, but not a substring. Feb 10 '12 at 8:57
• Also, this looks like a reasonable algorithms homework problem. Can you give some motivation/background? Feb 10 '12 at 8:59
• yes I meant substring. no it is not a hw problem. We use this in gapped-kmer kernel calculation of an SVM classification which has a lot of applications. The simplification I made here is that the sequences are not binary, they are usually in base four defined over an alphabet of {A,C,G,T} Feb 10 '12 at 12:12
• I assume you have at least $O(M^2 L)$ memory (at least for the answer), right? Is $L$ typically small? And on somewhat unrelated topic: was my answer on your previous question helpful? Feb 10 '12 at 14:34
• there appears to be a recursive simplication that if you have computed all the mismatches of substrings of size n-1 then you can use that to compute mismatches of strings of size n.
– vzn
Feb 10 '12 at 16:19

Upd: Space complexity improved to just $O(2^LML)$ which makes the solution much more practical.

The problem can be solved at least in $O(M^2NL + 2^L MNL)$ time and $O(2^L ML)$ additional memory as follows.

Let's denote by $S^j$ the prefix $s_1 s_2 \ldots s_j$ of $S$. First, for each binary string $s$ of length $L$ precompute mismatch profile $mp'(s, i, j, k) = mp(s, S_i^j, k)$ of that string with all prefixes of length $j$ of all strings $S_1, S_2, \ldots, S_M$. That will require $O(2^L MNL)$ memory to store and can be done in $O(2^L MNL)$ in straightforward way.

Then iteratively build mismatch profiles $mp_j(i,i',k)$ for each pair of $j$-prefixes of $S_1, S_2, \ldots, S_M$ for $j=L,L+1,\ldots,N$, so $mp_N(i,i',k)$ will be your desired answer. Let's denote by $end(S)$ the suffix of $S$ of the length $L$ (i.e. $S[(length(S)-L+1) \ldots length(S)]$). Notice that \begin{align*} mp_j(i,i',k) = &mp(S_i^{j-1},S_{i'}^{j-1},k) + mp(end(S_i^j), S_{i'}^{j-1}, k) + \\&mp(S_i^{j-1}, end(S_{i'}^j), k) + mp(end(S_i^j), end(S_{i'}^j),k) \end{align*} which is indeed \begin{align*} &mp_{j-1}(i,i',k) + mp'(end(S_i^j), i', j-1, k) + \\ &mp'(end(S_{i'}^j), i, j-1, k) + mp(end(S_i^j), end(S_{i'}^j),k) \end{align*}

The first 3 terms of that sum are already computed and the last one can be computed in amortized $O(1)$ time (it is $O(L)$ but you need to compute it only once for all values of $k$ with fixed $j,i,i'$). Thus you can compute $mp_j$ using $mp_{j-1}$ and $mp'$ in $O(M^2L)$ time, which leads to $O(M^2NL)$ term in the total complexity.

Notice that you need values of $mp'$ only for a fixed prefix length on each iteration. Thus you can compute $mp'$ not all of a piece at the beginning but rather iteratively for $j=L,L+1, \ldots, N$ in between of $mp_j$ computations. That can be easily done in $O(2^LML)$ time per iteration (so total complexity does not change) but requires only $O(2^LML)$ additional memory.

The main issue of this solution is exponential in $L$ time.

P.S. If $L$ is too small ($2^L < N$) then the naive method can be improved to $O(M^2 4^L)$ time by calculating the number of occurrences of each string $s$ of length $L$ in $S_i$. This allows you to process all different substrings of $S$ of length $L$ in $O(2^L)$ time instead of $O(N)$ (you will need to precompute pairwise differences between all such substrings first).

• Actually it can be done in $O(2^LML)$ additional memory! Feb 13 '12 at 9:50

the details would be slightly tricky to work out but (as a rough sketch) there appears to be a recursive simplification that if you have computed all the mismatches of substrings of size n-1 then you can use that to compute mismatches of strings of size n.

your problem also seems to have connections to boolean convolutions (ie using a sliding window and the mod 2 sum) & also matrix rigidity (originally proposed by Valiant)

if you are actually working with DNA bases as you indicate with your comment about ACGT, then why are you reformulating it in binary? it would probably be better to keep the ACGT formulation. the xor problem is not exactly equivalent and it would give wrong answers for subsequences off by one digit if you're really interested in a 4-symbol ACGT matching system.

you do not state this directly but from my knowledge of bioinformatics it would seem you might be interested in finding long subsequences such that there are low mismatches. if you add that criteria it significantly changes the optimization problem because one does not have to consider all mismatches, only those that are "better" than the current best record.

in that case there may be some similarity to this problem and the well studied problem of finding the longest matching substring in a string, with your case interested in the "longest near matching substring". similar mechanisms for that problem like the suffix tree may prove useful