Non-normal limiting distribution for optimal alignment scores of strings in binary alphabets

TL;DR

This paper analyzes the distribution of optimal alignment scores for binary strings, revealing a non-normal limiting distribution involving Tracy-Widom and normal components, with implications for string relatedness testing.

Contribution

It decomposes the alignment score into normal and Tracy-Widom parts, showing the non-normal limiting distribution under specific gap restrictions and alignment conditions.

Findings

Alignment scores have a non-normal limiting distribution involving Tracy-Widom.

The score decomposes into a normal part and a Tracy-Widom part under certain gap constraints.

The normal component is irrelevant for relatedness testing as it does not depend on alignment.

Abstract

We consider two independent binary i.i.d. random strings $X$ and $Y$ of equal length $n$ and the optimal alignments according to a symmetric scoring functions only. We decompose the space of scoring functions into five components. Two of these components add a part to the optimal score which does not depend on the alignment and which is asymptotically normal. We show that when we restrict the number of gaps sufficiently and add them only into one sequence, then the alignment score can be decomposed into a part which is normal and has order $O(\sqrt{n})$ and a part which is on a smaller order and tends to a Tracy-Widom distribution. Adding gaps only into one sequence is equivalent to aligning a string with its descendants in case of mutations and deletes. For testing relatedness of strings, the normal part is irrelevant, since it does not depend on the alignment hence it can be safely removed from the test statistic.

Figures

Equations

$$\begin{array}[]{c|c|c|c|c|c|c|c}Z&&A&A&T&C&U&A\\ \hline\cr X&&A&A&&C&&A\\ \hline\cr Y&&A&T&T&C&U&A\end{array}$$

$$\begin{array}[]{c|c|c|c|c|c|c|c}X&&A&A&&C&&A\\ \hline\cr Y&&A&T&T&C&U&A\end{array}$$

$$S_{\pi }(x,y)=S(A,A)+S(A,T)+S(G,T)+S(C,C)+S(G,U)+S(A,A)$$

$$S(x,y)=\underset{\pi }{\max }{S}_{\pi }(x,y),$$

$$\underset{0\le t_1<t_2<\dots <t_k\le 1}{\max }[B^1(t_1)-B^1(0)]+[B^2(t_2)-B^2(t_1)]+\dots +[B^k(t_k)-B^k(t_{k-1}))]$$

$$W^i(t)=B^i(t)+V(t),$$

$$\underset{0\le t_1<t_2<\dots <t_k\le 1}{\max }[W^1(t_1)-W^1(0)]+[W^2(t_2)-W^2(t_1)]+\dots +[W^k(t_k)-W^k(t_{k-1}))]=$$

$$=V(1)-V(0)+\underset{0\le t_1<t_2<\dots <t_k\le 1}{\max }[B^1(t_1)-B^1(0)]+[B^2(t_2)-B^2(t_1)]+\dots +[B^k(t_k)-B^k(t_{k-1}))].$$

$$VAR[L_n(S)]=O(n)$$

$$w((x_0,x_1,\dots ,x_i)):=\sum _{j=2}^{i}w((x_{j-1},x_j)).$$

$$\begin{array}[]{c|c|c|c|c|c}h&e&i&n&&i\\ \hline\cr h&e&&n&r&i\end{array}$$

$$\begin{array}[]{c||c|c|c|c|c}i&&&&&6\\ \hline\cr n&&&4&5&\\ \hline\cr i&&3&&&\\ \hline\cr e&&2&&&\\ \hline\cr h&1&&&&\\ \hline\cr\hline\cr 0&h&e&n&r&i\end{array}$$

$$1\le {\nu }_1<{\nu }_2<\dots <{\nu }_i\le n$$

$$1\le {\mu }_1<{\mu }_2<\dots <{\mu }_i\le n,$$

$$\sum _{j=1}^{i}S(x_{{\nu }_j},y_{{\mu }_j})-2(n-i)\cdot \delta ,$$

$$(({\nu }_j-1,{\mu }_j-1),({\nu }_j,{\mu }_j))$$

$$S:=\left(\begin{array}[]{ccc}S(a,a)&S(a,b)&S(a,G)\\ S(b,a)&S(b,b)&S(b,G)\\ S(G,a)&S(G,b)&S(G,G)\end{array}\right).$$

$$S_{0}:=\left(\begin{array}[]{rrr}1&1&1/2\\ 1&1&1/2\\ 1/2&1/2&0\end{array}\right),\,S_{1}=\left(\begin{array}[]{rrr}1&0&1/2\\ 0&-1&-1/2\\ 1/2&-1/2&0\end{array}\right),\,S_{2}=\left(\begin{array}[]{rrr}1&-1&0\\ -1&1&0\\ 0&0&0\end{array}\right)$$

$$S_{3}=\left(\begin{array}[]{rrr}0&0&1\\ 0&0&-1\\ 1&-1&0\end{array}\right),\,S_{4}=\left(\begin{array}[]{rrr}0&0&1\\ 0&0&1\\ 1&1&0\end{array}\right).$$

$$S=a_0{S}_0+a_1{S}_1+a_2{S}_2+a_3{S}_3+a_4{S}_4.$$

$$T(x,y):=h(x)+h(y)$$

$$T_{\pi }(X_1\dots X_n,Y_1\dots Y_n)=\sum _{i=1}^{n}h(X_i)+\sum _{i=1}^{n}h(Y_i).$$

$$h(a)=1/2,h(b)=-1/2,h(G)=0.$$

$$L_n(S_1)=(N_n^a-N_n^b)/2$$

$$L_n(S_1)=N_n^a-n.$$

$$VAR[L_n(S_1)]=VAR[N_n^a]=2n{p}_a(1-p_a).$$

$$L_n(S)=L_n(a_0{S}_0+a_1{S}_1+a_2{S}_2+a_3{S}_3+a_4{S}_4)=$$

$$=a_0{L}_n(S_0)+a_1{L}_n(S_1)+L_n(a_2{S}_2+a_3{S}_3+a_4{S}_4)=$$

$$(a_0-a_1)n+a_1\cdot N_n^a+L_n(a_2{S}_2+a_3{S}_3+a_4{S}_4).$$

$$\frac{N_n^a-2n{p}_a}{\sqrt{2n{p}_a(1-p_a)}}\approx \mathcal{N}(0,1),$$

$$L_n(a_0{S}_0+a_1{S}_1)=a_{0}n+a_1(N_n^a-n)$$

$$S=\left(\begin{array}[]{rrr}1&0&0\\ 0&-1&0\\ 0&0&0\end{array}\right)=S_{1}-S_{3}/2.$$

Full text

We show that when we restrict the number of gaps sufficiently and add them only into one sequence, then the alignment score can be decomposed into a part which is normal and has order $O(\sqrt{n})$ and a part which is on a smaller order and tends to a Tracy-Widom distribution. Adding gaps only into one sequence is equivalent to aligning a string with its descendants in case of mutations and deletes. For testing relatedness of strings, the normal part is irrelevant, since it does not depend on the alignment hence it can be safely removed from the test statistic.

1. Introduction

Optimal Alignments are among the main two components used by most DNA-alignment algorithms which try to find similar substrings in DNA-sequences. Because of its great practical importance, these problems have received a tremendous amount of attention.

Let us explain the key points from the DNA perspective. DNA sequencing is one of the most important areas of study in biology. This is still rather expensive and very complex. Suppose now that in a certain animal population, a certain gene is detected. This gene is probably likely to appear in the DNA of other animals too, but in a slightly different form. Thus, instead of repeating the same costly experiments on other animals to find the gene, they just search for a DNA-segment which looks similar to the gene in the first population. This is why often DNA-strings in different animals are compared in order to find places that are related. For more details on the analysis of DNA we refer the reader to [11] and in particular the iid assumption which appears also in this paper is explained on page 14 of the [11].

Here is a more concrete example. Assume $Z=AATCUA$ is a hypothetical DNA of an ancestor with two descendants, one with the DNA string $X=AACA$ (which lost the $T$ and the $G$) and another with the DNA equal to $Y=ATTCUA$ (which had an $A$ mutated into a $T$). We align these strings (with gaps), as follows:

