Semaine 8 - Annotation de gènes par chaînes de Markov Cachées
Exercice 1 : Rappels de biologie
Dans ce TME nous allons voir comment les modèles statistiques peuvent être utilisés pour extraire de l'information des données biologiques brutes. Le but sera de spécifier des modèles de Markov cachées qui permettent d'annoter les positions des gènes dans le génome.
Le génome, support de l'information génétique, peut être vu comme une longue séquence de caractères écrite dans un alphabet à 4 lettres: A , C , G et T . Chaque lettre du génome est aussi appelée pair de base (ou bp).
Il est maintenant relativement peu coûteux de séquencer un génome (quelques milliers d'euros pour un génome humain). Cependant on ne peut pas comprendre, simplement à partir de la suite de lettres, comment cette information est utilisée par la cellule (un peu comme avoir à disposition un manuel d'instructions écrit dans une langue inconnue).
Un élément essentiel est le gène , qui après transcription et traduction produira les protéines, les molécules responsables de la grande partie de l'activité biochimique des cellules.
La traduction en protéine est faite à l'aide du code génétique qui, à chaque groupe de 3 lettres (ou bp) transcrites fait correspondre un acide aminé. Ces groupes de 3 lettres sont appelés codon et il y en a {$4^3$}, soit {$64$}. Donc, en première approximation, un gène est défini par les propriétés suivantes (pour les organismes procaryotes):
- Le premier codon, appelé codon start est
ATGdans la majorité des cas (quelquefoisGTGouTTG). - Il y a 61 codons qui codent pour la séquence d'acides aminés.
- le dernier codon, appelé codon stop, marque la fin du gène et est l'une des trois séquences
TAA,TAGouTGA. Il n'apparaît pas dans le gène.
Voici un exemple coupé au milieu.
ATG.AGG.GAG.ATG.GCG.GATG.CCG.GAA-----CTG.CAT.GGT.AAG.AAT.CAT.TTG.GTC.CGT.GAA.TAA ### ### codon start codon stop
Nous allons intégrer ces différents éléments d'information pour prédire les positions des gènes. Notez que pour simplifier nous avons omis le fait au début que la molécule d'ADN est constituée de deux brins complémentaires, et donc que les gènes présent sur le brin complémentaire sont vus "à l'envers" sur notre séquence. Les régions entre les gènes sont appelées les régions intergéniques .
Exercice 2 : Modélisation de gènes
Question 1 : Téléchargement des données
Nous travaillerons sur le premier million de bp du génome de E. coli (souche 042). Les données sont préparées au format pickle à l'adresse suivante lien . Plutôt que de travailler avec les lettres A , C , G et T , nous allons les recoder avec des numéros (A =0, {$\ldots$}, T =3).
import pickle as pkl
data = pkl.load(file("genome_genes.pkl","rb"))
Xgenes = data.get("genes") #Les genes, une array de arrays
Genome = data.get("genome") #le premier million de bp de Coli
Annotation = data.get("annotation") ##l'annotation sur le genome
##0 = non codant, 1 = gene sur le brin positif
### Quelques constantes
DNA = ["A", "C", "G", "T"]
stop_codons = ["TAA", "TAG", "TGA"]
Chacune des séquences de gène commence par un codon start et fini par un des codons stop.
Out[31]: [0, 3, 2, 1, 2, 0]
#6 premieres bp: ATGCGA
Xgenes[0][-6:]
Out[32]: [2, 3, 1, 3, 0, 0]
#6 dernieres bp: GTCTAA
Question 2 : Un modèle séparant les codons et l'intergénique
Comme modèle le plus simple pour séparer les séquences de codons des séquences intergéniques, on va définir la chaîne de Markov caché dont le graphe de transition est donné ci dessous.
Un tel modèle se défini de la manière suivante. On peut rester dans les régions intergéniques, et quand on démarre un gène, la composition de chaque base du codon est différente. Il va falloir, afin de pouvoir utiliser ce modèle pour classifier, connaître les paramètres pour la matrice de transition (donc ici uniquement les probas {$a$} et {$b$}), et les lois {$(b_i, i=0,\ldots,3)$} des observations pour les quatre états.
Si {$a$} et {$b$} ont été fixés une telle chaîne de Markov se définit de la manière suivante en python.
Dans la suite du TME, deux options s'offrent à vous: tester la bibliothèque sklearn (bibliothèque renommée de machine learning en python) ou continuer dans l'univers des deux TME précédents (un code propre vous est fourni pour -re-partir d'une base saine). Les résultats intermédiaires attendus sont (évidemment) les mêmes quelque soit le modèle utilisé.
Dans l'univers sklearn (=scikit-learn)
| En restant dans les fonctions du TME précédent (fonctions mises au propre: lien) |
|
# import obligatoire from sklearn import hmm n_states_m1 = 4 # syntaxe objet python: créer un objet HMM model1 = hmm.MultinomialHMM(n_components = n_state_m1) Pi_m1 = np.array([1, 0, 0, 0]) ##on commence dans l'intergenique A_m1 = np.array([[1-a, a , 0, 0], [0 , 0 , 1, 0], [0 , 0 , 0, 1], [b , 1-b, 0, 0 ]]) # paramétrage de l'objet model1._set_startprob(Pi_m1) model1._set_transmat(A_m1) # [cf question d'après pour la détermination] model1._set_emissionprob(B_m1) |
from markov_tools import * Pi = np.array([1, 0, 0, 0]) A = np.array([[1-a, a , 0, 0], [0 , 0 , 1, 0], [0 , 0 , 0, 1], [b , 1-b, 0, 0 ]]) B = ... # pour utiliser le modèle plus loin: s, logp = viterbi(Genome,Pi,A,B) |
Ici B_m1 est une matrice à n_states_m1 lignes et 4 colonnes qui décrit toutes les lois d'observation (ligne 1: loi de l'intergénique, ligne 2: loi du codon position 0,...). On va voir dans la suite comment on l'estime.
Sous-question :
Estimez les probabilités {$a$} et {$b$} pour la matrice de transition. On utilise pour cela le fait que le temps de séjour dans un état suit une loi géométrique. Pour les régions intergénique, on fait l'hypothèse que la longueur moyenne est de 200 bp. Déduisez en la probabilité {$a$}
Sous-question :
Estimez la longueur moyenne des séquences des gène donnés dans Xgenes et déduisez en la probabilité {$b$}
Sous-question :
Estimez la distribution des quatre lettres sur tout le génome, pour la loi des observations de l'intergénique. Attention, on doit renvoyer une matrice de shape (1,4) .
Out[4]: array([[ 0.242654, 0.248799, 0.265139, 0.243408]])
Sous-question :
Estimez la distribution pour chacune des positions des codons et ranger le résultat dans une matrice à 3 lignes et 4 colonnes. Pour chaque gène on ignorerera le premier et le dernier codon.
Out[5]:
array([[ 0.25120572, 0.24063005, 0.35462469, 0.15353954],
[ 0.29302417, 0.23050332, 0.18223301, 0.2942395 ],
[ 0.18018394, 0.26191285, 0.29654511, 0.2613581 ]])
Sous-question :
On peut maintenant initialiser la matrice de transition pour notre modèle et l'utiliser pour prédire les positions des gènes.
model1._set_emissionprob(B_m1)
vsbce, pred = model1.decode(Genome)
#vsbce contient la log vsbce
#pred contient la sequence des etats predits (valeurs entieres entre 0 et 3)
#on peut regarder la proportion de positions bien predites
#en passant les etats codant a 1
sp = pred
sp[np.where(sp>=1)] = 1
percpred1 = float(np.sum(sp == Annotation) )/ len(Annotation)
percpred1
Out[10]: 0.636212
On classe un peu plus de {$60\%$} des positions correctement
Sous-question :
Afficher sur le même graphique la classification de la CMC avec l'annotation pour les 6000 premières positions du génome. On mettra en abscisse la position sur le génome, et 0 ou 1 suivant que la position appartient ou non à un gène. On utilisera une couleur pour l'annotation et une pour la prédiction.
Question 3 :
Evaluons maintenant si cela s'améliore en prenant en compte les frontières des gènes en construisant un modèle avec codon start et codon stop.
On veut maintenant d'intégrer l'information complémentaire qui dit qu'un gène commence "toujours" par un codon start et fini "toujours" par un codon stop avec le graphe de transition ci dessous.
Sous-question :
Ecrivez la matrice de transition correspondante, en mettant les probabilités de transition entre lettres pour les codons stop à 0.5.
Sous-question :
Adaptez la matrice des émissions pour tous les états du modèle. Les états correspondant au codons stop n'émettrons qu'une seule lettre avec une probabilité 1.
Pour le codon start, on sait que les proportions sont les suivantes:
ATG : 83% GTG: 14%, et TTG: 3%
Sous-question :
Evaluez les performances du modèle. On doit obtenir quelque chose du type.
Out[11]: 0.642634
L'amélioration n'est pas vraiment marquante (moins de 1%!), comment expliquez vous celà?
Sous-question :
Affichez à nouveau sur un graphique, la classification avec le modèle 1, le modèle 2, et l'annotation pour les 6000 premières positions. Comparez si les résultats sont différents sur certaines positions.
Question 4 :
(optionnel) Estimer les observations directement. Si on ne possède pas d'annotation, on peut aussi utiliser l'algorithme EM sur la structure de la matrice de transition pour estimer les lois d'émission.
params = "es") #l'estimation est faite pour les emission "e" et le start "s"
model3._set_transmat(A_m2)
model3.fit(G)
model3.transmat_ # afficher la matrice estimée
Comparez la matrice estimée pour les observations avec la matrice B_m1. Est ce que les deux matrices se ressemblent ? Comment pourrait on améliorer les résultat ?
Question 5 :
(optionnel) Un modèle avec le brin complémentaire. On veut maintenant de rajouter l'information pour la partie des gènes qui est codée sur le brin complémentaire (diagramme ci desssous).
Rappel brin complémentaire : l'information génétique est codé sur deux brins qui sont complémentaires au sens suivant :
- En face d'un
Ail y a unT - En face d'un
Cil y a unG
et réciproquement (T :A et G :T ).
Donc si on aurait le mini gène ATG.CTT.ATC.TAG sur le brin complémentaire, il sera inversé puis complémenté sur le brin positif et on verrait CTA.GAT.AAG.CAT . Explication :
CTA.GAT.AAG.CAT brin positif ||| ||| ||| ||| GAT.CTA.TTC.GTA brin complémentaire
Il faut donc transformer la matrice des observations sur les codons en une matrice pour les complémentaires inversés des codons, et des codons start et stop.
Sous-question :
Définissez le modèle correspondant et évaluez l'évolution des performances. Attention, l'annotation donnée est limitée, et les gènes sur le brin négatif sont annotés comme du non codant dans ce cas. Pour faire la comparaison vous passerez donc les régions prédites comme des gènes sur le brin négatif à la valeur 0.