Pattern statistics on Markov chains and sensitivity to parameter estimation

In order to study pattern occurrences in biological sequences, simple frequencies are not relevant in most cases because of pattern overlapping structure as well as composition bias in the sequences. A common workaround consists to compute the significance of an observation assuming the sequence X = X₁ ... X_ℓ over the finite alphabet $\mathcal{A}$. (size k) is generated according to an order m ≥ 1 homogeneous, stationary and ergodic Markov model. Let π (size k^m+1) defined by

π(w, a) = ℙ(X_m+1= a|X₁...X _m = w) ∀(w, a) ∈ $\mathcal{A}$ ^m × $\mathcal{A}$     (1)

be the parameter of this Markov model, Π its transition matrix (note that we have Π = π only if m = 1) and μ its stationary distribution (defined by μ × Π = μ).

We then introduce the pattern statistic defined by

$S=\{\begin{array}{rr}\hfill -{\log }_{10}\mathbb{P}(N\ge N_{\text{obs}})& \hfill \text{if }{N}_{\text{obs}}\ge \mathbb{E}[N]\\ \hfill {\log }_{10}\mathbb{P}(N\le N_{\text{obs}})& \hfill \text{if }{N}_{\text{obs}}<\mathbb{E}[N]\end{array}\left(2\right)$

where N is the random number of overlapping occurrences (i. e. X = aababaaba contains three overlapping occurrences of aba but only two non-overlapping ones) of a given fixed pattern on the random sequence X and N_obs is an observation.

When π is known (and hence μ), several statistical methods are available to compute S: exact computations [1–4], Gaussian [5, 6], binomial [7, 8], compound Poisson [9–11] or large deviations approximations [12]. But in general, the parameter π is not available and must be estimated. Let us denote by N₀ (resp. N₁) the (overlap) frequencies of all words of size m (resp. m + 1) in the sequence Y = Y₁ ... Y _n , then the Maximum-Likelihood Estimator (MLE) of π is given by

$\begin{array}{cc}\widehat{\pi }(w,a)=\frac{N_1(wa)}{{\displaystyle {\sum }_{b\in \mathcal{A}}{N}_1(wb)}}& \forall (w,a)\in {\mathcal{A}}^m\times \mathcal{A}\end{array}(3)$

and the MLE of μ (as a function of π) is therefore defined by $\widehat{\mu }$ × $\widehat{\Pi }$ = $\widehat{\mu }$ where $\widehat{\Pi }$ is the transition matrix associated to $\widehat{\pi }$

We introduce now the following estimators

$\begin{array}{cccc}{\mu }_N(w)=\frac{N_0(w)}{n-m+1}& and& {\pi }_N(w,a)=\frac{N_1(wa)}{N_0(w)}& \forall (w,a)\in {\mathcal{A}}^m\times \end{array}\mathcal{A}(4)$

which are known to be asymptotically equivalent with the MLE when n is large.

The quality of parameter estimation depends both on the number of parameters to estimate (k^m+1for an order m Markov model) and of the length (n) of the homogeneous sequence used for their estimation. When the same sequence (or set of sequences) is used both for observed frequencies and parameter estimation, m should not be greater than h – 2 for a pattern of length h (as else, the observed frequency of the pattern will be included in the model). As literature often suggests to use the highest possible order, it is hence common to consider m = 6 or more (for a DNA pattern of size h ≥ 8). Moreover, because of the homogeneity assumption of the model, the considered genomes have often to be segmented first. As a result, the sequences length used for parameter estimations are often dramatically reduced by such segmentation (e. g. n = 10⁵ to n = 10⁶ at the very best for DNA sequences). It is hence quite common to encounter high order Markov models estimated on rather short sequences which could result in high sensitivity to parameter estimation.

Considering that Y is generated through a Markov model of parameter π, the main goal of this paper is to study the distribution of S _N , the statistic S computed using the estimators μ _N and π _N , and the consequences of its variability in projects using pattern statistics. We first present in details how the delta-method can be used to get a Gaussian approximation for the distribution of S _N (using a binomial approximation to compute the pattern statistics). Then these approximations are validated through simulations and, at last, we consider a classical pattern study (finding the most over-represented patterns of a given size) and we evaluate the detrimental effect of parameter estimations both in terms of true positive rate and rank accordance.

The delta-method allows us to approximate the distribution of $\widehat{S}$ by a Gaussian distribution. This first requires to compute the expectation and covariance matrix of frequencies and then to study the derivative of a function which is specific of the method used to compute the pattern statistics. In the case of the binomial approximations, we have found an explicit expression of σ the standard deviation of $\widehat{S}$.

It is clear that our approximation of σ using the delta-method relies one two major assumptions: 1) the distribution of N is Gaussian; 2) F⁺ is regular enough (e.g. not too steep) around E. When m grows, E closes to the boundary of the definition range of F⁺ hence degrading assumption 2. Moreover, it is well known that Gaussian approximations for word frequencies become weaker when the expected numbers of their occurrences become smaller, thus degrading assumption 1. It is therefore obvious that our approximation of σ will get less and less reliable as m grows.