[TABLE]

In the analysis of the DNA one has access to the DNA of the descendants only. In terms of the alignments, one is interested in the historic alignment. In our example, we would have that the alignment $\pi$ is given by:

[TABLE]

and we call this the historic alignment. The letters which disappeared in one descendant are now represented by gaps. Which alignment with gaps corresponds to the evolution history of the strings is not known apriori. In other words, in our example, we know the strings $x=AACA$ and $Y=ATTCUA$, but not “the true historic” alignment given here by (1.1).

To guess which alignment is historic, the natural attempt is to find an alignment which maximizes similarity. For this, one starts with a scoring function $S:\mathcal{A}\times\mathcal{A}\rightarrow\mathbb{R}$ which measures similarity between letters. (Here $\mathcal{A}$ is the alphabet augmented by one symbol for the gap). The total score of an alignment is then the sum of the scores of the aligned symbol pairs. In alignment (1.1), the alignment score $S_{\pi}(x,y)$ is equal to

[TABLE]

where $G$ stands for a gap. An alignment $\nu$, which for given strings $x$ and $y$ maximizes the alignment score $S_{\nu}(x,y)$ is called an optimal alignment. The score of an optimal alignment is called optimal alignment score and we will denote it by

[TABLE]

where the maximum above is taken over all alignments with gaps of $x$ and $y$.

Which alignment is optimal depends ultimately on the scoring function used. Often, the scoring function biologists work with is a log likelihood. For this, one assumes models of letters evolving independently of their neighboring letters. One takes as scoring function $S(A,B)$, the logarithm of the probability that an ancestor letter evolves into an $A$ in one descendant and into a $B$ in the other descendant. Here $A$ and $B$ are any two letters from the alphabet under consideration, or also the gap symbol. With this choice of scoring function, the optimal alignment becomes the maximum-likelihood estimate for the historic alignment.

A high optimal alignment score indicates relatedness of the two DNA-strings. Thus in testing for relatedness of DNA-strings one usually takes as test statistics the optimal alignment score. The significance of such a test depends then on the order of magnitude of the fluctuation of the optimal alignment score.

It has been a long open question whether for two independent i.i.d. strings of length $n$, the fluctuation of the optimal alignment score is of order $\sqrt{n}$ or $\sqrt[3]{n}$. (See [29] and [9]). In this paper, we think we can see why these two different conjectures co-exist at the same time. We show that for two independent i.i.d strings of length $n$ the optimal alignment score (depending on the parameters) can contain a normal component with fluctuation of order $\sqrt{n}$. This component is not relevant for finding the optimal alignment score, because it does not depend on the alignment.

When testing for relatedness of DNA-strings one should remove this normal component, since it is only noise. This can be done by modifying the scoring function. Once we “cut out this normal component”, the question that remain is what is the asymptotic distribution of the remainder? Is it Tracy-Widom? We can’t prove it in full generality, however we check this for the special case where the number of gaps is a small power of the length of the strings and the gaps are allowed only into one string.

*We also exhibit a scoring function with an alignment score which is of the order $\sqrt{n}$ but not normal, nor Tracy-Widom. This, should then show that any claim that the order $VAR[L_{n}(S)]=O(n)$ automatically implies asymptotic normality is wrong! (Here $L_{n}(S)$ designates the optimal alignment score of two i.i.d. random strings of length $n$.) To prove normality, one needs a more detailed knowledge of the path structure of the optimal alignments! *

The optimal alignment can be viewed as a special case of Last Passage Percolation with correlations. To our knowledge, many of the proofs for models of First and Last Passage Percolation, where the asymptotic distribution is known, go back to proving asymptotic equivalence with the following functional of Brownian motions: Consider $B^{1},B^{2},\ldots$ a sequence of i.i.d. standard Brownian motions. Also consider when $k$ goes to infinity, the functional

[TABLE]

which properly rescaled becomes Tracy-Widom [5, 14, 25]. Now, in our case, we have a similar situation, where instead we have a sequence of dependent identically distributed standard Brownian motions $W^{1},W^{2},\ldots$. However, the correlation structure is rather simple. Precisely, there exists a Brownian motion $V(t)$ and i.i.d Brownian motions $B^{1},B^{2},\ldots$ so that

[TABLE]

for all $i=1,2,\ldots$. Hence, instead of independent Brownian motions we have correlated ones and the maximum given in (1.2) can be written as

[TABLE]

Hence, we get a normal component plus a term on a smaller scale which is indeed Tracy-Widom. This is the reason why we can prove in the current paper, under a very special situation (adding gaps only into one string and keeping their number a small power of $n$), that in the limit we obtain a normal plus Tracy-Widom on a smaller scale. However, the normal part is often irrelevant for the optimal alignment problem.

Note here that the normal component $V(1)-V(0)$ does not dependent on the choice of $t_{1}<t_{2}<\ldots<t_{k}\leq 1$ (which in what follows, correspond to alignments). Again, the normal component is in some sense irrelevant for finding the optimal alignment. We should leave it out because it is only added noise.

1.1. History of this problem

We consider two random i.i.d. strings $X=X_{1}X_{2}\ldots X_{n}$ and $Y=Y_{1}\ldots Y_{n}$. In the sequel, $L_{n}(S)$ denotes the optimal alignment score of $X$ and $Y$.

Chvatal-Sankoff proved [9] that $E[L_{n}(S)]/n$ converges to a constant $\gamma_{S}$, but $\gamma_{S}$ is not known in most cases. In [15], [16] we have good confidence intervals for $\gamma$. We have been able to prove the order

[TABLE]

in many situations [23],[7], [2], [17]. This proves the conjecture of Waterman [29] (for certain specific distributions of the letters and scoring functions), but is different from what Chavatal-Sankoff have conjectured. The order conjectured by Chavatal-Sankoff is the order which one would expect knowing that $L_{n}$ is also the weight of the heaviest path from $(0,0)$ to $(n,n)$ in an LPP-formulation of our problem. (In the physics literature, [18], the order for the standard deviation of the heaviest path from $(0,0)$ to $(n,n)$ in a vast class of LPP-models is conjectured to be $n^{1/3}$ not $n^{1/2}$. For more details see also [20]).

Let us explain what LPP is and how our optimal alignment can be formulated as an LPP-problem with correlated weights. In LPP, one considers an oriented graph $(E,V)$ with a random weight function $w:E\rightarrow\mathbb{R}^{+}$ of the edges. Let $x,y\in V$. A path from $x$ to $y$ in the graph $(E,V)$ is then defined as a sequence of vertices: $x=x_{0},x_{1},x_{2},\ldots,x_{i}=y$ where $x_{j}\in V$ and $(x_{j-1},x_{j})\in E$ for all $j=2,3,\ldots,i$. The total weight of the path $x=x_{0},x_{1},x_{2},\ldots,x_{i}=y$ is given by:

[TABLE]

an optimal path from $x$ to $y$ is then a path which maximizes the total weight among all paths from $x$ to $y$.

Now, our optimal alignment problem can be formulated as an LPP on the integer lattice with the optimal alignment score $L_{n}(S)$ being the weight of the path of maximal weight from $(0,0)$ to $(n,n)$. For this we take the set of vertices $V$ to be $\mathbb{N}\times\mathbb{N}$ and the edges to go always one to the right, one up or diagonally up to the next vertex. The weight for horizontal and vertical edges is minus the gap penalty. For the edge $((i-1,j-1),(i,j))$ the weight is $S(X_{i},Y_{j})$. With this setting when we align a letter with a gap, this corresponds to moving one unit vertically or horizontally. Aligning $X_{i}$ with $Y_{j}$ corresponds to moving along the edge $((i-1),(j-1)),(i,j))$. The optimal path then defines an optimal alignment where for every edge $((i-1,j-1),(i,j))$ contained in the optimal path, aligns $X_{i}$ with $Y_{j}$. The transversal fluctuation measures how far the optimal path from $(0,0)$ to $(n,n)$ deviates from the diagonal.

Let us give an example. Take the two strings $X=heini$ and $Y=henri$. Consider the following alignment with gaps:

[TABLE]

