Modéliser l’évolution des séquences d’ADN

Bien que je m’empêche de réagir à chaud, l’actualité médiatique est une source permanente d’inspiration. Alors quand les mathématiciens avertissent du desintérêt croissant pour leur discipline et que l’apparition de virus HINI mutants crée des frayeurs, l’envie me prend de dire quelques mots mêlant modèle mathématique et mutations dans des séquences d’ADN.

Attention le billet est long mais à la fin puisse votre persévérance être récompensée. Non seulement vous connaîtrez une formule estimant la distance entre deux séquences d’ADN, mais vous serez aussi passés par toutes les étapes du calcul (ce qui est généralement négligé partout ailleurs sur internet). Et puis le fait de rentrer dans le coeur d’un modèle, suivant pas à pas la démarche du modélisateur, pourrait réconcilier certains lecteurs avec les maths ? En tout cas moi ça m’a bien servi d’écrire ce billet.

On l’a vu précédemment, l’ADN est la molécule support de l’information génétique. Elle est formée de quatre « briques » appelées nucléotides, chacune symbolisée par une lettre: A, T, G et C. Chaque cellule possède une molécule d’ADN qui est copiée dans son intégralité au moment de la division cellulaire. Or, lors de cette réplication, des erreurs peuvent être commises, par exemple un A remplacé par un G. On parle alors de mutation (de substitution pour être plus précis). Par simplicité on ne traite pas le cas des insertions/délétions.

Depuis les années 1970, avec la possibilité de séquencer de l’ADN, on peut observer un fragment du génome d’un organisme et le comparer avec le même fragment mais provenant d’un autre organisme. En alignant ces deux séquences, on peut avoir une idée du nombre de match (:) et mismatch (*) les séparant, comme le montre l’image ci-dessous:

figure 1: un alignement global entre deux séquences d’ADN

En voyant cela, on se rappelle qu’une séquence d’ADN est une chaîne de caractères de longueur finie. Chaque position est appelée un site. Chaque site est dans un état particulier parmi quatre possibles (A, T, G ou C). Lorsque plusieurs séquences ont un ancêtre commun, on parle de séquence homologues. Sur la figure 1 la séquence S_1 et la séquence S_2 ont la même longueur L, 20 nucléotides, l’alignement résultant a donc 20 sites.

Pour mesurer la distance entre deux séquences, le plus simples est de calculer la proportion de sites différents. Sur la figure 1 une distance de 0.1 sépare la séquence S_1 de la séquence S_2 (2 mismatches sur 20 sites). Le problème c’est qu’en faisant cela on néglige les mutations cachées. Par exemple au site n°2 supposons que l’ancêtre commun ait été le nucléotide A. La séquence S_1 n’a pas eu de mutation mais la séquence S_2 aurait pu en avoir deux, d’abord de A vers C puis de C vers A. Or avec la distance décrite ci-dessus on ne compte pas de mutation puisque les deux séquences sont dans l’état A au 2e site. Dans la même veine, au 8e site supposons que la séquence ancestrale ait été dans l’état G. Aujourd’hui on observe que la séquence S_1 est toujours dans l’état G mais la séquence S_2 est dans l’état A, on compte donc une mutation, mais qui nous dit qu’il n’y en a pas eu plusieurs, par exemple de G vers T puis de T vers A ? Alors allons-y, écrivons un modèle mathématique prenant cela en compte !

On suppose que tous les sites suivent la même distribution de probabilités et que chaque site évolue de façon indépendante des autres. Ainsi la probabilité de passer de la séquence S_1 à la séquence S_2 (toutes les deux de longueur L) est donnée par:

\mathbb{P}( S_1 \rightarrow S_2 ) = \displaystyle \prod_{i = 1} ^L \mathbb{P}( S_1[i] \rightarrow S_2[i] )

Dire que deux probabilités sont indépendantes revient à les multiplier. Comme on suppose que tous les suites suivent la même loi de probabilité, on peut se concentrer sur l’évolution d’un seul site. Grâce à la formule ci-dessus, si on arrive à modéliser \mathbb{P}( S_1[i] \rightarrow S_2[i] ) on arrivera à modéliser \mathbb{P}( S_1 \rightarrow S_2 ).

Supposons maintenant que le temps avance en « tic », comme les aiguilles d’une horloge, et qu’à chaque « tic » une mutation peut subvenir ou non. Pour modéliser cela on va utiliser une chaîne de Markov. Une chaine de Markov est un processus stochastique (synonyme de probabiliste). Par exemple, si on note X l’évènement « le dé jeté affiche la valeur x« , on dit que X est une variable aléatoire prenant les valeurs 1, 2 ... 6. Et bien, si on étudie plusieurs lancés de dé au cours du temps, on se retrouve à étudier un processus probabiliste: pas compliqué…

Si maintenant on considère que la probabilité au temps t_{n+1} ne dépend que de l’état présent, c’est-à-dire du temps t_{n}, et non des états passés, les temps t_{n-1}, t_{n-2}, etc, on dit que le processus possède la propriété de Markov. Résumé d’une autre façon: « le futur ne dépend du passé qu’au travers de l’instant présent ». Et c’est tout naturellement ce qui arrive en génétique: une mutation à la génération des petit-enfants ne va pas dépendre du nucléotide en question à la générations des grand-parents mais uniquement du nucléotide en question à la génération des parents.