However, the approximation of σ has been validated through simulations and appears to be very reliable (even for m = 5 or 6). As pattern statistics computed through binomial approximations are close to the exact statistics [8], the value of σ should not differ a lot when another statistical method is used. We have compared our approximations to the empiric distribution obtained using compound Poisson and large deviations approximations and, as expected, our approximations remains quite reliable even for these statistical methods.

The variability due to parameter estimation is of course related to the Markov model order m and to the size k of the alphabet (as we have k^m+1parameters for this model) and to the length n of the sequence used for this estimation. For example, considering an order m = 6 model with n = 4639675 (Escherichia coli K12 complete genome) requires to estimate 3 × 4⁶ = 4096 free parameters which results roughly in 400 observation per free parameter. Although this situation seems quite comfortable, we have seen with our simulations that it leads an unacceptable variability for pattern statistics.

As literature often advices to use the highest possible Markov order for a given pattern problem (which means m = h - 2 for pattern of size h) it is easy to understand that such a practice could have very detrimental effects on the computed statistics unless huge data are available for estimation purpose. Even if we consider the more reasonable attitude to choose m using the classical framework of model selection (e.g. using the Akaike Information Criterion – AIC –) we get m = 5 for Mycoplasma genitalium and m = 6 for Escherichia coli K12 hence resulting in both cases in the same catastrophic results in terms of false positive and even worse ones in terms of ranking.

Moreover, we assumed here that our model was homogeneous all along the considered sequences. This is obviously completely false when complete genomes are considered. So it is more likely that the sample size n would be far smaller than a million on classical pattern studies (even of human genomes for example). As a result, the variability we pointed out in this paper will have a considerable detrimental effect on most studies unless the Markov order is carefully set.

In order to do so, we advice to compute our approximation of σ each time a pattern statistic is produced and then to evaluate, either by simulation (like in this paper) or by a theoretical work the impact of this variability on the considered study.

We give here the expression of the covariance matrix C introduced in section "distribution of N = (N₀, N₁)". The sequence Y (of length n) is generated by an homogeneous, stationary and ergodic order m Markov model of parameter π and stationary distribution μ. We want to compute the covariance of the vector N of random frequencies of size m and m + 1 words.

For any word w (of size h _w ), we introduce the following notation for h _w ≤ i ≤ n

$I_i(w)={\mathbb{I}}_{\left\{wendinpositioni\right\}}={\mathbb{I}}_{\{Y_{i-h_w+1}^i=w\}}(45)$

where $Y_i^j$ = Y _i ... Y _j for all i ≤ j. If h _w ≥ m, we denote by

$p(w)=\mu (w_1^m)\Pi (w_1^m,w_{m+1})\mathrm{...}\Pi (w_{h_w-m}^{h-1},w_h)(46)$

the probability to see one occurrence of w at a given position in the sequence. At last, if we consider another word v (of size h _v = m) and if h _w = m, we denote by

${\Pi }_{\delta }(v,w)={\displaystyle \sum _{x\in {\mathcal{A}}^{\delta }}p(vxw)}(47)$

the probability to see occurrences of v and w separated by a gap of length δ.

For any words v an w (to simplify, we suppose that h _v ≥ h _w ) then, for all δ ∈ $\mathbb{Z}$ and