We can represent this alignment as a path on the integer points in two dimension. We are always allowed to go to the left by one unit or up by one unit or diagonal to the next integer point. When we go horizontal or vertical, this means that we align a letter with a gap. When we go diagonal this means that we align two letters with each other. Below, we show how we represent the alignment given in (1.3) by a path on the two dimensional integer lattice:

[TABLE]

To visualize, the path represented by the last table above, draw short segments between the points: [math] to $1$, then $1$ to $2$, ….and up to a segment from $5$ to $6$. The path is then given by this sequence of segments. So, the first such segment from [math] to $1$ is diagonal, and corresponds to aligning the two $h$’s with each other. The score-contribution for this is $S(h,h)$. Then, the segment from $1$ to $2$ is diagonal and corresponds to aligning the $e$’s. Its score contribution is $S(e,e)$. The segment from $2$ to $3$ is vertical and hence corresponds to aligning a letter with a gap. In the current case, it is the first $i$ in $heini$ which gets aligned with a gap, and the score for this is the gap penalty $\delta<0$. Then $34$ is again diagonal and corresponds to aligning $n$ with $n$, with score $S(n,n)$. Then, the segment $45$ is horizontal and corresponds to aligning $r$ in $henri$ with a gap, which contributes to minus a gap penalty to the total score. The final segment $56$ is diagonal, and corresponds to aligning the last $i$ of $henri$ with the last $i$ of $heini$. Score is then $S(i,i)$. In other words, the horizontal or vertical moves contribute minus a gap penalty each, whilst the diagonal moves contribute to the score by aligning the corresponding letters. IN THIS EXAMPLE, WE ALIGNED ONLY IDENTICAL LETTERS. This is not always the case. In DNA-analysis, when letters can mutate we also align similar letters with each other, not just identical.

Formally we can view an alignment with gaps of the strings $x=x_{1}x_{2}\ldots x_{n}$ and $y=y_{1}y_{2}\ldots y_{n}$ as a couple consisting of two increasing integer sequences. The sequence $\nu$

[TABLE]

and the sequence $\mu$

[TABLE]

are non-negative integers. Here $i\leq n$ is any integer non-negative number. The alignment with gaps defined by $(\mu,\nu)$ is then the alignment which aligns $x_{\nu_{j}}$ with $y_{\nu_{j}}$ for all $j=1,2,\ldots,i$ and aligns all remaining letters with a gap. The alignment score is then

[TABLE]

where is $2(n-i)\delta$ represents the total gap penalty. (The penalty for one given gap being $-\delta$). In our current example, we have the number of aligned letter pairs is $j=4$. Furthermore, we find: $\nu_{1}=1,\nu_{2}=2,\nu_{3}=4,\nu_{4}=5$, whilst $\mu_{1}=1,\mu_{2}=2,\mu_{3}=3,\nu_{4}=5$. The path representing the alignment given by $(\nu,\mu)$ is the path containing all the segments:

[TABLE]

for $j=1,2,\ldots,i$ and additionally a minimum amount of vertical and horizontal unit length integer segments so as to make this a connected path.

First Passage Percolation (FPP) and LPP are part of a vast area of statistical physics [18] which is concerned with random growth models for which physicists expect some universality properties. More specifically, one considers growth of a cluster where material is being attached randomly on the surface of a nucleus. There are many fundamental questions open for decades, such as the universality of the fluctuation exponents. Many theoretical physicists, crystallographers and probabilists have worked on this and many conjectures have been formulated [18], as for example that the fluctuation should behave like $n^{1/3}$ and the transversal fluctuation should be or order $n^{2/3}$ (see KPZ-conjecture in [20]). These are difficult conjectures and rigorous results were obtained only under some special circumstances of LPP models, like Longest Increasing Subsequence [3], LPP on $\mathbb{Z}^{2}$ with exponential or geometric waiting times [19] and the case of percolation in a narrow tube [8]. To our knowledge, the best results show only a fluctuations of less than $n/\ln(n)$ [4] and above $\ln(n)$ [27].

We have been able to prove our results for several Longest Common Subsequence (LCS) and Optimal Alignments (OA) models.

LCS are a special case of OA. For two strings $x=x_{1}x_{2}\ldots x_{n}$ and $y=y_{1}y_{2}\ldots Y_{n}$ a common subsequence of $x$ and $y$ is a sequence which is a subsequence of $x$ and at the same time a subsequence of $y$. A LCS of $x$ and $y$ is a common subsequence of $x$ and $y$ of maximum length. The length of the LCS of $x$ and $y$ can be viewed as an optimal alignment score of $x$ and $y$. For this we take $S(a,b)$ to be $1$ if $a=b$ and $S(a,b)=$ otherwise. On top the gap penalty should be taken [math]. With this choice, the length of the LCS is also equal to the optimal alignment score.

Our results together with the results of Baik, Deift and Johansson [3], are among the few rigorous results determining the exact order of the fluctuation exponents for some non-trivial FPP or LPP related 2-D models. In their work, as in ours, the fluctuation is of order $n^{1/3}$ and the rescaled asymptotic limiting law is Tracy-Widom. Another interesting paper is for the TASEP model, by Borodin, Ferrari, Praehofer and Sasamoto [1]. There they do not use the same approach of reducing the problem to formula (3.11) and (3.12), and it is not the GUE which appears in the limit, but GOE. Other methods used to establish the solvability of LPP models are the vertex operator and fock space formalism put foward by Okounkov [26] in the study of Schur process.

Our paper is organized as follows. In Section 2, we introduce the main decomposition of the scoring function and show examples where the limiting distribution of the optimal alignment score is normal alongside with other examples which are non-normal. In the next section, namely Section 3, we restrict the number of gaps to a certain number $k$ which is finite and fixed and allow the gaps only in one sequence. Here we give a general description of the optimal alignment in terms of multiple dependent random walks and outline an argument which supports the conjecture that under some conditions, if $k$ is replaced by $n^{\alpha}$ with small $\alpha$, then the (rescaled) limiting distribution of the optimal alignment is Tracy-Widom. We can not prove this statement rigorously in general, but in the next section, namely Section 4, we prove this for the key scoring function $S(x,y)=xy$ under the assumption that the number of gaps is $k=n^{\alpha}$ for small $\alpha$ and the second sequence of letters is distributed normally rather than $\pm 1$. The main rigorous result is stated in Theorem 4.1 and appears in Section 4.

Section 5 is dedicated to technical results while the last one, Section 6 contains a few simulations.

2. Linearly decomposing the space of scoring functions

In this section we discuss the general decomposition of the scoring function in a natural basis. The point of this is that some of the scoring functions are in some sense artificial because the optimal score does not depend on the alignment. This will give a normal behavior in the limit for the optimal score.

We consider symmetric scoring for a binary alphabet. We do not allow alignments of gaps with gaps, so the scoring functions we consider are given by a $3\times 3$ matrix, with one irrelevant entry, namely

[TABLE]

Here $G$ stands for a gap. In our case we will put $S(G,G)=0$ and we will always assume that the above matrix is symmetric hence $S(a,b)=S(b,a)$ and $S(a,G)=S(G,a)$, whilst $S(b,G)=S(G,b)$. The linear space of such $3\times 3$ symmetric matrices is generated by the following five elements

[TABLE]

and

[TABLE]

Note that the matrices $S_{0}$, $S_{1}$ and $S_{2}$ are orthogonal to each other (with respect to the canonical inner product given by $\langle A,B\rangle=Tr(AB)$). Also, the matrices $S_{3}$ and $S_{4}$ concern only with the gap penalty. Thus, any scoring function we will consider can be written as

[TABLE]

First, it is important to note that the part $a_{0}S_{0}+a_{1}S_{1}$ gives the same alignment score for any alignment. This is so, because $S_{0}$ and $S_{1}$ can both be viewed as linear additive functions in the following way. For a given a function $h:\{a,b,G\}\rightarrow\mathbb{R}$, take the “linear “scoring function $T$ defined by:

[TABLE]

for all $x,y\in\{a,b,G\}$.

We also assume that $h(G)=0$. Clearly for such a “linear” scoring function, with value $h(G)=0$, the alignment score for any alignment $\pi$ does not depend $\pi$. Indeed, it is trivial to check that

