Cours

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 ATG dans la majorité des cas (quelquefois GTG ou TTG ).
  • 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 , TAG ou TGA . 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 numpy as np
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.

Xgenes[0][:6]
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) .

Binter
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.

Bgene
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.

B_m1 = np.vstack((Binter, Bgene))

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.

percpred2
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.

model3 = hmm.MultinomialHMM(n_s_m2, n_iter = 100, thresh = 1e-6,   #les paramètres de la  
                            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 A il y a un T
  • En face d'un C il y a un G

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.

Quelques références supplémentaires

  • SHoW un système de CMC pour la prédiction de gènes dans les génomes bactériens : lien
  • Une thèse célèbre sur la prédiction de gènes chez les eucaryotes lien