Exact Asymptotic Results for a Model of Sequence Alignment

Finding analytically the statistics of the longest common subsequence (LCS) of a pair of random sequences drawn from c alphabets is a challenging problem in computational evolutionary biology. We present exact asymptotic results for the distribution of the LCS in a simpler, yet nontrivial, variant of the original model called the Bernoulli matching (BM) model which reduces to the original model in the large c limit. We show that in the BM model, for all c, the distribution of the asymptotic length of the LCS, suitably scaled, is identical to the Tracy-Widom distribution of the largest eigenvalue of a random matrix whose entries are drawn from a Gaussian unitary ensemble. In particular, in the large c limit, this provides an exact expression for the asymptotic length distribution in the original LCS problem.

Finding analytically the statistics of the longest common subsequence (LCS) of a pair of random sequences drawn from c alphabets is a challenging problem in computational evolutionary biology. We present exact asymptotic results for the distribution of the LCS in a simpler, yet nontrivial, variant of the original model called the Bernoulli matching (BM) model which reduces to the original model in the c → ∞ limit. We show that in the BM model, for all c, the distribution of the asymptotic length of the LCS, suitably scaled, is identical to the Tracy-Widom distribution of the largest eigenvalue of a random matrix whose entries are drawn from a Gaussian unitary ensemble. In particular, in the c → ∞ limit, this provides an exact expression for the asymptotic length distribution in the original LCS problem. Sequence alignment is one of the most useful quantitative methods used in evolutionary molecular biology [1][2][3]. The goal of an alignment algorithm is to search for similarities in patterns in different sequences. A classic and much studied alignment problem is the so called 'longest common subsequence' (LCS) problem. The input to this problem is a pair of sequences α = {α 1 , α 2 , . . . , α i } (of length i) and β = {β 1 , β 2 , . . . , β j } (of length j). For example, α and β can be two random sequences of the 4 base pairs A, C, G, T of a DNA molecule, e.g., α = {A, C, G, C, T, A, C} and β = {C, T, G, A, C}. A subsequence of α is an ordered sublist of α (entries of which need not be consecutive in α), e.g, {C, G, T, C}, but not {T, G, C}. A common subsequence of two sequences α and β is a subsequence of both of them. For example, the subsequence {C, G, A, C} is a common subsequence of both α and β. There can be many possible common subsequences of a pair of sequences. The aim of the LCS problem is to find the longest of such common subsequences. This problem and its variants have been widely studied in biology [4][5][6][7], computer science [8][9][10]2], probability theory [11][12][13][14][15][16] and more recently in statistical physics [17][18][19]. A particularly important application of the LCS problem is to quantify the closeness between two DNA sequences. In evolutionary biology, the genes responsible for building specific proteins evolve with time and by finding the LCS of the same gene in different species, one can learn what has been conserved in time. Also, when a new DNA molecule is sequenced in vitro, it is important to know whether it is really new or it already exists. This is achieved quantitatively by measuring the LCS of the new molecule with another existing already in the database.
For a pair of fixed sequences of length i and j respectively, the length L i,j of their LCS is just a number. However, in the stochastic version of the LCS problem one compares two random sequences drawn from c alphabets and hence the length L i,j is a random variable.
A major challenge over the last three decades has been to determine the statistics of L i,j [11][12][13][14][15]. For equally long sequences (i = j = n), it has been proved that L n,n ≈ γ c n for n ≫ 1, where the averaging is performed over all realizations of the random sequences. The constant γ c is known as the Chvátal-Sankoff constant which, to date, remains undetermined though there exists several bounds [12,14,15], a conjecture due to Steele [13] that γ c = 2/(1 + √ c) and a recent proof [16] that γ c → 2/ √ c as c → ∞. Unfortunately, no exact results are available for the finite size corrections to the leading behavior of the average L n,n , for the variance, and also for the full probability distribution of L n,n . Thus, despite tremendous analytical and numerical efforts, exact solution of the random LCS problem has, so far, remained elusive. Therefore it is important to find other variants of this LCS problem that may be analytically tractable.
Computationally, the easiest way to determine the length L i,j of the LCS of two arbitrary sequences of lengths i and j (in polynomial time ∼ O(ij)) is via using the recursive algorithm [2,19] subject to the initial conditions L i,0 = L 0,j = L 0,0 = 0. The variable η i,j is either 1 when the characters at the positions i (in the sequence α) and j (in the sequence β) match each other, or 0 if they do not. Note that the variables η i,j 's are not independent of each other. To see this consider the simple example -matching of two strings α = AB and β = AA. One has by definition: η 1,1 = η 1,2 = 1 and η 2,1 = 0. The knowledge of these three variables is sufficient to predict that the last two letters will not match, i.e., η 2,2 = 0. Thus, η 2,2 can not take its value independently of η 1,1 , η 1,2 , η 2,1 . These residual correlations between the η i,j variables make the LCS problem rather complicated. Note however that for two random sequences drawn from c alphabets, these cor-relations between the η i,j variables vanish in the c → ∞ limit. A simpler but natural variant of this LCS problem is the Bernoulli matching (BM) model where one ignores the correlations between η i,j 's for all c [19]. The BM model reduces to the original LCS problem only in the c → ∞ limit. The length L BM i,j of the BM model satisfies the same recursion relation in Eq. (1) except that η i,j 's are now independent and each drawn from the bimodal distribution: p(η) = (1/c)δ η,1 + (1 − 1/c)δ η,0 . The BM model, though simpler than the original LCS problem, is still nontrivial due to the nonlinear recursion relation in Eq. (1). Using the cavity method of spin glass physics [20], the asymptotic behavior of the average length in the BM model was determined analytically [19], , same as the conjectured value of the Chvátal-Sankoff constant γ c for the original LCS model. However, other properties such as the variance or the distribution of L BM n,n remained untractable even in the BM model.
The purpose of this Letter is to present an exact asymptotic formula for the distribution of the length L BM n,n in the BM model for all c. Our main result is that for large n, where χ is a random variable with a n-independent distribution, Prob(χ ≤ x) = F TW (x) which is the well studied Tracy-Widom distribution for the largest eigenvalue of a random matrix with entries drawn from a Gaussian unitary ensemble [21]. For a detailed form of the function F TW (x), see [21]. We show that for all c, This allows us to calculate the average including the subleading finite size correction term and the variance of L BM n,n for large n, where one can use the known exact values [21], χ = −1.7711 . . . and χ 2 − χ 2 = 0.8132 . . .. In particular, we note that in the limit c → ∞, Eqs. In the BM model, the length L BM i,j can be interpreted as the height of a surface over the 2-d (i, j) plane constructed via the recursion relation in Eq. (1). A typical surface, shown in Fig. (1a), has terrace-like structures. It is useful to consider the projection of the level lines separating the adjacent terraces whose heights differ by 1 (see Fig.2) onto the 2-d (i, j) plane. Note that, by the rule Eq. (1), these level lines never overlap each other, i.e., no two paths have any common edge. The statistical weight of such a projected 2-d configuration is the product of weights associated with the vertices of the 2d plane. There are five types of possible vertices with nonzero weights as shown in Fig.2, where p = 1/c and q = 1 − p. Since the level lines never cross each other, the weight of the first vertex in Fig. (2) is 0.  Fig.(1a). The adjacent table shows the weights of all vertices on the 2-d plane.
Consider first the limit c → ∞ (i.e., p → 0). The weights of all allowed vertices are 1, except the ones shown by black dots in Fig.(2), whose associated weights are p → 0. The number N of these black dots inside a rectangle of area A = ij can be easily estimated. For large A and p → 0, this number is Poisson distributed with the mean N = pA. The Bethe ansatz analysis shows that BM corresponds to the sector of the 5-vertex model [22] where the density α of empty edges in a row of vertical edges is close to the boundary α ≈ 1 − . The careful examination of the free energy near this boundary allows one to conclude that the leading contribution in p (for p → 0) to N comes exactly from the line of phase transitions in a 5-vertex model. The subleading corrections to N are of order ∼ p 3/2 and are ensured by small deviations from the critical line being beyond the Poisson approximation [23].
The height L BM i,j is just the number of level lines N inside this rectangle of area A = ij. The problem of estimating N has recently appeared in a number of interface models such as a polynuclear growth model [24] and a ballistic deposition model [25]. By using a mapping to the longest increasing subsequence (LIS) of the equally likely permutations of a set of integers and then, by applying a celebrated result due to Baik, Deift and Johansson (BDJ) [26], it was shown [24,25] that the number of level lines N inside the rectangle (for large A), appropriately scaled, has a limiting behavior, where χ is a random variable with Tracy-Widom distribution. Using N = pA = ij/c, one then obtains in the limit p → 0, In particular, for large equal length sequences i = j = n, we get for c → ∞ Note that since the BM and the original LCS model are equivalent in the limit c → ∞, the exact results in Eqs. (6)-(7) also hold for the LCS model. Note that only the leading behavior of the average L n,n was known before [16] in the c → ∞ limit of the original LCS model. For finite c, while the above mapping to the LIS problem still works, the corresponding permutations of the LIS problem are not generated with equal probability and hence one can no longer use the BDJ results. To make progress for finite c, we map the BM model exactly to a 3-d anisotropic directed percolation (ADP) model first considered by Rajesh and Dhar [27]. This ADP model can further be mapped to a (1 + 1)-d directed polymer problem studied by Johansson [28]. For this specific directed polymer problem, Johansson derived exact asymptotic result for the distribution of the polymer energy. Translating these results back to the BM model, we derive our main results in Eqs. (3)- (5). Note that the recursion relation in Eq. (1) can also be viewed as a (1 + 1)-d directed polymer problem [18,19] and some asymptotic results (such as the O(n 2/3 ) behavior of the variance of L n,n for large n) can be obtained using the arguments of universality [18]. However, this does not provide precise results for the full distribution which are obtained here.
Let us consider a directed bond percolation on a simple cubic lattice. The bonds are occupied with probabilities p x , p y , and p z along the x, y and z axes and are all directed towards increasing coordinates. Imagine a source of fluid at the origin which spreads along the occupied directed bonds. The sites that get wet by the fluid form a 3-d cluster. In the ADP problem, the bond occupation probabilities are anisotropic, p x = p y = 1 (all bonds aligned along the x and y axes are occupied) and p z = p. Hence, if the point (x, y, z) gets wet by the fluid then all the points (x ′ , y ′ , z) on the same plane with x ′ ≥ x and y ′ ≥ y also get wet. Such a wet cluster is compact and can be characterized by its bounding surface height h(x, y) as shown in Fig.(1b). It is not difficult to see that the height h(x, y) satisfies the following recursion relation [27], where ξ i,j 's are i.i.d. random variables taking nonnegative integer values with Prob(ξ i,j = k) = (1 − p) p k for k = 0, 1, 2, . . .. One can also interpret the height h(x, y) in Eq. (8) as the energy of a directed polymer in the (x − y) plane. Precisely this particular version of the polymer problem was studied by Johansson [28] who obtained the asymptotic distribution of the height for large x and y, where q = 1 − p, χ is a random variable with a Tracy-Widom distribution. While the terrace-like structures of the ADP surface look similar to the BM surfaces (compare Figs.(1a) and (1b)), there is an important difference between the two. In the ADP model, the level lines separating two adjacent terraces can overlap with each other [27], which does not happen in the BM model. However, by making the following change of coordinates in the ADP model [27] ζ = x + h(x, y); η = y + h(x, y) one gets a configuration of the surface where the level lines no longer overlap. Moreover, it is not difficult to show that the projected 2-d configuration of level lines of this shifted ADP surface has exactly the same statistical weight as the projected 2-d configuration of the BM surface. Denoting the BM height byh(x, y) = L BM x,y , one then has the identity,h(ζ, η) = h(x, y), which holds for each configuration. Using Eq. (10), one can rewrite this identity as Thus, for any given height function h(x, y) of the ADP model, one can, in principle, obtain the corresponding height functionh(x, y) for all (x, y) of the BM model by solving the nonlinear equation (11). This is however very difficult in practice. Fortunately, one can make progress for large (x, y) where one can replace the integer valued discrete heights by continuous functions h(x, y) and h (x, y). Using the notation ∂ x ≡ ∂/∂x it is easy to derive from Eq. (10) the following pair of identities, In a similar way, one can show that We then observe that Eqs. (12) and (13) are invariant under the simultaneous transformations Since the height is built up by integrating the derivatives, this leads to a simple result for large ζ and η, χ, (16) where p = 1/c and q = 1 − 1/c. For equal length sequences x = y = n, Eq. (16) then reduces to Eq. (3).
Thus the discrepancy falls off as a power law for large ζ, indicating that indeed our solution is asymptotically exact. We have also performed numerical simulations of the BM model using the recursion relation in Eq. (1) for c = 2, 4, 9, 16, 100. Our preliminary results [23] for relatively small system sizes (up to n = 5000) are consistent with our exact results in Eqs. (3)-(5). The Tracy-Widom distribution of the random matrix theory has appeared recently in a number of problems [21,29,28,24,25]. In this Letter, we have shown that it also describes the asymptotic distribution of the length of the longest common subsequence in a sequence matching problem. While a possible link between the two problems was speculated before [29], a precise connection, so far, was missing and is provided here.