[TABLE]

Furthermore, the expression on the right side of the above equality is a linear function of a binomial random variable. Now, in our case, $S_{1}$ is of the above type with the function $h$ given by

[TABLE]

Thus, we find that

[TABLE]

where $N_{n}^{a}$ denotes the total number of $a$’s in both strings $X$ and $Y$ combined whilst $N_{n}^{b}$ denotes the number of $b$’s in both strings combined. Noting that $N_{n}^{a}+N_{n}^{b}=2n$, we finally find that

[TABLE]

Here $N_{n}^{a}$ is a binomial with parameter $2n$ and $p_{a}=P(X_{1}=a)$, consequently,

[TABLE]

Similarly, the alignment score according to $S_{0}$ does not depend on the alignment and is also equal to $n$. Since the alignment scores according to $S_{0}$ and $S_{1}$ do not depend on the alignment we find for any scoring function $S=a_{0}S_{0}+a_{1}S_{1}+a_{2}S_{2}+a_{3}S_{3}+a_{4}S_{4}$ that the optimal alignment score can be decomposed as follows:

[TABLE]

According to the central limit theorem, $N_{n}^{a}$ is asymptotically normal on the scale $\sqrt{n}$ and thus we can conclude that

[TABLE]

where $p_{a}=P(X_{1}=a)$.

Now, since the part $L_{n}(a_{0}S_{0}+a_{1}S_{1})$ does not depend on the alignment it has no affect on which alignment is optimal. Also, it is just noise added for testing if sequences are related or not. This noise can be removed, since it is known. Equivalently, we can just take $a_{2}S_{2}+a_{3}S_{3}+a_{4}S_{4}$ as scoring function instead of $S=a_{0}S_{0}+a_{1}S_{1}+a_{2}S_{2}+a_{3}S_{3}+a_{4}S_{4}$.

2.1. Normal and non-normal cases

In this section we present a simple example for which the optimal alignment score, properly rescaled, converges to a distribution which is not normal. Interestingly, the scaling is still the same as in the case of the previous section, namely, $\sqrt{n}$.

As mentioned, for scoring functions of the type $S=a_{0}S_{0}+a_{1}S_{1}$, the alignment score does not depend on the alignment $\pi$ of the strings $X_{1}\ldots X_{n}$ and $Y_{1}\ldots Y_{n}$. Furthermore, the limiting distribution of the optimal alignment score is normal on the scale $\sqrt{n}$ since

[TABLE]

where $N^{a}_{n}$ is the total number of $a$’s in the concatenated string $X_{1}\ldots X_{n}Y_{1}Y_{2}\ldots Y_{n}$. Since $N^{a}_{n}$ is binomial with parameter $2n$ and $p_{a}=p(X_{1}=a)$, the above optimal score, properly rescaled converges to the normal. Thus we can interpret the scoring function with $a_{2}=a_{3}=a_{4}=0$ as degenerate and hence not relevant for practical purposes.

Next we consider the following scoring function

[TABLE]

If $N_{a}^{X}$ designate the number of $a$’s in $X=X_{1}X_{2}\ldots X_{n}$ and $N_{a}^{Y}$ designate the number of $a$’s in $Y=Y_{1}\ldots Y_{n}$, then, the optimal alignment score according to the scoring function $S$ is simply equal to the minimum of $N_{a}^{X}$ and $N^{Y}_{a}$, namely,

[TABLE]

This is so because aligning a $b$ with a $b$ gives a negative score, whilst aligning something with a gap gives [math]. So, the optimal alignment is going to align all $b$’s with gaps. Since an $a$ aligned with an $a$ scores positively, the optimal alignment is going to align a maximum number of $a$’s with each other.

Note that $N_{a}^{X}$ and $N_{a}^{Y}$ have the same expectation and both tend asymptotically to be normal when re-scaled by $\sqrt{np(1-p)}$. Hence, we get

[TABLE]

where $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ are two standard normal independent of each other, which is certainly not normal. Thus, we have a non-normal limiting distribution, but still we have a fluctuation on a scale $\sqrt{n}$!

Notice also the following interesting fact, namely that there are many different alignments, which are very far away from each other. Indeed, say that $N_{a}^{X}>N_{a}^{Y}$. Then, the optimal alignments can be characterized as follows. Choose a subset of $a$’s in $X$ with cardinality equal to $N_{a}^{X}-N_{a}^{Y}$. These $a$’s are aligned with gaps. Then align all the remaining $a$’s in $X$ with $a$’s from $Y$. Finally, align all the $b$’s with gaps.

Now, real life scoring functions should not have $a_{4}$ too small, since in practical situations, like genetic evolution, there are usually not too many gaps. This reflects the fact that there are not too many deletes from an ancestors DNA. In the next section we show that if we look for the optimal alignment under the constrain of the total number of gaps being small, then we get that $L_{n}(a_{2}S_{2}+a_{3}S_{3}+a_{4}S_{4})$ can be decomposed into a normal part on the scale $\sqrt{n}$ plus a smaller part whose limiting distribution is non-normal.

3. Fixing the number of gaps and allowing them only in one

sequence

In this section, under the assumptions that the number of gaps is finite and fixed, we present a general outline of an argument which should be sufficiently convincing that the limiting behavior of the optimal alignment is normal in the scale $\sqrt{n}$ plus a Tracy-Widom on a smaller scale.

Let us assume that we only consider alignments with a fixed number of gaps equal to $k$, where $k$ is a fixed constant. To simplify the problem we will allow the gaps only in one sequence.

Let $k$ be a constant and let $X=X_{1}X_{2}\ldots X_{n-k}$ and $Y=Y_{1}Y_{2}\ldots Y_{n}$ be the two random strings of letters. Thus, we will align exactly $k$ letters from $Y$ with gaps, but all the letters of $X$ will be aligned with letters of $Y$. We will consider the optimal alignment of $X$ and $Y$ under the constrain that there are exactly $k$ gaps all inserted in $X_{1},\dots X_{n-k}$. The corresponding optimal score will be denoted by $L_{n}^{k}(S)$. To specify an alignment we simply need to specify which $k$ letters of $Y$ get aligned with gaps, thus we need to specify a sequence of $k$ increasing integers $0<c_{1}<c_{2}<\ldots<c_{k}\leq n$. The vectors consisting of these integers then define our alignment, namely, align $Y_{c_{i}}$ with a gap for all $i=1,2,3,\ldots,k$.

Let $R^{i}$ denote the random walk defined recurrently on $j$ by

[TABLE]

and then:

[TABLE]

for all $i=1,2,\ldots,k+1$ and for $j\in[1,n-1]$. (To simplify the notation we assume that the sequence of the $X_{i}$’s is a double infinite sequence of i.i.d. variables

[TABLE]

In the problem which we will consider only a finite number of $X_{i}$’s will have $i\leq 0$ and their contribution will be negligible).

For example, the first random walk $R^{1}$ is then defined by

[TABLE]

Similary, the second random walk $R^{2}$ is defined as follows:

[TABLE]

Consider the alignment score under the constrain that the $l$-th gap is aligned with $Y_{c_{l}}$, for all $l=1,2,\ldots,k$, where $0<c_{1}<\ldots<c_{k}\leq n$ is a given increasing sequence of integers. The alignment with exactly $k$ gaps is uniquely defined by this constrain and its alignment score is:

[TABLE]

Recall that we assume a zero gap penalty.

Let us give an example. Take the string $Y=Y_{1}Y_{2}Y_{3}Y_{4}Y_{5}Y_{6}Y_{7}Y_{8}$ and $X=X_{1}X_{2}X_{3}X_{4}X_{5}X_{6}$ we could align $Y_{3}$ and $Y_{6}$ with a gap. Then we have $k=2$ gaps and $c_{1}=3$ whilst $c_{2}=6$. So, this defines the $2$-gap-alignment $\pi$ given by:

[TABLE]

Hence, the alignment score with a [math]-gap penalty is equal to:

[TABLE]

Let us denote by $L_{n}^{k}(S)$ the optimal alignment score under the constrain that the alignment contains exactly $k$ gaps all aligned with letters from $Y$. Since the gap penalty is [math] the formula for $L_{n}^{k}(S)$ is then given by

