Création d'un HMM pour prédire les promoteurs de type sigmaA de Bacillus subtilis


Installation du programme

Télécharger l'archive compressée dans le répertoire Software :

- pour le système d'exploitation Linux Debian DEBIAN
- pour le système d'exploitation Linux Fedora FEDORA

Une fois décompressée et désarchivée (tar -xvf SHOW.tar.gz), un répertoire show-debian-etudiant est créé dans lequel vous trouverez les différents sous-répertoires nécessaires au programme, deux fichier sigma_em et sigma_vit qui vous seront utiles par la suite ainsi qu'un sous-répertoire seqdna contenant un ensemble de séquences de régions promotrices de Bacillus subtilis qui vous permettront d'apprendre les paramètres de votre modèle HMM et de le tester.

Pour tester si cela marche :
Se mettre dans le répertoire show-debian-etudiant et lancer la commande suivante :

bin/show_emfit
Si vous obtenez une réponse commençant par Usage :.... C'est que c'est bon

Si vous avez l'erreur suivante :
"error while loading shared libraries libgsl.so.0",  c'est que le chemin pour utiliser cette librairie n'est pas reconnu. Il faut donc expliciter ce chemin. La librairie se trouve dans show-debian-etudiant/lib.
Faire pwd pour avoir le chemin complet de votre directory show-debian-etudiant
Ensuite sur le terminal taper les deux commandes suivantes:
LD_LIBRARY_PATH= le résultat de pwd/lib
export LD_LIBRARY_PATH
Refaire bin/show_emfit. Maintenant cela doit être bon

Les séquences de travail

Comme dit ci-dessus, elles se trouve dans le sous-répertoire seqdna  dans show-debian-etudiant.  Chaque séquence se trouve dans un fichier individuel en format Fasta ayant l'extension .dna. Sur la première ligne à la suite du nom de la séquence sont indiquées les positions connues des régions -35 et -10 des promoteurs sigma A.


Création du modèle HMM

En vous basant sur les logos ci-dessous établir le schéma de votre HMM où chaque position de la séquence correspondra à un état. Deux logos différents ont été établis pour représenter les régions -35 et -10 ainsi que le +1 de transcription pour une région promotrice de type sigma A normal (le premier) ou une région promotrice étendue (le second). Chaque logo a été obtenu en combinant 3 logos obtenus indépendamment (petite barre verticale), un pour la région -35, un pour la région -10 et un pour la région +1. Un seul HMM sera établi pour tenter d'identifier simultanément ces deux types de régions promotrices dans les séquences de Bacillus subtilis.  Pour la modélisation, il faudra déterminer la taille de chacune de vos régions conservées ainsi que la distance entre ces régions.
La distance entre la boîte -35 et -10 peut varier de 16 à 22 nucléotides et la distance entre la boîte -10 et le +1 de transcription peut varier de 4 à 10 nucléotides.

Logo


Vous allez utiliser SHOW développé à l'INRA de Jouy en Josas pour construire votre HMM. Pour cela il faut prendre connaissance de sa syntaxe. Vous avez à votre disposition la documentation fournie avec le programme ainsi qu'un extrait commenté de cette syntaxe. 

Modélisation du HMM et implémentation dans la syntaxe de SHOW
Pour vous faciliter le travail un template pour construire le modèle vous ai fourni. Il faudra le compléter. Ce template est disponible ici.
Une fois votre fichier modèle construit, il va falloir réaliser l'estimation des paramètres du modèle (probabilités d'émission et de transition) en utilisant les séquences qui vous ont été fournies. Pour cela, en plus du fichiers modèle, deux autres fichiers sont requis

Préparation des fichiers pour SHOW :

- Fichier séquences :

Les exécutables de SHOW peuvent travailler soit sur une seule séquence, soit sur un ensemble de séquences. Les séquences à analyser sont référencer dans le fichier –seq<file>. Ce fichier doit suivre la structure suivante :
seq_identifier: sigma.dna
seq_type: dna
seq_files:
aadk.dna
abr.dna
...
seq_identifier fait référence à un identifiant choisi pour la (les) séquence(s). Il doit être le même que celui donné dans le fichier modèle (-model) pour le mot clef seq au niveau de la description des observations.
seq_type correspond à la nature des séquences (pour le moment seulement adn)
seq_files correspond aux noms des fichiers contenant les séquences à analyser

Pour créer ce fichier dans le répertoire seqdna contenant vos fichiers séquences, faire :
ls *.dna sigma_seq.txt (créera un fichier appelé sigma_seq.txt et contenant la lite des noms de vos fichiers)
Ajouter au début de ce fichier, les trois lignes seq_identifier, seq_type et seq_files

