# Une histoire de mots inattendus et de génomes

Sophie Schbath

► To cite this version:

Sophie Schbath. Une histoire de mots inattendus et de génomes. ALEA 2017, Mar 2017, Luminy, France. pp.63. hal-01570855

HAL Id: hal-01570855

https://hal.science/hal-01570855

Submitted on 5 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Une histoire de mots inattendus et de génomes

Sophie Schbath

<figure><img src="image_2.png" /></figure>

MaiAGE

MATHEMATIQUES ET INFORMATIQUE APPLIQUEES DU GENOME A L'ENVIRONNEMENT

<figure><img src="image_3.png" /></figure>

INRA
SCIENCE & IMPACT

ALEA 2017, Marseille, 22 mars 2017
Introduction
* **DNA** : Long molecule, sequence of nucleotides

* **Nucleotides** : A(denine), C(ytosine), G(uanine), T(hymine).

<figure><img src="image_4.png" /></figure>

...GTTCAATCGTAGGTAGGTACTGAATGGTAGGTATGTTGA...
* **DNA** : Long molecule, sequence of nucleotides

* **Nucleotides** : A(denine), C(ytosine), G(uanine), T(hymine).

* **Motif** (= oligonucleotides) : short sequence of nucleotides, e.g. **AGGTA**

<figure><img src="image_5.png" /></figure>

...GTTCAATCGTAGGTAGGTACTGAATGGTAGGTATGTTGA...
* **Functional motif** : recognized by proteins or enzymes to initiate a biological process

<figure><img src="image_6.png" /></figure>
Some functional motifs

* **Restriction sites**: recognized by specific bacterial restriction enzymes → double-strand DNA break.

E.g. **GAATTC** recognized by **EcoRI**

* **Chi motif**: recognized by an enzyme which processes along DNA sequence and degrades it → enzyme degradation activity stopped and DNA repair is stimulated by recombination.

E.g. **GCTGGTGG** recognized by **RecBCD** (*E. coli*)

* **parS**: recognized by the *Spo0J* protein → organization of *B. subtilis* genome into macro-domains.

TGTTAACACGTGAAACA
c t c t t

* **promoter**: structured motif recognized by the RNA polymerase to initiate gene transcription.

E.g. **TTGAC** (16;18) - - | **TATAAT** (*E. coli*).
* **Restriction sites**: recognized by specific bacterial restriction enzymes → double-strand DNA break.

E.g. **GAATTC** recognized by **EcoRI**

very rare along bacterial genomes

* **Chi motif**: recognized by an enzyme which processes along DNA sequence and degrades it → enzyme degradation activity stopped and DNA repair is stimulated by recombination.

E.g. **GCTGGTGG** recognized by **RecBCD (E. coli)**

very frequent along *E. coli* genome

* **parS**: recognized by the *SpoJ* protein → organization of *B. subtilis* genome into macro-domains.

TGTTAACACGTGAAACA
c t

very frequent into the ORI domain, rare elsewhere

* **promoter**: structured motif recognized by the RNA polymerase to initiate gene transcription.

E.g. **TTGAC** (16;18) - **TATAAT** (*E. coli*).

particularly located in front of genes
# Prediction of functional motifs

Most of the functional motifs are unknown in the different species.

For instance,

* which would be the **Chi motif** of *S. aureus* ? [Halpern *et al.* (07)]

* Is there an equivalent of **parS** in *E. coli* ? [Mercier *et al.* (08)]

**Statistical approach** : to identify candidate motifs based on their statistical properties.

The most over-represented
8-letter words under M1

*E. coli* ($\ell$ = 4.6 $10^6$)

<table><thead><tr><th>word</th><th>obs</th><th>exp</th><th>score</th></tr></thead><tbody><tr><td>gctggtgg</td><td>762</td><td>84.9</td><td>73.5</td></tr><tr><td>ggcgctgg</td><td>828</td><td>125.9</td><td>62.6</td></tr><tr><td>cgctggcg</td><td>870</td><td>150.8</td><td>58.6</td></tr><tr><td>gctggcgg</td><td>723</td><td>125.9</td><td>53.3</td></tr><tr><td>cgctggtg</td><td>619</td><td>101.7</td><td>51.3</td></tr></tbody></table>

The most over-represented families
*abcdefg* under M1

*H. influenzae* ($\ell$ = 1.8 $10^6$)

<table><thead><tr><th>motif</th><th>obs</th><th>exp</th><th>score</th></tr></thead><tbody><tr><td>gntggtgg</td><td>223</td><td>55.3</td><td>22.33</td></tr><tr><td>anntcatc</td><td>469</td><td>180.3</td><td>21.59</td></tr><tr><td>anatcgcc</td><td>288</td><td>87.8</td><td>21.38</td></tr><tr><td>tnatcgcc</td><td>279</td><td>84.5</td><td>21.18</td></tr><tr><td>gnagaaga</td><td>270</td><td>83.6</td><td>20.10</td></tr></tbody></table>
Here are some quantities of interest.

* **Number of occurrences** (overlapping or not) :

- Is $N^{\text{obs}}(\mathbf{w})$ significantly high ?

