3 Application 1: local score
We propose in this part to apply our results to the computation of exact p-values for local score. We first recall the definition of the local score of one sequence (section 3.1) and design a FMCI allowing to compute p-value in the particular case of an integer and i.i.d. score (section 3.2). We explain in sections 3.5 and 3.6 how to relax these two restrictive assumptions to consider rational or Markovian scores. The main result of this part is given in section 3.4 where we propose an algorithm improving the simple application of the general ones by using a specific asymptotic behaviour presented in section 3.3. As numerical application, we propose finally in section 3.7 to find significant hydrophobic segments in the Swissprot database using the Kyte-Doolittle hydrophobic scale. Our exact results are compared to the classical Gumble asymptotic approximations and discussed both in terms of numerical performance and reliability.

3.1 Definition
We consider S = S1,...,Sn a sequence of real scores and we define the local score Hn of this sequence by
which is exactly the highest partial sum score of a subsequence of S.
This local score can be computed in O(n) using the auxiliary process
U0 = 0     and for 1 ≤ j ≤ n      = max{0, Uj-1 + Sj}     (10)
because we then have Hn = maxj Uj.
Assuming the sequence S is random (Bernoulli or Markov model), we want to compute p-values relative to the event En = {Hn ≥ a} where a > 0.

3.2 Integer score
In order to simplify, we will first consider the case of integer scores (and hence a ∈ ≁) then we will extend the result to the case of rational scores.
In the Bernoulli case, [12] introduced the FMCI Z defined by
(resulting with a sequence of length n + 1) with 0 as the only starting state and a as the final absorbing state. The transition matrix ∏ is given by
where
p(i) = ℙ(S1 = i)     f(i) = ℙ(S1 ≤ i)     g(i) = ℙ(S1 ≥ i)     ∀i ∈ ℤ     (13)
It is possible to apply to this case the general algorithm 1 with L = a + 1 and k = 1 (please note that we have added Z0 to the sequence and n must then be replaced by n + 1 in the algorithm to get correct computations) to compute the p-value we are looking for. In the worst case, R has ζ = a2 non zero terms and the resulting complexity is O(a2) in memory and O(n × a2) in times. But in most cases, S1 support is reduced to a small number of values and the complexities decrease accordingly.

3.3 Asymptotic development
Is it possible to compute this p-value faster ? In the case where R admits a diagonal form, simple linear algebra could help to cut off the computations and answer yes to this question.
Proposition 4. If R admits a diagonal form we have
where []1 denotes the first component of a vector, with R∞ = limi→∞ Ri/λi, where 0 <λ < 1 is the largest eigenvalue of R and ν is the magnitude of the second largest eigenvalue. We also have v = [g(a),...,g(1)]'.
Proof. By using the corollary 15 (appendix A) we know that
Ri - λiR∞ = O(νi)     (15)
uniformly in i so we finally get for all α
uniformly for all n ≥ α and the proposition is then proved by considering the first component of equation (16).     □
Corollary 5. We have
and
Proof. Simply replace the terms in (17) and (18) with equation (14) to get the results.     □

3.4 Algorithm
The simplest way to compute ℙ(Hn ≥ a) is to use the algorithm 2 in our particular case. As the number of non zero terms in R is then a2, the resulting complexity is O(n × a2). Using the proposition 4, it possible to get the same result a bit faster on very long sequence by computing the first two largest eigenvalues magnitudes λ and ν (complexity in O(a2) with Arnoldi algorithms) and to use them to compute a p-value.
As the absolute error is in O(να) we obtain a require ε error level using a α proportional to log(ε)/log(ν) which results in a final complexity in O(log(ε)/log(ν) × a2). Unfortunately, this last method requires to use delicate linear algebra techniques and is therefore more difficult to implement. Another better possibility is to use the corollary 5 to get the following fast and easy to implement algorithm:
algorithm 3: local score p-value
x a real column vector of size a, (pi)i≥1 and (λi)i≥3 to sequences of real and i an integer
initialization x = [g(a),...,g(1)]', p1 = g(a), and i = 0
main loop while (i <n and (λi) has not yet converged towards λ)
    • i = i + 1
    • x = R × x (sparse product)
    • pi = pi-1 + x1
    • λi = (pi - pi-1)/(pi-1 - pi-2) (if defined)
