On se propose ici d'estimer la p-value P(N>=Nobs) d'un motif à l'aide de simulations. On commencer pour cela par compter le nombre d'occurrence de aaba dans la sequence ab1000.fasta:
sspatt ab1000.fasta -a "ab" -p "aaba"
On va ensuite obtenir une réalisation de la variable aléatoire N du nombre d'occurrences du motif aaba dans une chaîne de Markov (stationnaire) de longueur 1000 dont les paramètres sont donnés dans ab_model1.markov:
# on génére la sequence simulée simul_m -m ab_model1.markov -A "1:ab" -l 1000 -o sim_ab1000.fasta # on compte le nombre d'occurrences du motif sur cette séquence sspatt sim_ab1000.fasta -a "ab" -p "aaba"
Pour obtenir un échantillon de N il suffit de répéter les deux étapes précédentes à l'aide du script suivant (attention, ce script est une boucle sans fin, le laisser tourner une vingtaine de secondes puis l'interrompre avec CTRL-C):
rm -f sample.data while true ; do simul_m -m ab_model1.markov -A "1:ab" -l 1000 -o sim_ab1000.fasta ; sspatt sim_ab1000.fasta -a "ab" -p "aaba" >> sample.data ; done; ^C
Il ne reste alors plus qu'a charger l'échantillon dans R et à examiner quelque unes de ces caractéristiques:
s=read.table("sample.data");s=s[,2];mean(s>=16);
hist(s,nc=30,freq=F);
x=seq(range(s)[1],range(s)[2],length=100);
lines(x,dnorm(x,mean=mean(s),sd=sd(s)),lwd=2);
mean(s);sd(s);
mean(s>=19);mean(s>=20);mean(s>=21);
mean(s<=3);mean(s<=4);mean(s<=5);
Peut-on qualifier cette distribution de gaussienne ? Quelle p-value obtient-on à l'aide de cet échantillon ? Est-elle significative ? Quelle p-value aurait-on obtenue avec Nobs=3 ? Et avec Nobs=21 ?
Reprendre maintenant la même démarche mais en utilisant cette fois-ci les paramètres du fichier ab_model2.markov pour les simulations:
rm -f sample.data while true ; do simul_m -m ab_model2.markov -A "1:ab" -l 1000 -o sim_ab1000.fasta ; sspatt sim_ab1000.fasta -a "ab" -p "aaba" >> sample.data ; done; ^C
Peut-on qualifier cette distribution de Gaussienne ? Quels espérance et écart-type trouve-t-on pour N ? Quelle est la p-value de l'observation ? Conclusion ?
On va maintenant s'intéresser à la PMC associée au comptage d'un motif. Commençons par les occurrences du motif aaba sur l'alphabet binaire dans une chaîne de Markov d'ordre 1 dont les paramètres sont dans le fichier ab_model1_short.markov créé à partir de ab_model1.markov via:
head -n -1 > ab_model1_short.markov head -n -1 > ab_model2_short.markov
(On enlève en fait la distribution stationnaire du fichier ab_model1.markov pour se conformer au format de fichier modèle utilisé par spatt). Comme on l'a déjà vu dans le premier TP, il est possible d'obtenir les l'automate associé au motif via:
spatt -a "ab" -p "aaba" --dfa dfa.dot dot -Tps -Grankdir=LR dfa.dot -o dfa.eps gv dfa.eps
Mais spatt peut également générer un fichier pmc.sci contenant la matrice de transition de la PMC associée au format Scilab à l'aide de la commande:
spatt -a "ab" -p "aaba" -m 1 -M ab_model1_short.markov
Afin d'examiner et de manipuler cette matrice de transition, nous allons simplement la charger dans Scilab à l'aide de la commande:
scilab -nw -f pmc.sci
On entre alors dans le mode intéractif de Scilab (que l'on peut quitter à tout instant à l'aide de la commande quit). Dans ce mode, examiner et commenter le résultat des commandes suivantes:
L starts final full(P) full(Q)
(remarquons au passage que la commande full() est nécessaire pour afficher de manière lisible le contenu des matrices creuses P et Q).
Calculons la loi stationnaire de la PMC à l'aide d'une méthode aussi simple qu'efficace:
Pi=full(P+Q); Pi^1000
Avec quelle probabilité sommes-nous dans un état terminal sous la loi stationnaire ?
Comparer cette probabilité à mu(a)P(a|a)P(b|a)P(a|b) (où mu est la loi stationnaire et P(.|.) la matrice de transition du fichier ab_model1.markov). Conclusion ?
Utiliser une approche similaire pour calculer la probabilité d'occurrence du motif dans une chaîne de Markov stationnaire de paramètres ab_model2.markov.
On reprend la PMC du motif aaba dans une chaîne de Markov d'ordre 1 de paramètres ab_model1.markov. On va tout d'abord calculer l'espérance exacte du nombre d'occurrences dans une séquence de longueur 1000 dont le premier terme est x1=a à l'aide de Scilab (que l'on charge avec scilab -nw -f pmc.sci):
tic();
Pi=full(P+Q);
mu1=zeros(1,L);mu1(starts(1))=1;
eF=zeros(1,L);eF(final)=1;
esp=0;
for i=1:1000;
esp=esp+mu1*Pi^(i-1)*eF';
end;
printf("esperance=%.6f (in %f s)\n",esp,toc());
On fait maintenant le même travail mais en utilisant cette fois-ci un algorithme récursif et les matrices creuses:
tic();
Pi=P+Q;
mu1=zeros(1,L);mu1(starts(1))=1;
eF=zeros(1,L);eF(final)=1;
aux=mu1;
esp=0;
for i=1:1000;
esp=esp+aux*eF';
aux=aux*Pi;
end;
printf("esperance=%.6f (in %f s)\n",esp,toc());
Conclusion ?
Avec spatt on peut obtenir directement l'espérance et l'écart-type exacts avec la commande:
spatt -a "ab" -p "aaba" -m 1 -M ab_model1_short.markov -S ab1000.fasta --gaussian
Comparer les valeurs exactes de l'espérance et de la variance au valeurs empiriques obtenues au début du TP.
On peut égalemement utiliser spatt pour obtenir la distribution exacte de N dans la séquence ainsi que la p-value exacte:
spatt -a "ab" -p "aaba" -m 1 -M ab_model1_short.markov -S ab1000.fasta --over
On peut aussi obtenir la p-value exact (mais cette fois-ci sous l'hypothèse de stationnarité) via xspatt:
xspatt -a "ab" -p "aaba" -m 1 -U ab_model1_short.markov ab1000.fasta
Cette fois-ci c'est cependant le -log10 de la p-value (+log10 dans le cas d'un mot sous-représenté) qui est retourné (vérifier la proximité du -log10 de la p-value calculé par spatt et cette valeur).
Notons également qu'il est facile de (re-)calculer exactement la probabilité P(N=0) à l'aide de script Scilabsuivant:
tic();
P=full(P);
mu1=zeros(1,L);mu1(starts(1))=1;
p=sum(mu1*P^(1000-1));
printf("P(N=0)=%.6e (in %f s)\n",p,toc());
Utiliser spatt pour calculer la p-value du motif aaba dans la sequence ab1000.fasta avec le modèle de référence ab_model2.markov. Ce résultat est-il cohérent avec le résultat des simulations faites en début de TP ?
On va maintenant comparer les résultats exacts aux résultats approchés (approximations gaussiennes et grandes déviations). Pour cela, on commence par calculer les p-values de tous (en échelle logarithmique décimale) les mots de longueur 8:
xspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model1_short.markov > x.res gspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model1_short.markov > g.res ldspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model1_short.markov > ld.res ldspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model1_short.markov --precise > ldp.res
On peut alors combiner ces différents fichiers résultats en un seul fichier:
join x.res g.res | cut -d " " -f 1,2,4,8 > aux1.res; join aux1.res ld.res | cut -d " " -f 1-4,7 > aux2.res; join aux2.res ldp.res | cut -d " " -f 1-5,8 > all.res;
Et il ne reste alors plus qu'à lire ce fichier dans R pour pouvoir comparer les différentes approches
d=read.table("all.res");x=d[,3];g=d[,4];ld=d[,5];ldp=d[,6];
plot(x,g); abline(0,1);
plot(x,ld); abline(0,1);
plot(x,ldp); abline(0,1);
c(mean(abs(x-g)),mean(abs(x-ld)),mean(abs(x-ldp)));
c(max(abs(x-g)),max(abs(x-ld)),max(abs(x-ldp)));
Quels sont les qualité et les défauts de ces différentes approximations ?
On reprend maintenant la même démarche en utilisant cette fois-ci le deuxième modèle:
xspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model2_short.markov > x.res gspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model2_short.markov > g.res ldspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model2_short.markov > ld.res ldspatt ab1000.fasta -a "ab" -l 8 --all-words -m 1 -U ab_model2_short.markov --precise > ldp.res
Que peut-on alors dire des différentes approximations ?
Nous allons enfin aborder l'intéressante question de la dépendance des statistiques de motifs à l'estimation des paramètres du modèle. Pour cela, on commence par estimer un modèle d'ordre 1 sur U00093.fasta qui va constituer notre modèle de référence:
estim_m U00093.fasta -d 1 -o dna.markov
Il est alors possible de calculer des statistiques de motifs par rapport à ce modèle de référence:
head -n -1 dna.markov > dna_short.markov ldspatt U00093.fasta -a "acgt" -p "tcatcatc" -m 1 -U dna_short.markov
Mais que se passe-t-il si on ne dispose que d'un estimateur des paramètres de notre modèle ? Pour le découvrir, il suffit de simuler un échantillon de cet estimateur et d'utiliser les paramètres correspondant pour calculer la significativité de notre motif:
rm -f sample.data i=0; while [ $i -lt 100 ]; do let i=i+1; echo $i; simul_m -m dna.markov -l 562639 -o simu.fasta estim_m simu.fasta -d 1 -o simu.markov head -n -1 simu.markov > simu_short.markov ldspatt U00093.fasta -a "acgt" -p "tcatcatc" -m 1 -U simu_short.markov >> sample.data done ;
On peut alors facilement examiner la distribution de notre statistique de motif avec R:
s=read.table("sample.data");s=s[,4];
hist(s); mean(s); sd(s);
En tenant compte de la sensibilité à l'estimation des paramètres, quelle marge d'erreur avons nous sur la p-value du motif considéré ?
On va maintenant reprendre cette démarche en utilisant cette fois-ci un modèle d'ordre 4. Pour garder des statistiques d'ordres de grandeur comparables, on change cependant au passage de motif d'intérêt:
estim_m U00093.fasta -d 4 -o dna.markov head -n -1 dna.markov > dna_short.markov ldspatt U00093.fasta -a "acgt" -p "tatatata" -m 4 -U dna_short.markov rm -f sample.data i=0; while [ $i -lt 100 ]; do let i=i+1; echo $i; simul_m -m dna.markov -l 562639 -o simu.fasta estim_m simu.fasta -d 4 -o simu.markov head -n -1 simu.markov > simu_short.markov ldspatt U00093.fasta -a "acgt" -p "tatatata" -m 4 -U simu_short.markov >> sample.data done ;
Quelle marge d'erreur obtient-on alors ? Que peut-on conclure ?