Dans cette partie on va essentiellement travailler avec le logiciel R. Pour lancer le programme, il suffit de taper la commande suivante:
R
pour obtenir l'affichage (ce texte peut légèrement changer selon la version de R intallée):
R version 2.5.1 (2007-06-27) Copyright (C) 2007 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R est un logiciel libre livré sans AUCUNE GARANTIE. Vous pouvez le redistribuer sous certaines conditions. Tapez 'license()' ou 'licence()' pour plus de détails. R est un projet collaboratif avec de nombreux contributeurs. Tapez 'contributors()' pour plus d'information et 'citation()' pour la façon de le citer dans les publications. Tapez 'demo()' pour des démonstrations, 'help()' pour l'aide en ligne ou 'help.start()' pour obtenir l'aide au format HTML. Tapez 'q()' pour quitter R. >
qui se termine par l'invite >. Le programme est alors près a exécuter vos commandes.
Nous allons commencer par créer la matrice de transition d'une chaîne de Markov sur l'alphabet binaire {a,b}:
pi=matrix(c(0.4,0.6,0.2,0.8),ncol=2,nrow=2,byrow=T)
Pour voir la matrice pi ainsi construite il suffit de l'appeler pour obtenir:
[,1] [,2]
[1,] 0.4 0.6
[2,] 0.2 0.8
Il s'agit bien d'une matrice à deux lignes (ncol=2), deux colonnes (nrow=2) et dont les coefficient sont (par ligne car byrow=T): c(0.4,0.5,0.2,0.8). On peut facilement retrouver ces options dans l'aide de ligne de R qui est très bien faite. Il suffit pour cela de taper ?matrix (on quitte l'aide avec la touche Q)
On suppose que la chaîne de Markov commence toujours par un a c'est à dire que sa loi initiale est donnée par:
mu=c(1.0,0.0); print(mu);
Il est alors possible de calculer la loi marginale de la deuxième lettre de cette chaîne de Markov (à noter que l'opérateur %*% correspond au produit matriciel:
mu=mu%*%pi; print(mu);
En rappelant à plusieurs reprise cette ligne de commande (flêche vers le haut puis "Entrée"), on peut ainsi obtenir la loi marginale de la troisième lettre, de la quatrième, etc.
On constate que ce processus converge rapidement (en combien d'itérations ?) vers une valeurs de
Nous savons que cette loi stationnaire, ainsi que la vitesse de convergence de la loi marginale vers elle est déterminée par les valeurs propres et vecteurs propres de la matrice de transition. On peut facilement obtenir ces quantité via la commande (la fonction t() sert juste à transposer la matrice pi:
eigen(t(pi))
Si la plus grande valeur propre de la matrice est évidemment 1.0 (c'est une propriété des matrices stochastiques) c'est la seconde valeurs propre qui va nous parler de la vitesse de convergence:
nu=eigen(t(pi))$values[2]
Et pour obtenir la loi stationnaire, il nous suffit de normaliser un vecteur propre associé à la valeur propre 1.0:
v=eigen(t(pi))$vectors[,1]; v=v/sum(v);
Refaire le même travail (nombre d'itération avant convergence vers la loi stationnaire, loi stationnaire, seconde valeur propre) avec la matrice de transition définie par:
pi=matrix(c(0.90,0.10,0.05,0.95),ncol=2,nrow=2,byrow=T)
Proposer une matrice de transition convergeant au moins 10 fois moins vite vers sa loi stationnaire.
Compter tous les mots de longueur 2 dans la sequence ab1000.fasta en déduire une estimation de la matrice de transition d'une chaîne de Markov d'ordre 1. De manière similaire, donnée la première ligne d'une chaîne de Markov d'ordre 2 estimée sur cette séquence.
Faire le même travail sur la séquence d'ADN U00093.fasta.
On va maintenant utiliser la librairie seq++ pour vérifier nos calculs. La commande d'estimation markovienne est estim_m et s'utilise ainsi:
estim_m ab1000.fasta -A "1:ab" -d 1 -o ab_m1.markov estim_m ab1000.fasta -A "1:ab" -d 2 -o ab_m2.markov estim_m U00093.fasta -A "1:acgt" -d 1 -o dna_m1.markov estim_m U00093.fasta -A "1:acgt" -d 2 -o dna_m2.markov
Comparer le contenu des fichiers ainsi produits avec vos propres calculs.
On se pose maintenant la question de la variabilité de nos estimations. Pour étudier ce phénomène, on va générer un échantillon d'estimateurs et regarder se propriétés.
La commande simul_m permet de simuler facilement une chaîne de Markov dont les paramètres sont connues. Par exemple, la commande:
simul_m -A "1:ab" -m ab_m1.markov -l 1000 -o simu_ab.fasta
permet de simuler une chaîne de Markov de longueur 1000. Il nous est alors possible d'estimer les paramètres de notre chaîne de Markov à l'aide de cette séquence simulée pour obtenir une réalisation de l'estimateur:
estim_m simu_ab.fasta -A "1:ab" -d 1 -o simu_ab_m1.markov
Il ne reste plus qu'à répéter l'expérience un grand nombre de fois en stockant simplement les résultats (ici, uniquement la première ligne de la matrice de transition) pour obtenir l'échantillon souhaité:
#!/bin/bash rm -f sample.data i=0; while [ $i -lt 1000 ]; do let i=i+1; echo $i; simul_m -A "1:ab" -m ab_m1.markov -l 1000 -o simu_ab.fasta estim_m simu_ab.fasta -A "1:ab" -d 1 -o simu_ab_m1.markov grep -v "^#" simu_ab_m1.markov | head -n 1 >> sample.data done ;
On peut ensuite examiner les propriétés de cet échantillon avec les commandes R suivantes:
d=read.table("sample.data");
s=d[,1];hist(s,freq=F);
x=seq(range(s)[1],range(s)[2],length=200);
lines(x,dnorm(x,mean=mean(s),sd=sd(s)),lwd=2);
print(rbind(apply(d,2,mean),apply(d,2,sd)));
Quel est l'ordre de grandeur de l'écart-type de l'échantillon par rapport à la valeur de sa moyenne ? La variabilité d'estimation est-elle acceptable ?
Refaire le même travail en simulant des séquences de longueur 100, puis de longueur 10. Dans chaque cas, discuter les propriétés de l'estimateur (bais, normalité, écart-type).
Avec une approche similaire, comment pourrait-on faire pour savoir quel ordre de Markov il est raisonnable d'utiliser pour modéliser la séquence U00093.fasta ?
Plutôt que d'utiliser comme ci-dessus un échantillon d'estimateur pour choisir l'ordre de Markov d'un modèle, on va ici utiliser des critère d'ajustement de type vraisemblance pénalisée (ex: BIC ou AIC).
Les options -L, -I et -B de estim_m permettent justement de calculer (respectivement) la log-vraisemblance, le critère AIC et le critère BIC:
estim_m U00093.fasta -d 1 -L -I -B
Utiliser estim_m pour calculer les valeurs de ces trois quantités pour des modèles markovien d'ordre 1,2,3,4,5 et 6. Pour chaque critère quel est le meilleur modèle ?
On choisit arbitrairement d'utiliser le modèle sélectionné par le critère AIC. Mettre en oeuvre une simulation permettant d'obtenir la variabilité de la première transition de l'estimateur sachant que la longueur de la séquence U00093.fasta est 572017 (attention, pour éviter des temps de calculs trop élevés, on se limitera ici à une taille d'échantillon de 100).
Comparer à la variabilité d'un modèle d'ordre 6. Conclusion ?
On va maintenant utiliser les chaînes de Markov pour identifier les CDS du fichier L42023.fna qui contient le génome complet d'Haemophilus influenzae.
Pour cela, on commence par identifier toutes les sous-séquences comprises entre un codon start et un codon stop en phase. On va classer ses séquences en trois catégories en fonction de leurs longueurs L: non CDS si L <= 50, CDS si L >= 300 et indéterminées si 50 < L < 300.
Les commandes suivantes permettent d'effectuer cette tâche à l'aide du programme getorf du package EMBOSS:getorf -sequence L42023.fna -find 3 -minsize 50 -outseq noncds.fasta getorf -sequence L42023.fna -find 3 -minsize 300 -outseq cds.fasta getorf -sequence L42023.fna -find 3 -minsize 51 --maxsize 299 -outseq test.fasta
On va maintenant estimer des modèles markoviens (d'ordre 4) correspondant d'une part aux cds et d'autre part au non cds:
echo "> temp seq" > temp.fasta && grep "^>" -v cds.fasta >> temp.fasta estim_m temp.fasta -d 4 -o cds.markov echo "> temp seq" > temp.fasta && grep "^>" -v noncds.fasta >> temp.fasta estim_m temp.fasta -d 4 -o noncds.markov
Il est alors possible de calculer pour chaque sequence dont le statut est indéterminé, sa vraisemblance sous l'un et l'autre de ces modèles:
estim_m -l test.fasta -m cds.markov | sed "s/complete genome /=/" | cut -d "=" -f 2 > res_cds.txt estim_m -l test.fasta -m noncds.markov | sed "s/complete genome /=/" | cut -d "=" -f 2 > res_noncds.txt
Il ne reste alors plus qu'à lire les résultats et à les comparer pour obtenir le statut de chaque séquence indéterminée. Pour cela on se place sous R et on utilise les commandes suivantes:
cds=read.table("res_cds.txt"); cds=cds[[1]];
noncds=read.table("res_noncds.txt"); noncds=noncds[[1]];
sum(cds>noncds)
sum(cds<=noncds);
which(cds>noncds);