end • p = pi
    • if (i <n) then p = p + (pi - pi-1)
    • return p
At any step i of the main loop we have pi = ℙ(Hi ≥ a) and the final value taken by i is the α of proposition 4. One should note that only the last three terms of (pi)i≥1 and (for a simple convergence testing) the last two terms of (λi)i≥3 are required by the algorithm.

3.5 Rational scores
What if we consider now a rational score instead of an integer one ? If we denote by  ⊂ ℚ the support of S1, let us define M = mini∈ℕ{i ⊂ ℤ}. Changing the scale of the problem by the factor M allows us to get back to the integer case:
ℙ(Hn ≥ a) = ℙ(M Hn ≥ M a)     (19)
This scale factor will obviously increase the complexity of the problem, but as the support cardinal (denoted η) is not changed during the process, the resulting complexities are O(M × a × η) in memory and O(M × n × a × η) in time (n could vanish from the time complexity thanks to the faster algorithm presented above).
For example, if we consider the Kyte-Doolittle hydrophobicity score of the amino-acids (see [10] and table 1), it takes only η = 20 values and M = 10, the resulting complexity to compute ℙ(Hn ≥ a) is then O(200 × n × a). If we consider now the more refined Chothia score ([4]), the scale factor increases from M = 10 to M = 100 and the resulting complexities are multiplied by 10.
Table 1  Distribution of amino-acids estimated on Swissprot (release 47.8) database and Kyte-Doolittle hydrophobic scale. Mean score is -0.244.
a. a.  F  M  I  L  V  C  W  A  T  G
ℙ in %  4.0  2.4  5.9  9.6  6.7  1.5  1.2  7.9  5.4  6.9
score  2.8  1.9  4.5  3.8  4.2  2.5  -0.9  1.8  -0.7  -0.4
a. a.  S  P  Y  H  Q  N  E  K  D  R
ℙ in %  6.9  4.8  3.1  2.3  3.9  4.2  6.6  5.9  5.3  5.4
score  -0.8  -1.6  -1.3  -3.2  -3.5  -3.5  -3.5  -3.9  -3.5  -4.5

3.6 Markov case
All these results can be extended to the Markov case but this require to define a new FMCI allowing us to trace the last score (in the case of an order one Markov chain for the sequence S, if a higher order m is considered, we just have to add the corresponding number of preceding scores to Z instead of one):
Doing this now we get k = η (the cardinal of the score support) starting states instead of one so we need a starting distribution μ (which could be a Dirac) to compute the p-value.
We will not detail here the structure of the corresponding sparse transition matrix ∏ (see [13]) but we need to know its number ζ of non zero terms. If a is an integer value (we suppose here that the scale factor has been already included in it) then the order of R is M × a × η and ζ = O(M × a × η2) (and we get O(M × a × ηm+1) when an order m Markov model is considered).