[TABLE]

Hence, the maximum is obtained by choosing the increasing sequence $0<c_{1}<c_{2}<\ldots<c_{k}<n$ which maximizes the following sum:

take the increment of the first random walk on $[0,c_{1}]$, then we add the increment of the second random walk on $[c_{1},c_{2}-1]$ and keep on going until the increment of the last random walk on $[c_{k},n]$ has been added.

Next, we are going to put the different one-dimensional random walks $R^{i}$ together into one vector:

[TABLE]

In this manner $\vec{R}$ becomes a $k+1$-dimensional ”correlated” random walk. Thus we have $\vec{R}(0)=(0,0,0,\ldots,0)^{T}$ and

[TABLE]

If we subtract the drift and rescale by $\sqrt{n}$, the multidimensional correlated random walk $\vec{R}(t)$ converges to a $k+1$ dimensional Brownian motion and we refer to [21, 22] for more details. This means, that holding $k$ fixed, we can find a $\delta_{n}$ which goes to [math] as $n$ goes to infinity so that there exists a $k+1$-dimensional Brownian motion $\vec{B}(t)=(B^{1}(t),B^{2}(t),\ldots,B^{k+1}(t))^{T}$ for which

[TABLE]

holds with probability $1-\epsilon_{n}$ where $\epsilon_{n}$ goes to [math] as $n\rightarrow\infty$. Here the Brownian motion $\vec{B}(t)$ depends on $n$, but of course, its distribution does not. If $s\in[0,n]$ is not an integer then we take for $\vec{R}(s)$ the linear interpolation of $k\mapsto\vec{R}(k)$. Let us introduce a small modification of $L_{n}^{k}(S)$. This modification will be denoted by $\tilde{L}_{n}^{k}(S)$ and is equal to

[TABLE]

Clearly with this definition,

[TABLE]

where $C:=\max_{a,b\in\mathcal{A}}|S(a,b)|.$ Note that for all the random walks $R^{i}$, the drift is the same and equal to $E[S(X_{1},Y_{1})]$. Now, we can rewrite $\tilde{L}^{k}_{n}(S)$ using the unbiased random walk

[TABLE]

Assume next that inequality (3.2) holds. Then, we can approximate the unbiased random walks on the right side of (3.4) by components of a Brownian motion according to (3.2). Since the sum contains $2(k+1)$ terms, the approximation error will be at most $\delta_{n}\cdot 2(k+1)$ (if we assume (3.2) to hold) and thus

[TABLE]

where $L^{k}(\vec{B})$ refers to the maximum

[TABLE]

We can now use inequality (3.5) together with (3.3) to find

[TABLE]

The component of our Brownian motions $\vec{B}(t)$ could be correlated. To see this, let’s calculate the covariance between the first two random walks, that is between $R^{1}(n)-R^{1}(0)$ and $R^{2}(n)-R^{2}(0)$. The idea is that most terms of the sum corresponding to the random walk are not correlated, however the closer ones are correlated. To the point, take for example the first term $S(X_{1},Y_{1})$ in the sum

[TABLE]

Since, the $X_{i}$’s and $Y_{i}$’s are all independent of each other, an expression of the type $S(X_{i},Y_{j})$ can only be correlated with $S(X_{1},Y_{1})$ if either $i=1$ or $Y=1$. Therefore, in the sum

[TABLE]

there are only two terms which can have non-zero correlation to $R^{1}(1)-R^{1}(0)=S(X_{1},Y_{1})$. These terms are $S(X_{0},Y_{1})$ and $S(X_{1},Y_{2})$. This analysis holds true for every term in the sum for the first random walk. Only in the very last term in the sum on the right side of (3.7) there is a term for which there are not two terms in the sum on the right side of (3.8) correlated with it. Hence, the covariance between (3.7) and (3.8) is about $2n$ times constant, where the constant is given as the covariance between $S(X_{1},Y_{1})$ and $S(X_{0},Y_{1})$.

Let us expand on this calculation a little bit by looking at the covariance terms

[TABLE]

We can work out the same calculation for all pairs $i,j\leq k+1$ of random walks $R^{i}$ and $R^{j}$ provided $i\neq j$. The covariance is always about $2n\cdot COV(S(X_{1},Y_{1}),S(X_{0},Y_{1}))$. Hence, the multi-dimensional Brownian motion in our approximation (3.2) has a covariance matrix with all non-diagonal entries equal to $COV(S(X_{1},Y_{1}),S(X_{0},Y_{1}))$. Hence, we find

[TABLE]

for $i\neq j$ and

[TABLE]

This is a simple covariance structure since there is one common factor contained in all the Brownian motions $B^{i}(t)$ if $COV(S(X_{1},Y_{1}),S(X_{0},Y_{1}))$ is non-zero. In other words, we find that there exists a multidimensional Brownian motion $\vec{W}(t)=(W^{1}(t),W^{2}(t),\ldots,W^{k+1}(t))$ with all components independent and having variance

[TABLE]

so that

[TABLE]

where $I$ represents the $k+1$-dimensional identity matrix and $V(t)$ is a one dimensional Brownian motion which is independent of $\vec{W}(t)$. (To understand more in detail how we get the $V$ and the $W$ consider this: simply take a Brownian motion $V$ with [math] expectation starting at [math] and independently a collection of Brownian motions $W^{1},W^{2},\ldots,W^{k+1}$ which among themselves are i.i.d. Also, the $W_{j}$’s have all [math] expectation and start at [math]. You can now chose all these Brownian motion, so that $V+W^{1},V+W^{2},\ldots,V+W^{k}$ have the same correlation structure then the collection $B^{1},B^{2},\ldots,B^{k+1}$. Then, because you deal with gaussian process you have the same distribution. Then you can make a perfect coupling between the $V+W^{1},V+W^{2},\ldots,V+W^{k}$ and the $B^{1},B^{2},\ldots,B^{k+1}$ since they have the same distribution). We also have

[TABLE]

Hence, we find that

[TABLE]

where

[TABLE]

It is known (see [5, 14, 25]), that the rescaled $L^{k}(W)$ goes asymptotically to Tracy-Widom, more precisely,

[TABLE]

where convergence is in law as $k\rightarrow\infty$. Hence, there is a sequence of random variable $\mathcal{E}_{k}$ which goes to [math] in probability as $k\rightarrow\infty$, and a sequence of variables $F_{TW}^{k}$ each having same Tracy-Widom distribution, so that

[TABLE]

We can now rewrite expression (3.6) using (3.10) and get

[TABLE]

which, with the help of equation (3.13), becomes

[TABLE]

For the above inequality to be useful, we need the right side of it to be smaller than $\frac{1}{k^{1/6}}$. This will depend on the value of $\delta_{n}$. Now, we want to let $k$ go to infinity at a slower rate then $n$. Most papers on closeness of random walk and Brownian motion are for fixed dimension. Also, in our case, any $k$ consecutive components of the multi-dimensional random walk $\vec{R}_{i}$ are correlated. However, still assuming that the approximations hold for fixed dimension and correlated multidimensional random walk, we could take $\delta_{n}$ to be

[TABLE]

whilst

[TABLE]

for $a>0$ a constant not depending on $n$ or $k$. In this case for the inequality (3.15) to be useful, we would need

[TABLE]

to be of smaller order than $\frac{1}{k^{1/6}}$. This is the case when $k=n^{\alpha}$, with $\alpha$ small enough.

4. Letting the number of gaps $k=n^{\alpha}$ and normal letters

In this section, we are going to treat the case where the number of gaps $k$ increases to infinity with $n$ and the simplest such scenario is the case of

[TABLE]

for some small $\alpha>0$ which we assume is not depending on $n$. Again, we look for the best alignment score under all alignments with the prescribed number of gaps being $k$ and the gaps appear only in the sequence of $X$’s, thus we allow gaps to align with $Y$’s.

The highest alignment score among all alignments which align exactly $k$ letters of $Y=Y_{1}\ldots Y_{n}$ with gaps, but no letter of $X=X_{1}\ldots X_{n}$ with gaps, is again denoted by

[TABLE]

