ESTIMATING THE GUMBEL SCALE PARAMETER FOR LOCAL ALIGNMENT OF RANDOM SEQUENCES BY IMPORTANCE SAMPLING WITH STOPPING TIMES

Abstract

The gapped local alignment score of two random sequences follows a Gumbel distribution. If computers could estimate the parameters of the Gumbel distribution within one second, the use of arbitrary alignment scoring schemes could increase the sensitivity of searching biological sequence databases over the web. Accordingly, this article gives a novel equation for the scale parameter of the relevant Gumbel distribution. We speculate that the equation is exact, although present numerical evidence is limited. The equation involves ascending ladder variates in the global alignment of random sequences. In global alignment simulations, the ladder variates yield stopping times specifying random sequence lengths. Because of the random lengths, and because our trial distribution for importance sampling occurs on a different sample space from our target distribution, our study led to a mapping theorem, which led naturally in turn to an efficient dynamic programming algorithm for the importance sampling weights. Numerical studies using several popular alignment scoring schemes then examined the efficiency and accuracy of the resulting simulations.

Fig. 1

Gapped global alignment scores and the corresponding directed paths for two subsequences A[1, 10] = TACTAGCGCA and B[1, 9] = ACGGTAGAT, drawn from the nucleotide alphabet {A, C, G, T}. Figure 1 uses a nucleotide scoring matrix, where s(a, b) = 5 if a = b and –4 otherwise, and the affine gap penalty w_g = 3 + 2g. The vertex (i, j) is in the northeast corner of the cell (i, j), with the origin (0, 0) at the southwest corner of Figure 1. The cell (i, j) displays the global score S_{i, j}, calculated from (2.2). The optimal global path ending at the point (10, 8), for example, consists of 12 edges, in order: 1 east of length 1, 2 northeast, 1 north of length 2, 3 northeast, 1 east of length 3, and 1 northeast. The optimal global score S_{10, 8} = –5 + 5 + 5 – 7 + 5 + 5 + 5 – 9 + 5 = 9 is the sum of the corresponding edges and represents the path of greatest weight starting at (0, 0) and ending at (10, 8). The corresponding optimal global alignment of the subsequences A[1, 10] and B[1, 9] is $$\begin{array}{cc}\hfill & \mathrm{TAC}--\mathrm{TAGCGCA}\hfill \\ \hfill & -\mathrm{ACGGTAG}---\mathrm{A}.\hfill \end{array}$$ The edge maxima are M₁ = –4, M₂ = 0, M₃ = 5, M₄ = 1, M₅ = 3, M₆ = 8, M₇ = 13, M₈ = 9, M₉ = 6. The shading and the double lines indicate squares where a vertex (surrounded by double lines) generated an SALE β(k). The SALE scores are M_β(1) = M₃ = 5, M_β(2) = M₆ = 8, M_β(3) = M₇ = 13; and the global maximum M for A and B is no less than 13, the largest global score shown.

Fig. 2

Two examples of alignment path π generated by a Markov chain. As in Figure 1, the shading and the double lines indicate squares where a vertex (surrounded by double lines) generated an SALE. The SALEs determine the stopping time N = β(3). In Figure 2, the first SALE is determined by the score at the vertex (3, 3); the second SALE, the vertex (7, 6); the third SALE, the vertex (9, 10). Therefore, N = β(3) = 10. The vertical ray I_{N, N} and the horizontal ray D_{N, N} are indicated by double circles. The lower path π (solid line) ends at (N + 2, N) with a final transition to S; the upper path π (long-dashed line), at (N, N + 4) with a final transition to D. The closed vertices indicate intersection with the square corresponding to ${\omega }^{\prime }=({\omega }_{\mathbf{A}}^{\prime },{\omega }_{\mathbf{B}}^{\prime })=(\mathbf{A}[1,N],\mathbf{B}[1,N])$.

Fig. 3

Plot of relative errors against computation time (sec). Both axes are in logarithmic scale. Computation time was measured on a 2.99 GHz Pentium® DCPU. Relative errors for BLOSUM45 with Δ(g) = 14 + 2g are shown by ■; BLOSUM62 with Δ(g) 11 g, by ◆; BLOSUM80 with Δ(g) = 10 + g, by ▲; PAM70 with Δ(g) = 10 + g, by •; PAM30 with Δ(g) = 9 + g, by ○.