Utiliser les commandes suivantes pour produire un NFA associé au motif a.ab sur l'alphabet binaire {a,b}:
# spatt génère l'automate au format graphviz (.dot) spatt -a "ab" -p "a.ab" --nfa nfa.dot # graphviz produit un fichier image au format eps dot -Tps -Grankdir=LR nfa.dot -o nfa.eps # ghostview affiche le fichier image gv nfa.eps
A noter qu'il est possible d'enchainer ces trois commandes sur en une seule ligne:
# "&&" veut dire "et si la commande précédente a bien fonctionné" spatt -a "ab" -p "a.ab" --nfa nfa.dot && dot -Tps -Grankdir=LR nfa.dot -o nfa.eps && gv nfa.eps
De la même façon, il est évidemment possible de construire le DFA minimal associé à ce motif:
spatt -a "ab" -p "a.ab" --dfa dfa.dot && dot -Tps -Grankdir=LR dfa.dot -o dfa.eps && gv dfa.eps
Toujours sur l'aphabet binaire, donner pour chacun des quatre motifs suivants les caractéristiques (nombre d'états et nombre d'états terminaux) des NFA et DFA associés: a.(1)ab, a.(2)ab, a.(3)ab et a.(4)ab.
Avec ce dernier motif, on remarque que spatt n'affiche pas les détails de l'automate lorsque le nombre d'états dépasse 50. On peut néanmoins afficher les caractéristiques de l'automate même dans un tel cas à l'aide de la commande suivante:
spatt -a "ab" -p "a.(4)ab" -V | grep nstates | tail -n 1
Utiliser cette astuce pour poursuivre l'analyse sur les motifs
Créer un fichier seq_ab.fasta contenant la sequence suivante (au format FASTA):
> my_simple_binary_sequence abbabbabaaabbbababbbababaaabbababbbababaababbbaaaa babbabbbbababaaaabbabbababbabbaaabbabbbbabababbbab
Que fait la commande ci-après ?
spatt -a "ab" -p "abbba" -S seq_ab.fasta -m -1 --repartition | grep "1$" | cut -f 1
Et celle-ci ?
spatt -a "ab" -p "abbba" -S seq_ab.fasta -m -1 --repartition | grep "1$" | cut -f 1 | wc --lines
Utiliser des commandes similaires pour localiser et compter toutes les occurrences dans seq_ab.fasta des motifs suivants: aba et aba.(1-3)aba.
Créer un fichier pattern_count.pl contenant:
#!/usr/bin/perl
use Bio::Seq;
use Bio::SeqIO;
$pattern='aba.{1,3}aba';
$file='seq_ab.fasta';
$seqio_obj = Bio::SeqIO->new(-file => $file, -format => 'Fasta');
while ($seq_obj = $seqio_obj->next_seq) {
$seq = $seq_obj->seq;
$title=$seq_obj->display_name;
$count=0;
while ($seq=~/(?=$pattern)/g) {
print "$title: $pattern starts at position $-[0]\n";
$count++;
}
print "Total count of $pattern in $title is $count\n";
}
Rendre le script exécutable, puis le lancer:
chmod +x pattern_count.pl ./pattern_count.pl
Que fait ce script ? Quelles sont les différences avec spatt ?
Chercher sur le site de PROSITE la signature PDOC01004. La convertir au format des expressions régulières de perl à l'aide du script prosite2regex.pl via la commande:
prosite2regex.pl PROSITE_SIGNATURE
Modifier enfin
le script pattern_count.pl pour trouver les occurrences
de cette signature dans le fichier
sample_uniprot_sprot.fasta.
Refaire le même travail avec la signature
créer un fichier stree.pl contenant les lignes suivantes:
#!/usr/bin/perl
use Tree::Suffix;
$tree = Tree::Suffix->new('mississippi');
$tree->dump();
Exécuter ce script perl en redirigeant sa sortie dans le fichier mississippi.stree.
Si le script ne rend pas la main, appuyer sur
chmod +x stree.pl ./stree.pl > mississippi.stree ^C
Le fichier mississippi.stree contient maintenant le descriptif du suffix tree associé au texte 'mississipi' au format stree (vous pouvez examiner son contenu avec cat mississipi.stree par exemple). Comme dans le cas des automates, nous voulons utiliser le programme dot pour visualiser ce suffix tree. Pour cela, il nous faut donc convertir cette description en un format lisible par ce programme.
Enregistrer le fichier source C++ stree2dot.cc puis effectuer les commandes suivantes:
# compilation du fichier source g++ stree2dot.cc -o stree2dot # conversion .stree > .dot ./stree2dot mississippi.stree > mississippi.dot # production et affichage de la représentation graphique dot -Tps -Grankdir=LR mississippi.dot -o mississippi.eps && gv mississippi.eps
A l'aide de cet arbre, localiser les position de début des occurrences de 'ss' dans 'mississippi', puis déterminer la plus longue sous-chaine répétée.
Modifier le script stree.pl en:
#!/usr/bin/perl
use Tree::Suffix;
$tree = Tree::Suffix->new('mississippi');
@loc=$tree->find('ss');
foreach (@loc) {
$s=$$_[1];
$e=$$_[2];
print "\'ss\' starts in $s and ends in $e\n";
}
$rep=$tree->lrs;
print "longuest repeated substring = \'$rep\'\n";
Que fait la nouvelle version du script ?
S'inspirer des scripts pattern_count.pl et stree.pl pour produire un nouveau script localisant les occurrences de abba et la plus longue séquence répétée dans seq_ab.fasta à l'aide d'un arbre de suffixes.
On s'intéresse maintenant au problème suivant: comment obtenir le comptage de tous les mots d'une longueur donnée dans une séquence ? On va examiner ici plusieurs réponses à cette question et discuter des avantages et inconvénients de chacunes d'entre elles.
On commence par l'approche "force brute" consistant à compter successivement tous les mots d'ADN de longueur 5 dans la séquence U00093.fasta (il s'agit du chromosome VIII de la levure) à l'aide d'automates.
Créer un script bf.pl avec le contenu suivant:
#!/usr/bin/perl
use Bio::Seq;
use Bio::SeqIO;
$seqio_obj = Bio::SeqIO->new(-file => 'U00093.fasta', -format => 'Fasta');
$seq_obj = $seqio_obj->next_seq;
$seq = $seq_obj->seq;
foreach $token5 ('a','c','g','t') {
foreach $token4 ('a','c','g','t') {
foreach $token3 ('a','c','g','t') {
foreach $token2 ('a','c','g','t') {
foreach $token1 ('a','c','g','t') {
$word="$token5$token4$token3$token2$token1";
$count=0;
while ($seq=~/(?=$word)/g) {
$count++;
}
print "$word $count\n";
}
}
}
}
}
Exécuter ce script en mesurant sont temps d'exécution:
chmod +x bf.pl time ./bf.pl
On va maintenant s'intéresser au même problème en utilisant un arbre de suffixe. On commence par créer cet arbre, puis on s'en sert pour compter tous les mots. Mesurer le temps d'éxécution du script st.pl ayant le contenu suivant:
#!/usr/bin/perl
use Bio::Seq;
use Bio::SeqIO;
use Tree::Suffix;
$seqio_obj = Bio::SeqIO->new(-file => 'U00093.fasta', -format => 'Fasta');
$seq_obj = $seqio_obj->next_seq;
$seq = $seq_obj->seq;
$tree = Tree::Suffix->new($seq);
foreach $token5 ('a','c','g','t') {
foreach $token4 ('a','c','g','t') {
foreach $token3 ('a','c','g','t') {
foreach $token2 ('a','c','g','t') {
foreach $token1 ('a','c','g','t') {
$word="$token5$token4$token3$token2$token1";
$count=$tree->find($word);
print "$word $count\n";
}
}
}
}
}
On considère enfin deux approches fondée sur des arbres de préfixes: le programme wordcount du package EMBOSS (implémentation fondée sur des arbres) et le programme sspatt du package SPatt serie 1.x (implémentation utilisant de simples tableaux).
Voici la ligne de commande permettant d'utiliser wordcount et de mesurer son temps d'exécution (notons au passage que se programme n'affiche pas les mots de comptage nul et ordonne les mots par fréquences décroissantes):
time wordcount -sequence U00093.fasta -wordsize 5 -outfile stdout
Et voici enfin la commande ayant le même effet pour sspatt:
time sspatt U00093.fasta -a acgt -l 5 --all-words
Afin de pouvoir comparer de manière plus systématique les performances respectives de ces quatre approches, enregistrer localement les scripts perl: bfcounf.pl et stcounf.pl qui reprennent simplement les scripts bf.pl et st.pl en permettant de spécifier sur la ligne de commande le fichier séquence concerné, l'alphabet à utiliser et la longueur de mot. Exemples:
./bfcount.pl U00093.fasta acgt 2 ./stcount.pl U00093.fasta acgt 2
Pour h variant de 1 à 12, mesurer le temps d'exécution de chacune de les performances des ces quatres approches pour compter tous les mots d'ADN de longueur h dans U00093.fasta. A noter qu'on pourra éviter les affichages intempestifs à l'écran en redirigeant les résultats de ces commandes dans /dev/null. Que peut-on conclure sur les complexités respectives de ces différentes approches ?