- Fichier contenant les informations pour initialiser et faire tourner l’algorithme :
Fichier –em<file> Ce fichier sigma_em vous a été fourni dans l'archive. Il contient les lignes suivantes :
Les deux variables suivantes contiennent la taille des segments utilisée pour la mémoire sauvant les approximations lors de l’étape forward-backward de l’algorithme
estep-segment: 2000
estep_overlap: 100
nb_sel, niter_sel, eps_sel ne sont requises que si pobs vaut random
nb_sel: 3
niter_sel: 100
eps_sel: 0.01
niter: 1000
epsi: 0.001
niter et epsi: critère d’arrêt.
L’algorithme stop quand l’augmentation de la vraisemblance entre deux itérations est plus petite que la valeur de epsi ou quand le nombre maximal d’itérations défini par niter est atteint.
nb_sel correspond au nombre de point de départ aléatoire de l’algorithme EM
niter_sel se rapporte au nombre maximal d’itérations réalisé par point de départ
eps_selcorrespond au critère d’arrêt de l’algorithme pour chaque point de départ

Etape d'estimation des paramètres :

L'estimation des paramètres du modèle (probabilités d'émission et de transition) se fait en utilisant l'algorithme de Baum-Welch implémenté dans SHOW. Pour cela il va falloir lancer la commande suivante en se placant dans le répertoire contenant les fichiers séquences et votre fichier sigma_seq.txt donc dans seqdna.   :
../bin/show_emfit -model  ../sigma_model -em ../sigma_em -seq sigma_seq.txt > /dev/null

La commande > /dev/null va rediriger les warnings qui s'affichent à l'écran lorsque le programme tourne vers 'rien'.
Une fois le programme exécuté vous obtiendrez plusieurs fichiers. Celui contenant votre modèle avec les nouvelles probabilités estimées se nomme sigma_seq.model (si votre fichier séquence s'appelle sigma_seq.txt, sinon l'extension .model remplacera l'extension .txt). Vous le trouverez dans répertoire seqdna.

Etape de prédiction des positions des régions -35 et -10 en utilisant votre modèle :

Pour réaliser ces prédictions, l'algorithme de Viterbi implémenté dans SHOW va être utilisé. Nous allons ici tester notre modèle sur les mêmes séquences que celles utilisées pour l'apprentissage des paramètres. En général, d'autres séquences sont utilisées lors de l'étape test. Le fichier devra avoir le même format que sigma_seq.txt.
Il faut se mettre dans la directory contenant les fichiers séquences et votre fichier sigma_seq.txt. Il faudra lancer la commande suivante :
../bin/show_viterbi -model  ../sigma_seq.model -vit ../sigma_vit -seq sigma_seq.txt
le fichier sigma_vit vous a été donné dans l'archive pour l'installation de SHOW. Il contient les paramètres pour la gestion de la mémoire.

Le programme va créer un fichier .vit par séquence. L'entête du fichier récapitule les états de votre modèle HMM et attribue un numéro à chaque état, le premier recevant le chifrre 0, le second le chiffre 1 etc... Les lignes suivantes donnent la valeur numérique associée à l'état dans lequel la position de la séquence a été prédite.
Exemple d'un extrait d'un fichier .vit :
# viterbi reconstruction
# 0 : (bound)    1 : (background_1)    2 : (sigma-35+_1)    3 : (sigma-35+_2)    4 : (sigma-35+_3)    5 : (sigma-35+_4)....
1
1
1
1
2
3
4
5
6
7
8
9
...

La sortie n'est donc pas directement exploitable pour calculer les performances de prédiction du modèle, c'est à dire le pourcentage de fois où le modèle a correctement identifié les boîtes -35 et -10 des promoteurs sigma A (première position prédite de chacune de ces boîtes correspond bien aux positions connues dans les séquences analysées).

Pour pouvoir estimer le pouvoir prédictif du modèle, il va falloir parser ces ichiers.
Le programme vous permettant de parser les résultats du viterbi vous est fourni ici. Le sauvegarder dans la directory show_etudiant.  Si le fichier n'apparaît pas comme exécutable, changer son mode. Le programme prend comme argument en entrée le nom du répertoire contenant les fichiers séquences, donc ici seqdna. Un peu de souplesse a été acceptée pour considérer qu'une région avait été bien prédite. On accepte un écart < 3 entre position réelle et prédite. Le programme vous affichera les erreurs à l'écran, et donnera le pourcentage à la fin. Rediriger la sortie à l'écran dans un fichier pour garder les résultats.
Ce programme peut être largement amélioré et généralisé.