Une chaîne de Markov est caractérisée par sa matrice de transition P. Quand on modélise l’évolution d’une séquence d’ADN, la chaîne a 4 états (pour A, T, G et C), et la matrice aura 4 lignes et 4 colonnes. La valeur au croisement de la ligne i et de la colonne j est la probabilité p_{ij} d’être dans l’état i et de passer dans l’état j.

En 1969, Jukes et Cantor proposaient de modéliser l’évolution d’une séquence via une chaîne de Markov sous l’hypothèse que la probabilité f(t) de passer d’un nucléotide à un autre pendant la durée t était constante au cours du temps. Voici la matrice de transition correspondante:

P = \begin{pmatrix}1 - 3f(t) & f(t) & f(t) & f(t) \\ f(t) & 1 - 3f(t) & f(t) & f(t) \\ f(t) & f(t) & 1 - 3f(t) & f(t) \\ f(t) & f(t) & f(t) & 1 - 3f(t) \end{pmatrix}

Je précise de manière arbitraire que les 1e ligne et colonne correspondent au nucléotide A, les 2e au nucléotide T, les 3e à G et les 4e à C: p_{23} à l’intersection de la 2e ligne et de la 3e colonne est la probabilité que le site soit dans l’état T et mute vers l’état G.

i \ne j: \mathbb{P}( i \rightarrow j ) = p_{ij}(t) =f(t)

i = j: \mathbb{P}( i \rightarrow j ) = p_{ii}(t) = 1 - 3f(t)

Si on est à la génération t et que le nucléotide que l’on est en train d’analyser est un « A », deux scénarios sont possibles:

  • il y a une mutation, de A vers T, G ou C, chaque évènement ayant une probability f(t) d’arriver;
  • il n’y a pas de mutation, ceci avec une probabilité 1 - 3f(t) (la somme de tous les évènements doit faire 1).

Maintenant calculons f(t) = p_{ij}(t) et pour ceci commençons par différencier la matrice de transition P, c’est-à-dire regardons ce que vaut cette matrice P à l’instant t + \Delta t (c’est-à-dire très peu de temps après l’instant t) on a:

P( t + \Delta t ) = P(t) \times P(\Delta t)

Donc si maintenant on fait tendre \Delta t vers 0 (ça vous rappelle la définition de la dérivée n’est-ce pas ?):

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P( t + \Delta t ) - P(t)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P(t) \times P(\Delta t) - P(t + 0)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} \frac{P(t) \times P(\Delta t) - P(t) \times P(0)}{\Delta t}

P'(t) = \lim_{\Delta t \rightarrow 0} P(t) \frac{P(\Delta t) - P(0)}{\Delta t}

P'(t) = P(t) \lim_{\Delta t \rightarrow 0} \frac{P(\Delta t) - P(0)}{\Delta t}

Ainsi: P'(t) = P(t) \times P'(0)

P'(0) = \begin{pmatrix} - 3f'(0) & f'(0) & f'(0) & f'(0) \\ f'(0) & - 3f'(0) & f'(0) & f'(0) \\ f'(0) & f'(0) & - 3f'(0) & f'(0) \\ f'(0) & f'(0) & f'(0) & - 3f'(0) \end{pmatrix}

Posons f'(0) = \alpha, on a alors:

P'(0) = \begin{pmatrix} - 3\alpha & \alpha & \alpha & \alpha \\ \alpha & - 3\alpha & \alpha & \alpha \\ \alpha & \alpha & - 3\alpha & \alpha \\ \alpha & \alpha & \alpha & - 3\alpha \end{pmatrix}

En multipliant, par exemple, la 1e ligne de P(t) avec la 2e colonne de P'(0) on obtient:

p'_{12}(t) = \alpha - 3\alpha f(t) - 3\alpha f(t) + \alpha f(t) + \alpha f(t)

p'_{12}(t) = \alpha - 4 \alpha f(t)

Or on sait aussi que si i \ne j, on a p'_{12}(t) = p'_{ij}(t) = f'(t).

C’est-à-dire: f'(t) = \frac{df}{dt} = \alpha - 4 \alpha f(t)

Et maintenant on intègre cette équation différentielle:

\int \frac{df}{\alpha - 4 \alpha f(t)} = \int dt

\Rightarrow \frac{ln(\alpha - 4 \alpha f(t))}{-4\alpha} = t + c

\Rightarrow ln( \alpha - 4 \alpha f(t) ) = - 4 \alpha t + c

\Rightarrow \alpha - 4 \alpha f(t) = c \exp^{- 4 \alpha t}

\Rightarrow f(t) = \frac{1}{4} - \frac{c}{4 \alpha} \exp^{- 4 \alpha t}

Il nous faut maintenant calculer c qui est la constance d’intégration. Pour cela supposons que f(0) = 0 ce qui signifie qu’au temps t=0 on commence dans un état constant. Par exemple, si à t=0 on est dans l’état A alors la probabilité d’avoir une substitution de ce « A » à t=0 vaut 0. Ainsi:

