Fast algorithm for scoring a sequence

Question

Let us define:

• $U = \left\{ u_j \right\}, 1 \leq j \leq N = 2^{L}$, the set of all different binary sequences of length $L$.

• $V = \left\{ v_i \right\}, 1 \leq i \leq M = \binom{L}{k}2^{k}$, the set of all different gapped binary sequences with $k$ known bits and $L-k$ gaps.

We have a score function $f : V \rightarrow R$, that associates a real number $c_i$ to any of $v_i$'s; i.e. $c_i = f(v_i)$.

Given any pair of $u_j$ and $v_i$, we can count the number of mismatches, by counting the known bits of $v_i$ that are different in $u_j$. For example given $L=7, k=4$:

$\text{diff}(\text{"0001100"}, \text{"01...01"}) = 2$,

because the second bit and the last bit are different (gaps can match to any value).

We have a weight associated to the number of mismatches. In general, the number of mismatches can be any number between $0$ and $k$. so we have $k+1$ such weights given by $\left \{ w[0],w[1],\ldots,w[k] \right \}$.

For a given sequence $S$ of $T$ bits, $S = s_0s_1 \ldots s_{T-1}$, we define the total score as the weighted sum of all the $c_i$'s for all the the contiguous $L$-bits subsequences of $S$ weighted by their distance:

$g(S)=\sum_{j=0}^{T-L}\sum_{i=1}^{M} c_i \times w[\text{diff}(v_i, s_js_{j+1} \ldots s_{j+L-1})]$.

the obvious algorithm is as follows:

g:=0; for j:=0 to T-L for i:=0 to M g:=g+c[i]*calcDif(S,j,i)

where calcDif() takes $O(L)$ to count number of mismatches. so the overall time would be in the order of $T*M*L$.

What is a better approach to calculate the total score more efficiently? Typical values are: $L=20, k=10, T=2000$?

Answer 1

The problem can be solved at least in $O(M \cdot T \log T)$ which is theoretically better, but may be slower than $O(M\cdot T\cdot L)$ in practice if $L$ is comparable to $\log T$.

Let's show how to find a score of each $v_i$ independently in $O(T \log T)$. For any (possibly gapped) string $W = w_0 w_1 \ldots w_{n-1}$ define two polynomials $zeros_W(x)$ and $ones_W(x)$ as follows:

$$ zeros_W(x) = \sum_{i=0}^{n-1} [w_i = ''0'']x^i$$ $$ ones_W(x) = \sum_{i=0}^{n-1} [w_i = ''1'']x^i$$

where $[P]$ is Iverson notation.

For example, $zeros_{01...01}(x) = 1 + x^5$, $ones_{01...01}(x) = x + x^6$, $zeros_{0001100}(x) = 1+x+x^2+x^5+x^6$ and $ones_{0001100}(x) = x^3 + x^4$.

Then compute coefficients of the following polynomial $$c(x) = zeros_S(x) \cdot ones_{rev(v_i)}(x) + ones_{S}(x) \cdot zeros_{rev(v_i)}(x) = \sum_{j}c_j x^j$$ where $rev(W) = w_{n-1} \ldots w_1 w_0$ is just reversed $W$. It's easy to see that $c_{L-1+k}$ will be equals exactly to the difference between $s_k s_{k+1} \ldots s_{k+L-1}$ and $v_i$. For example, for $S=''0100''$ and $v_i = ''10''$ $c(x)$ will be equal to $(1+x^2+x^3)\cdot x + x \cdot 1 = 2x + x^3 + x^4$ which gives us all interesting differences (2, 0, 1) so we can calculate score in linear time. Thus the problem is reduced to polynomial multiplication, which can be done in $O(T \log T)$ using FFT.

The possible disadvantage of such approach is that FFT has essential constant hidden in its big-O, so using $O(T^{1.585})$ Karatsuba multiplication or something similar may work better in practice for not very large $T$.