When the number of gaps is fixed, in determining which alignment is optimal, the penalty for gaps does not matter much. Indeed, the gap penalty will affect the optimal alignment only by a term of order $n^{\alpha}$. Moreover, if we assume that the penalty of aligning a letter with a gap does not depend on the letter, we can simply assume that the gap penalty is actually equal to [math].

Therefore, the only part of the scoring functions which matters is the $2\times 2$ symmetric matrix

[TABLE]

The space of all such $2\times 2$ symmetric matrices constitute a $3$ dimensional space where a basis is formed by the matrices

[TABLE]

In addition, if we define two scoring matrices to be equivalent if the optimal alignments are the same, then it is easy to see that for any arbitrary scoring function $S$, the scoring itself $S$ and $S-aS_{0}-bS_{1}$ are equivalent. Thus the main task is to understand what happens under alignments with the scoring function $S_{2}$. Thus we assume for the rest of the section that the scoring function $S$ is equal to $S_{2}$, so that:

[TABLE]

Now, in the case considered here of a two dimensional alphabet, that scoring function can be viewed as a product scoring function by observing that if we introduce a numeric code for the letters, namely take $a$ to be $+1$ and $b$ to be $-1$. Then we find that

[TABLE]

Using this in (3.10), the key observation is that the part $V(t)$ has covariance given in formula (3.9) by

[TABLE]

which in the present case becomes

[TABLE]

due to independence and mean [math] of $X_{1},Y_{1},Y_{2}$. Hence, the term $V(t)$ disappears from inequality (3.15) for this choice of the scoring function.

Next, we notice that for a one dimensional “uncorrelated” random walk, the best coupling of the walk with a Brownian motion in (3.2) can be attained with $\delta_{n}=\frac{\ln(n)}{n^{1/4}}$ (see [22]).

With this choice of $\delta_{n}$, the term $\frac{2C\times k}{\sqrt{n}}$ in (3.15) is of order $O(n^{\alpha-1/2})$, thus is much smaller than $2\delta_{n}(k+1)$ which is of order $O(n^{\alpha-1/4}\ln(n))$. However, since we have to correlate a large number of random walks (exactly $k$ and not independent), we expect $\delta_{n}$ to be larger than $O(\ln(n)/n^{1/4})$.

A different look at the dominating term on the right hand side of (3.15), gives that

[TABLE]

which is enough (at least heuristically) to guarantee that the rescaled optimal score is asymptotically Tracy-Widom.

In order to deal with the approximation of the random walk by Brownian motions, we need to get good couplings of the random walk with Brownian motions. This is what we want to do next.

We should also quickly mention, that at this stage, it should be heuristically clear why we should expect Tracy-Widom as limiting distribution as soon as the constant $\alpha>0$ is small enough. Indeed, when $\alpha>0$ is very small, then $\delta_{n}$ should be as “close as we want” to the formula for the one dimensional random walk with independent increments. Now, for the one-dimensional random walk, it is known to allow a coupling to a Brownian motion with $\delta_{n}=\frac{\ln(n)}{n^{1/4}}$. Again, the only problem is to find a good coupling of the multidimensional random walk $\vec{R}(t)$ to a multidimensional Brownian motion. The random walk $\vec{R}(t)$ has $k+1$ components. Furthermore, there is dependency among the increments of the random walk. Each step (increment) of the random walk, can depend on $k$ steps before and after. Hence, $\vec{R}(t)$ is independent of $\vec{R}(t+k+1),\vec{R}(t+k+2),\vec{R}(t+k+3),\ldots$ for all $t\geq 0$.

The main question is thus how well can a random walk with dependencies be approximated by Brownian motion? (See [21, 22] for the one dimensional case, [12] for severeal dimensions and [24, 10] for overviews.)

Is it possible to find such a coupling when the dimension of the random walk grows to infinity? This is delicate as the constants involved depend on the dimension and thus are not very useful in our case.

Unfortunately we can not treat both problems simultaneously and it seems highly technical to deduce such a result in great generality. Instead, we present here a simple approach for a simplified situation.

First, we use the fact that our scoring function is a product. Furthermore, for the rest of the paper we also assume the variables $Y_{1},Y_{2},\ldots$ to be i.i.d. normal with

[TABLE]

instead of binary variables. The variables $X_{1},X_{2},\ldots$ are still binary variables in this section with $P(X_{i}=-1)=P(X_{i}=+1)=0.5.$

Again, since $E[X_{i}]=0$ the component $V(1)-V(0)$ in the equation (3.15) dissapears.

The normal assumption on $Y$’s is in place here to simplify the calculations. We believe that macroscopically it amounts to the same as taking these variables binary. The fact that the variables $Y_{1},Y_{2},\ldots$ are standard normal instead of binary has a very important consequence on the calculation side, namely, the optimal alignment still tries to align the maximum number of letter with the same sign and also will try to get that those $Y_{i}$’s which get aligned with an $X_{i}$ of opposite sign to be small in absolute value.

The key strategy of the proof is to condition on the sequence $X$ and exploit the normal structure of the $Y$. To set up the ground, let $X=(X_{1},X_{2},\ldots,X_{n})$ and let $x=(x_{1},x_{2},\ldots,x_{n})\in\{-1,1\}^{n}$ be a non-random binary vector. We are going to work conditional on $X=x$.

Now, consider our random walk $\vec{R}$ in the current case where $S(a,b)=a\cdot b$. We have

[TABLE]

where $\vec{R}(0)=(0,0,0,\ldots,0)^{t}$ and

[TABLE]

Looking at the above expression for the random walk $\vec{R}$, we see that the steps are dependent but once we condition the $X_{i}$’s to be constants, i.e. $X=x$ and $Y_{1},Y_{2},\ldots$ to be i.i.d $N(0,1)$, the $\Delta R(i)$’s become independent normal vectors. Hence, conditional on $X=x$, we get that the multidimensional random walk $\vec{R}$ has independent normal increment with expectation [math].

This random walk is not yet a Brownian motion restricted to discrete times but we are going to show that we can construct such a Brownian motion which is a good approximation. In order to do this we consider at first the random walk $\vec{R}(t)$ restricted to the times $t$ which are a multiple of

[TABLE]

with some small $\beta>0$ which will be determined later and $[x]$ denoting the integer part of $x$.

Thus, we will consider

[TABLE]

where to simplify notation we assume $n$ to be a multiple of $j$, though we can formally replace $n$ by $jn_{j}$ with $n_{j}=[n/j]$ and run essentially the same argument.

Using the independence of the increments conditional on $X=x$, we calculate the covariance matrix of the position of the random walk at time $j$ in the following way:

[TABLE]

The key fact here is that with very high probability the last sum in the above equation is approximately $jI$ where $I$ is the $k+1$ identity matrix.

Now, using this, and the fact that $jI$ is the covariance matrix of a standard $k+1$-dimensional Brownian motion at time $j$ we will show that there is a coupling between the random walk and a standard Brownian motion such that $\vec{R}(j)$ is close to the Brownian motion position at time $j$. Moreover the argument can be extended to include also the increments of the random walk. More precisely, we show that the coupling can be chosen such that the increments

[TABLE]

are also close to the increment of the coupled Brownian motion.

To construct such a coupling, we will build a Brownian motion

[TABLE]

by coupling at the times $j,2j,3j,\ldots,n$ at first. More precisely, we couple $\vec{R}(j)$ with $\vec{C}(j)$, then the increment $\vec{R}(2j)-\vec{R}(j)$ with the increment $\vec{C}(2j)-\vec{C}(j)$ and then we continue coupling the next increment $\vec{R}(3j)-\vec{R}(2j)$ with $\vec{C}(3j)-\vec{C}(2j)$ and so on and each step we do this independently of the construction at the previous steps. Recall that we do this conditionally on $X=x$.

The key is that the construction at each step is a consequence of Lemma 5.1 in the next Section 5.1 where we show how to couple a normal random variable having covariance close to a multiple of the identity. Essentially the main assertion here is that if $\vec{A}$ is a normal vector such that

[TABLE]