f(0)=0 \Rightarrow f(0) = \frac{1}{4} - \frac{c}{4\alpha} = 0 \Rightarrow c = \alpha

On a doncf(t) = \frac{1}{4} - \frac{\exp^{- 4 \alpha t}}{4}

Pour calculer \alpha on se place à t=0 puisque \alpha = f'(0), et on appelle \Pi_i la probabilité d’être dans l’état i au temps t=0 (\Pi_i = 1/4):

\sum_{i} \sum_{j \ne i} \Pi_{i} P'(0)_{ij} = 1

\Pi_{A} P'(0)_{AT} + \Pi_{A} P'(0)_{AG} + \Pi_{A} P'(0)_{AC} + \Pi_{T} P'(0)_{TA} + ... = 1

12 \frac{1}{4} \alpha = 1

\alpha = \frac{1}{3}

Finalement:

f(t) = \frac{1}{4} - \frac{exp(-4t/3)}{4}

C’est bien gentil vous allez me dire, on connaît maintenant f(t) mais ça ne résout pas notre problème initial qui était de prendre en compte les mutations cachées… En fait si mais il reste encore un peu de calcul à faire. Pour cela on doit estimer la valeur de la variable t, c’est-à-dire la distance qui sépare nos deux séquences S_1 et S_2. Grâce à notre modèle décrit ci-dessus cette distance est bien sûr reliée aux nombres de mutations observées entre les deux séquences tout en prenant en compte le fait que certaines mutations soient arrivées sans qu’on puisse les voir.

Afin d’estimer t on va appliquer la méthode du maximum de vraisemblance. Je rappelle que la vraisemblance (notée L pour likelihood) est la probabilité d’observer les données sachant le modèle: P(data/model). Dans notre cas on veut calculer la probabilité que la séquence S_1 (de taille n) ait pu évoluer en S_2, ce qui s’écrit:

L = P(S_{1}[1]) P(S_{1}[1] -> S_{2}[1]/t) P(S_{1}[2]) P(S_{1}[2] -> S_{2}[2]/t) ... P(S_{1}[n]) P(S_{1}[n] -> S_{2}[n]/t)

avec P(S_1[i]) la probabilité d’observer le nucléotide en question au i-ème site de la séquence S_1, et P(S_1[i] \rightarrow S_2[i]/t) la probabilité d’avoir muté au i-ème site du nucléotide de S_1 vers le nucléotide de S_2 pendant le temps t.

Quand on a un produit (multiplications) on aime bien le transformer en somme (additions). Pour cela on utilise la fonction logarithme:

ln L = ln P(S_1[1]) + ... + ln P(S_1[n]) + ln P(S_1[1] \rightarrow S_2[1]/t) + ... + ln P(S_1[n] \rightarrow S_2[n]/t)

Afin de trouver le maximum de cette vraisemblance on fait comme au lycée: on dérive la fonction ln L et on cherche les valeurs auxquelles la dérivée s’annule. Les n premiers termes de la formule valent une constante donc leur dérivée est nulle. Pour les n autres on peut poser que m_1 correspond aux mutations d’un nucléotide vers un autre (p_{ij}) et m_2 correspond aux mutations d’un nucléotide vers lui-même (p_{ii}).

\frac{d (ln L)}{dt} = \frac{m_{1}}{p_{ij}(t)} p_{ij}'(t) + \frac{m_{2}}{p_{ii}(t)} p_{ii}'(t) = 0

Comme on a calculé un peu plus haut f(t) et que f(t) = p_{ij}(t) on peut remplacer dans l’équation ci-dessus. Je vous épargne les calculs mais à la fin on obtient:

\hat{t} = -\frac{3}{4} ln ( 1 - \frac{4m_{1}}{3(m_{1}+m_{2})} )

Et pour résoudre enfin notre problème on peut définir p comme étant la proportion de sites différents entre nos deux séquences. Ainsi, alors qu’on a commencé par estimer la distance entre nos deux séquences par:

p = \frac{m_{1}}{m_{1}+m_{2}}

on estime maintenant cette distance par:

\hat{t} = -\frac{3}{4} ln ( 1 - \frac{4}{3} p )

Et c’est cela qu’on appelle la distance de Jukes-Cantor. Dans le cas de la figure 1, p=0.1 alors que \hat{t}=0.107. La distance de Jukes-Cantor est bien légèrement plus grande car elle prend en compte des mutations qui ont pu arriver mais qu’on ne voit pas.

Alors bien sûr, comme toujours en modélisation, on simplifie beaucoup, mais depuis l’article de Jukes Cantor en 1969 les modèles ont été perfectionnés et cela permet de bien mieux comprendre le génome des êtres vivants: vitesse d’apparition des mutations, importance fonctionnelle de certaines séquences, relations phylogénétiques entre les espèces… Mais ce serait trop pour ce billet !

ps: une bonne revue sur ces questions est disponible ici.

 

 

 

 

 

 

 

 

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion / Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion / Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion / Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion / Changer )

Connexion à %s

%d blogueurs aiment cette page :