Bases Théoriques :
Une chaîne de Markov , ça vous dit quelque chose ? Une classe de modèles de Markov, appelée Modèle de Markov Caché (Hidden Markov Model, HMM), est un modèle mathématique permettant de segmenter un signal observé en régions (états cachés) définis par le modèle. Lequel est composé de quatre éléments : N états, une matrice d’émissions, des conditions initiales et une matrice de transition. Le rôle de la matrice de transition est de définir la probabilité de passer d’un état à un autre ou de le maintenir inchangé au cours d’un certain nombre d’observations. Quant à la matrice d’émissions, elle définit la probabilité d’émission d'un certain signal en fonction de l'état caché.
Compliqué ? Pas du tout. Voyez par exemple le cas d’un modèle à deux états. Ce type d’algorithme a fait ses preuves dans différents domaines tels que la reconnaissance vocale (Juang Rabiner, Technometrics 1991), la détection de domaines topologiques issus de données de séquençage High-C (Dixon,Nature 2012), et dans le « foot printing » de données DNAse1-seq (Boyle, Genome Res 2011) entre autres.
Nous allons ici détecter des patterns dans des séquences d’ADN à l’aide d’une librairie R nommée « HMM » et apprendre par l’exemple certaines des fonctions de cette librairie.
ISOCHORE dans la ligne de mire
Objectif : détecter automatiquement une région ISOCHORE riche en GC de 10 paires de bases (pb) dans une séquence d’ADN (ATGC) d’une longueur de 100 pb. Un ISOCHORE est défini par un contenu en GC de 50% et un contenu en AT de 50%. Dans le reste du génome (ou de notre séquence à analyser) le contenu en GC est de 20% et le contenu en AT est de 80%.
Première étape : définir le modèle. À savoir, le nombre d’états, la matrice d’émission et la matrice de transition. Dans le cas présent, deux états : ISOCHORE et NON-ISOCHORE. Sachant que la longueur d’un isochore est de 10 pb environ et que la séquence analysée est de 100 pb, voici une définition de la matrice (de dimensions 2x2) de transition entre les deux états :
| matrice de transition | ISOCHORE | NON-ISOCHORE |
| ISOCHORE | 9/10 (1-p) | 1/10 (p) |
| NON-ISOCHORE | 1/90 (q) | 89/90 (1-q) |
Pour cet exemple, la probabilité de rester dans l’état ISOCHORE est de 9/10 ; en effet, un isochore fait environ 10 pb de long et la probabilité de passer de l’état ISOCHORE à l’état NON- ISOCHORE est de 1-9/10. De plus, sachant que la longueur totale de la séquence est de 100 pb, 90 observations sur 100 sont dans l’état NON-ISOCHORE. Ainsi, la probabilité de passer de l’état NON-ISOCHORE à l’état ISOCHORE est de 1/90= 1/(100-10) et la probabilité de rester dans l’état NON-ISOCHORE de 1-1/90 =89/90.
La matrice d’émissions du modèle est donc définie comme suit :
| Matrice d’emmission |
ISOCHORE |
NON-ISOCHORE |
|
A |
0.25 |
0.4 |
|
T |
0.25 |
0.4 |
|
G |
0.25 |
0.1 |
|
C |
0.25 |
0.1 |
Elle permet de déterminer la probabilité d’une observation donnée selon l’état caché. En ce qui concerne ISOCHORE, nous avons 50% de GC et 50% de AT, donc une probabilité égale de 25% pour chaque observation A, T ,G,C. Dans le cas de l’état NON-ISOCHORE, un biais apparaît, avec 80% de AT donc 40% de A, 40% de T, 10% de G et 10% de C.
Enfin, dans ce modèle, la probabilité initiale d’être dans l’un ou l’autre des états a été définie à 90% pour l’état NON-ISOCHORE et à 10% pour l’état ISOCHORE.
Le modèle complet est représenté dans la figure suivante.
Les miracles de la librairie HMM pour R
Deuxième étape : installer la librairie HMM pour R. Pour cela, utilisez les commandes suivantes :
|
1 2 |
install.packages("HMM") library(HMM) |
Puis définissez le modèle décrit plus haut :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
states < ;- c("N","I") emissions < ;- c("C","G","A","T") initProb < ;- c(0.9, 0.1) transProb < ;- t(matrix(c(89/90, 0.1, 1/90, 0.9), 2)) emissionProb < ;- matrix(c(0.1, 0.25, 0.1, 0.25, 0.4, 0.25, 0.4, 0.25), 2) hmm < ;- initHMM(states, emissions, initProb, transProb, emissionProb) hmm # imprime un résumé du modèle |
Ensuite, il faut simuler. Soit créer une simulation de la séquence observée à partir du modèle ou alors analyser une séquence d’intérêt puis détecter la région ISOCHORE. Marche à suivre pour créer une séquence d’ADN de longueur 100 bp contenant une ISOCHORE :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
# Séquence à analyser observation < ;- c("A", "T", "G", "A", "A", "A", "C", "A", "A", "T", "T", "T", "T", "T", "A", "T", "A", "A", "C", "T", "A", "T", "A", "T", "A", "T", "G", "T", "A", "C", "A", "T", "A", "T", "T", "T", "T", "T", "A", "T", "T", "A", "A", "T", "T", "A", "G", "C", "G", "A", "C", "G", "C", "C", "G", "C", "A", "C", "T", "A", "G", "A", "A", "T", "A", "T", "A", "A", "T", "T", "C", "A", "T", "T", "A", "T", "G", "T", "T", "T", "A", "T", "T", "G", "A", "A", "A", "A", "C", "A", "T", "A", "C", "T", "T", "G", "A", "A", "A", "A") # Simulation à partir du modèle simulation < ;- simHMM(hmm,100) simulation # imprime le résultat de la simulation # notamment la séquence d’observations générée par le modèle # et les états cachés pour chaque observation. |
Pour analyser une séquence d’observations non simulée par le modèle, vous pouvez détecter les états cachées avec la fonction Viterbi de la librairie HMM.
|
1 2 3 |
viterbi < ;- viterbi(hmm, observation) viterbi # affiche les états cachés les plus probables selon la séquence analysée et le modèle. |
Si vous n’y voyez pas clair, c’est normal. Pour vous aider :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
obsToVal < ;- function(obs){ if((obs == "A") | (obs=="T")){return(2)} else {return(4)} } stateToVal < ;- function(st){ if(st == "N"){return (1)} if(st == "I"){return (3)} else {return ("NA")} } par(mfrow=c(2,1)) plot(sapply(observation, obsToVal),pch=19,ylim=c(0,8), col=sapply(observation, obsToVal),cex=.5) legend("topright",c("A/T","G/C"),pch=19,col=c("red","blue")) plot( sapply(viterbi, stateToVal),type="l",ylim=c(0,5)) |
Voici ce que vous devez obtenir :
Le panel supérieur montre la localisation des G/C ou des A/T le long de la séquence de 100 pb. Quant au panel inférieur, il indique l’analyse du modèle par rapport à la séquence. Le HMM a détecté un ISOCHORE entre les bases 44 et 58.
Pour les probabilités en marche avant et en marche arrière, c’est par ici :
|
1 2 3 4 5 6 7 |
logForwardProb < ;- forward(hmm, simulation$observation) exp(logForwardProb) logBackwardProb < ;- backward(hmm, simulation$observation) exp(logBackwardProb) |
Enfin, cette librairie permet d’utiliser l’algorithme de Baum-Welch pour améliorer les performances du modèle et ajuster les paramètres à un test set. Voici comment :
|
1 2 3 |
bw = baumWelch(hmm, observation, 10) # 10 iterations print(bw$hmm) |
Ami-e-s en quête des mystères des HMM, j’espère avoir éclairé votre lanterne avec ce bref tutoriel dont j’ai volontairement retiré les détails algorithmiques pour ne pas aggraver votre mal de tête ! Et si le Graal vous hante, ne renoncez pas à vous en emparer en fouillant dans la littérature existante sur les algorithmes de Viterbi, « Forward » , « Backward » , et autres joyeusetés telles que Baum-Welch.
Que la force soit avec vous !
Un grand merci à nallias, Estel, muraveill, Aline Jaccottet et Hatice Akarsu pour leur relecture.
Laisser un commentaire