- Is $N^{\text{obs}}(\mathbf{w})$ significantly higher than $N^{\text{obs}}(\mathbf{w}')$ ?

- Is $N_1^{\text{obs}}(\mathbf{w})$ significantly more unexpected than $N_2^{\text{obs}}(\mathbf{w})$ ?

* **Distance** between motif occurrences :

- Are there significantly rich regions with motif **w**

- Are two motifs significantly correlated ?

* **Waiting time** till the first occurrence :

- Is the presence of a motif **w** significant ?
# A model to define what to expect

Assessing the significance of an observed value (count, distance, occurrence, etc.) requires to define a **null model** to set what to expect.

A model for random sequences :

* **Markov chain models** : a Markov chain of order *m* (M*m*) fits the *h*-mers frequencies for *h* = 1, ..., (*m* + 1).

* Hidden Markov models allow to integrate heterogeneity.
Assessing the significance of an observed value (count, distance, occurrence, etc.) requires to define a **null model** to set what to expect.

## A model for random sequences :

* **Markov chain models** : a Markov chain of order *m* (*Mm*) fits the *h*-mers frequencies for *h* = 1, ..., (*m* + 1).

* Hidden Markov models allow to integrate heterogeneity.

## A model for the occurrence processes :

* **(compound) Poisson processes** allow to fit the number of occurrences and then to study the significance of inter-arrival times ([Robin (02)], or to compare the exceptionality of a word in two sequences ([Robin et al. (07)]).

* **Hawkes processes** allow to estimate the dependence between occurrence processes ([Gusto and S. (05)], [Reynaud and S. (10)])
Let $X_1, X_2, X_3, \dots, X_\ell, \dots$ be a stationary Markov chain of order $m$ on $\mathcal{A} = \{\text{a}, \text{c}, \text{g}, \text{t}\}$, i.e.

$$ \mathbb{P}(X_i = b | X_1, X_2, \dots, X_{i-1}) = \mathbb{P}(X_i = b | X_{i-m}, \dots, X_{i-1}). $$

Transition probabilities are denoted by

$$ \pi(a_1 \cdots a_m, b) = \mathbb{P}(X_i = b | X_{i-m} \cdots X_{i-1} = a_1 \cdots a_m), $$

whereas the stationary distribution is given by

$$ \mu(a_1 a_2 \cdots a_m) := \mathbb{P}(X_i = a_1, \dots, X_{i+m-1} = a_m), \quad \forall i. $$
Let $X_1, X_2, X_3, \dots, X_\ell, \dots$ be a stationary Markov chain of order $m$ on $\mathcal{A} = \{\text{a}, \text{c}, \text{g}, \text{t}\}$, i.e.

$$ \mathbb{P}(X_i = b | X_1, X_2, \dots, X_{i-1}) = \mathbb{P}(X_i = b | X_{i-m}, \dots, X_{i-1}). $$

Transition probabilities are denoted by

$$ \pi(a_1 \cdots a_m, b) = \mathbb{P}(X_i = b | X_{i-m} \cdots X_{i-1} = a_1 \cdots a_m), $$

whereas the stationary distribution is given by

$$ \mu(a_1 a_2 \cdots a_m) := \mathbb{P}(X_i = a_1, \dots, X_{i+m-1} = a_m), \quad \forall i. $$

The MLE are

$$ \begin{gathered} \hat{\pi}(a_1 \cdots a_m, a_{m+1}) = \frac{N^{\text{obs}}(a_1 \cdots a_m a_{m+1})}{N^{\text{obs}}(a_1 a_2 \cdots a_m)}, \qquad \hat{\mu}(a_1 \cdots a_m) = \frac{N^{\text{obs}}(a_1 \cdots a_m)}{\ell - m + 1} \\ \rightarrow \hat{\mathbb{E}} N(a_1 \cdots a_m a_{m+1}) \approx N^{\text{obs}}(a_1 \cdots a_m a_{m+1}) \end{gathered} $$
# Overlapping occurrences

Occurrences of words may overlap in DNA sequences (no space between words).

⇒ occurrences are not independent.

* Occurrences of overlapping words will tend to occur in clumps. For instance, they are 3 overlapping occurrences of CAGCAG below :

TAGACAGATAGACGAT CAGCAGCAGCAG ACAGTAGGCATGA...

* On the contrary, occurrences of non-overlapping words will never overlap.
# Overlapping occurrences (2)

All results on word occurrences will depend on the overlapping structure of the words.

Classically, this structure is described thanks to the **periods** of a word :

$p$ is a period of $\mathbf{w} := w_1 w_2 \cdots w_h$ iff $w_i = w_{i+p}, \quad \forall i$

meaning that $2$ occurrences of $\mathbf{w}$ can overlap on $h - p$ letters.

<figure><img src="image_7.png" /></figure>
Overlapping occurrences (2)

All results on word occurrences will depend on the overlapping
structure of the words.

Classically, this structure is described thanks to the **periods** of a word :

$p$ is a period of $\mathbf{w} := w_1 w_2 \cdots w_h$ iff $w_i = w_{i+p}, \quad \forall i$

meaning that 2 occurrences of **w** can overlap on $h - p$ letters.

<figure><img src="image_8.png" /></figure>

We also define the overlapping indicator :

$\varepsilon_{h-p}(\mathbf{w}) = 1$ if $p$ is a period of $\mathbf{w}$, and 0 otherwise
Detecting words with significantly unexpected counts
Let $N(\mathbf{w})$ be the number of occurrences of the word $\mathbf{w} := w_1 w_2 \cdots w_h$ in the sequence $X_1 X_2 X_3 \cdots X_\ell$ (model M1) :

$$N(\mathbf{w}) = \sum_{i=1}^{\ell-h+1} Y_i$$

where

$$Y_i = \mathbf{1}\{\mathbf{w} \text{ starts at position } i\} \sim \mathcal{B}(\mu(\mathbf{w}))$$

and

$$\mu(\mathbf{w}) = \mu(w_1) \prod_{j=1}^{h-1} \pi(w_j, w_{j+1}).$$
Problem

Let $\mathcal{N}(\mathbf{w})$ be the number of occurrences of the word $\mathbf{w} := w_1 w_2 \cdots w_h$
in the sequence $X_1 X_2 X_3 \cdots X_\ell$ (model M1) :

$$
\mathcal{N}(\mathbf{w}) = \sum_{i=1}^{\ell-h+1} Y_i
$$

where

$$
Y_i = \mathbf{1}\{\mathbf{w} \text{ starts at position } i\} \sim \mathcal{B}(\mu(\mathbf{w}))
$$

and

$$
\mu(\mathbf{w}) = \mu(w_1) \prod_{j=1}^{h-1} \pi(w_j, w_{j+1}).
$$

**Question** : how to decide if $\mathcal{N}^{\text{obs}}(\mathbf{w})$ is significantly unexpected (under model M1) ?

**Ideally** : one should compute the *p*-value $\mathbb{P}(\mathcal{N}(\mathbf{w}) \ge \mathcal{N}^{\text{obs}}(\mathbf{w}))$ or at least compare $\mathcal{N}^{\text{obs}}(\mathbf{w})$ with the expected count $\mathbb{E}\mathcal{N}(\mathbf{w})$
If the $Y_i$'s were independent, we would have

$$ \sum_{i=1}^{\ell} Y_i \sim \mathcal{B}(\ell, \mu(\mathbf{w})) \text{ approx by } \begin{cases} \mathcal{P}(\ell\mu(\mathbf{w})) & \text{if } \ell\mu(\mathbf{w}) \text{ small,} \\ \mathcal{N}(\ell\mu(\mathbf{w}), \ell\mu(\mathbf{w})(1-\mu(\mathbf{w}))) & \text{if } \ell\mu(\mathbf{w}) \sim \infty. \end{cases} $$

But the $Y_i$'s are not independent (overlaps) :

* For non-overlapping words, such as **ATGAC**, $Y_i = 1 \Rightarrow Y_{i+1} = 0$.

* For overlapping words, such as **ATGAT**,

$$ \mathbb{P}(Y_{i+3} = 1 | Y_i = 1) > \mathbb{P}(Y_{i+3} = 1). $$
* In the 80's, the ratio $\frac{N^{\text{obs}}(\mathbf{w})}{\mathbb{E}N(\mathbf{w})}$ was used with

$$ \mathbb{E}N(\mathbf{w}) = (\ell - h + 1)\mu(\mathbf{w}) = (\ell - h + 1)\mu(w_1) \prod_{j=1}^{h-1} \pi(w_j, w_{j+1}) $$

→ problem with the variability around 1 : $\text{Var}(N)$?
* In the 80's, the ratio $\frac{\mathcal{N}^{\text{obs}}(\mathbf{w})}{\mathbb{E}\mathcal{N}(\mathbf{w})}$ was used with

$$ \mathbb{E}\mathcal{N}(\mathbf{w}) = (\ell - h + 1)\mu(\mathbf{w}) = (\ell - h + 1)\mu(w_1) \prod_{j=1}^{h-1} \pi(w_j, w_{j+1}) $$

$\rightarrow$ problem with the variability around 1 : $\text{Var}(N)$ ?

* Normalization by $\mathbb{E}\mathcal{N}(\mathbf{w})$ like for a Poisson variable [Brendel et al. (86)]

$$ \frac{\mathcal{N}^{\text{obs}}(\mathbf{w}) - \mathbb{E}\mathcal{N}(\mathbf{w})}{\sqrt{\mathbb{E}\mathcal{N}(\mathbf{w})}} $$

$\rightarrow$ problem with the variability around 0 : $\text{Var}(N) \neq \mathbb{E}(N)$.
## Scores (2)

* The variance formula was published in 1992 by Kleffe & Borodowsky. The overlapping structure of the word clearly appears.

$$ 
\begin{aligned}
\mathrm{Var}[N(\mathbf{w})] = & (\ell - h + 1)\mu(\mathbf{w})[1 - \mu(\mathbf{w})] \\
& + 2\mu(\mathbf{w}) \sum_{d=1}^{h-1} (\ell - h - d + 1) \left[ \varepsilon_{h-d}(\mathbf{w}) \prod_{j=h-d+1}^{h} \pi(w_{j-1}, w_j) - \mu(\mathbf{w}) \right] \\
& + 2\mu^2(\mathbf{w}) \sum_{t=1}^{\ell-2h+1} (\ell - 2h - t + 2) \left[ \frac{1}{\mu(w_1)} \pi^t(w_h, w_1) - 1 \right]
\end{aligned}
$$

One then uses the z-score and the Central Limit Theorem :

$$ \frac{N(\mathbf{w}) - E[N(\mathbf{w})]}{\sqrt{\mathrm{Var}(N(\mathbf{w}))}} \rightarrow \mathcal{N}(0, 1) \text{ as } \ell \rightarrow \infty. $$
* The variance formula was published in 1992 by Kleffe & Borodowsky. The overlapping structure of the word clearly appears.

One then uses the z-score and the Central Limit Theorem :

$$ \frac{N(\mathbf{w}) - \mathbb{E}N(\mathbf{w})}{\sqrt{\mathrm{Var}(N(\mathbf{w}))}} \rightarrow \mathcal{N}(0, 1) \text{ as } \ell \rightarrow \infty. $$

$\rightarrow$ problem when parameters $(\pi, \mu)$ are unknown and have to be estimated by their MLE $(\hat{\pi}, \hat{\mu})$

Indeed,

$$ \mathrm{Var}(N - \mathbb{E}(\hat{N})) \neq \mathrm{Var}(\hat{N}) $$
* Prum, Rodolphe, de Turckheim (95) proposed an appropriate normalizing factor $\hat{\sigma}$ for $N - \widehat{\mathbb{E}N}$ which depends on the overlapping structure of the word.

It leads to the following score

$$ \frac{\mathcal{N}(\mathbf{w}) - \widehat{\mathbb{E}N}(\mathbf{w})}{\widehat{\sigma}(\mathbf{w})} \rightarrow \mathcal{N}(0, 1) \text{ as } \ell \rightarrow \infty. $$

and an approximation of the p-value :

$$ \mathbb{P}(N \geq \mathcal{N}^{\text{obs}}) \approx \mathbb{P} \left( \mathcal{N}(0, 1) \geq \frac{\mathcal{N}^{\text{obs}} - \widehat{\mathbb{E}N}}{\widehat{\sigma}} \right) $$

A similar score has been derived under model Mm.
Which model to use ?

<figure><img src="image_9.png" /></figure>

Scores of exceptionality for the 65,536 **8-letter words** in the *E.coli* backbone.
Example : gctggtgg occurs 762 times in the E. coli's genome (leading strands, $ℓ = 4.610^6$).

<table><thead><tr><th>model</th><th>fit</th><th>expected</th><th>score</th><th>p-value</th><th>rank</th></tr></thead><tbody><tr><td>M00</td><td>length</td><td>70.783</td><td></td><td></td><td></td></tr><tr><td>M0</td><td>bases</td><td>85.944</td><td>72.9</td><td>&lt; 10<sup>-323</sup></td><td>3</td></tr><tr><td>M1</td><td>dinucl.</td><td>84.943</td><td>73.5</td><td>&lt; 10<sup>-323</sup></td><td>1</td></tr><tr><td>M2</td><td>trinucl.</td><td>206.791</td><td>38.8</td><td>&lt; 10<sup>-323</sup></td><td>1</td></tr><tr><td>M3</td><td>tetranucl.</td><td>355.508</td><td>22.0</td><td>1.4 10<sup>-107</sup></td><td>5</td></tr><tr><td>M4</td><td>pentanucl.</td><td>355.312</td><td>22.9</td><td>2.3 10<sup>-116</sup></td><td>2</td></tr><tr><td>M5</td><td>hexanucl.</td><td>420.867</td><td>19.7</td><td>1.0 10<sup>-86</sup></td><td>1</td></tr><tr><td>M6</td><td>heptanucl.</td><td>610.114</td><td>10.6</td><td>1.5 10<sup>-26</sup></td><td>3</td></tr></tbody></table>
Influence of the model (2)

<table><thead><tr><th></th><th colspan="2">gctggtgg</th><th colspan="2">ggcgctgg</th><th colspan="2">ccggccta</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>M0</td><td>85.944</td><td>&lt; 10<sup>-323</sup> (3)</td><td>85.524</td><td>&lt; 10<sup>-323</sup> (2)</td><td>70.445</td><td>0.47 (25608)</td></tr><tr><td>M1</td><td>84.943</td><td>&lt; 10<sup>-323</sup> (1)</td><td>125.919</td><td>&lt; 10<sup>-323</sup> (2)</td><td>48.173</td><td>10<sup>-3</sup> (13081)</td></tr><tr><td>M2</td><td>206.791</td><td>&lt; 10<sup>-323</sup> (1)</td><td>255.638</td><td>10<sup>-283</sup> (3)</td><td>35.830</td><td>10<sup>-8</sup> (4436)</td></tr><tr><td>M3</td><td>355.508</td><td>1.4 10<sup>-107</sup> (5)</td><td>441.226</td><td>10<sup>-78</sup> (15)</td><td>14.697</td><td>10<sup>-49</sup> (47)</td></tr><tr><td>M4</td><td>355.312</td><td>2.3 10<sup>-116</sup> (2)</td><td>392.252</td><td>10<sup>-120</sup> (1)</td><td>15.341</td><td>10<sup>-46</sup> (21)</td></tr><tr><td>M5</td><td>420.867</td><td>1.0 10<sup>-86</sup> (1)</td><td>633.453</td><td>10<sup>-22</sup> (24)</td><td>27.761</td><td>10<sup>-18</sup> (36)</td></tr><tr><td>M6</td><td>610.114</td><td>1.5 10<sup>-26</sup> (3)</td><td>812.339</td><td>0.16 (14686)</td><td>25.777</td><td>10<sup>-26</sup> (4)</td></tr></tbody></table>

Expected counts and *p*-values (rank) under models *Mm*,

$m = 0, 1, \dots, 6$,

estimated from the *E. coli*'s genome (4 638 858 bps, leading strands).
Another approximation for rare words

The Gaussian approximation appeared to be not accurate for
expectedly rare words ($\mathbb{E}(N(\mathbf{w})) = O(1)$ as $\ell \to +\infty$).

Here “**w**” is rare along the sequence.

* If **w** is not self-overlapping : $N(\mathbf{w}) \sim \text{Pois}(\mathbb{E}[N(\mathbf{w})])$ (Chen-Stein method).

* In the general case, $N(\mathbf{w})$ is approximated by a Geometric-Poisson distribution with parameter $((1 - a(\mathbf{w})))\mathbb{E}[N(\mathbf{w})]$; $a(\mathbf{w})$ (S. (95)).

Both (compound) Poisson approximations are still valid when plugging
the estimated parameters of the Markov model.
In the general case : clump decomposition

<figure><img src="image_10.png" /></figure>

$$N(\mathbf{w}) = \sum_{c=1}^{\tilde{N}(\mathbf{w})} K_c$$

where

* $\tilde{N}(\mathbf{w})$ is the number of clumps

* $K_c$ is the size of the c-th clump
# Compound Poisson approximation ($\mathbb{E}[N(\mathbf{w})] = O(1)$)

In the general case : clump decomposition

<figure><img src="image_11.png" /></figure>

$$N(\mathbf{w}) = \sum_{c=1}^{\tilde{N}(\mathbf{w})} K_c \qquad a(\mathbf{w}) = \sum_{\substack{\text{main periods} \\ p \mid j=1}} \prod_{j=1}^{p} \pi(w_j, w_{j+1}).$$

where

is the overlapping proba.

* $\tilde{N}(\mathbf{w})$ is the number of clumps
  can be approximated by $\text{Pois}((1 - a(\mathbf{w}))\mathbb{E}[N(\mathbf{w})])$ (Chen-Stein)

* $K_c$ is the size of the c-th clump
# Compound Poisson approximation ($\mathbb{E}[N(\mathbf{w})] = O(1)$)

In the general case : clump decomposition

<figure><img src="image_12.png" /></figure>

$$N(\mathbf{w}) = \sum_{c=1}^{\tilde{N}(\mathbf{w})} K_c \qquad a(\mathbf{w}) = \sum_{\text{main periods } p} \prod_{j=1}^{p} \pi(w_j, w_{j+1}).$$

where

is the overlapping proba.

* $\tilde{N}(\mathbf{w})$ is the number of clumps can be approximated by $\text{Pois}((1 - a(\mathbf{w}))\mathbb{E}[N(\mathbf{w})])$ (Chen-Stein)

* $K_c$ is the size of the c-th clump follows a geometric distribution $G(a(\mathbf{w}))$
Chen(75), Stein(71), Arratia et al. (89)

$$Y_i \sim \mathcal{B}(p_i) \qquad N = \sum_{i \in I} Y_i$$

$$Z_i \sim \mathcal{P}o(p_i) \text{ indep.} \qquad Z = \sum_{i \in I} Z_i$$

$$d_{\text{TV}}(\mathcal{L}(N) - \mathcal{L}(Z)) \le 2(b_1 + b_2 + b_3)$$

$$\text{where} \quad b_1 = \sum_{i \in I} \sum_{j \in B_i} \mathbb{E} Y_i \mathbb{E} Y_j, \quad B_i \text{ is any neighborhood of } i \text{ in } I$$

$$b_2 = \sum_{i \in I} \sum_{j \in B_i \setminus \{i\}} \mathbb{E}(Y_i Y_j)$$

$$b_3 = \sum_{i \in I} \mathbb{E} \left| \mathbb{E}(Y_i - p_i \mid \sigma(Y_j, j \notin B_i)) \right|.$$
* Probability for a clump of **w** to start at a given position

<figure><img src="image_13.png" /></figure>

$$ (1 - a) \quad \mu(\mathbf{w}) $$
* Probability for a clump of **w** to start at a given position

(1 - a) $\mu(\mathbf{w})$

<figure><img src="image_14.png" /></figure>

* Probability for a *k*-clump of **w** to start at a given position

(1 - a) $\underbrace{\underline{a \ a \ a \ a}}_{a^{k-1}} \mu(\mathbf{w}) (1 - a)$

<figure><img src="image_15.png" /></figure>
Under model M1 (with known parameters)

* the exact distribution of the count $N(\mathbf{w})$ can be computed

* via its generating function [Régnier (00)],

* via the duality equation $\mathbb{P}(N(\mathbf{w}) \ge x) = \mathbb{P}(T_x \le \ell)$ where $T_x$ is the position of the $x$-th occurrence; The distribution of $T_x$ can be obtained by recursion or via its generating function (Robin & Daudin (99), Stefanov (03)).

* large deviation technique can be used to directly approximate the $p$-value [Nuel (04)] :

$$ \mathbb{P}(N \ge N^{\text{obs}}) \approx \exp(-\ell I(N^{\text{obs}})). $$

It is a very accurate (but numerically costly) method for very exceptional words.
Prediction and identification of functional DNA motifs
* Motif involved in the repair of double-strand DNA breaks. Chi needs to be frequent along bacterial genomes.

* Chi motifs have been identified for few bacterial species. They are not conserved through species.

* Known Chi motifs are 5 to 8 nucleotides long and can be degenerated.

* Moreover, Chi activity is strongly orientation-dependent (direction of DNA replication). It is present preferentially on the leading strands (high skew).

*The skew of a motif w is defined by $N^{obs}(w)/N^{obs}(\bar{w})$ where $\bar{w}$ is the reverse complementary of w.*
E. coli as a learning case

* 8-letter word GCTGGTGG

* 762 occurrences on the leading strands ($l = 4.6 \times 10^6$)

* Among the most over-represented 8-letter words (whatever the model $Mm$)
  ⇒ its frequency cannot be explained by the genome composition.

* Its rank is improved if one analyzes only the backbone genome (genome conserved in several strains of the species).

* Its skew equals 3.20 ($p$-value of $3.310^{-11}$).

*The skew significance can be evaluated thanks to the Gaussian approximation of word counts.*
Halpern *et al.* (2007)

* Analysis of the *S. aureus* backbone (ℓ = 2.44 10⁶).

* 8-letter words : none of the most over-represented and skewed motifs were frequent enough.

* 7-letter words :

<figure><img src="image_16.png" /></figure>

A=gaaaatg (1067),

B=ggattag (266),

C=gaagcgg (272),

D=gaattag (614)
Toward more complex motifs
# Signature Motif of the Ter Macrodomain of *E. coli*

## The MatP/matS Site-Specific System Organizes the Terminus Region of the *E. coli* Chromosome into a Macrodomain

Romain Mercier,¹ Marie-Agnès Petit,² Sophie Schbath,³ Stéphane Robin,⁴ Meriem El Karoui,² Frédéric Boccard,¹* and Olivier Espaillie²*

### Cell (2008)

#### Use of R'MES software :

* exceptional frequency

* exceptional contrast

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>
**Definition** : a word whose one or more positions may tolerate different nucleotides. The IUPAC alphabet maybe used (R=A or G, Y= C or T, N=A or C or G or T, etc.)

**Examples** :

* the Chi motif of *H. influenzae* is gNtggtgg

* the matS motif of *E. coli* is gtgacRNYgtcac

→ one will consider them like a family *W* of words :

$$N(\mathcal{W}) = \sum_{\mathbf{w} \in \mathcal{W}} N(\mathbf{w})$$
# Degenerated motifs (2)

Gaussian approximation [S. (95)] :

* $\mathbb{E}(N(\mathcal{W})) = \sum_{\mathbf{w} \in \mathcal{W}} \mathbb{E}(N(\mathbf{w}))$

* $\text{Var}(N(\mathcal{W})) : \text{need for } \text{Cov}(N(\mathbf{w}), N(\mathbf{w}'))$ 
  → one needs to know all possible overlaps between $\mathbf{w}$ and $\mathbf{w}'$
# Degenerated motifs (2)

Gaussian approximation [S. (95)] :

* $\mathbb{E}(N(\mathcal{W})) = \sum_{\mathbf{w} \in \mathcal{W}} \mathbb{E}(N(\mathbf{w}))$

* $\text{Var}(N(\mathcal{W})) : \text{need for } \text{Cov}(N(\mathbf{w}), N(\mathbf{w}'))$
  → one needs to know all possible overlaps between $\mathbf{w}$ and $\mathbf{w}'$

Compound Poisson approximation [Roquain and S. (07)] :

* mixed clumps need to be considered (again it requires all possible overlaps between any $\mathbf{w}$ and $\mathbf{w}'$ in the $\mathcal{W}$ family)

* the clump size is no more geometric, the overlap probability $a(\mathbf{w})$ is replaced by a matrix $A = (a(\mathbf{w}, \mathbf{w}'))$

* we still get a compound Poisson distribution for $N(\mathcal{W})$
Here is an example of a PWM of length *h* = 5 :

$$
\mathbf{m} = \begin{pmatrix} 1 & 0.25 & 0 & 0.25 & 0.3 \\ 0 & 0 & 0 & 0.25 & 0.1 \\ 0 & 0.75 & 1 & 0.25 & 0 \\ 0 & 0 & 0 & 0.25 & 0.6 \end{pmatrix} \begin{array}{l} A \\ C \\ G \\ T \end{array}
$$

$m_{a,j}$ = probability of letter *a* at motif position *j*

Such representation induces a set of “compatible” words having
different probabilities (or “weights”)
How to count occurrences of a PWM?

$$ \mathbf{m} = \begin{pmatrix} 1 & 0.25 & 0 & 0.25 & 0.3 \\ 0 & 0 & 0 & 0.25 & 0.1 \\ 0 & 0.75 & 1 & 0.25 & 0 \\ 0 & 0 & 0 & 0.25 & 0.6 \end{pmatrix} \begin{array}{c} A \\ C \\ G \\ T \end{array} $$

Weights :

$$ \nu_i = \prod_{j=1}^{h} m(X_{i+j-1}, j). $$

...GTTCGT**AGGTACGGTACTGATGGT**AAGTATGAGGCT... weights 0.05 0 0.02 0.1
## How to count occurrences of a PWM?

$$ \mathbf{m} = \begin{pmatrix} 1 & 0.25 & 0 & 0.25 & 0.3 \\ 0 & 0 & 0 & 0.25 & 0.1 \\ 0 & 0.75 & 1 & 0.25 & 0 \\ 0 & 0 & 0 & 0.25 & 0.6 \end{pmatrix} \begin{array}{l} A \\ C \\ G \\ T \end{array} $$

Weights :

$$ \nu_i = \prod_{j=1}^{h} m(X_{i+j-1}, j). $$

...GTTCGTAGGTACGGTACTGATGGTAAGTATGAGGCT... weights 0.05 0 0.02 0.1

**Classical approach** : to count the number of “hits” i.e. $\sum_i 1\{\nu_i \ge \alpha\}$

→ If the set of words $W = \{w, \nu(w) \ge \alpha\}$ is not too large, one can use previous results for a word family
## How to count occurrences of a PWM?

$$ \mathbf{m} = \begin{pmatrix} 1 & 0.25 & 0 & 0.25 & 0.3 \\ 0 & 0 & 0 & 0.25 & 0.1 \\ 0 & 0.75 & 1 & 0.25 & 0 \\ 0 & 0 & 0 & 0.25 & 0.6 \end{pmatrix} \begin{array}{c} A \\ C \\ G \\ T \end{array} $$

Weights :

$$ \nu_i = \prod_{j=1}^{h} m(X_{i+j-1}, j). $$

...GTTCGTAGGTACGGTACTGATGGTAAGTATGAGGCT...

**Classical approach** : to count the number of “hits” i.e. $\sum_i 1\{\nu_i \ge \alpha\}$

→ If the set of words $W = \{w, \nu(w) \ge \alpha\}$ is not too large, one can use previous results for a word family

Otherwise, there exists dedicated results [Touzet & Varré (07)], [Pape et al. (08)], [Turatsinze et al. (08)].
Drawbacks of the classical approach :

* choice of the threshold $\alpha$

* the hits are not weighted anymore

→ New approach : to directly study the distribution of the weighted count defined by

$$T(\mathbf{m}) = \sum_{i=1}^{\ell-h+1} \nu_i.$$
Drawbacks of the classical approach :

* choice of the threshold $\alpha$

* the hits are not weighted anymore

→ New approach : to directly study the distribution of the weighted count defined by

$$T(\mathbf{m}) = \sum_{i=1}^{\ell-h+1} \nu_i.$$

Note : if $\mathbf{m}$ is a word, there is a unique compatible word and both counts are equal to the total word count
Expectation and variance of $T(m)$ can be analytically derived

A Gaussian approximation can be performed

A compound Poisson approximation has been derived for

$$T = \sum_{c=1}^{C} K_{c} :$$

* the number $C$ of clumps of compatible words can be approximated by a Poisson variable with explicit parameter (Chen-Stein method)

* the distribution of $K_c$, the total weight of the $c$th clump can be simulated

A compound Poisson approximation is better than a Gaussian approximation as soon as occurrences of compatible words are rare ($h$ large enough and $\text{card(compatible words)} \ll 4^h$).
# NBS motif

<table><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th></tr></thead><tbody><tr><th>A</th><td>.88</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>.88</td></tr><tr><th>C</th><td>0</td><td>0</td><td>0</td><td>.11</td><td>1</td><td>.88</td><td>.73</td><td>.73</td><td>.88</td><td>1</td><td>.11</td><td>0</td><td>0</td><td>0</td></tr><tr><th>G</th><td>0</td><td>0</td><td>0</td><td>.11</td><td>0</td><td>0</td><td>.12</td><td>.12</td><td>0</td><td>0</td><td>.11</td><td>0</td><td>0</td><td>0</td></tr><tr><th>T</th><td>.12</td><td>1</td><td>1</td><td>.78</td><td>0</td><td>.12</td><td>.15</td><td>.15</td><td>.12</td><td>0</td><td>.78</td><td>1</td><td>1</td><td>.12</td></tr></tbody></table>
To be done :

* check the Gaussian approximation is good when few 0's in the PWM

* study the influence of the parameter estimation

* generalize to Markovian sequences
Structured motifs

<figure><img src="image_19.png" /></figure>

What is the probability for a structured motif to occur in a given sequence ?

**Difficulty** : even for 2 boxes, previous results on word counts cannot be used because the overlapping structure is too complicated. The structure of the motif need to be considered.
The 2 next approaches rely on the exact distribution of the following intersite distances in Markovian sequences (recursive formula or probability generating function) [Robin and Daudin (01)], [Stefanov (03)]:

* $T_{\alpha,w}$, the waiting time to reach pattern $w$ from state $\alpha$

* $T_{w,w'}$, the waiting time to reach pattern $w'$ from pattern $w$
# Structured motifs (3)

A first order approximation ([Robin et al. (02)])

* The probability $P(N(\mathbf{m}) = 0) = 1 - P(N(\mathbf{m}) \ge 1)$ is approximated by

$$ (1 - \mu(\mathbf{m})) (1 - \gamma(\mathbf{m}))^{\ell-|\mathbf{m}|} $$

where $\mu(\mathbf{m}) = P(\mathbf{m} \text{ occurs at position } i)$

$\gamma(\mathbf{m}) = P(\mathbf{m} \text{ at } i | \mathbf{m} \text{ not at } i-1)$

* The occurrence probability of $\mathbf{m}$ is calculated like

<figure><img src="image_20.png" /></figure>

where $T_{s,w_2}$ is the waiting time to reach pattern $w_2$ from state $s$
Structured motifs (4)

An exact approach via random sums ([Stefanov et al. (06)])

* Assumption : $\mathbf{w}_2$ should not occur in between the 2 boxes.

* The explicit formula for the pgf of $\tau_m$ is given thanks to the following decomposition :

$$ \tau_m \stackrel{\mathcal{D}}{=} T_{\alpha, \mathbf{w}_1} + \sum_{b=1}^{L'} \left( \underbrace{\sum_{a=1}^{L_1} X_{\mathbf{w}_1, \mathbf{w}_1}^{(ab)} + F_{\mathbf{w}_1, \mathbf{w}_2}^{(b)} + T_{\mathbf{w}_2, \mathbf{w}_1}^{(b)}}_{\stackrel{\mathcal{D}}{=} T_{\mathbf{w}_1, \mathbf{w}_2} \mid \text{failure}} + \underbrace{\sum_{c=1}^{L_2} X_{\mathbf{w}_1, \mathbf{w}_1}^{(c)} + S_{\mathbf{w}_1, \mathbf{w}_2, \mathbf{w}_2}}_{\stackrel{\mathcal{D}}{=} T_{\mathbf{w}_1, \mathbf{w}_2} \mid \text{success}} \right) $$

<figure><img src="image_21.png" /></figure>

$L_1$, $L_2$ and $L'$ are independent geometric variables.
## Overviews :

* REINERT, G., SCHBATH, S. and WATERMAN, M. (2005). *Applied Combinatorics on Words*. volume 105 of *Encyclopedia of Mathematics and its Applications*, chapter Statistics on Words with Applications to Biological Sequences. Cambridge University Press.

* REINERT, G., SCHBATH, S. and WATERMAN, M. (2000). Probabilistic and statistical properties of words : an overview. *J. Comp. Biol.* **7** 1–46.

* ROBIN, S., RODOLPHE, F. and SCHBATH, S. (2005). *DNA, Words and Models*. Cambridge University Press.

* ROBIN, S., RODOLPHE, F. and SCHBATH, S. (2003). *ADN, mots et modèles*. BELIN.

* SCHBATH, S. and ROBIN, S. (2008). *How pattern statistics can be useful for DNA motif discovery?*. To appear in *Scan Statistics - Methods and Applications*, Glaz, J., Pozdnyakov, I. and Wallenstein, S. Eds., Statistics for Industry and Technology series, Birkhauser.
## Word count (Poisson approximations) :

* ERHARDSSON, T. (2000). Compound Poisson approximation for counts of rare patterns in Markov chains and extreme sojourns in birth-death chains *Ann. Appl. Prob.*, **10**, 573–591.

* GESKE, M. X., GODBOLE, A. P., SCHAFFNER, A. A., SKOLNICK, A. M. and WALLSTROM, G. L. (1995). Compound Poisson approximations for word patterns under Markovian hypotheses. *J. Appl. Prob.*, **32**, 877–892.

* REINERT, G. and SCHBATH, S. (1998). Compound Poisson and Poisson process approximations for occurrences of multiple words in Markov chains. *J. Comp. Biol.*, **5**, 223–253.

* ROQUAIN, E. and SCHBATH, S. (2007). Efficient compound Poisson approximation for the number of occurrences of multiple words in Markov chains. *Adv. Appl. Prob.*, **39**, 128–140.

* SCHBATH, S. (1995). Compound Poisson approximation of word counts in DNA sequences. *ESAIM : Probability and Statistics*. **1**, 1–16.
## Word count (others) :

* KLEFFE, J. and BORODOVSKY, M. (1992). First and second moment of counts of words in random texts generated by Markov chains *Comp. Applic. Biosci.*, **8**, 433–441.

* NUEL, G. (2004). LD-SPatt : Large Deviations Statistics for Patterns on Markov chains. *J. Comp. Biol.*

* PRUM, B. RODOLPHE, F., TURCKHEIM, É. (1995). Finding words with unexpected frequencies in DNA sequences *J. R. Statist. Soc. B*, **57**, 205–220.

* PUDLO, P. (2004). Estimations précises de grandes déviations et applications à la statistique des séquences biologiques. PhD thesis, Université Lyon I.

* RÉGNIER, M. (2000). A unified approach to word occurrence probabilities *Discrete Applied Mathematics*, **104**, 259–280.

* RÉGNIER, M. and DENISE, A. (2004). Rare events and conditional events on random strings. *Discrete Mathematics and Theoretical Computer Science* **6** 191–214.

* ROBIN, S. and SCHBATH, S. (2001). Numerical comparison of several approximations of the word count distribution in random sequences. *J. Comp. Biol.* **8** 349–359.

* ROBIN, S., SCHBATH, S. and VANDEWALLE, V. (2007). Statistical tests to compare motif count exceptionalities. *BMC Bioinformatics* **8**:84 1–20.
Distances and waiting times :

* ROBIN, S. (2002). A compound Poisson model for words occurrences in DNA sequences
  *J. Royal Statist. Soc., C series, 51*, 437–451.

* ROBIN, S. and DAUDIN, J.-J. (1999). Exact distribution of word occurrences in a random sequence of letters *J. Appl. Prob.*, **36**, 179–193.

* ROBIN, S. and DAUDIN, J.-J. (2001). Exact distribution of the distances between any occurrences of a set of words *Ann. Inst. Statist. Math.*, **36**, 895–905.

* ROBIN, S., DAUDIN, J.-J., RICHARD, H., SAGOT, M.-F. and SCHBATH, S. (2002).
  Occurrence probability of structured motifs in random sequences. *J. Comp. Biol.* **9**
  761–773.

* STEFANOV, V. (2003). The intersite distances between pattern occurrences in strings generated by general discrete - and continuous- time models: an algorithmic approach. *J. Appl. Prob.*, **40**.

* STEFANOV, V., ROBIN, S. and SCHBATH, S. (2007). Waiting times for clumps of patterns and for structured motifs in random sequences. *Discrete Applied Mathematics*. **155** 868–880.
## Prediction of functional motifs :

* HALPERN, D., CHIAPELLO, H., SCHBATH, S., ROBIN, S., HENNEQUET-ANTIER, C., GRUSS, A. and EL KAROUI, M. (2007). Identification of DNA motifs implicated in maintenance of bacterial core genomes by predictive modelling. *PLoS Genetics.* **3**(9) e153.

* MERCIER, R., PETIT, M.-A., SCHBATH, S., ROBIN, S., EL KAROUI, M., BOCCARD, F. and ESPELI, O. (2008). The MatP/matS site specific system organizes the Terminus region of the E. coli chromosome into a Macrodomain. *Cell*.

* TOUZAIN, F., SCHBATH, S., DEBLED-RENNESSON, I., AIGLE, B., LEBLOND, P. and KUCHEROV, G. (2008). SI.GffRid : a tool to search for σ factor binding sites in bacterial genomes using comparative approach and biologically driven statistics. *BMC Bioinformatics.* **9**:73 1–23.

## Others :

* GUSTO, G. and SCHBATH, S. (2005). FADO : a statistical method to detect favored or avoided distances between motif occurrences using the Hawkes' model. *Statistical Applications in Genetics and Molecular Biology*.

* REYNAUD-BOURET, P. and SCHBATH, S. (2010). Adaptive estimation for Hawkes' processes; Application to genome analysis. *Annals of Statistics.* **38**(5) 2781–2822.