then we can find a normal vector $\vec{B}$ of i.i.d normal components $N(0,\sigma^{2})$ such that $|Cov(\vec{A}-\vec{B})|\leq\epsilon$. Here we use the operator norm $|A|$ for a matrix and measure the closeness in the operator norm, in particular all components of $\vec{A}-\vec{B}$ have variance at most $\epsilon$.

Once all the increments $\vec{C}((i+1)j)-\vec{C}(ij)$ are defined, we define $\vec{C}(t)$ for $t=j,2j,3j,\ldots,n$ as the sum

[TABLE]

Conditional on $X=x$, the pieces

[TABLE]

for $i=1,2,\ldots,n/j$ are independent. Since

[TABLE]

depend only on (4.9), this implies that conditional on $X=x$, the sequence of random vectors

[TABLE]

constitutes a multidimensional standard Brownian motion restricted to the times $0,j,2j,3j,\ldots,n$. We can for instance use the construction of Levy for the Brownian motion to extend these vectors to a Brownian motion such that the increments are given exactly by $\vec{C}((i+1)j)-\vec{C}(ij)$ for $i=1,2,\dots,n/j$. We do not look for other details about the Brownian motion during the times $[ij,(i+1)j]$, we will simply use some basic estimates for the maximum of Brownian motions between these time intervals. For instance using the reflection principle we can show that

[TABLE]

In particular what this yields is that for some constant $c$ independent of $n$ and $\alpha$

[TABLE]

thus the deviation of the Brownian motion on each of these intervals is very well controlled.

Now define a few events related to our coupling which will take care of the main estimates.

• •

Let $\bf{F^{n}}$ be the event that none of the components of the standard multidimensional Brownian motions $\vec{C}(t)$ goes further than $\log(n)\sqrt{j}$ during any of the time intervals $[ij,(i+1)j]\subset[0,n]$ for $i=0,1,2,\ldots,n/j$ That is let $F^{n}_{l}$ be the event that

[TABLE]

for all $i=1,2,\ldots,n/j$. Set

[TABLE]

• •

Let $\bf{G^{n}}$ be the event that the multidimensional random walk $\vec{R}(t)$ does not move further than $\log(n)\sqrt{j}$ during any of the time intervals $[ij,(i+1)j]\subset[0,n]$. Hence, let $G^{n}_{l}$ be the event

[TABLE]

for all $i=1,2,\ldots,n/j$. Let

[TABLE]

• •

The event $\bf{A^{n}}$ takes care of the fact that the covariance matrix of $\vec{R}(ij)-\vec{R}((i-1)j)$ for all $i=1,2,\ldots,n/j$ is close to $jI$,were $I$ denotes the $(k+1)\times(k+1)$ identity matrix. More specifically, let thus $A^{n}_{1}$ be the event that

[TABLE]

where $C>0$ is a constant which will be defined later. Similarly for any integer $i$ so that

[TABLE]

holds, let $A^{n}_{i}$ be the event that

[TABLE]

Finally let $A^{n}$ be the event

[TABLE]

• •

The next event $\bf{D^{n}}$ assures that the random walk $\vec{R}(t)$ and the Brownian motion $\vec{C}(t)$ are close to each other at times which are multiples of $j$. More specifically, let $D^{n}_{l}$ be the event that the $l$-th component of $\vec{C}(t)-\vec{R}(t)$ at times which are multiples of $j$ is bounded as follows:

[TABLE]

Let $D^{n}$ be intersection of these for all $l$, precisely:

[TABLE]

Now, we summarize the picture we have obtained so far. Rescale the random walk and the Brownian motion $\vec{C}_{t}$ by a factor $\sqrt{n}$ and the time by a factor $n$ by setting for $t\in[0,1]$

[TABLE]

By scaling the Brownian motion, we know that $\vec{B}(t)$ is again a Brownian motion. When the events $D^{n}$, $F^{n}$ and $G^{n}$ all hold, inequality (3.2) is satisfied for the multi-dimensional Brownian motion $\vec{B}_{t}:[0,1]\rightarrow\mathbb{R}^{k}$ with $\delta_{n}$ given as

[TABLE]

where again, the number of gaps $k$ and the time step are given by

[TABLE]

We can now go back to inequality 3.15, to find that:

[TABLE]

holds with probability at least $1-\epsilon_{n}$ where

[TABLE]

Recall also, that $\mathcal{E}_{1},\mathcal{E}_{2},\ldots$ designates a sequence of random variables which go to [math] in probability. The sequence $F_{TW}^{1},F^{2}_{TW},\ldots$ designates on the other hand a sequence of random variables which are all identically distributed according to Tracy-Widom.

We eliminated the term $\frac{2C\cdot k}{\sqrt{n}}$ by increasing the coefficient in front of $\max\left\{\frac{k\log(n)\sqrt{j}}{\sqrt{n}},\frac{k\log(n)\sqrt{k}}{j^{1/4}}\right\}$.

Furthermore, from inequality 4.12, we obtain the important conclusion that in order to show that the re-scaled optimal alignment score with $k=n^{\alpha}$ imposed gaps goes to Tracy Widom we only need to prove prove two things:

I: $\epsilon_{n}$ goes to zero, hence

[TABLE]

as $n\rightarrow\infty$, which we solve in Section 5.2.

II: The term

[TABLE]

in the bound in 4.12 is of order less than $\frac{1}{k^{1/6}}$. This is the place where we choose $\alpha$ and $\beta$.

We have to verify that there exists a choice $\alpha>0$ and $\beta>0$ making expression 4.15 less than order $1/k^{1/6}$. In other words, we want

[TABLE]

which is satisfied for $n$ large enough as soon as

[TABLE]

This is of course possible if and only if

[TABLE]

which leads to $\alpha<\frac{1}{9}$. Thus, for any $\alpha<\frac{1}{9}$, we have convergence of the re-scaled optimal alignment score to Tracy-Widom. Hence, we are able to prove our result provided there are less than order $n^{1/9}$ gaps in the optimal alignment! Let us formulate this fact in the form of a Theorem.

Theorem 4.1**.**

Let $\alpha\in[0,1/9)$ be a constant. Let $X_{1},X_{2},\ldots$ be an i.i.d. sequence of binary variables so that

[TABLE]

and let $Y_{1},Y_{2},\ldots$ be an i.i.d. sequence of standard normal random variables. Let $S_{2}$ denote the scoring function defined by $S_{2}(c,d)=c\cdot d$ for any two letters $c,d$ of our binary-alphabet. Let $k=[n^{\alpha}]$. Then, the optimal alignment score with constrain to have $k$ gaps and properly rescaled as

[TABLE]

converges in law to Tracy-Widom.

Proof.

Multiplying (4.12) by $k^{1/6}$, we find that that

[TABLE]

holds with probability at least $1-\epsilon_{n}$, In the next section it is proven that $\epsilon_{n}$ goes to [math] as $n\rightarrow\infty$. (Recall also inequality (4.14)). Furthermore, since $\alpha<1/9$, we have from (4.16) that

[TABLE]

goes to [math] as $n\rightarrow\infty$. Hence, since also $\mathcal{E}_{k}$ go to [math] in probability, we have that

[TABLE]

goes to [math] in probability, when $n\rightarrow\infty$. Hence, since $F_{TW}^{k}$ always has Tracy-Widom distribution, it follows, that

[TABLE]

converges in Law to Tracy-Widom as $n\rightarrow\infty$. But, when we replace $k$ by $n^{\alpha}$ in expression 4.20, we obtain expression (4.17) from our theorem. So, we just finished proving that (4.17) converges to Tracy-Widom in Law. ∎

Remark. We believe that taking the variables $Y_{1},Y_{2},\ldots$ binary instead of standard normal our theorem should holds for $\alpha>0$ small enough. However, the proof for approximating the multidimensional random walk $\vec{R}^{j}$ by a multidimensional Brownian motion will be much more difficult. The intricate part is very similar to passing from the CLT to Berry-Essen type CLT theorem is when there are non-finite correlations or dimension goes to infinity at the same time as $n$ (see [6], [13]).

Now we got Tracy-Widom limiting distribution for the binary case, by using the scoring function $S_{2}$ and fixing the number of gaps. If instead we use a scoring function

[TABLE]