max(h _v , h _w - δ) ≤ i ≤ min(n, n - δ) we have

which do not depend on i.

It is therefore easy to show that

$\begin{array}{cc}\mathbb{E}[N(v)N(w)]={\displaystyle \sum _{i=h_v}^n{\displaystyle \sum _{\delta =h_w-i}^{n-i}{D}_{\delta }(v,w)}}& \left(49\right)\\ ={\displaystyle \sum _{\delta =h_w-n}^{n-h_v}{N}_{\delta }{D}_{\delta }(v,w)}& \left(50\right)\\ =M(v,w)+O(v,m)& \left(51\right)\end{array}$

where the main part (2n - h _v - h _w + 2 terms) is given by

$M(v,w)={\displaystyle \sum _{\delta =h_v}^{n-h_w}{N}_{-\delta }{D}_{-\delta }}(v,w)+{\displaystyle \sum _{\delta =h_w}^{n-h_v}{N}_{\delta }{D}_{\delta }}(v,w)(52)$

and the overlapping part (h _v + h _w - 1 terms) by

and with

$N_{\delta }=\{\begin{array}{ll}n-h_w+1+\delta \hfill & \delta \in [h_w-n,h_w-h_v[\hfill \\ n-h_v+1\hfill & \delta \in [h_w-h_v,0]\hfill \\ n-h_v+1-\delta \hfill & \delta \in ]0,n-h_v]\hfill \end{array}(54)$

As we have

C(v, w) = M (v, w) + O (v, w) - E (v) E (w)     (55)

the problem is hence to compute M and O for all pairs of size m or m + 1 words. In order to simplify, we will just treat here the case of a pair of size m words (other cases can be derived from this special case).

For the main part we obtain

$\begin{array}{c}M(v,w)={\displaystyle \sum _{\delta =m}^{n-m}{N}_{-\delta }\mu (w){\Pi }_{\delta -m+1}(w,v)}\\ +{\displaystyle \sum _{\delta =m}^{n-m}{N}_{\delta }\mu (v){\Pi }_{\delta -m+1}(v,w)}\end{array}(56)$

(2n - 2m + 2 terms). As P _k (v, w) quickly converges toward μ(w) when k grows (convergence speed is given by λ ^k where λ is the magnitude of the second eigenvalue of the transition matrix Π). So there exists a rank r ≥ m such as

$\begin{array}{l}M(v,w)\simeq \mu (v)\mu (w){\displaystyle \sum _{\delta =r}^{n-m}(N_{-\delta }+N_{\delta })}\\ +{\displaystyle \sum _{\delta =m}^{r-1}{N}_{-\delta }\mu (w){\Pi }_{\delta -m+1}}(w,v)\\ +{\displaystyle \sum _{\delta =m}^{r-1}{N}_{\delta }\mu (v){\Pi }_{\delta -m+1}}(v,w)(57)\end{array}$

which has only 2r - 2m + 1 terms.

And for the overlapping part we get

$\begin{array}{l}O(v,w)=N_0\times \mu (v)\times {\mathbb{I}}_{\{v=w\}}\\ +{\displaystyle \sum _{\delta =1}^{m-1}{N}_{-\delta }}\times p(w{v}_{m-\delta +1}^m)\times {\mathbb{I}}_{\left\{v_1^{m-\delta }=w_{1+\delta }^m\right\}}\\ +{\displaystyle \sum _{\delta =1}^{m-1}{N}_{\delta }\times p(v{w}_{m-\delta +1}^m)}\times {\mathbb{I}}_{\left\{v_{1+\delta }^m=w_1^{m-\delta }\right\}}(58)\end{array}$

which has 2m + 1 terms.

So the overall complexity for the computation of one term of C is hence O(r) where the value of r is directly connected to the magnitude λ of the second eigenvalue of the transition matrix.

In the particular case of an order one Markov model (m = 1), we give here the complete expressions of M and O.

For all a, b, c, d ∈ $\mathcal{A}$, we have

$\begin{array}{l}M(a,b)\simeq (n-r+1)(n-r)\mu (a)\mu (b)\\ +{\displaystyle \sum _{\delta =1}^{r-1}(n-\delta )\left(\mu (b){\Pi }^{\delta }(b,a)+\mu (a){\Pi }^{\delta }(a,b)\right)(59)}\end{array}$

O (a, b) = nμ(a) $\mathbb{I}$_{{a = b}}    (60)

$\begin{array}{l}\frac{M(ab,c)}{\Pi (a,b)}\simeq (n-r)(n-r-1)\mu (a)\mu (c)\\ +{\displaystyle \sum _{\delta =1}^{r-1}(n-\delta -1)\left(\mu (c){\Pi }^{\delta }(c,a)+\mu (a){\Pi }^{\delta }(b,c)\right)(61)}\end{array}$

$\frac{O(ab,c)}{\Pi (a,b)}=(n-1)\mu (a)({\mathbb{I}}_{\{a=c\}}+{\mathbb{I}}_{\{b=c\}})(62)$

$\begin{array}{l}\frac{M(ab,cd)}{\Pi (a,b)\Pi (c,d)}\simeq (n-r-1)(n-r-2)\mu (a)\mu (c)\\ +{\displaystyle \sum _{\delta =1}^{r-1}(n-\delta -2)\left(\mu (c){\Pi }^{\delta }(d,a)+\mu (a){\Pi }^{\delta }(b,c)\right)(63)}\end{array}$

$\begin{array}{l}\frac{O(ab,cd)}{\Pi (a,b)}=(n-1)\mu (a){\mathbb{I}}_{\{ab=cd\}}\\ +(n-2)\Pi (c,d)\left(\mu (c){\mathbb{I}}_{\{a=d\}}+\mu (a){\mathbb{I}}_{\{b=c\}}\right)(64)\end{array}$

With the example given in section "validation" we get for the expectation

and

The magnitude of the second eigenvalue of Π is λ = 0.3, then rank r = 19 give a relative error < 10 ^-10 and we get for the covariance

and

$C_{1,1}=\left[\begin{array}{cccc}1536.8& -390.0& -390.0& -756.9\\ -390.0& 581.2& 581.0& -772.2\\ -390.0& 581.0& 581.2& -772.2\\ -756.9& -772.2& -772.2& 2301.4\end{array}\right](69)$

The beta function is defined by

for all a, b > 0. The incomplete beta function for all x ∈ [0,1] is then defined by

$\beta (x,a,b)={\displaystyle {\int }_0^x{t}^{a-1}}{(1-t)}^{b-1}dt(71)$

and

$\begin{array}{cc}{\beta }^{-}(x,a,b)=\beta (a,b)-\beta (x,a,b)& \left(72\right)\\ ={\displaystyle {\int }_x^1{t}^{a-1}}{(1-t)}^{b-1}dt& \left(73\right)\end{array}$

Using a continued fraction representation, these functions can be quickly numerically evaluated in O($\sqrt{\max (a,b)}$) in the worst case [15, Chapter 6].

A great interest of this function is that it is connected to the cumulative distribution function of a binomial distribution by the following relation:

$\mathbb{P}(ℬ(n,p)\ge k)=\frac{\beta (p,k,n-k+1)}{\beta (k,n-k+1)}(74)$

with (n, k) ∈ ℕ* × ℕ, 0 ≤ k ≤ n and p ∈ [0,1].

Finally, let us remark that the incomplete beta function is differentiable in x and that

$\frac{\partial \beta (x,a,b)}{\partial x}=x^{a-1}{(1-x)}^{b-1}(75)$

We give here the complete expression of σ for a single pattern in the special case of an order m = 0 homogeneous Markov model of parameter μ.

The MLE of μ is given by

${\mu }_N=\frac{N_1}{n}(76)$

where N₁ is the frequency of all letters.

A Gaussian approximation gives

with E₁ = nμ and, for all a, b ∈ $\mathcal{A}$,

C_1,1 (a, b) = nμ (a) $\mathbb{I}$_{a = b}- nμ(a) × nμ(b)     (78)

We have also

$P(N)=\frac{1}{n^h}{\displaystyle \prod _{a\in \mathcal{A}}{N}_1{(a)}^{A_1(a)}}(79)$

which implies for all a ∈ $\mathcal{A}$ that

So finally we get

where Q is either defined by equation if the pattern is over-represented or by equation if under-represented.