3.7 Numerical results
In this section, we apply the results presented above to a practical local score study. We consider the complete protein database of Swissprot release 47.8 and the classical amino acid hydrophobic scale of Kyte-Doolittle given in table 1 ([10]). The database contains roughly 200 000 sequences of various lengths (empiric distribution given in figure 1).
Figure 1  Empiric distribution of Swissprot (release 47.8) protein lengths. In order to improve readability, 0.5% of sequences with length ∈ [2 000, 9 000] have been removed from this histogram. Once the best scoring segment has been determined for each of these sequences, we need to compute the corresponding p-values. According to [9], the asymptotic distribution of Hn is given (if mean score is < 0, which is precisely the case here) by the following conservative approximation:
ℙ(Hn ≥ a) ≃ 1 - exp (-nKe-aλ)     (21)
where constants λ and K depend on the scoring distribution.
With our hydrophobic scale and a distribution of amino-acids estimated on the entire database we get
λ = 5.144775 × 10-3     and     K = 1.614858 × 10-2
(computation performed with a C function implemented by Altschul). Once the constants are computed we could get all the approximated p-values very quickly (a few seconds for the 200 000 p-values).
On the other hand, our new algorithm allows to compute (for the very first time) the exact p-values for this example. As the chosen scoring function has a one digit precision level, we need to use a scale factor of M = 10 to fall back to the integer case. A C++ implementation (available on request) performed all the computations in roughly three hours on a Pentium 4 CPU 2.8 GHz (this means approximately 20 p-values computed by second).
We can see on figure 2 the comparison between exact values and Karlin's approximations. The conservative design of the approximations seems to be successful except for very short unsignificant sequences. While the approximations are rather close to perfection for sequences with more than 2 000 amino-acids, the smaller the sequence is, the worse the approximations get. This is obviously consistent with the asymptotic nature of Karlin's formula but seems to indicate that these approximations are not reliable for 99.5% of the sequence in the database (protein of length < 2 000).
Figure 2  Exact p-value against Karlin ones (in log scale). Color refers to a range of sequence lengths: smaller than 100 in black (≃ 20 000 sequences), between 100 and 200 in red (≃ 40 000 sequences), between 200 and 500 in orange (≃ 90 000 sequences), between 500 and 1000 in yellow (≃ 30 000 sequences), between 1000 and 2 000 in blue (≃ 6 000 sequences) and greater than 2 000 in green (≃ 1 000 sequences). The solid line represents y = x. Range have been chosen for readability and few dots with exact p-value smaller than 10-30 are hence missing. One should object that it exists ([1,2]) a well known finite size correction to formula (21) that might be useful, especially when considering short sequences. Unfortunately in our case, this correction does not seems to improve the quality of the approximations (data not shown) and we hence make the choice to ignore it.
In table 2 we compare the number of sequences predicted to have a significant hydrophobic segment at a certain e-value level by the two approaches. If the Karlin's approximations are used, many proteins are considered unsignificant while they are. For example, with the classical database threshold of 10-5, only few sequences (6%) are correctly identified by Karlin's approximations.
Table 2  Number of e-value smaller than a threshold are given for exact computations (exact) and asymptotic Karlin's approximations (Karlin). The last row gives the accuracy of asymptotic predictions (accuracy = Karlin/exact).
e-value  10-1  10-2  10-3  10-4  10-5  10-6
exact  9473  7772  6271  4563  3232  2348
Karlin  3417  2047  1056  439  195  96
accuracy  34%  26%  17%  10%  6%  4% We have seen that Karlin's approximations are often far too conservative to give accurate results, but what about the ranking ? Table 3 proposes the Kendall's tau rank correlation (see [16] chapter 14.6 for more details) which is equal to 1.0 for a complete rank agreement and equal to -1.0 for a complete inverse rank agreement. As we will certainly be interested in the most significant sequences produced by our study, we compute our Kendall's tau only on these sequences. When all sequence lengths are considered, Karlin's approximations show their total irrelevance to give correct ranking for the first 10 or 50 most significant p-values. Even when the 100 first p-values are taken into account, relative ranks given by Karlin's approximations are wrong in 63% of the cases, which is huge. However, in the case where the approximations values are close to the exact ones (sequence lengths greater than 2 000, which correspond only to 0.5% of the database), p-values obtained with both methods are highly correlated.
Table 3  Kendall's tau (rank correlation) comparing the most significant exact p-values (the reference) to the Karlin's approximations. The column "all" gives the result for all sequences while the Ri give the results for a certain range of sequence lengths: smaller than 100 for R1, between 100 and 200 for R2, between 200 and 500 for R3, between 500 and 1 000 for R4, between 1 000 and 2 000 for R5 and greater than 2 000 for R6.
number of p-values  all  R1  R2  R3  R4  R5  R6
10  0.30  0.64  0.24  -0.20  0.58  0.64  0.97
50  0.14  0.73  0.50  0.46  0.56  0.78  0.97
100  0.37  0.70  0.67  0.62  0.61  0.80  0.98