then we get an additional normal term on the scale $\sqrt{n}$. Indeed, $S_{0}$ and $S_{1}$ do not affect which alignment is optimal, thus,

[TABLE]

As argued before, $L_{n}^{k}(a_{0}S_{0}+a_{1}S_{1})$ depends only on the number of $a$’s and $b$’s in each of the sequences $X$ and $Y$. It is even a linear functional of these numbers of occurrences and hence asymptotically normal on the scale $\sqrt{n}$. On, the other hand, according to our theorem, we have that after rescaling $L_{n}^{k}(S_{2})$ is asymptotically Tracy-Widom on the scale $n^{1/2-\alpha/6}$, (provided $\alpha\in[0,1/6)$ and $k=n^{\alpha}$). Thus, there is a larger normal part and Tracy-Widom on a lower scale. But, for optimal alignment of DNA sequences, one should always set $a_{0}=a_{1}=0$, since the part $a_{0}S_{0}+a_{1}S_{1}$ does not influence which alignment is optimal. In that case when $a_{0}=a_{1}=0$, the limiting distribution is purely Tracy-Widom.

*The main question to ask here is what happens if we allow the number of gaps to be of order $O(n)$ which we leave as an open problem here. *

5. Technical results

5.1. How to couple normal random vectors with close

covariance matrix

The main coupling result is contained in the following Lemma.

Lemma 5.1**.**

Let $\vec{v}$ be a $d$-dimensional normal random vector with expectation [math]. Let $\epsilon>0$ and $j>0$. Assume that

[TABLE]

Then there exists a $d$-dimensional normal random vector $\vec{N}$ with covariance matrix $jI$ and expectation [math] so that

[TABLE]

In particular all the components of $\vec{v}-\vec{N}$ have variance at most $\epsilon$.

Proof.

The main idea is to construct a normal vector simply by rescaling the vector $\vec{v}$. The construction uses the fact that if $A$ is a linear transformation from $\mathbb{R}^{d}$ into itself, then $A\vec{v}$ is again a normal random variable of covariance matrix $A\,COV(\vec{v})A^{\prime}$. In this framework we take $\Sigma=COV(\vec{v})$ and

[TABLE]

Clearly, we now have

[TABLE]

Furthermore,

[TABLE]

Taking the eigenvalues of $COV(\vec{v})$ to be $\lambda_{1},\lambda_{2},\dots,\lambda_{d}$, the condition (5.1) reads as

[TABLE]

On the other hand since $\lambda_{i}\geq 0$, we also know have that $|\sqrt{\lambda_{i}}-\sqrt{j}|\leq|\sqrt{\lambda_{i}}+\sqrt{j}|$. Thus the conclusion is reached from the fact that the eigenvalues of $\Sigma^{1/2}$ are $\lambda_{i}^{1/2}$ and the following chain of inequalities

[TABLE]

∎

5.2. Proof that the events $A^{n}$, $D^{n}$, $F^{n}$ and $G^{n}$ all

have high probability

Lemma 5.2**.**

We have that $A^{n}$ has probability close to one. More precisely, assuming $\beta>\alpha$, then for any constant $\gamma<\beta-\alpha$, the following holds

[TABLE]

Proof.

Consider the sum on the right side of (4.2). Instead of taking the full sum, we are going to split the sum in the following form

[TABLE]

where for simplicity we assume that $j$ is a lot larger than $k$ and that $j$ is a multiple of $k$.

Since for each $1\leq l\leq k+1$, the sequence

[TABLE]

are iid when $i$ runs from $1$ to $j/k$. The covariance matrix of these vectors is the $k+1$-dimensional identity matrix $I$. So,

[TABLE]

is an estimate of the covariance matrix of $\vec{X}_{1,1}$. Hence, expression 5.2 is approximately equal to $I$ and

[TABLE]

Now, according to [28] the precision of such an approximation in operator norm is the square root of the dimension times the square root of the sample size. More precisely, we have

[TABLE]

Here, $C$ is a constant which does not depend on $j$. Hence, for the full sum, this means

[TABLE]

So, we have proven that

[TABLE]

Hence, since $P((A^{n}_{i})^{c})=P((A^{n}_{1})^{c})$ and since

[TABLE]

we have

[TABLE]

Since, $j=n^{\beta}$ and $k=n^{\alpha}$ with $\alpha>\beta$ it follows that the dominant term in the bound on the right side of 5.4 is $e^{-j/(k+1)}$ which is an exponential decaying function in a fractional power of $n$.

∎

The next lemma shows that when the event $A^{n}$ holds, then $D^{n}$ has high probability

Lemma 5.3**.**

We have that

[TABLE]

In particular, $P((D^{n})^{c}|A^{n})\xrightarrow[n\to\infty]{}0$.

Proof.

We work conditional on $X=x$. We have that when the event $A^{n}$ holds, then

[TABLE]

where $I$ is a $k+1$ dimensional identity matrix. So, our coupling of $\vec{R}((i+1)j)-\vec{R}(ij)$ to $\vec{C}((i+1)j)-\vec{C}(ij))$ is such that the covariance matrix of the difference has norm at most $Ck\sqrt{j}$. Thus, conditional on $X=x$, and assuming $A^{n}$ holds, we get

[TABLE]

for all $i=1,2,\ldots,\frac{n}{j}$. This implies that the $l$-th component of the difference

[TABLE]

has a variance less or equal to $Ck\sqrt{j}$. But conditional on $X=x$, we have that the differences (5.5) are normal, independent with expectation $\vec{0}$. Taking the $l$-th component of

[TABLE]

at the times $j,2j,3j,\dots,n$ is thus like taking a standard Brownian motion at the times $t_{1}<t_{2}<\ldots<t_{n}$ where

[TABLE]

Since the variance increases at each step by at most $Ck\sqrt{j}$, we get that

[TABLE]

So, the maximum of the $l$-th component of the difference $\vec{R}(t)-\vec{C}(t)$ restricted to the times $t=j,2j,3j,\ldots,\frac{n}{j}$ is thus stochastically bounded by the maximum of a standard Brownian motion during the time interval

[TABLE]

The standard Brownian motion starting at the origin. By the reflection principle for Brownian motion, this distribution function can be bounded by the distribution function of the standard Brownian motion at time

[TABLE]

We have just shown that conditional on $X=x$ and when $A^{n}$ holds, then

[TABLE]

where the last inequality was obtained by using the lemma 5.6 below. We have just proven that

[TABLE]

∎

Lemma 5.4**.**

We have that $P((F^{n})^{c})\leq e^{-\log^{2}(n)/3}$ for all $n$ large enough.

Proof.

Let $F^{n}_{li}$ be the event that:

[TABLE]

So, we have

[TABLE]

where $B_{t}$ designates a standard Brownian motion starting at the origin in a one dimensional space. By the reflection principle, we get that the right side of 5.6 is equal to

[TABLE]

where $\mathcal{N}(0,1)$ designates a standard normal. Using Lemma 5.6, we find that the right side of 5.7 can be bounded by $e^{-\log^{2}(n)/2}/2$ so that

[TABLE]

Now

[TABLE]

where in the last interesection above $l$ ranges over $1,2,\ldots,k+1$ and $i$ ranges over $1,2,\ldots,n/j$. Hence, we get

[TABLE]

where the dominant term is $e^{-\log^{2}(n)/2}$ which beats any power of $n$. Note that $\frac{n(k+1)}{j}$ is just a power of $n$. ∎

Lemma 5.5**.**

We have that $P((G^{n})^{c})\leq e^{-\log^{2}(n)/3}$ for all $n$ large enough.

Proof.

Same as for $F^{n}$. So we leave it to the reader.

∎

Lemma 5.6**.**

Let $\mathcal{N}(0,1)$ be a standard normal. Let $s>0$. Then, we have

[TABLE]

6. Simulations

In the pictures below we have simulated the optimal alignment with length $n$ running from $1$ to $10000$ and each time we compute the estimated standard deviation (of 100 simulations from a $0,1$ sequence). The gaps are allowed only in one sequence and the number of gaps is restricted to be either a fixed number, a number $k=n^{\alpha}$ or a proportional number to $n$.

7. Acknowledgements

We want to thank the reviewers of this paper for the comments which improved the paper.