Introduction
Alignement de deux séquences
Alignement d'une séquence sur une base de données
Alignement multiple
Phylogénie
Phylogénie : méthode de Parsimonie
Calculs de la distance entre deux séquences
Méthodes de distances
Notes sur l'implémentation des algorithmes de reconstruction phylogénétique
Maximum de vraissembalnce et phylogénie
Simulation de Monte-Carlo à l'aide de chaines de Markov (MCMC)
Recherche de mot dans un texte
Recherche de mot dans un texte : indexation d'un texte
Chaines de Markov cachées
Strucure secondaire de l'ARN et grammaires stochastiques (SCFG)
Bibliographie (livres)
Dans ce document, nous présentons les principaux algorithmes utilisés en bioinformatique : il s'agit surtout des manipulations de séquences (comparaison (ou alignement) de deux séquences, comparaison d'une famille de séquences), de reconstruction phylogénétique (construction d'un arbre qui décrive le mieux les données dont on dispose), de manipulation de chaines de caractères (recherche dans un texte) ou de texte (indexation d'un texte ou d'un génome). Au cours de notre voyage dans le domaine de la bioinformatique, nous constaterons que la notion de graphe est omniprésente (arbres phylogénétiques, chaines de Markov, automates probabilistes, chaines de Markov cachées, SCFG, etc.) ; nous découvrirons aussi un peu de statistiques (tests statistiques, p-valeur, estimateur du maximum de vraissemblance, méthodes bayesiennes, simulation de Monte-Carlo à l'aide de chaines de Markov, etc.)
(Il manque quelques illustrations -- je ne suis pas vraiment doué en dessin --, mais j'espère que c'est quand-même compréhensible.)
(Il y a aussi quelques passages qui mériteraient d'être (ré)écrits (en particulier celui sur les E-valeurs) -- mais comme les journées n'ont que 24 heures, ils resteront comme ça.)
Alignement simple Dotplot Programmation dynamique : alignement global Variantes : Alignement local, semi-global Matrices BLOSUM et PAM Heuristiques : fasta, blast Interprétation statistique du résultat Alignement multiple Scores : SP, entropie minimale, arbre Programmation dynamique, MSA Étoile UPGMA (phylogénétique rudimentaire) Multalign Consensus, PSSM Chaines de Markov Introduction HMM (Modèles de Markov cachés) Algorithme forward Algorithme de Viterbi Algorithme Forward-Backward, aka Baum-Welch, aka EM HMM pour décrire un profil HMM pour repérer introns et exons Phylogénétique Méthodes de distance : UPGMA, NJ Parcimonie Maximum de vraissemblance Recherche de motif dans un texte Automate fini, expressions régulières Recherche de mot dans un texte Algorithme Naïf Algorithme de KMP Algorithme de Boyer-Moore Trie Arbre de suffixes Liens préfixes Algorithme de McCreight Algorithme d'Ukkonen Structure secondaire des ARN Algorithme de Nussinov-Jacobson SCFG (HMM=Stochastic Regular Grammar) Graphes Dans un autre document... Application : contigs, cartes de restriction Optimisation Bientôt -- peut-être
(Pour les lecteurs qui auraient à enseigner l'algorithmique pour la bioinformatique.)
Vérifier qu'on a bien compris chacun de ces algorithmes, qu'on est capable de les dérouler sur des exemple, qu'on est capable de les implémenter -- on a choisi de les implémenter en C (ou parfois, selon mon humeur, en R ou en Python -- pour les algorithmes gourmants en calcul, ça n'est pas une très bonne idée), mais n'importe quel langage propre ferait l'affaire.
Voici quelques exemples de sujets de TP -- pour ceux qui auraient à enseigner la bioinformatique.
1. Implémenter l'algorithme d'alignement de Needleman-Wunsch.
S'il reste du temps : le modifier pour qu'il utilise une valeur différente pour la création et l'extension d'une brèche. Le modifier pour qu'il calcule un alignement local ou semi-local. Le modifier pour qu'il affiche tous les alignements optimaux, même s'il y en a plusieurs. Le modifier pour qu'il calcule la p-valeur du test "l'alignement est dû au hasard" (pour cela, mélanger les caractères d'une séquence, calculer le score de l'alignement, recommencer, après quelques centaines d'itérations, comparer le score initialement obtenu avec la distribution des scores après permutation).
2. Implémenter l'algorithme UPGMA
S'il reste du temps : modifier le programme pour qu'il calcule un alignement multiple à l'aide de l'arbre obtenu. Bootstrapper l'arbre.
3. Chaines de Markov
Prendre quelques gènes humains pour lesquels on connait la décomposition en introns et exons. En déduire une chaine de Markov cachée, à deux états (intron et exon). Implémenter l'algorithme de Viterbi qui utilise cette chaine pour prédire les introns et les exons. Quelle est la qualité du résultat ? (Penser à la valisation croisée.)
Variante : Prendre quelques protéines appartenant à une famille. Calculer une chaine de Markov (non cachée : les états sont les 20 acides aminés) qui reconnait cette famille. Ecrire un programme qui dit si une nouvelle protéine appartient à la famille ou pas (on donnera une p-valeur). Quelle est la qualité de cette prédiction (penser à la valisation croisée) ? Peut-on l'améliorer avec des chaines de Markov d'ordre supérieur ?
4. Implémenter l'algorithme de Knuth-Moris-Pratt.
S'il reste du temps : implémenter l'algorithme de Boyer-Moore. Implémenter l'algorithme d'Aho-Corasick.
5. Calculer naïvement l'arbre des suffixes d'un texte.
S'il reste du temps : implémenter l'algorithme de McCreight. A l'aide de simulations, vérifier qu'il est linéaire (en temps) en la taille du texte. Y a-t-il une différence entre la taille de l'arbre des suffixes des parties codantes et non-codantes du gènome ? Comment peut-on l'expliquer ?
Voici un exemple de sujet d'examen.
1. Donner un algorithme naïf permettant de trouver les répétitions dans le mot CGCTTACGCT (Indication : penser au dotplot). Quelle est sa complexité ?
2. Donner l'arbre des suffixes du mot CGCTTACGCT$. Pourquoi a-t-on ajouté un $ à la fin du mot ?
3. Ajouter les liens suffixes de McCreight.
4. Comment utiliser cet arbre des suffixes pour trouver les répétitions dans ce mot ? (Indication : comment peut-on interpréter les noeuds internes de l'arbre des suffixes ?) Donner tous les facteurs répétés. Quelle est la complexité cet l'algorithme de recherche de répétitions ?
5. Donner un algorithme de recherche de répétitions approchées dans un mot. (Indication : penser à la programmation dynamique) (Autre indication : il y a aussi un algorithme, sans programmation dynamique, qui repose sur un simple parcours en profondeur.)
Solution :
1. Faire un dotplot de la chaine contre elle-même et chercher la diagonale la plus longue (en dehors de la diagonale principale).
Remplissage du tableau : O(n^2)
Recherche de la diagonale la plus longue : O(n^2)
2.
--ACGCT$ --C |--T | |--$ | |--TACGCT$ |--GCT | |--$ | |--TACGCT$ --$ --T |--ACGCT$ |--$ |--TACGCT$ --GCT |--$ |--TACGCT$
On ajoute un $ (n'importe quel caractère qui n'apparait pas dans la chaine convient) afin qu'aucun suffixe ne soit préfixe d'un autre suffixe, i.e., pour que les feuilles de l'arbre soient en bijection avec les suffixes. Si on ne fait pas ça, on obtient ce qu'on appelle un "arbre des suffixes implicites".
3.
4. sous-chaines répétées = préfixes des noeuds internes.
C CT CGCT, CGC, CG T GCT, GC, G
Construction de l'arbre : O(n) (Ukkonen ou McCreight)
Recherche dans l'arbre (parcours en profondeur) : O(n) (Je rappelle que l'arbre est de taille O(n).)
5a. Procéder comme pour l'alignement local, mais en initialisant le tableau de programmation dynamique avec des zéros sur le première ligne, la première colonne et la diagonale.
5b. Pour des répétitions de taille m avec 1 erreur, chercher les chaines répétées de longueur au moins m/2 et essayer de les prolonger.
On dispose de deux séquences (d'ADN génomique, d'ARNm ou de protéines) et on essaye de voir à quel point elles sont semblables (alignement global), si elles contiennent des motifs communs (alignement local), si l'une est incluse dans l'autre (alignement semi-global) ou si on peut les recoller (c'est ce dont on a besoin quand on séquence un génome : on le découpe en morceaux, on séquence les morceaux et on essaye de les recoller -- c'est encore de l'alignement semi-global).
Etant données deux protéines, définies par leur structure primaire (la suite des acides aminés qui les composent), ou deux séquences nucléiques (génes, ARNm, etc.) on cherche si ces protéines sont semblables ou pas. Pour cela, on tente de les aligner, i.e., d'apparier chaque acide aminé de la première avec un acide aminé de la seconde, de manière à maximiser le nombre d'acides aminés en commun.
ACGTAC || || CGGATA
On peut aussi ajouter des trous (on parle de "brèche" ou, en franglais, de "gap") dans ces alignements.
-ACGTAC | || | TA-GTTC
D'un point de vue biologique, les gènes ou les séquences d'acides aminés évoluent : il arrive qu'un acide aminé soit remplacé par un autre, il arrive que des morceaux de séquence disparaissent, il arrive que de nouveaux morceaux de séquence (qui viennent d'une autre partie du génome, d'un autre individu (c'est courrant chez les bactéries), ou d'une autre espèce (c'est comme ça que les virus se reproduisent)) soient introduits. Ce sont ces phénomènes que les algorithmes d'alignement essayent de modéliser.
On dit que deux séquences sont similaires quand elles se ressemblent -- elles peuvent donc être plus ou moins similaires.
On dit que deux séquences sont homologues si elles ont une origine commune -- elles sont homologues ou elles ne le sont pas.
(Des séquences homologues dans une même espèce sont dites paralogues, des séquences homologues dans deux espèces différentes sont dites orthologues -- ces définitions simplistes sont provisoires, nous verrons les vraies définitions quand nous aborderons la phylogénie.)
La manière la plus simple de comparer deux séquences est de considérer une matrice, dont les lignes correspondent aux caractères de la première séquence, dont les colonnes correspondent aux lettres de la seconde, et dans laquelle on met un 1 aux cases pour lesquelles la lettre de la ligne coïncide avec la lettre de la colonne (on parle d'appariement (match, en anglais)) et 0 aux autres (mésappariement (mismatch)).
En général, on ne regarde pas les correspondances caractère par caractère, mais k caractères à la fois, afin d'éliminer le bruit (par exemple, k=10 pour l'ADN ou l'ARN, k=4 pour les protéines).
Voici quelques exemples.
Les deux séquences suivantes ont visiblement la même origine, et la même longueur : les segments correspondent aux parties communes.
Dans le graphique suivant, on voit un ARNm et l'ADN qui a servi à le produire : les deux sont identiques, mais l'ARNm a été obtenu en retirant certains morceaux de l'ADN (les introns -- les morceaux qui restent sont les exons).
Voici un dotplot d'une séquence par rapport à elle-même (les paramètres ne sont pas les mêmes que pour les graphiques précédents, car on désire plus de détails -- ici, on n'a que du bruit).
Voici un dotplot d'une séquence qui comporte des répétitions par rapport à elle-même.
Dans un dotplot, on ne considère que deux types d'appariement de deux nucléotides (ou acides aminés) : soit ils sont identiques, sont ils sont différents. On peut être un peu plus fin en tenant compte de la proximité des propriétés physico-chimiques des acides aminés ou des nucléotides : ils peuvent être identiques, très semblables, moyennement semblables, plutôt différents, complètement différents.
Par exemple, pour des séquences nucléiques, on peut affecter un score de 2 aux appariements, 1 aux appariements entre deux pyrimidines ou deux purines (des transitions) et 0 aux autres (des transversions). La matrice de similarité suivante résume ces score.
A C T G A 2 0 0 1 C 0 2 1 0 T 0 1 2 0 G 1 0 0 2
Le score d'un alignement de séquences nucléiques sera alors la somme des scores des appariements qui le constituent.
On peut faire la même chose pour les séquences protéiques, mais c'est un peu plus compliqué, car il y a 20 acides aminés et on ne sait pas vraiment quelles sont les propriétés physico-chimiques à retenir pour avoir un alignement qui reflète l'évolution des protéines. Les matrices de similarité pour les protéines ne seront pas empiriques, comme pour les nucléotides, mais seront estimées à partir de familles de protéines (i.e., des ensembles de protéines dont on sait qu'elles ont la même origine). Il y a essentiellement deux méthodes et deux familles de matrices de similarité pour les protéines : BLOSUM et PAM.
Pour affecter un score à l'appariement de deux acides aminés, on peut prendre des ensembles de séquences déjà alignées (par d'autres méthodes, la plus simple étant : les protéines ont la même fonction, la même longueur et en les mettant brutalement les unes au dessus des autres, on constate que la plupart des acines aminés s'alignent) et regarder quels sont les appariements les plus fréquents : on leur donnera le score le plus élevé.
Plus précisément, on part d'un ensemble de séquences identiques à plus de 62%, d'où l'appellation "BLOSUM62". Le nombre qui suit est le pourcentage de similitude entre les séquences utilisées pour construire la matrice -- il ne s'agit pas du nombre de mutations, mais d'un pourcentage de similarité : plus il est bas, plus les séquences sont éloignées ; de plus, il est toujours compris entre 0 et 100.
De plus, on part d'alignements dans brèche -- sinon la construction de l'alignement serait délicate.
Il reste encore à calculer le "score" : comment faire ? On part d'une matrice de taille 20*20, dont la case (i,j) contient le nombre d'appariements entre l'acide aminé i et l'acide aminé j. Cette matrice est symétrique.
Le score sera en fait un logarithme de rapports de vraissemblance :
Probabilité que i et j soient là pour des raisons évolutives
Score(i,j) = log_10 --------------------------------------------------------------
Probabilité que i et j soient là par hasard
Probabilité que i et j soient là pour des raisons évolutives
= log_10 --------------------------------------------------------------
(Probabilité d'observer i) * (Probabilité d'observer j)
A_ij/n
= log_10 -----------------
A_i./n * A_.j/n(Le lecteur est sensé avoir fait un peu de statistiques, en particulier connaître la notion de vraissemblance, la notion d'estimateur du maximum de vraissemblance et le test des rapport de vraissemblance.)
Les matrices PAM sont calculées presque de la même manière : on part d'un ensemble de séquences assez proches dont on sait qu'elles ont une origine commune. On calcule une phylogénie pour expliquer leur évolution et on reconstruit toutes les séquences ancestrales. Au lieu de compter le nombre d'occurrences de chaque type d'appariement entre les séquences actuelles, comme pour les matrices BLOSUM, on compte le nombre d'occurrences de chaque type d'appariement dans tous les couples (père,fils).
On compte donc le nombre de mutations et pas le nombre de différences : ainsi, le nombre qui suit PAM (par exemple, PAM120) désigne (en pourcents) le nombre moyen de mutations par site entre deux séquences actuelles (en suivant un chemin qui passe par les ancètres des séquences).
Une autre différence avec les matrices BLOSUM, c'est qu'il y a un modèle probabiliste derrière, qui tient compte de l'évolution (celui des matrices BLOSUM est très rudimentaire). Ses hypothèses sont les suivantes : les mutations sont indépendantes (indépendantes entre elles, indépendantes des mutations intervenues auparavant, indépendantes de la position, indépendantes des bases alentour), la probabilité d'une mutation dans un sens est la même que la probabilité de la mutation dans l'autre (P(i-->j)=P(j-->i)).
L'hypothèse d'indépendance nous dit qu'on peut multiplier toutes ces probabilités ; comme on a pris le logarithme, il revient au même d'ajouter ces logarithmes -- c'est numeriquement plus stable.
Remarque : les matrices PAM (construites à partir d'un ensemble de séquences très proches : au moins 85% de similarité) sont un bon choix si on veut comparer des séquences très proches, mais pas pour comparer des séquences éloignées.
BLOSUM PAM
-----------------------------------------------------------------
Séquences de même taille Séquences de même origine,
et fonction, donc alignées pour lesquelles on calcule une
phylogénie et on reconstruit
les séquences ancestrales.
-----------------------------------------------------------------
On n'utilise que les séquences On utilise à la fois les séquences
actuelles actuelles et les séquences
ancestrales reconstruites
-----------------------------------------------------------------
Nombre de différences par site Nombre de mutations par site
-----------------------------------------------------------------
BLOSUM 80 PAM 1 (séquences proches)
BLOSUM 62 PAM 120
BLOSUM 45 PAM 250 (séquences éloignées)
-----------------------------------------------------------------
Modèle probabiliste sous-jacent Modèle probabiliste sous-jacent
rudimentaire, qui ne connait qui tient compte de l'évolution
pas la notion d'évolution
-----------------------------------------------------------------
Plus facile Plus difficile
(pas de reconstruction
phylogénétique ni de calcul
des séquences ancestrales)
-----------------------------------------------------------------
Pour des séquences pas trop Pour des séquences proches
proches (car alors, le
chemin suivi par l'évolution
est difficile à estimer, cette
estimation est peu fiable, et
ne sert plus à grand-chose).Pour plus de détails, lire :
http://lectures.molgen.mpg.de/statistik/docs/Kapitel_12.pdf http://www.esil.univ-mrs.fr/%7Edgaut/Cours/pam.html http://www.genopole-lille.fr/fr/formation/cib/doc/anaseq/matrices.pdf
JTT (et peut-être aussi WAG) est une variante de PAM, un peu plus récente. Lire l'introduction de
http://www.bioinf.org/molsys/data/WAG_MATRICES.pdf
Mentionnons aussi STR et metREV.
Etant données deux séquences, on peut calculer un score de similarité (d'autant plus élevé que les séquences sont proches) ou un score de distance (positif ou nul, d'autant plus proche de zéro que les séquences sont proches).
Pour l'alignement global, on a le choix, mais pour l'alignement local, les distances poseraient des problèmes d'interprétation : on prendra des scores de similarité.
La programmation dynamique est une généralisation du principe "diviser pour régner", que l'on retrouve dans beaucoup d'algorithmes. Quand ce principe s'applique, on peut représenter le déroulement de l'algorithme par un arbre, i.e., un graphe sans cycle : on a besoin de chaque résultat intermédiaire une seule fois.
A FAIRE Donner un exemple Dessiner l'arbre
Ces algorithmes "diviser pour régner" se programment très naturellement par récurrence (la récurrence sert en fait à explorer des arbres).
Parfois, ce n'est pas un arbre, i.e., on peut avoir besoin de réutiliser plusieurs fois des résultats intermédiaires : dans une telle situation, on garde ces résultats intermédiaires en mémoire (très souvent dans un tableau). On parle alors de programmation dynamique.
C'est par exemple dans ce cadre que rentre le calcul des coefficients binomiaux, définis par
C(n,0) = 1 C(n,n) = 1 C(n+1,k) = C(n,k) + C(n,k-1)
On a l'impression que ce calcul peut s'effectuer par récurrence (d'ailleur beaucoup de livres d'algorithmiques commettent cette erreur -- et après on va s'étonner que les étudiants trouvent la notion de récurrence compliquée), mais c'est une très mauvaise idée. En effet, comme
C(n+1,k) = C(n,k) + C(n,k-1)
= C(n-1,k) + C(n-1,k-1) + C(n-1,k-1) + C(n-1,k-2)on serait amené à calculer deux fois C(n-1,k-1). Pire, plus on descend dans l'"arbre", plus le nombre de calculs redondants augmente. La raison de ces problèmes, c'est que l'"arbre" qu'on explore n'est pas un arbre.
A FAIRE : dessin de ce "non-arbre" (c'est le triangle de Pascal).
Ici, il suffit de calculer progressivement les valeurs de C(n,k), pour des valeurs croissantes de n.
C[0,0] <- 1
Pour n de 1 à N:
C[n,0] <- 1
Pour k de 1 à n-1:
C[n,k] <- C[n-1,k] + C[n-1,k-1]
C[n,n] <- 1En général, un algorithme de programmation dynamique se déroule ainsi : on commence par calculer tous les résultats intermédiaires en les stockant dans un tableau, puis, à l'aide de ce tableau, on construit l'objet demandé. Sauf dans des situations triviales, il y a deux étapes : remplissage du tableau (qui ne donne pas directement la solution), puis construction de la solution, de proche en proche, en s'aidant du tableau.
Voici quelques exemples de problèmes se résolvant à l'aide de la programmation dynamique.
1. On veut multiplier entre elles n matrices, de tailles différentes : dans quel ordre effectuer ces multiplications (i.e., comment trouver un parenthésage de ce produit -- on peut représenter un parenthésage par un arbre binaire) pour faire le moins d'opérations possible ?
2. Trouver un alignement de deux protéines (ce sont les algorithmes de Needleman & Wunsch et de Smith & Waterman, que nous allons détailler ci-après).
Ces algorithmes ont en fait des applications en dehors de la bioinformatique. Ainsi, pour extraire les sigles et leurs définitions dans un texte (par exemple, un corpus de résumés d'articles de recherche), on peut commencer par repérer les mots en majuscules entre parenthèses et utiliser un algorithme d'alignement (basé sur de la programmation dynamique) pour trouver la correspondance entre les lettres des mots précédents et les lettres du sigle.
http://biotext.berkeley.edu/papers/psb03.pdf
3. On dispose d'un automate fini dont les transitions sont chacune affectées d'une probabilité et dont les états sont coloriés. On connait la suite des couleurs et on essaye de trouver la suite d'états dont elle provient le plus probablement. C'est l'algorithme de Viterbi, que nous rencontrerons un peu plus loin, quand nous parlerons de Modèles de Markov Cachés (HMM).
4. Nous verrons d'autres exemples de programmation dynamique au cours de ce document : par exemple l'algorithme de Nussinov (pour calculer la structure secondaire des ARN), l'algorithme de Fitch-Hartigan (pour calculer les séquences ancestrales dans un alignement phylogénétique -- ce dont nous avons besoin pour calculer les matrices PAM), l'algorithme de Sankoff (pour calculer le score de parsimonie d'un arbre).
Pour plus de détails, lire le chapitre correspondant du Cormen-Leisserson-Rivest et faire les exercices.
A FAIRE : l'URL du livre (chez Amazon)
Pour aligner deux séquences, on utilise un tableau, dont les lignes représentent les lettres de la première séquence et les colonnes les lettres de la deuxième (comme pour le dotplot). Un alignement correspond alors à un chemin dans le tableau.
-ACGTAC | || | TA-GTTC ACGTAC T * A ** G * T * T * C *
On cherche l'alignement qui présente la plus grande similarité ; la similarité est définie à partir du choix d'une matrice de similarité et d'une pénalité pour les brèches (gaps) (generalement, on prend une pénalité pour l'ouverture de la brèche et une autre proportionnelle à sa longueur ou au logarithme de sa longueur).
Pour chaque case du tableau, on cherche la similarité du meilleur chemin vers l'origine, et on met cette similarité dans la case. Une fois le tableau rempli, on connait la similarité du chemin le meilleur, on sait où il aboutit et on peut le reconstruire de proche en proche.
Pour l'algorithme de Needleman et Wunsch, on initialise la première ligne et la première colonne du tableau (qui correspondent à des brèches au début d'une des deux séquences) et on commence la reconstruction de l'alignement à partir de la case en bas à droite.
http://bioportal.weizmann.ac.il/course/ATIB02/lecture1/lecture1.pdf http://www.maths.tcd.ie/~lily/pres2/
Remarque : le résultat va dépendre du choix de la matrice de similarité et de la pénalité choisie pour les brèches. On n'hésitera pas à les faire varier, pour voir comment l'alignement change.
Remarque : il peut y avoir plusieurs solutions.
Remarque : la solution que l'on veut (celle qui correspond effectivement à l'évolution des séquences) est peut-être sub-optimale (elle a un score très proche des solutions optimales) ou locale (avec des brèches au bord).
Pour pouvoir comprendre cette valeur, on peut prendre plein de protéines au hasard et regarder la distribution des similarités. On peut aussi prendre des séquences au hasard, construites à partir de morceaux de protéines (comme pour le bootstrap des séries temporelles).
A FAIRE : Donner l'URL du numéro de Rnews qui parle du bootstrap des séries temporelles.
Voici l'algorithme
(Non testé)
Soient
s1, s2 les deux séquences à aligner, de longueurs m et n
M une matrice m*n, que l'on remplira de sorte que
M[i,j] = score du meilleur alignement de s1[1..i] sur s2[1..j]
a : une pénalité pour les brèches et les mesappariements
0 : le score (la "distance") d'un appariement
Initialisation :
M[i,1] <- a*(i-1)
M[1,j] <- a*(j-1)
(ici, les indices des tableaux commencent à 1)
Remplissage :
D[i,j] <- M[i-1,j-1] + 0 si s1[i]==s2[j] appariement
M[i-1,j-1] + a sinon mesappariement
P[i,j] <- M[i-1,j] + a brèche
Q[i,j] <- M[i,j-1] + a brèche
M[i,j] <- min( D[i,j], P[i,j], Q[i,j] )Exemple d'implémentation :
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <string.h>
#define MAX_LENGTH 1000
#define TRANSITION 1
#define TRANSVERSION 0
#define IDENTITY 3
char chaine1[] = "atggccctgtggatgcgcctcctgcccctgctggcgctgctggccctctggggacctgacccagccgcagcctttgtgaaccaacacctgtgcggctcacacctggtggaagctctctacctagtgtgcggggaacgaggcttcttctacacacccaagacccgccgggaggcagaggacctgcaggtggggcaggtggagctgggcgggggccctggtgcaggcagcctgcagcccttggccctggaggggtccctgcagaagcgtggcattgtggaacaatgctgtaccagcatctgctccctctaccagctggagaactactgcaactag";
char chaine2[] = "atggccctgttggtgcacttcctacccctgctggccctgcttgccctctgggagcccaaacccacccaggcttttgtcaaacagcatctttgtggtccccacctggtagaggctctctacctggtgtgtggggagcgtggcttcttctacacacccaagtcccgccgtgaagtggaggacccacaagtggaacaactggagctgggaggaagccccggggaccttcagaccttggcgttggaggtggcccggcagaagcgtggcattgtggatcagtgctgcaccagcatctgctccctctaccagctggagaactactgcaactaa";
int Gamma(char i, char j) {
if(i==j)
return IDENTITY;
if ( (i=='a' && j=='g') || (i=='c' && j=='t') )
return TRANSITION;
if ( i=='-' || j=='-' )
assert(0==1); // GAP non utilisé
return TRANSVERSION;
}
int max(int a, int b, int c) {
if( a >= b && a >= c )
return a;
if( b >= a && b >= c )
return b;
return c;
}
int min(int a, int b) {
return (a>b) ? b : a;
}
void strrev(char * s) {
int i;
int n = strlen(s);
for (i=0; i<n/2; i++) {
char tmp = s[i];
s[i] = s[n-1-i];
s[n-i-1] = tmp;
}
}
void affiche_scores (int ** score, int n) {
int i,j;
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
printf("%3d ", score[i][j]);
}
printf("\n");
}
}
#define MODULO 100
void affiche_alignement (char * a, char * b) {
int n = strlen(a);
int i,j;
assert(strlen(a) == strlen(b));
for (i=0; i< (int) ceil(((double)n)/MODULO) ; i++) {
for (j=MODULO*i; j<min(n, MODULO*(i+1)); j++) {
printf("%c",a[j]);
if (j % 10 == 9) printf(" ");
}
puts("");
for (j=MODULO*i; j<min(n, MODULO*(i+1)); j++) {
if( a[j] == b[j] )
printf("|");
else
printf(" ");
if (j % 10 == 9) printf(" ");
}
puts("");
for (j=MODULO*i; j<min(n, MODULO*(i+1)); j++) {
printf("%c",b[j]);
if (j % 10 == 9) printf(" ");
}
puts("\n");
}
}
int main (int argc, char ** argv) {
int i,j,k;
int GAP_START = -5;
int GAP_EXT = -2;
if(argc>1)
GAP_START = atoi( argv[1] );
if(argc>2)
GAP_EXT = atoi( argv[2] );
int score_d [MAX_LENGTH][MAX_LENGTH];
int score_p [MAX_LENGTH][MAX_LENGTH];
int score_q [MAX_LENGTH][MAX_LENGTH];
assert( strlen(chaine1) < MAX_LENGTH );
assert( strlen(chaine2) < MAX_LENGTH );
for (i=0; i<strlen(chaine1); i++)
for (j=0; j<strlen(chaine2); j++) {
score_d[i][j] = 0;
score_p[i][j] = 0;
score_q[i][j] = 0;
}
for (i=0; i<strlen(chaine1); i++) {
score_d[i][0] = GAP_START + (i-1)*GAP_EXT;
score_p[i][0] = GAP_START + (i-1)*GAP_EXT;
score_q[i][0] = GAP_START + (i-1)*GAP_EXT;
}
for (j=0; j<strlen(chaine2); j++) {
score_d[0][j] = GAP_START + (j-1)*GAP_EXT;
score_p[0][j] = GAP_START + (j-1)*GAP_EXT;
score_q[0][j] = GAP_START + (j-1)*GAP_EXT;
}
score_d[0][0] = score_p[0][0] = score_q[0][0] = 0;
for (i=0; i<strlen(chaine1); i++) {
for (j=0; j<strlen(chaine2); j++) {
int d = score_d[i][j] + Gamma(chaine1[i],chaine2[j]) ;
int p = score_p[i][j] + GAP_START;
int q = score_q[i][j] + GAP_START;
score_d[i+1][j+1] = max(d,p,q);
d = score_d[i][j+1] + GAP_START;
p = score_p[i][j+1] + GAP_EXT;
q = score_q[i][j+1] + GAP_START;
score_p[i+1][j+1] = max(d,p,q);
d = score_d[i+1][j] + GAP_START;
p = score_p[i+1][j] + GAP_START;
q = score_q[i+1][j] + GAP_EXT;
score_q[i+1][j+1] = max(d,p,q);
}
}
// affiche_scores(score_d, 20);
{
int i,j;
int n=20;
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
printf("%3d ", score_d[i][j]);
}
printf("\n");
}
}
i=strlen(chaine1)-1;
j=strlen(chaine2)-1;
k=0;
char solution1[MAX_LENGTH];
char solution2[MAX_LENGTH];
while(i>=0 && j>=0) {
int d = score_d[i+1][j+1];
int p = score_p[i+1][j+1];
int q = score_q[i+1][j+1];
int m = max(d,p,q);
if( m==d ) {
solution1[k]=chaine1[i];
solution2[k]=chaine2[j];
// printf("d%d%c%c ", k, solution1[k], solution2[k]);
i--;
j--;
} else if ( m==p ) {
solution1[k] = chaine1[i];
solution2[k] = '-';
// printf("p%d%c%c ", k, solution1[k], solution2[k]);
i--;
} else {
solution2[k] = chaine2[j];
solution1[k] = '-';
// printf("q%d%c%c ", k, solution1[k], solution2[k]);
j--;
}
k++;
}
puts("\n");
// printf("1: %d, %d, %s\n", strlen(solution1), solution1[2], solution1);
solution1[k] = '\0';
solution2[k] = '\0';
strrev(solution1);
strrev(solution2);
//affiche_scores(score, 20);
affiche_alignement(solution1, solution2);
return EXIT_SUCCESS;
}
Dans l'exemple précédent, nous avons attribué une pénalité, toutjours la même, pour chaque position contenant une brèche. Ce n'est pas très réaliste : on peut séparer cette pénalité en deux, d'une part une pénalité de création de brèche, d'autre part une pénalité d'extension de brèche.
pénalité = pénalité de création + longueur * pénalité d'extension
Exercice : modifier l'algorithme précédent en conséquence. Quelle est sa complexité ? (Réponse : toujours O(n^2).)
On peut aussi considérer une pénalité non linéaire, par exemple
pénalité = pénalité de création + pénalité d'extension * log(longueur)
Dans le cas d'une pénalité quelconque,
pénalité = f(longueur) (où f est une fonction quelconque)
l'algorithme est en O(n^3).
On recherche des alignements locaux : l'algorithme est le même, avec quelques différences.
1. On initialise la matrice avec des zéros sur la première ligne et la première colonne ;
2. On utilise des scores de similarités et pas des distances (avec l'alignement global, on a le choix, mais pour l'alignement local, la notion de "distance" poserait des problèmes d'interprétation) ;
3. Lors des calculs de scores, on remplace les nombres négatifs par zéro ;
4. On construit le chemin à partir de la case contenant la valeur la plus élevée (et pas en partant de la case en bas à droite) et en s'arrétant quand on arrive à une case contenant zéro -- on peut d'ailleurs avoir plusieurs chemins.
Si x et y sont deux séquences, on note SW(x,y) le score de leur alignement obtenu par l'algorithme de Smith-Watermann. C'est un noyau. Cela signifie qu'in existe un espace vectoriel euclidien E et, pour toute séquence x, un vecteur f(x) dans E, tel que SW(x,y) soit le produit scalaire de f(x) et f(y).
On peut par exemple calculer la distance entre deux séquences :
d(f(x),f(y)) = ( f(y) - f(x) ) * ( f(y) - f(x) )
= f(y) * f(y) + f(x) * f(x) - 2 f(x) * f(y)
= SW(x,x) + SW(y,y) - 2 SW(x,y)Exercice : Soient x1, x2,...,xn des séquences. Calculer la distance entre f(x) et le centre de gravité de f(x1),...,f(xn).
(On pourra utiliser ce résultat un peu plus loin, pour trouver la séquence "la plus centrale", dans l'algorithme de l'étoile, qui calcule un alignement multiple.)
Voir :
http://www.cg.ensmp.fr/%7Evert/talks/021205inria/inria.pdf
Au lieu de travailler directement avec les séquences, on se place dans un espace vectoriel qui décrit un peu mieux les relations entre elles. Cette espace est peut-être de dimension démesurée : ça n'est pas grave, car on ne manipule jamais des coordonnées dans cet espace, on passe toujours par le noyau.
Cette idée se généralise : pour étudier des points dans un espace vectoriel, on peut plonger ces points dans un espace vectoriel plus gros, suceptible de mieux décrire le phénomène. Par exemple, en remplaçant un point de coordonnées (x,y) par un point de coordonnées (x,y,x^2,y^2,x*y), on peut linéariser certains phénomènes. Si l'espace de départ est déjà de dimension élevée, cela peut sembler une mauvaise idée, car sa dimension sera encore pire : il n'en est rien, car on ne manipule pas les coordonnées dans ce nouvel espace, on se contente d'utiliser son produit scalaire.
Dans ce nouvel espace, on peut appliquer des méthodes classiques telles l'énalyse en composantes principales ou la régression. C'est d'ailleurs l'une des idées derrière les SVM (Support Vector Machine) -- j'ai déjà dit que c'était important, les SVM, non ? -- je le répèterai encore.
http://www.kernel-machines.org/papers/talk-dagm99_4.ps.gz
On dispose de deux séquences, de longueurs très différentes, et on soupçonne que l'une des deux est contenue dans l'autre (à quelques mutations et indels près). L'alignement local ne convient pas, car il nous donnerait des morceaux des deux séquences ; l'alignement global non plus, car il tiendrait compte d'une pénalité pour les brèches aux extrémités.
Nous allons donc simplement reprendre l'algorithme de l'alignement global sans mettre de pénalité pour les brèches aux extrémités.
On initialise la matrice avec des zéros sur la première ligne et la première colonne, on cherche un chemin qui aille d'une case de la premièer ligne ou de la première colonne à une case de la dernière ligne et de la dernière colonne -- on part donc de la case contenant la valeur maximale parmi les cases de la dernière ligne et de la dernière colonne.
Autre application : les extrémités des deux séquences sont semblables (par exemple, deux contigs dont on pense qu'ils se chevauchent).
C'est une modification des algorithmes précédents, linéaire en espace (dans les algorithmes précédents, on devait garder en mémoire un tableau contenant toutes les similarités : c'était quadratique).
A FAIRE : à défaut de détailler l'algorithme, mettre au moins un URL...
On dispose maintenant d'un ensemble de séquences (une "base de données", compter plusieurs Go) et on cherche, étant donnée une nouvelle séquence, s'il y a des séquences similaires dans la base.
Les algorithmes suivants sont des heuristiques pour calculer un tel alignement local : ils donnent un résultat approché mais rapide.
1. Rechercher des correspondances de mots de 10 lettres (on n'utilise pas de matrice, comme dans l'algorithme de SW, mais une table de hachage, qui met en correspondance les mots de dix lettres, les séquences dans lesquelles ils apparaissent et l'endroit exact où ils apparaissent) entre la séquence à étudier et chacune des séquences de la base de données. Dans cette étape, on demande des correspondances exactes.
2. Trouver des segments dans la matrice de SW, i.e., des suites contiguës de correspondances. On peut par exemple regarder s'il n'y a pas une diagonale comportant beaucoup de segments (le numéro de la diagonale passant par (i,j) est i-j), que l'on peut ensuite essayer de recoller, avec une matrice de similitude : on accepte des correspondances approchées. On essaye ensuite de prolonger, de part et d'autre, les dix segments les plus longs.
3. Essayer de relier ces segments (on peut par exemple tracer tous ces segments sur un dotplot, les relier de toutes les manières possibles et trouver un chemin de score maximal dans le graphe ainsi obtenu).
A FAIRE : dessin
4. Pour les séquences de la base de donnée pour lesquelles on a trouvé des correspondances intéressantes, chercher un alignement avec l'algorithme de SW. On ne regarde pas la matrice entière, mais on se limite à une bande autour de l'alignement trouvé au point précédent.
Vocabulaire de FastA : un hotspot est un mot de dix lettres commun entre les séquences, un diagonal run est un ensemble de hotspots proches et sur la même diagonale.
Il y a en fait 3 algorithmes (programmes) différents :
Blast1 (search of local similarities without gaps) NCBI-Blast2 WU-Blast2
On utilisera généralement le dernier, car en changeant ses options, on peut se ramener aux deux premiers.
A FAIRE : URL blast_fasta-uk.pdf http://acer.gen.tcd.ie/~amclysag/blast.html http://bioweb.pasteur.fr/seqanal/blast/intro-uk.html
1. Recherche des correspondances approchées de mots de 4 lettres (pour les séquences d'acides aminés, on prendra 4 lettres, pour les séquences de bases, on en prendra 10). Avec FastA, on utilisait une table de hachage, mais ici on utilise un automate fini -- mais c'est pareil.
http://www.dcs.kcl.ac.uk/teaching/units/csmacmb/DOC/03-acmb-sols.ps
A la différence de FastA, on ne se limite pas aux correspondances exactes : la table de hachage contient les mots de 4 lettres des séquences connues, et on l'interroge avec les mots de 4 lettres de la séquence à étudier et ceux qui leurs sont proches.
2. On étend (sans brèche (gap)) chacune des ces correspondances. On ne les étend pas jusqu'au bout, on fixe une valeur minimale en dessous de laquelle la similarité ne doit pas descendre.
Le résultat est la liste de ces correspondances, chacune avec son score (similarité) et sa p-valeur, i.e., la probabilité qu'une correspondance avec un score aussi élevé soit dur au hasard. On a aussi une p-valeur générale, probabilité d'observer un ensemble de correspondances dont la somme des scores soit aussi élevée.
Certains documents sur internet disent que le BLAST utilise un "automate fini", en laissant croire qu'il s'agit de quelque chose de compliqué et en passant les détails sous silence. En fait, c'est juste une manière simple de décrire "un mot de quatre lettres et tous ses voisins".
Pour les lecteurs qui ne connaissent pas encore la notion d'automate fini (c'est important, nous y reviendrons beaucoup plus en détail), c'est juste un graphe dont les arrêtes sont étiquetées par des lettres et que l'on cherche à parcourrir (d'un noeud "départ" jusqu'à un noeud "arrivée") en utilisant une à une les lettres d'un mot, dans l'ordre. Si on y arrive, on dit que le mot est reconnu par l'automate.
Prenons l'exemple du mot "ACTT". Il est reconnu par l'automate
A C T T
départ -----> * -----> * -----> * -----> arrivée,qui ne reconnait d'ailleurs aucun autre mot.
L'automate suivant reconnait le mot ACTT et tous les mots qui n'en diffèrent que par une seule lettre.
Le lecteur qui connait la notion d'expression régulière remarquera que cet automate fini est équivalent à
ACT. | AC.T | A.TT | .CTT
où "." désigne un caractère quelconque et | la disjonction (le "ou"). En fait, les notions d'automate fini et d'expression régulière sont équivalentes.
1. La première étape est la même : recherche des correspondances approchées de mots de quatre lettres entre la séquence à étudier et la base de données.
2. On ne tente pas de prolonger toutes les correspondances : on se limite à celles pour lesquelles il y a une autre correspondance sur la même diagonale, pas trop loin. Il s'agit toujours de prolongement sans brèche (gap).
3. Pour la (les) meilleur correspondance, on fait un prolongement avec des brèches (gaps), par un algorithme de SW modifié, en partant du milieu de la correspondance.
Idem, mais le dernier prolongement peut se faire par un algorithme de SW limité à une bande ou avec un similitude minorée.
En jouant avec les options, on peut lui demander de se comporter comme FastA, Blast1 ou NCBI-Blast2.
Par ailleurs, WU-Blast2 tient compte de la présence éventuelle de plusieurs alignements pour un même couple de séquences.
C'est un BLAST optimisé à l'aide de "profils" (un profil, c'est quelque chose qui permet de reconnaitre une famille de gènes ou de protéines : par exemple une séquence consensus, une expression régulière, une matrice PSSM, une chaine de Markov cachée, etc. -- nous en reparlerons plus loin, avec les alignements multiples)
1. Faire un blast, puis un alignement multiple et en déduire un profil.
2. Interroger la base de données sur ce profil, faire un alignement multiple, en déduire un profil. Recommencer.
C'est un Blast limité aux séquences qui contiennent un motif. L'alignement multiple devra respecter ce motif.
Attention : la notion de E-valeur est assez délicate à (bien) comprendre et à (bien) utiliser -- je ne suis d'ailleurs pas persuadé d'avoir tout compris. Cette partie est d'ailleurs trop courte, trop allusive, trop superficielle...
Généralement, le résultat du blast n'est pas accompagné d'une p-valeur (probabilité d'observer autant d'alignements dans des séquences prises au hasards), mais plutôt une E-valeur, qui est l'espérance du nombre d'alignements avec un score au moins aussi bon dans une séquence aléatoire de même longueur et composition.
Cette valeur dépend de la taille de la base de donnée, car plus la base de donnée est grande, plus elle est complète, plus on risque de trouver une séquence "proche" de celle qui nous intéresse -- mais cela risque d'être révélateur de la densité des séquences de la base de donnée plutôt que d'une homologie...
De plus, comme les bases de données grossissent sans cesse, ces E-valeurs croissent sans cesse.
Mais pourquoi prendre des E-valeurs, plus délicates à interpréter que les p-valeurs ? En effet, on a l'habitude d'interpréter des p-valeurs, pas des E-valeurs ; les p-valeurs sont dans [0,1] alors que les E-valeurs peuvent être supérieures à 1 ; on sait comment modifier des p-valeurs quand on fait des tests très nombreux, mais quid des E-valeurs ; etc.
Quel est l'avantage des E-valeurs sur les p-valeurs ?
(En fait, au moment où j'écris ces lignes, je ne suis moi-même pas convaincu -- aussi mes arguments paraitront-ils peut-être un peu légers.)
Voici une raison. On parle de p-valeur par rapport à un test. Lors d'une recherche dans une base de donnée, ce test serait
H0: "similitude fortuite"
contre
H1: "homologie, i.e., similitude résultant d'une origine commune".
le problème, c'est que ces hypothèses sont plus biologiques que statistiques.
Comment calculer ces E-valeurs ?
Notons p1...pn la proportion des différents AA dans la base et dans la séquence à étudier et S(i,j) la matrice de similarité (on suppose que l'un des S(i,j) est positif et que l'espérance du score,
Somme( S(i,j) * p_i * p_j ),
est négative). C'est le modèle de Karlin-Altschul.
Le score du meilleur alignement local sans brèche (gap) suit alors une distribution evd (extreme value distribution).
A FAIRE : URL d'un article expliquant ce que sont les EVD...
C'est ainsi qu'on calcule la p-valeur du blast.
En présence de brèche, on n'a pas de théorie : on utilise des simulations. Les simulations sont longues, elles sont faites avant, pour une matrice de similarité donnée.
Il faut faire très attention à l'interprétation de ces valeurs :
http://conferences.oreillynet.com/presentations/bio2003/dwan_chris.ppt
La E-valeur dépend de la taille de la requête et de la taille de la base. Ainsi, si on fixe une E-valeur seuil et si on augmente ensuite la taille de la base de données, on verra certains résultats disparaître -- car leur E-valeur aura augmenté et sera passée au dessus du seuil.
On prendra aussi garde à ne pas confondre "statistiquement significatif" et "biologiquement significatif".
Dans l'alignement de deux séquences, on cherchait à savoir si deux séquences étaient proches ou pas -- et si cette similarité avait une origine biologique (homologie). Ici, la démarche est inverse : on a un ensemble de séquences dont ont sait qu'elles sont homologues et on essaye de voir comment cette homologie se traduit au niveau des séquences.
http://www.infobiogen.fr/doc/tutoriel/PHYLO/phylogenie.html http://www.infobiogen.fr/services/deambulum/fr/bioinfo3.html http://www.unige.ch/sciences/biologie/biani/msg/teaching/evolution.htm http://tolweb.org/tree/phylogeny.html http://rdp.cme.msu.edu/html/ http://www3.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html http://www.treebase.org/treebase/
On sait définir le score de l'alignement de deux séquences. Il y a plusieurs manières de le généraliser à l'alignement de plus de deux séquences.
On peut tout d'abord faire la somme des score de tous les alignements de deux séquences définis par notre alignement multiple : c'est le score des "sommes de paires".
SP = Somme score(A_i, A_j)
i<jOn peut aussi faire jouer un rôle central à une séquence et considérer la somme des scores des alignements de chaque séquence avec cette séquence centrale.
S = Somme score(A_i0, A_i)
i!=i0Si on dispose d'un arbre phylogénétique, on peut calculer un score SP pondéré (le score d'un arbre est la moyenne du score de sa branche droite et de sa branche gauche, qui se calculent de manière récursive -- Exercice : modifier cette définition pour tenir compte de la longueur des branches).
Si on dispose d'un arbre phylogénétique, on peut faire la somme des scores des alignement de deux séquences voisines sur l'arbre (mais cela suppose qu'on a réussi à reconstruire les séquences ancestrales).
S = Somme( score(A_i,A_j) où (i,j) est une branche de l'arbre )
On peut aussi utiliser la notion d'entropie. En statistique, on définit l'edntropie d'un message de probabilité p comme
Entropie = - log2 p
Pour un alignement, on prendrait :
Entropie = Somme Somme c(a,i) log2 p(a,i)
i aoù
i = numéro de la colonne a = une lettre (A, C, T, G pour les séquences nucléiques, etc.) p(a,i) = probabilité d'observer la lettre a dans la colonne i de l'alignement c(a,i) = nombre d'occurrences de la lettre a dans la colonne i
Dans la suite, pour simplifier, nous nous limiterons essentiellement au score des sommes de paires (SP).
On peut toujours utiliser la programmation dynamique pour aligner trois séquences ou plus : au lieu de remplir un tableau à deux dimensions, l'algorithme remplit un tableau de dimension trois ou plus. Le problème, c'est que la complexité (en mémoire et en temps) est exponentielle par rapport au nombre de séquences.
Il existe néanmoins certaines heuristiques qui permettent de se limiter à certains morceaux de ce tableau multidimentionnel. Mentionnons MSA, DCA, Prrp.
http://www.genopole-lille.fr/fr/formation/cib/doc/anaseq/multiple.pdf
Les algorithmes d'alignement multiple suivants construisent l'alignement progressivement, en ajoutant une séquence à la fois. Cela peut poser un gros problème au niveau des brèches : si l'un des premiers alignements crée une brèche, cette brèche sera toujours présente à la fin de l'alignement multiple, quelles que soient les autres séquences.
En anglais, on dit "once a gap, always a gap" ("brèche un jour, brèche toujours").
Après avoir calculé la distance entre tous les couples de séquences, on essaye de trouver celle qui est la plus centrale : par exemple celle pour laquelle la somme des distances aux autres suites est la plus faible. On aligne ensuite chaque séquence avec cette séquence centrale.
Voici un exemple d'implémentation (si vous lancez le programme, vous verrez que python n'est peut-être pas un très bon choix pour ce genre d'algorithme : il faut attendre plus de dix secondes pour le calcul des premiers alignements deux à deux).
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
from Numeric import *
# Les séquences à aligner
sequences = [
"atggcagtgtggcttcaggctggtgctctgttggtcctgttggtcgtgtccagtgtaagcactaacccaggcacaccgcagcacctgtgtggatctcatctggtcgatgccctttatctggtctgtggcccaacaggcttcttctacaaccccaagagagacgttgagccccttctgggtttccttcctcctaaatctgcccaggaaactgaggtggctgactttgcatttaaagatcatgccgagctgataaggaagagaggcattgtagagcagtgctgccacaaaccctgcagcatctttgagctgcagaactactgtaactga",
"atggcagcgctctggctccaggccttctccctgctcgtcttaatgatggtttcgtggccgggctcccaggccgtcggtgggccacagcacctgtgcggctcccacctggtggatgccctgtacctggtctgtggggacagaggcttcttctacaaccccaggagagatgtggaccctctgcttggtttcctccctccaaaggcaggtggtgctgtggtgcaaggtggtgagaatgaagtgaccttcaaagaccagatggaaatgatggtgaagcgaggcattgtggaggaatgctgtcacaaaccctgtaccatcttcgacctgcagaactactgcaactga",
"atggccctgtggatgcacctcctgcccctgctggcgctgctggccctctggggacccgagccagccccggcctttgtgaaccagcacctgtgcggcccccacctggtggaagccctctacctggtgtgcggggagcgaggtttcttctacgcacccaagacccgccgggaggcggaggacctgcaggtggggcaggtggagctgggtgggggctctatcacgggcagcctgccacccttggagggtcccatgcagaagcgtggcgtcgtggatcagtgctgcaccagcatctgctccctctaccagctgcagaactactgcaactag",
"atggccctgtggacacgcctgcggcccctgctggccctgctggcgctctggcccccccccccggcccgcgccttcgtcaaccagcatctgtgtggctcccacctggtggaggcgctgtacctggtgtgcggagagcgcggcttcttctacacgcccaaggcccgccgggaggtggagggcccgcaggtgggggcgctggagctggccggaggcccgggcgcgggcggcctggaggggcccccgcagaagcgtggcatcgtggagcagtgctgtgccagcgtctgctcgctctaccagctggagaactactgtaactag",
"atggcagtttggatccaggctggtgctctgttgttccttttggccgtctccagtgtgaacgctaacgcaggggccccgcagcacctgtgtggatctcatctagtcgatgccctctacctggtctgtggtccaacaggtttcttctataaccccaagagagatgttgaccctcctctgggtttccttcctccaaaatctgcccaggaaactgaggtggctgactttgcatttaaagatcatgccgaggtgataaggaagagaggcattgtagagcaatgttgccacaaaccctgcagtatctttgagctgcagaattactgtaactaa",
"atggccctgtggatgcgcctcctgcccctgctggcgctgctggccctctggggacctgacccggtcccggcctttgtgaaccagcacctgtgcggctcccacctggtggaagccctctacctggtgtgcggggagcgaggcttcttctacacgcccaagacccgccgggaggcagaggacccgcaggtggggcaggtagagctgggcgggggccctggcgcaggcagcctgcagcccttggcgctggaggggtccctgcagaagcgcggcatcgtggagcagtgctgtaccagcatctgctccctctaccagctggagaactactgcaactag",
"atggctctctggatccgatcactgcctcttctggctctccttgtcttttctggccctggaaccagctatgcagctgccaaccagcacctctgtggctcccacttggtggaggctctctacctggtgtgtggagagcgtggcttcttctactcccccaaagcccgacgggatgtcgagcagcccctagtgagcagtcccttgcgtggcgaggcaggagtgctgcctttccagcaggaggaatacgagaaagtcaagcgagggattgttgagcaatgctgccataacacgtgttccctctaccaactggagaactactgcaactag"
]
for i in range(len(sequences)):
sequences[i] = sequences[i][1:50]
# Paramètres
gap=-3
identity=5
transition=3
transversion=0
# Fonctions diverses (qui existent peut-être sous d'autre noms, mais que je n'ai pas trouvées)
def order(l, desc=False):
# Beware: this assumes no ties!
n = len(l)
if desc:
s = -sort(-l)
else:
s = sort(l)
m = {}
for i in range(n):
m[s[i]]=i
r = zeros(n, Int)
for i in range(n):
r[i] = m[l[i]]
return r
def join(l):
a = "";
for i in l:
a = a + str(i)
return a
# Calcul du score de deux séquences
def gamma(u,v):
u = u.lower()
v = v.lower()
if u==v :
return identity
if u in ("a", "g") and v in ("a", "g") :
return transition
if u in ("c", "t") and v in ("c", "t") :
return transition
if u == "-" or v == "-":
return gap
return transversion
def score (a, b, g=gamma, affiche=False) :
na = len(a)
nb = len(b)
gap_start = gap
gap_extend = gap
# Initialisation du tableau de programmation dynamique
d = zeros((na+1,nb+1), Float)
for i in range(na+1):
d[i][0] = gap_start + i*gap_extend
for j in range(nb+1):
d[0][j] = gap_start + j*gap_extend
# Remplissage du tableau de programmation dynamique
for i in range(1,na+1):
for j in range(1,nb+1):
d1 = d[i-1][j-1] + gamma(a[i-1],b[j-1])
d2 = d[i-1][j] + gap
d3 = d[i][j-1] + gap
d[i][j] = max(d1, d2, d3)
if affiche:
# Affichage du début de la matrice
for i in range(15):
for j in range(15):
print string.rjust(repr(d[i][j]),6),
print
print
# Affichage du score
print "Score : " + str(d[na][nb]) + "\n"
# Construction de l'alignement
seqa = list()
seqb = list()
i = na
j = nb
while i>0 and j>0:
d1 = d[i-1][j-1] + gamma(a[i-1],b[j-1])
d2 = d[i-1][j] + gap
d3 = d[i][j-1] + gap
if d[i][j] == d1:
seqa = list(a[i-1]) + seqa
seqb = list(b[j-1]) + seqb
i = i-1
j = j-1
elif d[i][j] == d2:
seqa = list(a[i-1]) + seqa
seqb = list('-') + seqb
i = i-1
elif d[i][j] == d3:
seqa = list('-') + seqa
seqb = list(a[j-1]) + seqb
j = j-1
if i>0:
seqa = list(a[0:(i-1)]) + seqa
if j>0:
seqb = list(b[0:(j-1)]) + seqb
if affiche:
# Affichage de l'alignement
modulo = 50
for i in range(int(len(seqa)/modulo-.001)+2):
for j in range(modulo):
if modulo*i+j == len(seqa):
break
print seqa[modulo*i+j],
print
for j in range(modulo):
if modulo*i+j == len(seqa):
break
if seqa[modulo*i+j] == seqb[modulo*i+j]:
print "|",
else:
print " ",
print
for j in range(modulo):
if modulo*i+j == len(seqa):
break
print seqb[modulo*i+j],
print
print
return {"score": d[na][nb], "a": seqa, "b": seqb}
#score( sequences[0], sequences[1], affiche=True )
# Calcul des scores des alignements deux à deux
n = len(sequences)
d = zeros((n,n),Float)
for i in range(n-1):
for j in range(i+1,n):
print "Calcul de la similarité entre les séquences %d et %d" % (i,j)
d[i,j] = score( sequences[i], sequences[j] )["score"]
for i in range(n):
for j in range(n):
if i>j:
d[i,j] = d[j,i]
# Affichage des scores
print "Matrice des scores :"
for j in range (1,n):
print string.ljust(repr(j),3),
for i in range(j):
print string.rjust(repr(d[i][j]),6),
print
# Recherche de la séquence la plus centrale
dc = resize(array(0,Float),(n,))
for i in range(n):
for j in range(n):
if i != j:
dc[i] = dc[i] + d[i,j]
print "Scores de centralité :"
o = order(dc,desc=True)
for i in range(n):
print string.rjust(repr(i),6),
print
for i in range(n):
print string.rjust(repr(dc[i]),6),
print
for i in range(n):
print string.rjust(repr(o[i]),6),
print
print "Séquence la plus centrale :",
for i in range(n):
if o[i] == 0:
i0=i
break
print i0
print "Ordre des séquences"
o = d[i0]
o[i0] = 1e100
o = order(o,desc=True)
for i in range(n):
print string.rjust(repr(i),6),
print
for i in range(n):
print string.rjust(repr(o[i]),6),
print
for i in range(n):
print string.rjust(repr(d[i0][i]),6),
print
# Il faut calculer l'ordre des séquences...
# Pour cela, on regarde la distance avec
# Alignement successifs par rapport à cette séquence
n = len(sequences)
I = zeros(n)
for i in range(n):
I[o[i]] = i
for i in range(1,n):
print "\nAlignement des séquences %d (centrale) et %d" % (I[0],I[i])
r = score(sequences[I[0]], sequences[I[i]])
print " Début de l'alignement :"
print " " + join(r["a"][1:100])
print " " + join(r["b"][1:100])
print "Propagation des brèches"
# On ajoute des "-" dans toutes les séquences déjà considérées
# Pour la dernière, on se contente de prendre le résultat
sequences[I[i]] = join(r["b"])
# Pour les autres, on compare I[0] avec r["a"] et on ajoute des "-" quand il faut
if i>1:
s = {} # Pas propre...
for j in range(i):
s[j] = {}
k1 = 0
for k2 in range(len(r["a"])):
if r["a"][k2] != sequences[I[0]][k1]:
for j in range(i):
s[j][k2] = '-'
else:
for j in range(i):
if k1<len(sequences[I[j]]):
s[j][k2] = sequences[I[j]][k1]
else:
s[j][k2] = '-'
k1 = k1 + 1
for j in range(i):
sequences[I[j]] = join(s[j].values())
print " Résultat de la propagation :"
for j in range(i+1):
print " " + string.ljust(repr(j),4) + sequences[I[j]][1:100]
# Affichage de l'alignement multiple
print "\nAlignement multiple"
modulo = 120
N = len(sequences[I[0]])
for i in range(int(N/modulo-.001)+2):
for k in range(n):
print sequences[k][ (modulo*i) : (modulo*(i+1)) ]
print
L'algorithme de l'étoile, que nous venons de présenter, fait jouer un rôle central à une séquence. Au lieu de cela, on peut remplacer cette séquence centrale par la "moyenne" des séquences déjà alignées. L'ordre dans lequel on ajoute les séquences n'est pas donné par la distance à une séquence centrale mais par un arbre phylogénétique (obtenu, par exemple, par des méthodes de distance (UPGMA, NJ) -- ce sont des algorithmes de reconstruction phylogénétique qui ne partent pas d'un alignement multiple mais simplement du score des alignements deux par deux). Par conséquent, à chaque étape, on aligne deux alignements multiples (qui correspondent aux deux sous-arbres du noeud considéré).
A FAIRE : URL de clustalw.html
ClustalW construit des alignements multiples de la manière suivante.
1. Faire des alignements deux par deux, en déduire une matrice de distance entre les séquences, puis un arbre, par aggrégation (UPGMA) ou par l'algorithme NJ (Neighbour Joining : voir plus loin).
2. Construire progressivement l'alignement multiple, en commençant par les séquences les plus "proches". Au début, on aligne bien des séquences, mais ensuite, on doit aligner des alignements : au lieu d'avoir un caractère à chaque position, on a un vecteur de probabilités, donnant la répartition des caractères.
Ainsi, la séquence ACTTCAG serait représentée par
Position : 1 2 3 4 5 5 7 ------------------------ A 1 0 0 0 0 1 0 C 0 1 0 0 1 0 0 G 0 0 0 0 0 0 1 T 0 0 1 1 0 0 0
et l'alignement
ACTTCAG | || | AGTT-AC
serait représenté par
A 1 0 0 0 0 1 .5 C 0 .5 0 0 .5 0 0 G 0 .5 0 0 0 0 .5 T 0 0 1 1 0 0 0 - 0 0 0 0 .5 0 0
Les séquences sont pondérées, de sorte que plusieurs séquences très proches ne biaisent pas l'alignement final.
La pénalité pour les brèches (gaps) dépend de leur position -- elles n'apparaissent pas au hasard : on tient compte de la présence d'une brèche voisine et de la nature hydrophile ou hydrophobe des acides aminés.
La matrice de similarité n'est pas la même tout au long de l'algorithme : on s'efforce de choisir la "plus pertinente" selon l'éloignement des séquences (mais on reste dans la même famille, soit PAM, soit BLOSUM).
Il est possible que certaines des séquences, très éloignées, faussent le résultat. Pour contrer cet effet, on peut commencer par faire l'alignement sans la séquence problématique et ne la rajouter qu'à la fin. (Cela revient à modifier l'arbre à la main).
C'est un algorithme d'alignement multiple local -- les précédents étaient globaux.
On commence par faire des alignements locaux deux par deux, en cherchant de longs morceaux de diagonale dans les dotplots, et en les combinant progressivement.
http://www.genopole-lille.fr/fr/formation/cib/doc/anaseq/multiple.pdf http://bioinformatics.org/software/mollin_entry.php3?name=DIALIGN.html http://www.gsf.de/biodv/dialign.html ftp://ariane.gsf.de/pub/unix/dialign2
Je ne les détaille pas.
http://igs-server.cnrs-mrs.fr/%7Ecnotred/Publications/Ps_pdf/coffee.pdf
Pour construire l'alignement multiple, on utilise réellement les alignements deux par deux (qui ne servent donc pas uniquement à construire l'arbre) -- de plus, à chaque étape, on utilise tous ces alignements deux par deux, et pas uniquement ceux concernant les séquences que l'on est en train d'aligner. En outre, on peut utiliser des alignements d'origines diverses (ClustalW, Dialign, alignement structurel).
L'algorithme est le suivant :
1. Par défaut, on calcule deux alignements pour chaque couple de séquences, un global (avec ClustalW), un local (avec lalign).
A FAIRE : il faut avoir parlé de "lalign" avant... (et rappeler ici le nom de l'algorithme).
2. On présente le résultat sous forme de paires << le caractère numéro n1 de la séquence A correspond au caractère numéro n2 de la séquence B >>. On qualifie de "bibliothèque" (library) cet ensemble de paires.
3. On donne des poids à chacune de ces paires (je n'ai pas compris comment) et on combine les alignements de différentes origines (en faisant la somme des poids).
On peut voir ces alignements de caractères et leurs poids comme des contraintes.
4. On étend la bibliothèque, i.e., on modifie les poids des diférentes paires : on augmente le poids d'une paire consistante avec les autres alignements. On peut par exemple prendre les séquences trois par trois : quand les trois alignements coïncident (n1,A,n2,B) (n2,B,n3,C) (n1,A,n3,C), on augmente leur poids, quand ils ne coïncident pas, (n1,A,n2,B) (n2,B,n3,C) (n1,A,n4,C) (avec n3 différent de n4), on diminue leur poids (ou on n'y touche pas).
(On peut aussi présenter cette extension de la bibliothèque ainsi : on regarde les alignements direct de A et C, et les alignements indirects de A et C, qui passent par une autre séquence. Le poids final est la somme des poids.)
5. On construit ensuite un alignement multiple (à partir d'un arbre approché puis par programmation dynamique, comme avec ClustalW, mais ici on n'a même pas besoin de pénalité pour les brèches) qui tient compte de la bibliothèque : pour cela, on utilise une matrice de similarité qui dépend de la position dans les séquences et des poids correspondants.
http://igs-server.cnrs-mrs.fr/%7Ecnotred/Publications/Ps_pdf/tcoffee.pdf
On cherche un motif commun, de longueur L, dans un ensemble de séquences, sans préalablement les aligner. On procède comme suit.
0. Prendre un segment de longueur L sur chacune de ces séquences.
1. Choisir l'une des séquences, déplacer la fenêtre de longueur L et, pour chaque position, calculer le nouveau score de l'alignement avec les segments de longueur L des autres séquences. Choisir la meilleure position (variante : en choisir une au hasard, avec une probabilité proportionnelle à son score ; variante : recuit simulé).
2. Recommencer en 1.
Voir :
http://www.msci.memphis.edu/%7Egiri/compbio/papers/gibbs2.pdf
Exercice : Comment faire pour tenir compte des dégénérescences (nucléotides ou acides aminés inconnus) ?
Exercice : Comment faire si on ne veut pas fixer la largeur de la fenêtre a priori ?
Exercice : Si on lance l'algorithme plusieurs fois on obtient des résultats différents. Comment voir s'il s'agit de motifs réellement différent ou semblables ? Et si ces motifs se chevauchent ?
Nous reverrons l'échantilloneur de Gibbs dans un cadre un peu plus général, comme exemple de simulation de Monte Carlo à l'aide de chaines de Markov.
HMMT PIMA A FAIRE : URL
Plan prévisionnel de cette partie
Motivations, vocabulaire, arbres Parsimonie Petite parsimonie : Algorithme de Fitch Petite parsimonie pondérée : algorithme de Sankoff Branch and Bound, Heuristiques -- A FAIRE Compatibilité -- A FAIRE Méthodes de distances Distance ultramétrique, UPGMA Distance additive, NJ, moindres carrés A FAIRE : calculs de distances... Maximum de vraissemblance A FAIRE MCMC (fusionner avec la partie précédente ? -- non, c'est trop long) Motivations : calculs d'intégrales (=espérance, moyenne), forêt d'arbres phylogénétiques Intégration de Monte-Carlo Intégration de monte-Carlo à l'aide d'une chaine de Markov (parler d'ergodicité) Metropolis Metropolis-Hastings Gibbs Sampler Recuit Simulé Applications Choses plus récentes Utilisation de l'ordre des gènes (réarrangement des gènes)
Le but de la phylogénétique est de retrouver l'arbre de la vie (l'évolution des espèces) : si trois espèces A, B et C ont un ancètre commun (on postule que nous avons tous un ancètre commun), dans quel ordre se sont effectuées les différenciations ? Est-ce d'abord C qui s'est éloigné, B et C ne se différenciant que bien plus tard ? Ou bien B ? Ou A ?
Trouver l'origine d'une nouvelle souche bactérienne ou virale (SRAS, grippe aviaire, listéria, etc. : d'où l'infection est-elle partie ?)
Etudier les migrations humaines.
Rechercher des fraudes dans l'industrie alimentaire (on commence d'ailleurs à trouver des tests génétiques pour analyser l'origine de la viande de porc -- pour le "korobuta", au Japon).
A FAIRE : ajouter l'URL de la page HTML http://www.bbc.co.uk/radio4/science/rams/leadingedge_20040520.ram
Distinguer les gènes orthologues et les gènes paralogues, ce qui permet d'avoir une idée de leur fonction : les gènes orthologues ont généralement la même fonction, les gènes paralogues moins souvent.
Les arbres phylogénétiques sont une première étape d'autres algorithmes, par exemple les alignements multiples.
On peut peut aussi utiliser ces arborescences pour organiser de grandes bases de données de séquences, soit pour y faciliter les recherches, soit pour utiliser moins de place (on commence par reconnaître les familles de séquences puis, pour chaque famille, on calcule une séquence consensus et on se contente de stocker ce consensus et les mutations qui conduisent aux séquences -- on peut gagner jusqu'à 95% de place...).
On peut aussi appliquer les méthodes que nous allons présenter à des données non biologiques, pour étudier, par exemple, l'histoire des langues.
http://www.cis.upenn.edu/~histling/
On peut même combiner des données génétiques et linguistiques.
A FAIRE : l'URL exact (mettre les 2 URLs : la page HTML et le fichier *.ram) http://www.bbc.co.uk/radio4/materialworld/20031127
http://www.iro.umontreal.ca/%7Emabrouk/IFT6291/phylogenie.pdf http://evolution.gs.washington.edu/gs541/2003/ http://www.cs.tufts.edu/comp/150B/lectures/MolEvolPhyloLecture.ppt
Taxon (pluriel anglais : taxa) = les séquences qui apparaissent sur l'arbre, à la fois les séquences dont on part mais aussi les séquences ancestrales reconstruites.
Hypothèse de l'horloge moléculaire : "les vitesses d'évolution ne dépendent pas de l'espèce". Si on fait cette hypothèse, on se retrouve avec un arbre dont toutes les branches sont au même niveau (i.e., un arbre décrivant une distance "ultramétrique").
Gènes homologues : ils ont une origine commune.
Génes orthologues : gènes homologues, si la divergence qui les a créés est une spéciation.
Gènes paralogues : gènes homologues, si la divergence qui les a créés est une duplication.
Gènes xénologues : gènes homologues, créés par un transfert horizontal de gènes (i.e., entre espèces différentes -- les bactéries peuvent ainsi s'échanger des gènes, mais c'est aussi le cas des virus, qui insèrent leur matériel génétique dans celui de leur hôté).
Homoplasie : survenue de mutations identiques dans des branches différentes de l'arbre. Une similarité entre deux séquences peut être due à une homologie (une origine commune) ou à une homoplasie (au hasard).
Distance d'édition entre deux séquences = nombre de différences, i.e., nombre de mésappariements.
Distance ultramétrique : distance mesurée sur un arbre dont toutes les feuilles sont au même niveau. Si l'hypothèse de l'horloge moléculaire est vraie, la distance entre les taxons doit être ultramétrique.
Remarque : étant donnée une distance ultramétrique, il existe un unique arbre qui la décrive.
Distance additive : distance mesurée sur un arbre dont les feuilles ne sont pas nécessairement toutes au même niveau. Que l'hypothèse de l'horloge moléculaire soit vérifiée ou non, la distance entre les taxons doit être additive (sinon, ça n'est pas un arbre). Dans la pratique, elle l'est rarement (mais si elle l'est, c'est signe que les données sont cohérentes) et les algorithmes vont construire un arbre avec une distance additive qui s'approche "le plus" (en un sens souvent vague, qui varie d'un algorithme à l'autre).
Le fait que les distances soient additives a une conséquence un peu surprenante : si on remplace un ensemble de points A1, A2,... par leur "moyenne" A, alors la distance entre A et B est la moyenne des distances entre les Ai et B. Si on regarde ce que ça signifie sur un dessin, on constatera que c'est vrai pour la distance euclidienne (i.e., usuelle) en dimension un (sur une droite) mais faux en dimensions supérieures (dans le plan, dans l'espace).
Arbre enraciné : graphe sans cycle, dans lequel on a choisi un noeud, qu'on appellera la racine. On peut alors dessiner l'arbre avec la racine tout en haut, les feuilles en bas et les branches qui descendent. La racine représentera l'ancètre le plus lointain.
Arbre non enraciné : graphe sans cycle. Comme on n'a pas choisi de racine, on ne sait pas où est l'ancètre le plus ancien. Quand on le dessine, on mettra généralement les noeuds au milieu et les feuilles autour (Il n'y a pas de haut ni de bas).
Outgroup : Pour enraciner un arbre non enraciné, on peut prendre une taxon un peu plus éloigné et le rattacher à l'arbre, une fois celui-ci construit : comme il est plus éloigné, on peut choisir l'endroit où il est attaché à l'arbre comme ancètre le plus ancien, i.e., comme racine. On dit que ce taxon ajouté à la fin est un outgroup. (On peut aussi considérer des outgroups constitués de plusieurs taxons.)
Midpoint rooting : on peut aussi enraciner un arbre en prenant les deux feuilles les plus éloignées et en mettant la racine au milieu. (Sous l'hypothèse de l'horloge moléculaire, ça marche, mais en général, c'est un peu risqué...)
Remarque : une autre manière d'enraciner les arbres consiste à placer la racine au milieu de l'arrête la plus longue.
Phénétique : analyse phylogénétique basée sur les distances entre les taxons (NJ, UPGMA). En considérant simplement les distances, on perd beaucoup d'information.
Phénogramme ou dendogramme : résultat d'une analyse phylogénétique basée sur les distances.
Cladistique : analyse phylogénétique basée sur les caractères observés, généralement les séquences elles-mêmes (méthode de parsimonie ou maximum de vraissemblance).
Cladogramme : résultat d'une analyse phylogénétique cladistique (parsimonie ou du maximum de vraissemblance) ; on considère les taxons dans leur intégralité et pas seulement les distances qui les séparent.
Clade : groupe monophylétique, i.e., sous-arbre, i.e., groupe contenant un ancètre et tous ses descendants.
OTU (Operational Taxonomic Unit) : feuille de l'arbre
Polytomie ou multifurcation : noeud de l'arbre avec plus de trois branches.
Synapomorphie : des caractères homologues provenant d'un ancètre commun immédiat
Plésiomorphie : des caractères homologues provenant d'un ancètre commun lointain
A FAIRE :
D'autres mots qui ressemblent à "monophylétique".
Voir aussi :
http://www.infobiogen.fr/services/deambulum/fr/bioinfo3.html#PHYLOGENIE
On utilisera essentiellement Phylip.
http://evolution.genetics.washington.edu/phylip.html
Ça n'est pas un seul programme, mais plusieurs. Mentionnons :
dnapars Parcimonie pour des données nucléiques
protpars Parcimonie pour des données protéiques
dollop Parcimonie (Dollo) pour des données quelconques
penny Parcimonie (Wagner, Camin-Sokal) pour des données quelconques
dnadist Calcul d'une matrice de distance pour des données nucléiques
protdist Calcul d'une matrice de distance pour des données protéiques
neighbor Calcul d'un arbre par neighbor-joining ou UPGMA
fitch Fitch-Margoliash
kitsch FM avec horloge moléculaire
seqboot Construction d'échantillons de bootstrap
consense Construire un arbre à l'aide des morceaux qui
apparaissent le plus souvent dans uine liste d'arbres
(on a par exemple obtenu cette liste d'arbres à l'aide
de la méthode de parcimonie avec du bootstrap)
dnaml Maximum de vraissemblance, données nucléaires
dnamlk Idem, avec une horloge moléculaire
protml Maximum de vraissemblance, données protéiques
protmlk Idem, avec une horloge moléculaire
restdist Calcul de distances à partir de données concernant des sites de restriction
restml Maximum de vraissemblance à partir de sites de restrictiona
drawtree Tracé (Postscript) d'un arbre non enraciné précédemment obtenu
drawgram Tracé d'un arbre enracinéet aussi :
Pour avoir des résultats fiables, on utilisera les options "Jumble" (qui change l'ordre initial des séquences, pour les programmes dont le résultat est suceptible d'en dépendre) et "global rearrangement" (qui, une fois un arbre complet trouvé, essaye tous les bouturages possibles pour l'améliorer).
Il existe un autre logiciel : PAUP (commercial, essentiellement sous Mac).
http://www.bioinf.org/molsys/data/paupv1.pdf http://www.bioinf.org/molsys/data/paupv2.pdf
Les données à partir desquelles on fait une analyse phylogénétique peuvent être de natures très différentes.
On peut prendre des caractères visibles : présence d'ailes, d'écailles, de sabots, nombre d'estomacs, régime alimentaire, etc.
Très souvent, il s'agit de séquences d'ADN ou de protéines -- en particulier l'ADN mitochondrial. (On prend souvent le cytochrome C ou l'hémoglobine.) C'est un bon choix, car ces séquences évoluent de manière aléatoire et sont peu soumises à la sélection : on peut donc facilement modéliser leur évolution.
Il peut aussi s'agir de cartes génétiques, i.e., de marqueurs. Par exemple, la longueur des fragments obtenus après digestion par une enzyme de restriction (RFLP).
Il peut aussi s'agir de distances, estimées expérimentalement par la température d'hybridation de brins d'ADN d'un même gène chez des espèces différentes.
Il peut aussi s'agir, carrément, de l'ordre des gènes. En effet, on constate que la souris et l'homme ont à peu près les mêmes gènes, mais dans un ordre différent. Si on plus de deux espèces, on peut essayer de retrouver l'ordre des réarrangements du génome.
Exercice : lire la documentation des logiciels mix, restml, gendist de Phylip : quel type de données manipulent-ils ?
Dans une analyse phylogénétique basée sur des séquences, comment traiter les délétions/insertions ? Faut-il considérer le caractère "-" (indiquant une délétion) comme un cinquième nucléotide ? Faut-il simplement les oublier ? Aucune de ces solutions n'est satisfaisante.
Les vitesses d'évolution peuvent différer selon les espèces : l'hypothèse de l'horloge moléculaire, qui simplifie grandement l'analyse, n'est pas vérifiée.
Il arrive que certains gènes soient dupliqués -- ça complique l'analyse.
Il y a parfois des réarrangements au sein même des gènes : comment en tenir compte ?
La représentation d'une analyse phylogénétique sous forme d'arbre n'est pas toujours possible : ainsi, il arrive qu'un gène A provienne de la concaténation de deux gènes B et C. (Dans certains cas, on peut quand-même représenter la phylogénie sous forme d'un graphe -- on accepte la présence de cycles).
Autre exemple : le matériel génétique est parfois transféré horizontalement, i.e., d'une espèce à l'autre. C'est le cas par exemple des virus, qui intègrent leur patrimoine génétique à leur victime ; c'est aussi le cas de certaines symbioses -- la plus connue étant celle (putative) donnant naissance aux mitochondries, chez les eucaryotes. Dans ce cas, l'évolution n'est plus représentée par un arbre mais par un graphe.
A moins que les séquences soient très proches, il y a plusieurs évolutions possibles. (Pour ce problème, on verra qu'il est possible de demander non pas un seul arbre, mais plusieurs (en bootstrappant les données de départ), à l'aide desquels on pourra construire un nouvel arbre, qui ne sera plus binaire.) Ainsi, quand les séquences ont beaucoup évolué, il a pu y avoir des réversions.
On est amené à supposer que les mutations en deux endroits différents sont indépendantes -- pour l'ARN, c'est complètement faux, car sa structure secondaire est primordiale ; il en va probablement de même pour les gènes codant des protéines, car la structure secondaire de l'ARNm peut aussi avoir son importance.
La représentation naïve (chaque noeud est un couple de pointeurs vers ses noeuds fils) n'est pas toujours très pratique : en particulier, on aimerait pouvoir considérer des arbres non enracinés ou, quand on considère des arbres enracinés, pouvoir facilement changer la place de la racine. Or cette représentation imposerait une réécvriture complète de l'arbre.
Voici une autre représentation des arbres binaires, qui ne présente pas ces problèmes.
Chaque noeud est constitué de trois couples de pointeurs, un pointeur vers le couple suivant, un autre vers un autre noeud. Cette représentation n'est pas celle qui vient à l'esprit en premier, mais elle présente deux avantages : d'une part, elle se généralise très facilement aux arbres non binaires (un noeud n-aire est constitué de n+1 couples de pointeurs), d'autre part, l'arbre n'est pas enraciné.
(Personnellement, je préfère représenter un arbre comme un ensemble de noeuds, chaque noeud étant une liste de pointeurs sur les noeuds voisins).
Pour stocker un arbre dans un fichier texte, on peut utiliser la notation suivante :
((A,B),C)
qui représente l'arbre
---- | | | | --- | | | | A B C
On peut aussi ajouter la longueur de chaque branche :
((A:1,B:1):2,C:4)
La méthode de parsimonie essaye de trouver l'arbre le plus parcimonieux, i.e., celui qui explique l'état actuel avec le moins de mutations possible.
Les méthodes de distances commencent par calculer la distance entre les taxons et essaye de trouver l'arbre (ultramétrique ou non) qui approche le mieux cette distance.
La méthode du maximum de vraissemblance est basée sur un modèle évolutif et cherche l'évolution la plus probable.
On peut résumer cela dans un tableau :
Méthode Proximité des séquences Hypothèse de l'horloge moléculaire ---------------------------------------------------------------------------------------- Parsimonie Très proches Oui Méthodes de distance Proches Oui (UPGMA) ou Non (FM, NJ) Maximum de vraissemblance Peu importe Non
Etant donné un ensemble de séquences ou de caractères mesurés chez différentes espèces, on cherche l'arbre qui explique ces changements en faisant intervenir le moins de mutations possible -- l'arbre le plus parcimonieux.
On remarquera que cela exclut certaines évolutions possibles, en particulier, les réversions. On réservera donc les méthodes de parsimonie aux situations où il est raisonnable de penser qu'il n'y a que peu de mutations.
Un site est informatif s'il y a au moins deux individus dans un état et deux autres dans un autre état. Dans un alignement multiple, il y a en fait très peu de sites informatifs.
Par exemple, dans l'alignement suivant, nous avons mis une étoile au dessous des sites informatifs. Il y en a peu...
BABBBEAAABCBCBEAABAE CDDEABBECDCEDCABCCDD ECDCDACCCABEDBBACBBB CECDEDAABDBCACADBDDC DECCBCEDECCBDDCCEABE * ** * *
Remarque : Cette notion de site informatif est propre aux méthodes de parcimonie ; un site non informatif peut quand-même nous apporter quelque-chose si on utilise des méthodes de distance ou le maximum de vraissemblance.
Remarque : l'exemple ci-dessus a été obtenu ainsi :
n<-5
k<-20
m<-matrix(sample(LETTERS[1:5],n*k,replace=T),nr=n)
for (i in apply(m,1,paste, collapse='')) {
cat(i)
cat("\n")
}
v <- apply(m,2,function(x){sum(table(x)>1)>1})
cat(paste(ifelse(v,"*"," "),collapse='')); cat("\n")
Il s'agit de compter le nombre minimal d'évènements (mutations) dans _un_ arbre pour expliquer les séquences ou caractères observés, i.e., de calculer son score de parsimonie. (Ensuite, il faudrait faire ça pour tous les arbres : c'est la grande parsimonie.)
On considère un site à la fois.
On parcourt l'arbre en commençant par les feuilles et en remontant progressivement vers la racine (on peut le programmer ainsi : parcourrir l'arbre en largeur, se souvenir de l'ordre dans lequel on a visité les noeuds (i.e., d'abord la racine, puis les noeuds les plus internes, et enfin les feuilles), renverser cet ordre).
On décore chaque feuille par la valeur correspondante du site.
Chaque noeud sera ensuite décoré à l'aide des décorations de ses deux fils : soit par l'intersection de ces deux décorations, si elle est non vide, soit par leur union.
C A C A G
| | | | |
--CA-- | --AG--
| | |
| ---ACG---
| |
-------C-------
|
|On compte ensuite le nombre d'union, i.e., le nombre d'intersections vides : c'est le score de parsimonie recherché, i.e., le nombre minimal de mutations.
C A C A G
| | | | |
--*--- | --*---
| | |
| ----*----
| |
---------------
|
|
Nombre minimal de mutations = 3
Il s'agit toujours de calculer le score de parsimonie d'un arbre, mais cette fois-ci, les différentes mutations ont chacune un "poids" différent : on a recours à de la programmation dynamique.
On calcule :
S(j,i) = score du sous-arbre partant du noeud j s'il a la valeur i
en commençant par les feuilles :
| 0 si la feuille j a la valeur i
S(j,i) = |
| infini sinonen descentant progressivement :
S(j,i) = Min( c(i,k) + S(gauche,k) , k) + Min( c(i,k) + S(droite,k), k )
où gauche et droite sont les deux noeuds fils de j et c(i,k) est le coût d'une mutation de i en k (par exemple, 0 si i=k et 1 sinon).
Les deux algorithmes précédents se contentaient de donner le score de parsimonie d'_un_ arbre : il reste à trouver celui qui donnera le meilleur score.
On pourrait penser énumérer _tous_ les arbres, calculer le score de parsominie de chacun d'entre eux et retenir le meilleur. Mais il y en a beaucoup trop.
En fait, c'est un problème NP-complet : on aurra donc recours, le plus souvent, à des heuristiques.
Si le nombre d'espèces (i.e., de feuilles de l'arbre) n'est pas trop élevé, on peut énumérer tous les arbres.
On peut par exemple procéder récursivement : il n'y a qu'un seul arbre (non enraciné) à trois feuilles ; si on connait tous les arbres à n feuilles, on en déduit tous les arbres à n+1 feuilles : pour chaque branche de chaque arbre à n feuilles, on ajoute une (n+1)-ième feuille sur cette branche.
On peut représenter l'ensemble des arbres à examiner par un arbre (oui... un arbre dont les feuilles sont labellées par des arbres...) : la racine représente l'unique arbre à trois feuilles, ses fils sont les arbres (à 4 feuilles) obtenus en ajoutant une feuille, ses petits fils sont les arbres obtenus en ajoutant encore une feuille, etc.
La recherche exhaustive consiste alors à parcourir cet arbre (en profondeur).
L'arbre précédent, présentant toutes les topologies possibles, est immense. Mais en fait, il n'est pas nécessaire d'en explorer toutes les ramifications : en effet, pour chaque noeud de l'arbre, on peut calculer un score de parsimonie, que l'on peut voir comme une "parsimonie partielle", une étape dans le calcul de la parsimonie de tous les arbres labellant les feuilles descendant du neoud en question. Si cette parsimonie partielle est plus mauvaise que la parsimonie d'arbres (complets) déjà examinés, il n'est pas nécessaire d'aller jusqu'aux feuilles : on peut élaguer l'arbre des topologies à examiner à partir du noeud en question.
C'est toujours un algorithme exhaustif et pas une heuristique : si l'algorithme se termine, on aura bien l'arbre de meilleure parsimonie, car les topologies qui n'ont pas été examinée étaient plus mauvaises. Le Branch and Bound a en fait permi d'examiner (et de rejeter) plusieurs topologies en même temps.
On peut aussi construire l'arbre progressivement, en ajoutant les feuilles une par une (dans l'ordre dans lequel elles se présentent), en choisissant à chaque fois l'endroit le plus parcimonieux où les mettre.
C'est un algorithme glouton, car on ne revient jamais sur les choix qu'on a faits.
Pour améliorer le résultat, on peut essayer de réorganier l'arbre, par exemple en détachant un sous-arbre et en essayant de le mettre ailleurs (SPR : Subtree Pruning and Regrafting), par exemple en retirant une arrête (ce qui produit deux composantes connexes) et en en ajoutant une ailleurs (pour relier ces deux composantes) (TBR : Tree Bissection and Reconnection), par exemple en effectuant des mouvements de "fusion" (NNI : Nearest Neighbor Interchange).
On peut aussi recommencer en changeant l'ordre dans lequel on ajoute les feuilles.
Construire un arbre avec 3 espèces (au hasard), puis ajouter les autres, une par une, à l'endroit qui donne l'arbre le plus parcimonieux. On obtient un arbre dont le score de parsimonie est proche du score optimal.
La parsimonie de Wagner cherche l'arbre le plus parcimonieux parmi tous les arbres possibles -- c'est la parsomonie que nous avons vue.
La parsimonie de Camin-Sokal limite la recherche aux arbres ne comportant pas de réversion (i.e., pas de mutations de la forme A-->B-->A).
La parsimonie de Dollo limite la recherche aux arbres sans convergence.
La méthode de compatibilité est une simplification de la parsimonie, utile surtout quand les caractères considérés sont binaires.
On dit qu'un caractère (ou site) qui ne prend que deux valeurs (i.e., que deux allèles) est compatible avec un arbre s'il existe une branche de l'arbre le coupant en deux, avec d'un côté l'une des valeurs et de l'autre l'autre. Il revient au même de dire qu'on peut expliquer ce caractère à l'aide d'une unique mutation sur cet arbre.
On généralise cette définition à des caractères qui prennent plus de deux valeurs.
Le reconstruction phylogénétique par compatibilité consiste à chercher l'arbre compatible avec le plus de sites possibles.
On peut adapter n'importe quel algorithme de parsimonie pour rechercher un arbre satisfaisant au critère de compatibilité : il suffit, lors du calcul du score d'un arbre pour un caractère, de remplacer le score par 2 s'il est supérieur à 2.
Mais c'est un peu long : le critère de compatibilité est une simplification de la parsimonie, on voudrait un algorithme plus rapide. Nous allons présenter l'algorithme de Wilson.
On peut rechercher si deux caractères A et B sont compatibles avec un même arbre : il suffit de remplir le tableau 2x2 suivant. On met un 0 dans la case (i,j) s'il n'y a pas d'individu avec A=i et B=j et on met 1 sinon. Si le tableau contient au moins un zéro, les caractères sont compatibles, sinon, ils ne le sont pas.
Exemple :
A B A=0 A=1 1 1 ----------- 1 0 B=0 | 1 1 0 0 B=1 | 0 1 0 0 1 0 Ils sont compatibles
Théorème : un ensemble de caractères binaires est compatible avec un arbre si et seulement si ils sont compatibles deux à deux.
Pour trouver le plus grand ensemble de caractères compatibles, il suffit donc de considérer le graphe dont les sommets sont les caractères et comportant une arrête entre deux sommets si les caractères sont compatibles. Un ensemble de caractères compatible est une clique : on s'est donc ramené à chercher une clique maximale dans un graphe.
Une fois notre clique maximale trouvée, il reste à construire l'arbre. On commence par un graphe en forme d'étoile, avec tous les individus attachés à un seul noeud. Puis on considère les caractères de la clique maximale un par un : à chaque caractère, on est amené à séparer un noeud en deux.
http://www.cs.technion.ac.il/Labs/cbl/teaching/uw/lect16.pdf
Rappelons-le encore, l'algorithme que nous venons de présenter ne marche que pour des caractères binaires.
On peut aussi reprendre l'idée de la parsimonie pour aligner des séquences : dans la discussion précédente, nous nous étions limité à un seul type d'évènement pour passer d'une séquence à une autre, la mutation -- on peut aussi ajouter des insertions et des délétions.
La parsimonie produit alors à la fois un alignement multiple et une phylogénie.
On a fait l'hypothèse que les taux de mutations sont les mêmes dans toutes les branches de l'arbre. Si cette hypothèse n'est pas vérifiée, les résultats obtenus sont faux.
La parsimonie ne tient pas compte des réversions : un évènement de la forme A --> B --> A a tendance à être rejeté par les méthodes de parsimonie. (Si les séquences sont vraiment très proches, ce n'est pas trop grave).
La parsimonie est très sensible au biais des codons.
C'est lent...
dnapars, dnapenny, dnacomp, dnamove, dnainvar protpars dnacomp
Après la parsimonie, les méthodes de distance constituent l'une des grandes familles d'algorithmes de reconstruction phylogénétique. Ici, le problème est le suivant : étant donnée la matrice des distances entre les différentes feuilles d'un arbre (que l'on ne connait pas, dont chaque arrête a une longueur et dans lequel on mesure les distances entre les feuilles le long des arrêtes), peut-on reconstituer cet arbre ?
Avant d'utiliser ces algorithmes, il nous faut donc calculer la distance entre les séquences.
L'idée générale derrière tous ces calculs de distance est la suivante : on choisit un modèle décrivant l'évolution des séquences avec le temps (par exemple, « les mutations sont des processus poissonniens, dont la probabilité peut être donnée a priori ou être estimée à partir des données » -- ces probabilités peuvent aussi dépendre du type de mutation) et on calcule le temps "le plus probable", "le plus vraissemblable" -- on parle d'ailleurs de maximum de vraissemblance.
(Le lecteur statisticien pourra remplacer l'"estimateur du maximum de vraissemblance" par n'importe quel autre estimateur.)
Nous allons donc présenter ces différents modèles d'évolution des séquences.
http://conferences.computer.org/bioinformatics/CSB2003/NOTES/Corey.ppt
Très souvent, toutes les mutations ne sont pas considérées de la même manière : ainsi, l'apparition d'un A à la place d'un C est très différente de l'apparition d'un T à l' place d'un C. Le modèle est alors une chaîne de Markov. (Ai-je déjà dit que c'était important, les chaines de Markov ?)
David Bryant, Estimating Evolutionary Distances b etween Sequences http://www.crm.umontreal.ca/Bio2003/pdf/Bryant.pdf
Dans la discussion précédente, nous avons utilisé la "distance" entre les séquences : mais on a différentes manières, plus ou moins fines, pour la calculer.
On peut tout simplement compter le nombre de mesappariements entre les deux séquences dans l'alignement multiple : c'est la distance d'édition.
On peut aussi utiliser un modèle stochastique : on suppose que les bases subissent des substitutions de manière indépendante, suivant un processus de Markov. Dans le modèle de Jukes-Cantor, la probabilité d'un changement ne dépend pas des bases ; dans le modèle de Kimura, elle en dépend : A<->G et C<->T sont plus fréquentes, on a donc deux probabilités selon que la base reste une purine ou une pyrimidine, ou si sa nature change.
On en déduit la distance évolutive, i.e., l'estimation du nombre de substitutions (en fonction de la fréquence des différences effectivement observées).
On peut imaginer des modèles plus fins : par exemple, le modèle de Jin et Nei (le temps est continu (en fait, dans les modèles précédents, il l'était déja) et les modifications suivent une (des) lois gamma).
Notons u la probabilité d'une mutation (par unité de temps) et supposons que les trois mutations possibles en chaque site ont la même probabilité a=u/3. Lors du calcul de la distance, on tient compte de la possibilité de mutations multiples.
A G C T A a a a G a a a C a a a T a a a
Après calculs, on trouve :
d = -3/4 * ln( 1 - 4D/3 )
où D est la proportion de nucléotides qui diffèrent dans les deux séquences.
C'est comme le modèle de Jukes-Cantor, mais les transitions et les transversions ont une probabilité différente.
A G C T A b a a G b a a C a a b T a a b
Ici, la formule finale ne fait pas intervenir que la proportion de nucléotides qui diffèrent, mais la proportion de transitions et la proportion de transversions entre les deux séquences.
(je n'ai pas la formule sous la main)
F81 (Felsenstein) : toutes les substitutions ont la même probabilité (comme avec le modèle de Jukes-Cantor), mais les bases ont des fréquences différentes.
Hasegawa-Kishino-Yano (HKY) : comme le modèle de Kimura, mais on précise en plus la fréquence limite des quatre bases.
A G C T
A (a+b) p_G a p_C a p_T
G (a+b) p_A a p_C a p_T
C a p_A a p_G (a+b) p_T
T a p_A a p_G (a+b) p_CF84 : Idem, mais on précise simplement la fréquence limite des purines et des pyrimidines.
Tamura-Nei : généralisation des deux modèles précédents.
GTR (General Time-Reversible model) : la probabilité d'observer une (suite de) mutation(s) de A vers T est la même que de T vers A.
A G C T
A a p_G b p_C c p_T
G a p_A d p_C e p_T
C b p_A d p_G f p_T
T c p_A e p_G f p_CModèle général : on précise la probabilité de mutation de chaque acide aminé en chaque autre (12 paramètres).
Les modèles précédents sont en fait des chaines de Markov : quand un site est dans un état donné (disons, A), on connait les probabilités pour qu'il mute en C, G, T à l'instant suivant.
A FAIRE : dessin.
C'est néanmoins plus compliqué que les chaines de Markov auxquelles vous êtes êtes peut-être habitués, pour deux raisons.
D'une part, il y a une chaine de Markov pour chaque site. La séquence complète est aussi décrite par une chaine de Markov (le produit cartésien des chaines de Markov des sites), mais elle est énorme : on travaillera généralement avec les chaines de chacun des sites.
D'autre part, avec les chaînes de Markov, le temps est discret -- ce n'est pas le cas ici, le temps entre deux mutations n'est pas un multiple entier d'une certaine unité. On peut néanmoins garder le même formalisme, en prenant comme unité de temps "dt" (si l'analyse non standard était plus répandue dans l'enseignement secondaire, on pourrait dire que "dt" est un "infiniment petit") et en remplaçant une probabilité de transition u par "u dt". Le résultat ne s'appelle plus "chaine de Markov", mais processus de Markov.
[L'analyse non standard, ce sont les mathématiques avec des nombres réels infiniment grands et infiniment petit -- les physiciens les utilisent courremment, en croyant être peu rigoureux -- ils ont tors, c'est tout à fait correct.]
Dans les discussions précédentes, on supposait que les probabilités de mutation étaient les mêmes en chaque point de la séquence : ce n'est pas le cas. On constate que certains endroits des séquences tendent à accumuler les mutations (on les qualifie de "hot spots") alors que d'autres en contiennent très peu (ce sont probablement des zônes très importantes pour la fonction de la protéine).
On peut imaginer des modèles dans lesquels la probabilité des mutations n'est pas la même d'un site à l'autre : on suppose alors que ces taux de mutation sont des variables aléatoires, qui suivent une certaine loi. On utilise généralement une distribution Gamma pour modéliser ces taux (quand on veut modéliser une quantité positive comme un taux, on pense souvent à la loi exponentielle ou à ses généralisations -- comme, justement, la loi Gamma).
(Je ne sais pas d'où vient cette formule...)
d = - ln det F
où F est la matrice de divergence des deux séquences, de taille 4x4, qui donne la fréquence des appariements de bases dans l'alignement des deux séquences.
Dans la plupart des modèles précédents, on supposait que les deux séquences avaient la même composition, i.e., que la matrice F était symétrique.
Les modèles précédents permettaient de calculer des distances entre des séquences nucléiques, en tenant compte des propriétés physico-chimiques des nucléotides, mais en obérant complètement le code génétique. Pour des séquences non codantes, c'est effectivement ce qu'il convient de faire.
Par contre, pour des séquences codantes, la survenue d'une mutation dépend des acides aminés correspondants : une mutation qui ne change pas l'acide aminé est plus probable qu'une mutation qui le change ; une mutation qui change un acide aminé en un autre avec des propriétés physico-chimiques semblables est plus probable qu'une mutation qui en fait apparaître un complètement différent.
Ces remarques sont valables tant pour les séquences nucléiques codantes que (plus simplement), pour les séquences protéiques.
Dans le cas du calcul de la distance entre deux séquences protéiques prendre le score de leur alignement (obtenu par un algorithme de programmation dynamique, par exemple avec une matrice PAM ou JTT).
Les méthodes de distance tentent de reconstruire un arbre phylogénétique (un arbre, avec une longueur sur chacune des branches) en connaissant simplement la distance entre les feuilles (on mesure la distance le long des branches de l'arbre).
Ces méthodes varient quant aux propriétés de l'arbre reconstruit (la méthode UPGMA suppose que l'arbre est ultramétrique, les autres méthodes sont capable de retrouver un arbre quelconque) et quant à la robustesse face à des données bruitées.
Un arbre dont chaque branche est munie d'une longueur définit une matrice de distances entre ses feuilles.
Dire que la distance provient d'un arbre revient à dire que cette distance est additive, i.e., si A, B, C, D sont comme dans la figure suivante
A FAIRE : DESSIN A et B ont un ancètre commun, C et D ont un ancètre commun, ces ancètres communs sont reliés.
alors
d(A,B) + d(C,D) <= d(A,C) + d(B,D) = d(A,D) + d(B,C)
(C'est une propriété de la distance sur la droite réelle (en dimension un), qui ne se généralise pas aux dimensions supérieures (plan, etc.).)
(En d'autres termes, si on a simplement la matrice des distances et si elle est additive, alors elle provient d'un arbre -- on ne nous dit pas, pratiquement, comment le construire (voir plus loin), mais on sait qu'il existe.)
Si l'hypothèse de l'horloge moléculaire est vérifiée, toutes les feuilles sont à la même distance de la racine. Cette propriété se lit sur la matrice de distance : elle est "ultramétrique", i.e., pour tous A, B, C,
d(A,C) <= max(d(A,B),d(B,C))
L'algorithme UPGMA permet de reconstruire l'arbre correspondant à une distance ultramétrique -- on doit supposer qu'elle est ultramétrique.
Le principe est simple :
1. Prendre les deux feuilles les plus proches, appelons-les A et B. 2. Les enlever de l'ensemble des feuilles et les remplacer par une seule feuille, AB, dont les distances avec les autres feuilles sont obtenues ainsi : d(AB,C) = ( d(A,C) + d(B,C) ) / 2.
Le fait que l'arbre soit ultramétrique (i.e., on a fait l'hypothèse de l'horloge moléculaire) entraine que A et B forment bien un groupe monophylétique (dont AB est l'ancètre).
Le fait que la distance provienne d'un arbre fait qu'elle est additive et que la distance de AB aux autres feuilles est bien la moyenne des distances à A et à B.
Cet algorithme n'est pas très satisfaisant, pour deux raisons : d'une part, il fait l'hypothèse de l'horloge moléculaire, qui est fausse, d'autre part, il est très sensible au bruit. Il y a de bien meilleurs algorithmes. Néanmoins, la première étape dans un alignement multiple est souvent le calcul d'un arbre phylogénétique approché : quand un résultat approximatif suffit, l'algorithme UPGMA peut être utile.
Exemple d'implémentation :
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#define NOMBRE_DE_SEQUENCES 7
#define MAX_LENGTH 1000
char * sequences [] = {
"atggcagtgtggcttcaggctggtgctctgttggtcctgttggtcgtgtccagtgtaagcactaacccaggcacaccgcagcacctgtgtggatctcatctggtcgatgccctttatctggtctgtggcccaacaggcttcttctacaaccccaagagagacgttgagccccttctgggtttccttcctcctaaatctgcccaggaaactgaggtggctgactttgcatttaaagatcatgccgagctgataaggaagagaggcattgtagagcagtgctgccacaaaccctgcagcatctttgagctgcagaactactgtaactga",
"atggcagcgctctggctccaggccttctccctgctcgtcttaatgatggtttcgtggccgggctcccaggccgtcggtgggccacagcacctgtgcggctcccacctggtggatgccctgtacctggtctgtggggacagaggcttcttctacaaccccaggagagatgtggaccctctgcttggtttcctccctccaaaggcaggtggtgctgtggtgcaaggtggtgagaatgaagtgaccttcaaagaccagatggaaatgatggtgaagcgaggcattgtggaggaatgctgtcacaaaccctgtaccatcttcgacctgcagaactactgcaactga",
"atggccctgtggatgcacctcctgcccctgctggcgctgctggccctctggggacccgagccagccccggcctttgtgaaccagcacctgtgcggcccccacctggtggaagccctctacctggtgtgcggggagcgaggtttcttctacgcacccaagacccgccgggaggcggaggacctgcaggtggggcaggtggagctgggtgggggctctatcacgggcagcctgccacccttggagggtcccatgcagaagcgtggcgtcgtggatcagtgctgcaccagcatctgctccctctaccagctgcagaactactgcaactag",
"atggccctgtggacacgcctgcggcccctgctggccctgctggcgctctggcccccccccccggcccgcgccttcgtcaaccagcatctgtgtggctcccacctggtggaggcgctgtacctggtgtgcggagagcgcggcttcttctacacgcccaaggcccgccgggaggtggagggcccgcaggtgggggcgctggagctggccggaggcccgggcgcgggcggcctggaggggcccccgcagaagcgtggcatcgtggagcagtgctgtgccagcgtctgctcgctctaccagctggagaactactgtaactag",
"atggcagtttggatccaggctggtgctctgttgttccttttggccgtctccagtgtgaacgctaacgcaggggccccgcagcacctgtgtggatctcatctagtcgatgccctctacctggtctgtggtccaacaggtttcttctataaccccaagagagatgttgaccctcctctgggtttccttcctccaaaatctgcccaggaaactgaggtggctgactttgcatttaaagatcatgccgaggtgataaggaagagaggcattgtagagcaatgttgccacaaaccctgcagtatctttgagctgcagaattactgtaactaa",
"atggccctgtggatgcgcctcctgcccctgctggcgctgctggccctctggggacctgacccggtcccggcctttgtgaaccagcacctgtgcggctcccacctggtggaagccctctacctggtgtgcggggagcgaggcttcttctacacgcccaagacccgccgggaggcagaggacccgcaggtggggcaggtagagctgggcgggggccctggcgcaggcagcctgcagcccttggcgctggaggggtccctgcagaagcgcggcatcgtggagcagtgctgtaccagcatctgctccctctaccagctggagaactactgcaactag",
"atggctctctggatccgatcactgcctcttctggctctccttgtcttttctggccctggaaccagctatgcagctgccaaccagcacctctgtggctcccacttggtggaggctctctacctggtgtgtggagagcgtggcttcttctactcccccaaagcccgacgggatgtcgagcagcccctagtgagcagtcccttgcgtggcgaggcaggagtgctgcctttccagcaggaggaatacgagaaagtcaagcgagggattgttgagcaatgctgccataacacgtgttccctctaccaactggagaactactgcaactag"
};
char * nom_sequences [] = {
"AF036326.INS",
"AF038123.PE1",
"ATRINS.PE1",
"BOVPPI.PE1",
"CCINS01.PE1",
"CEPPINS.PE1",
"CHKINS3.PE1"
};
////////////////////////////////////////////////////////////
// Code de la dernière fois
#define TRANSITION 1
#define TRANSVERSION 0
#define IDENTITY 3
#define GAP -5
int gamma(char i, char j) {
if(i==j)
return IDENTITY;
if ( (i=='a' && j=='g') || (i=='c' && j=='t') )
return TRANSITION;
if ( i=='-' || j=='-' )
return GAP;
return TRANSVERSION;
}
int max(int a, int b, int c) {
if( a >= b && a >= c )
return a;
if( b >= a && b >= c )
return b;
return c;
}
int min(int a, int b) {
return (a>b) ? b : a;
}
double score(char * chaine1, char * chaine2) {
unsigned int i,j;
double score [MAX_LENGTH][MAX_LENGTH];
assert( strlen(chaine1) < MAX_LENGTH );
assert( strlen(chaine2) < MAX_LENGTH );
for (i=0; i<strlen(chaine1); i++)
for (j=0; j<strlen(chaine2); j++)
score[i][j] = 0;
for (i=0; i<strlen(chaine1); i++)
score[i][0] = (1.0*i)*GAP;
for (j=0; j<strlen(chaine2); j++)
score[0][j] = (1.0*j)*GAP;
for (i=0; i<strlen(chaine1); i++) {
for (j=0; j<strlen(chaine2); j++) {
double d = score[i][j] + gamma(chaine1[i],chaine2[j]) ;
double p = score[i][j+1] + GAP;
double q = score[i+1][j] + GAP;
score[i+1][j+1] = max(d,p,q);
}
}
unsigned int n=10;
for (i=strlen(chaine1)-n; i<=strlen(chaine1); i++) {
for (j=strlen(chaine2)-n; j<=strlen(chaine2); j++) {
printf("%3.0f ", score[i][j]);
}
printf("\n");
}
return score[strlen(chaine1)][strlen(chaine2)];
}
////////////////////////////////////////////////////////////
int min2(int a, int b) {
return (a>b) ? b : a;
}
int max2(int a, int b) {
return (a<b) ? b : a;
}
typedef struct _Arbre {
int numero_de_la_sequence;
struct _Arbre * droite;
double distance_gauche;
struct _Arbre * gauche;
double distance_droite;
} Arbre ;
Arbre * new_arbre_nomme (int nom) {
Arbre * a = (Arbre *) malloc(sizeof(Arbre));
a->gauche = (Arbre *) NULL;
a->droite = (Arbre *) NULL;
a->numero_de_la_sequence = nom;
return a;
}
Arbre * fusionne_arbres(Arbre * a, Arbre * b) {
Arbre * c = (Arbre *) malloc(sizeof(Arbre));
c->gauche = a;
c->droite = b;
return c;
}
void affiche_arbre(Arbre * a) {
printf("(");
if (a->gauche != NULL)
affiche_arbre(a->gauche);
if (a->droite != NULL)
affiche_arbre(a->droite);
if (a->droite == NULL && a->gauche == NULL)
printf("%d", a->numero_de_la_sequence);
printf(")");
}
/*
typedef struct _ListeArbres {
Arbre * valeur;
struct _ListeArbres * suivant;
} ListeArbres;
int length_liste_arbres(ListeArbres * l) {
int i=1;
while (l->suivant != NULL) {
l = l->suivant;
i++;
}
return i;
}
*/
int main () {
int i,j;
double scores[NOMBRE_DE_SEQUENCES][NOMBRE_DE_SEQUENCES];
for(i=0; i<NOMBRE_DE_SEQUENCES; i++) {
for(j=0; j<NOMBRE_DE_SEQUENCES; j++) {
scores[i][j] = 0;
}
}
// Calcul de la matrice des scores
for(i=0; i<NOMBRE_DE_SEQUENCES; i++) {
for(j=i+1; j<NOMBRE_DE_SEQUENCES; j++) {
scores[i][j] = score(sequences[i], sequences[j]);
printf("%d,%d: %f\n", i, j, score(sequences[i], sequences[j]));
}
}
/*
double scores [NOMBRE_DE_SEQUENCES][NOMBRE_DE_SEQUENCES] = {
{-0,-22.2,-28,-71.1,-69.4,-70.9},
{-22.2,-0,-19.5,-77.2,-75,-72.2},
{-28,-19.5,-0,-75.7,-74.4,-70.7},
{-71.1,-77.2,-75.7,-0,-11,-22.7},
{-69.4,-75,-74.4,-11,-0,-20.7},
{-70.9,-72.2,-70.7,-22.7,-20.7,-0}
};
*/
// Affichage de la matrice des scores
for(i=0; i<NOMBRE_DE_SEQUENCES; i++) {
for(j=0; j<NOMBRE_DE_SEQUENCES; j++) {
printf("%d, %d: %-3.3f\n", i,j,scores[i][j]);
}
printf("\n");
}
// Construction de l'arbre
// On part d'une forêt, i.e., une liste d'arbres
/*
ListeArbres * foret = (ListeArbres *) malloc(sizeof(ListeArbres));
ListeArbres * a;
a = foret;
for (i=0; i<NOMBRE_DE_SEQUENCES; i++) {
a->valeur = new_arbre_nomme(i);
if (i<NOMBRE_DE_SEQUENCES-1) {
a->suivant = (ListeArbres *) malloc(sizeof(ListeArbres));
} else { // Dernière séquence
a->suivant=NULL;
}
}
*/
Arbre * foret[NOMBRE_DE_SEQUENCES];
for (i=0; i<NOMBRE_DE_SEQUENCES; i++) {
foret[i] = new_arbre_nomme(i);
}
int n = NOMBRE_DE_SEQUENCES;
while (n > 1) {
// Affichage de la matrice des scores
printf("Il reste %d séquences\n", n);
for(i=0; i<n; i++) {
for(j=0; j<n; j++) {
printf("%-3.0f ", scores[i][j]);
}
printf("\n");
}
// On cherche le maximum dans la matrice.
int imax=0;
int jmax=1;
double m = scores[imax][jmax];
for (i=0; i<n; i++) {
for (j=i+1; j<n; j++) {
if( scores[i][j] > m ) {
m = scores[i][j];
imax = i;
jmax = j;
}
}
}
// On fusionne les arbres imax et jmax.
i = min2(imax,jmax);
j = max2(imax,jmax);
printf("On fusionne les arbres %d et %d\n", i,j);
foret[i] = fusionne_arbres(foret[i], foret[j]);
for (i=j+1; i<n; i++) {
foret[i-1] = foret[i];
}
// On fusionne les colonnes imax et jmax
double nouveaux_scores[NOMBRE_DE_SEQUENCES][NOMBRE_DE_SEQUENCES];
int k,l;
// On commence par copier les parties du tableau qui ne changent pas
// Ici, k et l représentent des indices du tableau scores
// (et pas nouveaux_scores)
// 1 (première ligne de la matrice en blocs)
for(k=0; k<i; k++)
for(l=k+1; l<j; l++)
nouveaux_scores[k][l] = scores[k][l];
// 2
for(k=0; k<i; k++)
for(l=i+1; l<j; l++)
nouveaux_scores[k][l] = scores[k][l];
// 3
for(k=0; k<i; k++)
for(l=j+1; l<n; l++)
nouveaux_scores[k][l-1] = scores[k][l];
// 4 (deuxième ligne de la matrice en blocs, au milieu du tableau)
for(k=i+1; k<j; k++)
for(l=k+1; l<j; l++)
nouveaux_scores[k][l] = scores[k][l];
// 5
for(k=i+1; k<j; k++)
for(l=j+1;j<n;j++)
nouveaux_scores[k][l-1] = scores[k][l];
// 6 (troisième et dernière ligne de la matrice par blocs)
for(k=j+1;j<n;j++)
for(l=k+1;l<n;l++)
nouveaux_scores[k-1][l-1] = scores[k][l];
// On met les moyennes des colonnes imax et jmax dans la colonne d'indice i
// Maintenant, k et l sont des indices de nouveaux_scores
// 1 (colonne i)
for(k=0;k<i;k++)
nouveaux_scores[k][i] = (scores[k][i] + scores[k][j])/2;
// 2 (ligne i, entre les colonnes i et j)
for(l=i+1;l<j;l++)
nouveaux_scores[i][l] = (scores[i][l] + scores[l][j])/2;
// 3 (ligne i, entre les colonnes j et n)
for(l=j;l<n-1;l++)
nouveaux_scores[i][l] = (scores[i][l+1] + scores[j][l+1])/2;
n--;
// Le nouveau tableau est complet : on le copie
for (k=0;k<n;k++)
for (l=k+1;l<n;l++)
scores[k][l] = nouveaux_scores[k][l];
}
affiche_arbre(foret[0]);
printf("\n");
return EXIT_SUCCESS;
}
Considérons un instant le problème à moitié résolu : on connait un arbre décrivant l'évolution d'un ensemble d'espèces (youpi !), on connait la distance entre ces espèces, et on veut assigner des longueurs à chacune des branches de l'arbre.
En général, il ne sera pas possible de donner des longueurs aux arrêtes qui redonnent exactement les distances entre les espèces.
On peut néanmoins trouver les longueurs "les plus conformes" à ces distances, au sens des moindres carrés : on cherche à minimiser
Somme ( w_ij * (D_ij - d_ij)^2 ) i,j
où D_ij est la distance observée entre les individus i et j et d_ij est la distance prédite par l'arbre.
Pour les poids, on peut prendre
w_ij = 1 (Cavalli-Sforza and Edwards)
ou
w_ij = D_ij (Margoliash)
Cela peut aussi s'écrire :
D_ij = Somme( x_ijk d_k ) + bruit_ij
k
où x_ijk = 1 si l'arrête k est sur le chemin entre les noeuds i et j
0 sinon
d_k = longueur de l'arrête kou, en termes matriciels :
D = X d + bruit
Si on connait D et X, on cherche d qui minimise la norme L2 du bruit -- c'est simplement un problème de régression.
Le problème, bien sûr, c'est qu'il faut connaître l'arbre, i.e., il faut connaître X. On cherche bien sûr l'arbre qui minimise la norme L2 du bruit ; énumérer tous les arbres serait trop long : nous aurons recours à des algorithmes heuristiques.
Supposons tout d'abord qu'il n'y a que trois séquences, A, B et C. On cherche a, b, c tels que
d(A,B) = a + b d(A,C) = a + c d(B,C) = b + c
On peut alors relier A, B et C à un ancètre commun D avec
a = d(A,D) b = d(B,D) c = d(C,D)
(Si on voulair que l'hypothèse de l'horloge moléculaire soit vérifiée, on imposera a=b=c, i.e., on remplacera ces trois nombres par leur moyenne.)
Considérons maintenant le cas général.
1. Prendre les deux séquences les plus proches, disons A et B. 2. Appliquer l'algorithme pour trois séquences à ces deux séquences et à la moyenne des autres séquences -- comme les distances sont supposées additives, la distance entre A et la moyenne des autres séquences est la moyenne des distances entre A et les autres séquences. 3. Remplacer A et B par le noeud obtenu à l'étape précédente. 4. Recommencer tant qu'il y a au moins trois points.
Graphiquement, on part d'une étoile, toutes les séquences sont reliées à un noeud central, on sélectionne les deux séquences les plus proches et on les détache de ce noeud central ; ensuite, on itère.
A FAIRE : DESSIN (En fait, ce dessin est général)
L'algorithme de Fitch-Margoliash souffre d'un gros défaut : on regroupe les deux taxons les plus proches, comme si on savait que ces deux taxons étaient voisins dans l'arbre. En fait, rien ne dit qu'il sont voisins !
Ainsi, si l'arbre réel était le suivant, il commencerait par regroupper C et D (au lieu de regrouper A et C ou B et D).
A
\ B
\ /
\ /
\ /
\ /
\ /
---
/ \
C DIl nous faut donc une autre méthode pour trouver deux taxons qui sont voisins dans l'arbre...
Le principe du Neighbor Joining est très semblable à l'algorithme de Fitch-Margoliash, c'est juste le choix des séquences à regrouper qui diffère : pour chaque paire de séquences, on regarde la somme des longueurs des branches de l'arbre obtenu si on les fusionne, et on choisit de fusionner la paire conduisant à la valeur la plus faible.
En termes techniques, on utilise le fait que deux taxons i et j qui minimisent
d(i,j) - moyenne( d(i,k), où k != i, j ) - moyenne( d(j,k), où k != i, j )
sont voisins.
En des termes encore différents, on cherche un arbre de longueur minimale. En fait, c'est général : minimiser la longueur de l'arbre revient à minimiser la somme des carrés des différences entre les distances observées et celles prédites par l'arbre (souvenez-vous des histoires de moindres carrés, un peu plus haut) -- c'est un théorème : je ne sais pas si on peut comprendre, intuitivement, pourquoi c'est effectivement équivalent. On parle parfois de "Minimum Evolution Method".
Algorithme de Sattath et Tversky : on considère cette fois-ci tous les ensembles de quatre sommets, disons A, B, C, D et...
BIONJ : une variante du NJ un peu moins sensible au bruit.
A FAIRE The Kishino-Hasegawa test to compare trees has been in use for a dozen years. A similar but better test is the recent Shimodaira-Hasegawa test
Phylip :
dnadist, protdist : calcul de matrices de distance
fitch : pas d'horloge moléculaire
kitsch : horloge moléculaire
neighbor : UPGMA (horloge moléculaire, arbre enraciné) ou NJ (pas
d'horloge, arbre non enraciné)Quelques exercices :
http://www.genopole-lille.fr/fr/formation/cib/doc/phylo/tp.html
Certains algorithmes partent d'une étoile (graphe avec un seul noeud interne, contenant tous les taxons), qu'ils séparent progressivement jusqu'à obtenir un arbre binaire. D'autres algorithmes ont une approche inverse : ils ajoutent les taxons un par un (on commence par les deux premières séquences, qui forment un arbre à deux noeuds, puis on ajoute la troisième, au meilleur endroit possible, etc.).
Cet algorithme glouton est souvent trop rudimentaire : on peut l'améliorer en tentant de réarranger l'arbre à chaque étape (uniquement des réarrangements "locaux" (par exemple, l'analogue des mouvements de fusion que j'utilise pour construire des chemins dans les espaces de modules de courbes dans ma thèse)).
On peut encore l'améliorer en essayant tous les bouturages possibles : pour chaque sous-arbre, on le coupe et on essaye de le greffer partout ailleurs.
Mais tout cela dépend de l'ordre des séquences : on peut changer cet ordre et recommencer, plusieurs fois (au moins dix), pour enfin choisir l'arbre le meilleur.
Avec des données peu nombreuses, on peut aussi énumérer tous les arbres : on les construit un par un, en ajoutant les feuilles une par une. Mais dès que le score devient trop mauvais (pire que le meilleur score obtenu jusqu'à présent), on arrête et on passe aux arbres suivants.
Ainsi, on n'a pas à examiner tous les arbres : tous ceux qui commençaient de la même manière que celui qu'on vient de laisser de côté sont aussi éliminés.
On peut construire un arbre progressivement, en ajoutant les branches une à une. Lors de chaque ajout, on doit choisir où ajouter la branche -- le nombre de choix peut être assez grand.
Une fois un arbre construit, on peut se promener dans l'espace des arbres ayant les mêmes branches à l'aide de mouvements de "fusion" sur les branches internes -- on parle de NNI (Nearest Neighbor Interchange).
Une autre manière de se promener dans cet espace consiste à prendre un sous-arbre, à le détacher et à le greffer ailleur.
Une variante consiste à retirer une arrête de l'arbre (on se retrouve avec deux composantes connexes) puis a rejoindre deux arrêtes.
On ajoute les espèces les unes après les autres (dans un ordre abritraire, mais le résultat dépendra de l'ordre...) à l'endroit produisant l'arbre le plus parcimonieux.
Les méthodes utilisées ne sont pas toujours très stables : un très faible changement dans les données de départ peut changer l'arbre obtenu. Le bootstrap consiste à simuler ces faibles changements et à regarder ce qui change sur l'arbre obtenu. On fait ensuite une sorte de "moyenne" ou de "consensus" sur ces arbres : les noeuds qui apparaissent stables resteront, les autres seront fusionnés ; l'arbre obtenu n'est donc plus binaire.
Le bootstrap consiste à ré colonnes, recalculer l'arbre, répéter au moins 1000 fois, regarder quels groupes apparaissent dans plus de 95% des cas, tenter de construire un arbre (qui ne sera plus binaire) à l'aide de ces groupes.
Il y a des variantes. Par exemple : enlever la moitié des colonnes, recalculer l'arbre, répéter, construire un arbre avec les groupes apparaissant dans plus de 95% des cas.
Voici une autre utilisation du bootstrap. Afin de regarder s'il y a une structure dans les données qu'on considère, on peut prendre la matrice des données et, pour chaque colonne, mettre les données dans le desordre. On obtient alors des données aléatoires et on peut comparer (par exemple) le nombre d'étapes de l'algorithme sur ces données aléatoires avec le nombre d'étapes sur les données initiales. En répétant cela plusieurs fois, on a une idée de la distribution de ce nombre d'étapes : on en déduit si le nombre observé avec les données initiales est significativement différent (on peut alors affirmer que nos données ne sont pas aléatoires et que rechercher une phylogénie a une sens) ou pas (on ne peut rien conclure).
Côté logithèque, on peut utiliser les programmes seqboot (pour faire les tirages avec remise) et consense (pour calculer l'arbre consensus) de Phylip (entre les deux, il convient d'appliquer la méthode de construction de l'arbre qu'on a choisie, par exemple NJ avec "dnadist" suivi de "neighbor" -- on n'oubliera pas d'utiliser l'option "M" pour dire qu'il y a plusieurs jeux de données).
Voici un exemple :
A FAIRE
Dans la discussion précédente, j'ai passé sous silence un "gros" détail : une fois qu'on a un ensemble d'arbres, i.e., une forêt, comment en déduit-on un "arbre consensus" ?
Commençons par reconstruire un arbre non enraciné. Une branche d'un tel arbre, c'est une partition de l'ensemble des feuilles en deux ensembles disjoints. Il suffit donc de construire un arbre à l'aide des branches les plus fréquentes.
A FAIRE : détailler l'algorithme et regarder ce que ça donne avec : AE | BCDF 3 ACE | BDF 3 ACEF | BD 1 AC | BDEF 1 AEF | BCD 1 ADEF | BC 2 ABDF | EC 1 ABCE | DF 3 A FAIRE (Les documents qu'on trouve sur internet sont remarquablement silencieux/confus sur ce point...)
Mettre de l'ordre dans ce qui suit.
Le problème est le suivant : on a choisi un arbre (par exemple, construit de manière approché, par une méthode de distances ; s'il y a peu de séquences, on peut carrément les examiner tous) dont les feuilles sont des séquences et on cherche les séquences correspondant à chacun des noeuds. (C'est par exemple ce dont on a besoin pour calculer les matrices PAM.)
C'est de la programmation dynamique. Pour tout noeud interne v, on note C(v) le score optimal du sous-arbre correspondant.
http:///www.lacim.uqam.ca/%7Echauve/Enseignement/BIF7001/Demonstrations/Phylogenie/ http://www.lifl.fr/%7Evarre/enseignement/dess/bioinfo/phylogenie/tp.html
La phénétique (taxinomie numérique) est la comparaison des phénotypes (i.e, des caractères) des espèces actuelles : on ne se préoccupe pas de la manière dont les espèces ont évolué pour en arriver là. Il est possible que des caractères semblables aient la même origine (on parle d'homologie) ou que la ressemblance soit purement fortuite (homoplasie).
Le cladisme s'efforce au contraire de regarder comment les espèces ont évolué et tente de se limiter homologie et d'éliminer les homoplasies. Pour décider si des espèces sont proches, on ne considère que des caractères évolués (des caractères moins évolués, comme l'épaisseur de la peau ou la présence d'ailes peuvent être dues à des convergences des espèces sans ancètre commun -- pour des caractères évolués, c'est moins probable). Il reste bien sûr à définir ce qu'on entend par "caractère évolué" (caractère que l'espèce a mis longtemps à acquérir).
(Je serais tenté de voir l'opposition phénétique/cladisme comme un analogue de l'opposition synchronie/diachronie en linguistique.)
Considérons un exemple simple : on tire à pile ou face une dizaine de fois et on obtient 7 piles et 3 faces ; on sait que la pièce est truquée et on cherche quelle est la probabilité p d'obtenir pile.
Le principe du maximum de vraissemblance consiste à calculer la probabilité d'obtenir les résultats effectivement obtenus pour chaque valeur du paramètre à estimer (ici, c'est la probabilité p -- plus loin, ce sera la topologie d'un arbre). Cette probabilité s'appelle "vraissemblance". Ici, il vient :
V(p) = p^7 * (1-p)^3.
L'estimateur du paramètre p au sens du maximum de vraissemblance est la valeur du paramètre qui maximise la vraissemblance. En d'autres termes, on cherche la valeur qui maimise la probabilité d'observer les résultats effectivement observés.
Ici :
p maximise V(p) <=> p maximise log V(p)
<=> p maximise 7 log p + 3 log(1-p)
=> 7/p - 3/(1-p) = 0
<=> 7(1-p) = 3p
<=> p = 7/10Le résultat n'est guère surprenant (on aurait obtenu le même résultat sans réfléchir ni faire le moindre calcul), mais l'intérêt, c'est que la méthode est généralisable : si le paramètre est plus complexe (par exemple, la topologie d'un arbre) ou si le phénomène est moins trivial (par exemple, une chaine de Markov), on peut toujours raisonner ainsi.
Remarque : Les estimateurs du maximum de vraissemblance présentent un autre avantage : ils sont asymptotiquement non biaisé (en d'autres termes : sur des données suffisemment volumineuses, ils donnent en moyenne le bon résultat -- sur des données trop peu nombreuses, ça peut être faux...).
Le maximum de vraissemblance présente un autre avantage : il permet d'effectuer des tests. Ainsi, dans notre exemple de calcul de probabilité d'obtention de "pile" avec la pièce truquée, nous pouvons regarder si la probabilité est significativement différente de 1/2. Il suffit pour cela de considérer deux modèles, le premier avec p quelconque, le second avec p=1/2, de calculer leur vraissemblance, puis de faire le quotient. Ce quotient suit (asymptotiquement) une loi du Chi^2 sont le nombre de degrés de liberté est la différence entre le nombre de paramètres des modèles (ici : un modèle à 1 paramètre (p), un modèle à 0 paramètre, différence = 1).
V(p) = p^7 * (1-p)^3
= 0.002223566
V(1/2) = 0.0009765625
V(p) 0.002223566
LR = -------- = -------------- = 2.276932
V(1/2) 0.0009765625
p-valeur = 0.13Pour que ce test ait un sens, il faut que l'un des modèles soit inclus dans l'autre.
(Si on a des doutes sur la pertinence de l'approximation du Chi2, on peut faire des simulations -- ça s'appelle un bootstrap paramétrique.)
C'est assez complexe, je ne suis pas persuadé d'avoir tout compris et je ne me sens pas capable de l'expliquer simplement, sans vous assommer de formule : je ne détaille donc pas.
http://evolution.gs.washington.edu/gs541/2003/lecture33.pdf
La méthode du maximum de vraissemblance permet de faire des tests sur les arbres obtenus : voici quelques idées de tels tests.
Exemple : on peut calculer la vraissemblance avec et sans l'hypothèse de l'horloge moléculaire, faire le quotient des deux et comparer avec une variable du Chi2 avec autant de degrés de liberté que d'équations imposées par l'hypothèse. Cela permet de savoir si cette hypothèse est vérifiée ou non.
Exemple : on veut comparer deux arbres. On n'utilise pas le quotient des vraissemblances, car on ne sait pas du tout comment il évolue (on n'a pas deux modèles imbriqués l'un dans l'autre, mais simplement deux modèles, très semblables, avec le même nombre de paramètres). Par contre, on peut calculer le quotient des vraissemblances en chaque site. On peut ensuite faire un test du signe ou un test de Wilcoxon (les tests paramétriques, Z ou T, sont aussi utilisés -- à tors). On peut aussi estimer la p-valeur à l'aide de simulations (RELL test).
A FAIRE : le test du signe ? google: sign test google: winning-sites test
On prendra garde à ne pas accorder une confiance indue aux résultats : même si différentes méthodes donnent le même résultat, ce résultat qu'on obtient n'est pas nécessairement significativement différent d'autres résultats qu'on n'obtient pas. En termes plus simples : ces programmes ne donnent pas d'"intervalle de confiance"...
A FAIRE : il manque l'attraction des longues branches.
Motivations Tirages au hasard Calculs d'intégrales (espérance : moyenne, variance, etc.) Méthodes bayesiennes Exemples d'espaces Algorithme de Métropolis, implémentation Généralisations et variantes Correction de Hastings Echantilloneur de Gibbs Métropolis-Hastings-Green Applications Phylogénie Alignement multiple
Pour en savoir plus :
http://students.washington.edu/~fkrogsta/bayes/stat538.pdf
ou demander "MCMC metropolis hastings" à Google.
On cherche à tirer au hasard un élément parmi la population (x1,x2,...,xn) selon la distribution (p1,p2,...,pn). Une manière compliquée de procéder consiste à construire une chaine de Markov dont c'est la distribution limite, i.e., construire une suite d'éléments (aléatoires) u1,u2,...,uk,... telle que pour k très grand, uk ait la bonne distribution.
Pour des exemples simplistes, ça semble _vraiment_ tordu, mais on se retrouve parfois face à des espaces (ou "univers") finis de taille gigantesque : par exemple, lors d'une étude phylogénétique, on peut vouloir faire des tirages dans l'"ensemble des arbres à 1000 feuilles".
Si on voulait un tirage selon une distribution équiprobable, on pourrait s'en tirer plus facilement : mais, afin de faire des simulations, on aimerait que les arbres qui décrivent assez bien l'évolution des espèces étudiées aient une probabilité plus importante.
Mais il y a encore pire : en général, on ne sait pas calculer cette probabilité. Pour chaque arbre a, on sait calculer un score S(a) et la probabilité d'un arbre a0 est
score(a0)
P(a0) = ----------------
Somme score(a)
aOr, la somme du dénominateur est beaucoup trop grande pour être calculée.
Les chaines de Markov vont nous permettre de nous promener dans l'espace des arbres, en passant d'un arbre à un arbre "voisin" (pour une notion de voisinage qui reste à définir). Les échantillons qu'on obtiendra ainsi ne seront bien évidemment pas indépendants, mais ils suffiront néanmoins pour effectuer quelques simulations et calculer certaines quantités associées à notre espace d'arbres.
A FAIRE : reprendre le paragraphe suivant -- c'est bien, mais hors sujet. Il faudra le mettre ailleurs.
On rencontre ces mesures de probabilité non équiprobables quand on applique des méthodes bayesiennes : quand on cherche à calculer la valeur d'un paramètre d'un modèle qui décrive le mieux les données observées, on veut souvent un seul résultat. Mais ça n'est pas très naturel, car ce résultat sera, par nature, imprécis : on le complète donc par un intervalle de confiance ou une variance. Mais cette variance n'est directement interprétable que si la distribution du paramètre estimé est normale -- ce qui n'est probablement pas le cas. Une manière de préciser le résultat consiste à donner, non pas une valeur du paramètre, mais toute sa distribution : si on veut une seule valeur, on prendra le mode (MLE) ou la moyenne de cette distribution ; si on veut des intervalles de confiance, on pourra facilement les calculer, quel que soit le risque. Cette distribution de probabilité est dite "a posteriori".
Si l'univers est complexe (par exemple, un ensemble d'arbres ou de graphes), on peut difficilement décrire cette distribution de probabilité. Mais ça n'est pas grave, car on n'a pas besoin de la décrire : si on peut faire des tirages aléatoires selon cette distribution, on pourra en déduire ce qui nous intéresse (le mode, la moyenne, des "intervalles" de confiance, etc.). Les méthodes que nous allons présenter (simulation de Monte Carlo à l'aide de chaines de Markov, MCMC) permettent justement d'effectuer ces tirages.
La méthode de Monte-Carlo permet de calculer une intégrale
Intégrale( f(x), x \in I )
de la manière suivante : prendre x1, x2, ..., xn au hasard dans I et calculer
1 --- Somme f(X_k) N k=1..N
Si l'ensemble d'intégration I est un intervalle fermé de la droite réelle, on prendra les x1 au hasard suivant une loi uniforme. Par contre, s'il s'agit de la droite réelle tout entière, ça n'est plus possible (il n'y a pas de loi de probabilité uniforme sur R tout entier) ; on essaiera d'écrire la fonction à intégrer sous la forme f(x) = g(x) p(x) où p est la densité d'une mesure de probabilité sur R. Il suffira alors de prendre les xi au hasard suivant la loi de p et de calculer
1 --- Somme g(X_k) N k=1..N
Si on compare avec les méthodes d'intégration numérique auxquelles on est habitué (rectangles, trapèzes, etc.), ça n'a pas l'air très efficace -- en particulier en dimension un, i.e., sur la droite réelle. C'est en dimension élevée que l'intégration de Monte-Carlo devient plus intéressante. Imaginons que l'on veuille calculer l'intégrale d'une fonction f de 1000 variables dans l'hypercube [0,1]^1000. La généralisation de la méthode des rectangles suggère d'évaluer la fonction f en chacun des coins de cet hypercube et de faire la moyenne. Le problème, c'est que cet hypercube a 2^1000 coins. C'est trop. L'intégration de Monte-Carlo se contentera d'évaluer f en quelques centaines de points au hasard dans cet hypercube et de faire la moyenne : c'est un peu moins précis, mais on a un résultat avant que le Soleil ne s'éteigne...
Cet exemple peut sembler tiré par les cheveux : rencontre-t-on vraiment des fonction de 1000 variables dans la vie de tous les jours ? Est-on réellement amené à calculer des intégrales sur de tels espaces dans des problèmes concrets ?
Eh bien oui. Face à un arbre phylogénétique retraçant l'évolution d'un millier d'espèces, on peut chercher le nombre de mutation sur chacune des branches, i.e., la longueur de ces branches. L'ensemble de ces longueurs (il y en a quelques milliers) est un point dans un espace de dimension élevée. En fait, l'espace est bien pire que ça, car on voudra considérer tous les arbres...
(Il y a des exemples dans d'autres domaines : par exemple, la valuation des "path-dependent options", en finance.)
Mais veut-on réellement calculer des intégrales sur de tels espaces ? Eh bien oui : moyenne (esp s'écrire sous forme d'intégrale. Par exemple, la probabilité que l'on soit dans un certain ensemble H0 d'hypothèses, sachant qu'on a observé O s'écrit
P(H \in H0) = intégrale( P(H|O), H \in H0 )
P(O|H) P(H)
= intégrale( ------------- , H \in H0 )
P(O)Autre exemple, la probabilité d'observer O' sachant qu'on a déjà observé O s'écrit
P(O'|O) = Intégrale( P(O'|H) * P(H|O), H )
Voici un algorithme un peu trop simpliste (c'est juste l'intégration de Monte-Carlo) :
Pour i = 1 à N: Choisir Xi au hasard suivant p
Le problème, c'est que p peut être arbitraire et, en particulier, difficile à calculer. Dans les exemples qui vont nous intéresser, l'univers sera fini mais très grand et les probabilités seront de la forme
w(a)
P(X=a) = ------------------
Somme( w(b), b )On peut calculer le numérateur, mais le dénominateur, somme sur tous les élements de l'univers, ne se calcule pas (la somme est trop grande).
Au lieu de cela, on peut essayer :
Choisir X0 au hasard. Pour i = 1 à N: Choisir Xi parmi les voisins de X(i-1)
Cela revient en fait à mettre une chaine de Markov sur notre espace, qui permet de passer d'un élément x(i) à un élément voisin x(i-1). Cette manière de choisir un point voisin s'appelle un noyau de transition.
Il y a un "léger" problème, c'est qu'on n'utilise plus p : réintroduisons-le.
Choisir X0 au hasard. Pour i = 1 à N: Choisir X' parmi les voisins de X(i-1) Poser Xi = X(i-1) ou Xi = X' (choisir selon p).
Mais comment fait-on pour "choisir selon p" ? D'après la forme de p mentionnée plus haut, on ne sait pas calculer la probabilité d'un élément, mais par contre, on sait calculer le quotient de deux probabilités. On peut donc procéder ainsi.
# Algorithme de Métropolis
Choisir X0 au hasard.
Pour i = 1 à N:
Choisir X' parmi les voisins de X(i-1)
Choisir un nombre u au hasard dans [0,1]
P(X')
Si u < -----------
P(X(i-1))
Alors
X(i) = X'
Sinon
X(i) = X(i-1)Cet algorithme convient (i.e., la chaine de Markov a bien la distribution limite souhaitée -- modulo quelques hypothèses sur le noyau de transition) si la chaine de Markov est symétrique, i.e., si la probabilité de passer de X à Y est la même que celle de passer de Y à X. Si ce n'est pas le cas, il convient de corriger l'algorithme.
# Algorithme de Métropolis-Hastings
Choisir X0 au hasard.
Pour i = 1 à N:
Choisir X' parmi les voisins de X(i-1)
Choisir un nombre u au hasard dans [0,1]
P( X(i-1) | X' ) P(X')
Si u < ------------------ -----------
P( X' | X(i-1) ) P(X(i-1))
Alors
X(i) = X'
Sinon
X(i) = X(i-1)On peut aussi avoir quelques scrupules, car les Xi ne sont pas indépendants -- en fait, ça n'est pas grave, du moins dans les applications qui nous intéressent -- des calculs d'intégrale.
P(X | X') ne dépend que de la distance entre X et X'.
Dans l'algorithme, on peut changer
Choisir X' parmi les voisins de X(i-1)
en
Choisir X' au hasard.
Si la distribution selon laquelle on choisit X' a des queues plus épaisses que la distribution cible, ça marchera très bien -- mais en cas contraire, ça marchera vraiment très mal...
(Exemple : on n'arrivera pas à simuler (de cette manière) une distribution de Cauchy à l'aide d'une distribution uniforme ou gaussienne.)
Dans un espace de dimension n, on peut traiter les coordonnées une à une : à chaque étape, on ne chage qu'une seule des coordonnées.
La littérature affirme souvent que c'est plus simple, plus rudimentaire que l'algorithme de Metropolis-Hastings : ça me semble plus complexe, et ça suppose qu'on connait somme toute beaucoup de choses sur la distribution.
L'échantilloneur de Gibbs est un cas particulier de la méthode précédente, qui traitait les coondonnées une à une : on impose que le noyau de transition soit la probabilité conditionnelle
P(xi | xj, j!=i)
Ca veut dire qu'il faut connaitre ces probabilités conditionnelles...
Considérons par exemple deux variables X1, X2, normales, de moyenne nulle, de variance 1, de covariance r. On a alors
X1|X2 ~ N(r*X2, 1-r^2)
et de même pour X2. A chaque étape de l'algorithme, on modifie soit X1, soit X2. Le chemin est donc "manhattanien".
my.gibbs <- function (x0=c(0,0), N=20, r=.5, plot.it=T) {
res <- matrix(NA, nr=N, nc=2)
res[1,] <- c(0,0)
for (i in seq(2,N-1,by=2)) {
x1 <- res[i-1,1]
x2 <- res[i-1,2]
res[i,2] <- x2
# Normalement, on devrait prendre un nouveau point n'importe comment
# et l'accepter ou le rejeter -- pour simplifier, j'en prend directement
# un suivant la probabilité marginale.
x1 <- rnorm(1, r*x2, sd=sqrt(1-r^2))
res[i,1] <- x1
res[i+1,1] <- x1
x2 <- rnorm(1, r*x1, sd=sqrt(1-r^2))
res[i+1,2] <- x2
}
if (plot.it) {
plot(res, type='l', xlab="x1", ylab="x2", main="Gibbs sampler")
points(res, pch=16)
invisible(res)
} else {
res
}
}
my.gibbs()
my.gibbs(N=1000)Avec l'algorithme de Métropolis-Hastings, une trajectoire ressemblerait plutôt à
r <- my.gibbs(N=40,plot.it=F)[2*(1:20),]
plot(r, type='l', xlab="x1", ylab="x2",
main="Trajectoire d'une simulation de Métropolis (presque)")
points(r, pch=16)
r <- my.gibbs(N=2000,plot.it=F)[2*(1:1000),]
plot(r, type='l', xlab="x1", ylab="x2",
main="Trajectoire d'une simulation de Métropolis (presque)")
points(r, pch=16)Remarque : On a déjà parlé de l'échantilloneur de Gibbs, à propos de l'alignement local multiple -- voyez-vous le lien ?
C'est une généralisation de Métropolis-Hastings dans laquelle on remplace les densités par des mesures de probabilité : cela généralise toutes les variantes précédentes, y compris l'échantilloneur de Gibbs.
Pour calculer
I = Intégrale( f(x) p(x), x \in I )
quand on ne sait pas prendre des éléments au hasard en suivant la distribution définie par p, on peut essayer d'écrire
p(x)
I = Intégrale( f(x) * ------ * q(x), x \in I )
q(x)où q définit une distribution pour laquelle on peut facilement faire des tirages.
On utilise cette méthode même quand on connait p, car en choisissant bien q, on peut faire baisser la variance (i.e., l'erreur) de I.
Le choix du noyau de transition (i.e., du "voisinage") est délicat.
D'une part, il faut que le taux d'acceptation soit suffisemment élevé pour ne pas rester coincé en un point. Pour cela, il fait que les "voisins" soient assez "proches".
D'autre part, il faut pouvoir atteindre assez rapidement des états assez éloignés : pour cela, il faut que les "voisins" ne soient pas trop "proches".
Pour voir si ces problèmes apparaissent ou pas, pour choisir ou paramétrer le noyau de transition, il convient de faire des dessins (en dimension un : série temporelle et ACF (ou PACF)). Un coup d'oeil à l'ACF permet aussi de construire des variables indépendantes, en prenant X(n), X(n+k), X(n+2k), X(n+3k), etc., avec k suffisemment grand pour que l'autocorrélation cor(X(n),X(n+k)) soit suffisemment faible.
Il convient aussi de faire plusieurs simulations, avec des points de départ différents, afin de voir si ça converge ou pas, et surtout au bout de combien de temps ça converge. On peut aussi faire une simulation beaucoup plus longue, afin de vérifier que ça a bien convergé.
N <- 200
K <- 5
res <- matrix(NA, nr=N, nc=K)
for (i in 1:K) {
r <- my.gibbs(x0=runif(2,-1,1), N=N, plot.it=F)
res[,i] <- apply(r,1,cumsum)[1,] / 1:N
}
matplot(res, type='l', lty=1,
ylab="Estimateur", xlab="Longueur de la simulation",
main="Gibbs sampler")Généralement, on retire le début de la série temporelle, jusqu'à ce que la chaine de Markov ait atteint sa distribution stationnaire (en anglais, on appelle "burn-in" cette phase de transition).
En R :
library(ts)
MH <- function (N = 1000,
voisin = function (x) { rnorm(1,x,1) },
p = dnorm, # Distribution de probabilité à simuler
q = function (x,y) { 1 } # Correction de Hastings
) {
res <- rep(NA,N)
x <- 0
for (i in 1:N) {
y <- voisin(x)
u <- runif(1)
if ( u < q(x,y)/q(y,x) * p(y)/p(x) ) {
x <- y
}
res[i] <- x
}
ts(res)
}
x1 <- ts(rnorm(1000))
x2 <- MH()
x3 <- MH(voisin = function (x) { rnorm(1,x,.01) })
x4 <- MH(voisin = function (x) { rnorm(1,x,50) })
op <- par(mfrow=c(4,1))
plot(x1, main="Hasard")
plot(x2, main="MCMC")
plot(x3, main="MCMC, voisins trop proches")
plot(x4, main="MCMC, voisins trop lointains")
par(op)
op <- par(mfrow=c(2,2))
qqnorm(x1, main="Hasard"); qqline(x1);
qqnorm(x2, main="MCMC"); qqline(x2);
qqnorm(x3, main="MCMC, voisins trop proches"); qqline(x3);
qqnorm(x4, main="MCMC, voisins trop lointains"); qqline(x4)
par(op)
op <- par(mfrow=c(2,2))
hist(x1, col='light blue', probability=T, main="Hasard")
curve(dnorm(x),col='red',lwd=3, add=T)
hist(x2, col='light blue', probability=T, main="MCMC")
curve(dnorm(x),col='red',lwd=3, add=T)
hist(x3, col='light blue', probability=T, main="MCMC, voisins trop proches")
curve(dnorm(x),col='red',lwd=3, add=T)
hist(x4, col='light blue', probability=T, main="MCMC, voisins trop lointains")
curve(dnorm(x),col='red',lwd=3, add=T)
par(op)
op <- par(mfrow=c(4,1))
acf(x1, main="Hasard")
acf(x2, main="MCMC")
acf(x3, main="MCMC, voisins trop proches")
acf(x4, main="MCMC, voisins trop lointains")
par(op)
op <- par(mfrow=c(4,1))
pacf(x1, main="Hasard")
pacf(x2, main="MCMC")
pacf(x3, main="MCMC, voisins trop proches")
pacf(x4, main="MCMC, voisins trop lointains")
par(op)(Pour une chaine de Markov d'ordre 1, c'est normal : il n'y a que la première autocorrélation à regarder)
Exemple d'application à des calculs d'intégrales (on devrait trouver 1 : c'est la variance d'une loi normale standard) :
> mean(x1^2) [1] 0.960904 > mean(x2^2) [1] 1.004407 > mean(x3^2) [1] 0.01575253 > mean(x4^2) [1] 0.9143307
Dans les algorithmes MCMC, on peut modifier la probabilité d'acceptation d'un nouvel état au cours du temps : au début, on est très tolérant et on accepte des état qu'on ne devrait pas accepter et, progressivement, on se rapproche de la probabilité d'acceptation de Métropolis-Hastings.
Cela permet de réduire la durée de la période de transition ("burn-in") et aussi d'éviter certaines pathologies de la chaine de Markov (irréductibilité, etc.).
On peut détourner cette méthode pour un but complètement différent : maximiser (ou minimiser) une certaine quantité, par exemple une log-vraissemblance. L'algorithme (qu'on appelle le "recuit simulé") est le suivant :
1. Prendre une valeur au hasard.
2. Prendre une nouvelle valeur "proche" de l'ancienne
3. Si c'est mieux
Alors
garder la nouvelle valeur
Sinon
garder l'ancienne valeur, avec une certaine probabilité
4. Recommencer en 2.Il y a très, très longtemps, j'ai écrit un mémoire sur quelques méthodes d'optimisation, dont le recuit simulé :
http://zoonek.free.fr/Ecrits/1994_recuit_simulé.ps.gz
Pour estimer un paramètre d'un modèle statistique par la méthode du maximum de vraissemblance, on est amené à calculer une log-vraissemblance, qui s'écrit comme une somme (ou une intégrale, c'est pareil). Si cette somme est trop grande, on peut utiliser une simulation de Monte-Carlo pour l'estimer.
A FAIRE
Souvent, une unique chaine de Markov ne fournit pas de résultats suffisemment fiables : elle peut passer à côté de régions importantes de l'univers (par exemple si cette partie importante de l'univers est séparée de notre point de départ par une zône de moindre probabilité).
Pour y remédier, on peut simplement calculer plusieurs chaines de Markov, en parallèle -- en général, cela produit des résultats plus fiables qu'avec une seule chaine, même plus longue.
On peut aussi généraliser cette idée en permettant aux chaines de s'"échanger" de l'information -- on se rapproche un peu des algorithmes génétiques...
A FAIRE : mettre une référence pour cet article :
A single chain is often not very successful to explore all
possibilities and can miss important areas of the search space
because the walk to these regions might lead through regions that
have very low probabilities. Charles Geyer and Elisabeth Thompson
came up with a scheme that allows a better exploration: use multiple
chains that exchange information about their location and
probability with each other, and that have different capabilities to
move around the space. They coined the term heated chains because
the hot chains accept more often and therefore walk around faster
than the chain of interest (the cold chain). The analogy comes from
physics where gas molecules bounce around much faster when hot than
when cold. To achieve this in the MCMC framework one is powering the
acceptance-rejection formula: the cold chains accepts with r^{1/1} ,
a twice as hot chain with r^{1/2} , a very very hot chain with
r^{1/1000}. if r is bigger than one each chain will always accept,
but if it's smaller than 1 the hot chain will accept more than the
cold chain (for example: r=0.25, the cold chain will accept wit
proability 0.25, the twice as hot chain will accept with 0.5 and the
very hot chain will accept with 0.999). Every step a pair of chains
compares their probabilities and uses and acceptance rejection
scheme, when accepted they swap their place. This allows the cold
chain to reach regions that are difficult to reach for it alone.
On cherche un mot dans un texte en déplaçant une fenêtre le long du texte :
Pour i variant de 1 à longueur(texte)-longueur(mot)+1 :
// On regarde si le mot se trouve en position i du texte
j=1
tant que j <= longueur(mot) et mot[j] = texte[i+i-1]
incrémenter j
si j = longueur(mot)
alors dire qu'on a une occurrence du mot en position i du texteExemple : recherche de "ACT" dans la chaine "ACGTAC". J'ai indiqué par une barre verticale les caractères qui correspondent et par un point les caractères qui ne correspondent pas.
ACGTAC ACGTAC ACGTAC ACGTAC ||. . . . ACT ACT ACT ACT
Dans l'algorithme précédent, on déplace la fenêtre de gauche à droite et, à l'intérieur de celle-ci, on compare les caractères de gauche à droite. Au lieu de cela, on peut déplacer la fenêtre de gauche à droite et, à l'intérieur de celle-ci, comparer les caractères de droite à gauche.
Pour i variant de 1 à longueur(texte)-longueur(mot)+1 :
j=longueur(mot)
tant que j >=1 et mot[j] = texte[i+i-1]
décrémenter j
si j = longueur(mot)
alors dire qu'on a une occurrence du mot en position i du texteExemple :
ACGTAC ACGTAC ACGTAC ACGTAC . .| . . ACT ACT ACT ACT
Dans l'algorithme naïf, si on a déjà comparé avec succès une dizaine de lettres du mot avec le texte quand on s'apperçoit que ça ne marche pas, on peut parfois utiliser cette connaissance que l'on a des lettres suivantes pour décaler la fenêtre de plus d'une case.
Reprenons l'exemple précédent.
ACGTAC ACGTAC ACGTAC ACGTAC ||. . . . ACT ACT ACT ACT
Lors de la première étape, on a regardé les trois premiers caractères ; or on sait que si les deux premiers caractères coïncident et pas le troisième, on peut avancer la fenêtre de deux cases (si on avance d'une seule case, le premier caractère sera C, ce qui ne convient pas). La recherche devient donc :
ACGTAC ACGTAC ACGTAC ||. . . ACT ACT ACT
Pour ce faire, il a suffit de remplir le tableau suivant avant de commencer la recherche.
Position du premier caractère Décalage de la fenêtre qui ne coïncide pas --------------------------------------------------------- 1 1 2 1 3 2
Pour d'autres chaînes de caractères, le tableau est un peu moins trivial. Voici ce qu'on obtient pour ACACGA :
Position du premier caractère qui ne coïncide pas Caractères Décalage ------------------------------------------------ 1 nonA 1 2 A nonC 1 3 A C nonA 3 4 A C A nonC 4 5 A C A C nonG 2 6 A C A C G nonA 6
Pour le remplir, il suffit de chercher, pour tout préfixe M[1..c] du mot à chercher M, le plus long suffixe de M[1..c] qui soit aussi préfixe de M.
Plus long suffixe
strict de M[1..c]
c M[1..c] qui soit pr
--------------------------------------------------------------------------
A?
2 AC vide AC? 2
3 ACA A ACA? 2
4 ACAC AC ACAC? 2
5 ACACG vide ACACG? 5
6 ACACGA A ACACGA? 5Ca n'est pas le même tableau : ici, on ne tient pas compte de l'information supplémentaire que le dernier caractère ne correspondait pas. Mais ça s'adapte facilement. (En fait, on peut aussi tenir compte de la valeur du dernier caractère, celui qui ne correspond pas : on obtient alors un automate fini.)
De plus, il n'est pas nécessaire de regarder à nouveau les premiers caractères de la fenêtre, qui ont déjà été comparés avec succès.
Avant de présenter l'algorithme, nous devons expliquer comment calculer la longueur de ces "suffixes qui sont aussi des préfixes" -- ça s'appelle des "bords". Voici quelques exemples de bords, avec la chaine sur la ligne du milieu, le préfixe au dessus, le suffixe en dessous (on pourrait aussi les encadrer ou les souligner). On remarquera que préfixes et suffixes peuvent se chevaucher.
A AC ACGTA ATA AGAGA ACA ACAC ACGTACGTA ATACCCGATA AGAGAGA A AC ACGTA ATA AGAGA
Exercice :
A FAIRE Calculer, récursivement, la longueur des bords. Ecrire et implémenter l'algorithme.
La manière la plus simple de rechercher un mot dans un texte consiste à construire un automate fini représentant ce mot. (Si vous ne savez pas ce qu'est un automate fini, reportez-vous au livre de P. Wolper, Introduction à la calculabilité.)
Ainsi, pour reconnaître la chaine ACACGA, on peut utiliser l'automate suivant.
A FAIRE : dessin
C'est une variance de l'algorithme de Knuth-Morris-Pratt dans lequel on parcourt la fenêtre de droite à gauche.
C'est une variante de l'algorithme de Knuth-Morris-Pratt dans lequel on cherche plusieurs mots.
http://courses.cs.vt.edu/~algnbio/algnbio_2001/lectures/AhoCorasick.html
Dans les algorithmes précédents, on effectuait des prétraitements sur le mot pour pouvoir ensuite lui soumettre des textes. On peut faire le contraire et faire subir un prétraitement au texte afin de pouvoir effectuer rapidement des recherches sur des mots.
Pour indexer un texte dans une langue humaine, on peut le découper en mot et mettre ces mots dans une table de hachage ou un B-arbre (les INDEX d'une base de donnée) qui nous renverra une liste chainée indiquant la position des différentes occurrences des mots.
En génomique, le problème est un peu différent : on veut indexer un texte, mais il n'y a pas de mots. On veut toutes les sous-chaines : pour cela, il nous faudra une structure de donnée un peu plus efficace -- une énumération de toutes les sous-chaines n'est pas raisonnable : il y en a beaucoup trop et elles sont très longues.
A FAIRE : détailler un peu plus (l'exemple que je prends n'a rien à voir avec la problématique qui nous intéresse : il s'agit d'indexer des mots, pas un texte).
Quand on indexe un texte, on constate que beaucoup de mots commencent de la même manière : au lieu de les stocker dans un tableau (table de hachage), on peut les mettre dans un arbre.
A FAIRE : exemple (dessin) On prendra un mot avec ses "conjugaisons" ex: normal, normale, normales, normaux.
L'exemple suivant montre le trie des mots normal, normale, normales, normaux.
n--o--r--m--a--l
| `--e
| `--s
`--u--xLe problème avec cet arbre, c'est qu'on peut aussi l'écrire
n--o--r--m--a--l--e--s
`--u--xPour bien voir quels sont les mots présents dans le trie, nous allons ajouter un caractère indiquant la fin des mots -- cela permet d'avoir une feuille pour chaque mot.
n--o--r--m--a--l--$
| |
| `--e--$
| |
| `--s--$
|
`--u--x--$ou :
n o r m a l $
*-----*-----*-----*-----*-----*-----*-----*
| |
| | e *
| `-----*-----*
| |
| | s $
| `-----*-----*
|
| u x $
`-----*-----*-----*Un tel arbre codant un ensemble de mots s'appelle un trie (prononcez "traille"). Les arrêtes sont labellées par des lettres.
Pour avoir des arbres plus compacts, on peut enlever les noeuds qui ne donnent pas naissance à plusieurs embranchements. Les arrêtes sont alors labellées par des chaines et non plus des lettres.
normaa l $
*----------*-----*-----*
| |
| | e *
| `-----*-----*
| |
| | s $
| `-----*-----*
|
| u x $
`-----*-----*-----*
Pour indexer tous les mots (toutes les sous-chaines) d'un texte, on peut mettre tous les suffixes du texte dans un trie. Si on prend un trie compact, sa taille n'est plus quadratique mais simplement linéaire.
Pour savoir si un mot est contenu dans le texte, il suffit de voir si le mot rentre dans le trie.
Il suffit simplement d'ajouter les suffixes dans l'arbre, les uns après les autres (par exemple, en commençant par le plus grand).
Voici un exemple d'implémentation :
Exercice laissé au lecteur : proposer une implémentation simpliste, avec un arbre non compressé.
Voici un autre exemple, cette fois-ci avec un arbre compressé.
// Construction de l'arbre des suffixes
// (c) 2003 Vincent Zoonekynd <zoonek@math.jussieu.fr>
// Distributed freely -- use at your own risk
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
// J'ai besoin de la fonction strndup dans ma routine d'affichage...
#define __USE_GNU
#include <string.h>
// Structure de donnée pour manipuler l'arbre des suffixes
#define ALPHABET_SIZE 255
typedef struct _Arbre {
// Le chemin vers ce noeud est indexé par mot[a..b]
char * mot;
int a;
int b;
struct _Arbre * child[ALPHABET_SIZE];
} Arbre;
Arbre * newArbre (char * mot, int a, int b) {
int i;
Arbre * arbre = (Arbre *) malloc(sizeof(Arbre));
arbre->mot = mot;
arbre->a = a;
arbre->b = b;
for (i=0; i<ALPHABET_SIZE; i++) {
arbre->child[i] = NULL;
}
return arbre;
}
void ajouteMot (Arbre * arbre, int k) {
char * mot = arbre->mot;
int n = strlen(mot);
assert(k<n);
printf("Ajout de %s\n", mot + k);
// On regarde si la première lettre du mot est reconnue par l'arbre.
if( arbre->child[ (int) mot[k] ] != NULL ) {
// On descend l'arbre et on regarde le plus grand préfixe du mot à
// insérer qui est d
Arbre * branche = arbre->child[ (int) mot[k] ];
int a = branche->a;
int b = branche->b;
int i=0;
while (a+i<b && k+i<n && mot[a+i] == mot[k+i])
i++;
assert(k+i<n); // Ca voudrait dire que le mot est déjà dans l'arbre
if (a+i == b) {
printf(" Descente\n");
ajouteMot(branche, k+i);
} else {
// Il faut couper le noeud courrant en deux.
Arbre * nouvelle_branche = newArbre(mot, k+i, n);
Arbre * nouveau_noeud = newArbre(mot, k, k+i);
arbre->child[ (int) mot[k] ] = nouveau_noeud;
nouveau_noeud->child[ (int) mot[a+i] ] = branche;
nouveau_noeud->child[ (int) mot[k+i] ] = nouvelle_branche;
branche->a = a+i;
printf(" Insertion d'un nouveau noeud\n");
}
} else {
// On crée une nouvelle branche
printf(" nouvelle branche\n");
arbre->child[ (int) arbre->mot[k] ] = newArbre(arbre->mot,k,n);
}
}
Arbre * arbreDesSuffixes (char * mot) {
Arbre * arbre = newArbre(mot,0,0);
int i;
int n = strlen(mot);
for (i=0; i<n; i++) {
ajouteMot(arbre, i);
afficheArbre(arbre);
}
return arbre;
}
void _afficheArbre(Arbre * arbre, int n) {
int i;
for (i=0; i<n; i++)
printf(" ");
printf("%s\n", strndup(arbre->mot + arbre->a, arbre->b - arbre->a));
for (i=0; i<ALPHABET_SIZE; i++) {
if( arbre->child[i] != NULL )
_afficheArbre( arbre->child[i], n+2 );
}
}
void afficheArbre(Arbre * a) {
_afficheArbre(a, 0);
printf("\n");
}
int main () {
char * chaine = "AAATAT$";
Arbre * arbre = arbreDesSuffixes(chaine);
afficheArbre(arbre);
return EXIT_SUCCESS;
}
Les sous-chaines répétées d'un texte sont les préfixes des noeuds internes (i.e., les noeuds qui ont au moins deux fils).
On en déduit un algorithme de recherche de toutes les sous-chaines répétées avec au plus k erreurs : il suffit de parcourrir l'arbre en profondeur, en gardant en mémoire le nombre d'erreurs commises.
Exercice laissé au lecteur : implémenter cet algorithme
Cet algorithme est efficace en espace (O(n), si n est la taille du texte), mais pas en temps (O(n^2)). Nous allons présenter un algorithme linéaire en temps.
On procède séquentiellement : on commence par construire l'arbre des suffixes du texte constitué par le premier caractère du texte, puis par les deux premiers, les trois premiers, etc. A chaque étape, on doit donc ajouter une lettre à chacune des feuilles (justement, l'un des problèmes, c'est que ce sont pas toujours des feuilles).
Voici par exemple le détail de la construction itérative de l'arbre des suffixes de la chaine CCGTATATACG$. (Faites le calcul à la main, en dessinant les arbres normalement -- je suis très mauvais en ASCII-art.)
Iteration 0 --C Iteration 1 --CC Iteration 2 --C |--CG |--G --G Iteration 3 --C |--GT |--CGT --GT --T Iteration 4 --A --GTA --C |--GTA |--CGTA --TA Iteration 5 --TAT --C |--CGTAT |--GTAT --GTAT --AT Iteration 6 --C |--GTATA |--CGTATA --ATA --TATA --GTATA Iteration 7 --C |--GTATAT |--CGTATAT --TATAT --GTATAT --ATAT Iteration 8 --C |--CGTATATA |--GTATATA --TATATA --ATATA --GTATATA Iteration 9 --A |--C |--TA | |--C | |--TAC --C |--CGTATATAC |--GTATATAC --GTATATAC --TA |--C |--TA | |--C | |--TAC Iteration 10 --A |--CG |--TA | |--CG | |--TACG --C |--GTATATACG |--CGTATATACG --GTATATACG --TA |--CG |--TA | |--CG | |--TACG Iteration 11 --A |--CG$ |--TA | |--TACG$ | |--CG$ --C |--G | |--TATATACG$ | |--$ |--CGTATATACG$ --$ --G |--TATATACG$ |--$ --TA |--CG$ |--TA | |--TACG$ | |--CG$
Pour bien distinguer les différentes étapes de l'algorithme, on peut les numéroter ainsi :
(1,1) : ajout du suffixe à 1 lettre du préfixe à 1 lettre (1,2) : ajout du suffixe à 1 lettre du préfixe à 2 lettres (2,2) : ajout du suffixe à 2 lettres du préfixe à 2 lettres (1,3) : ajout du suffixe à 1 lettre du préfixe à 3 lettres (2,3) : ajout du suffixe à 2 lettres du préfixe à 3 lettres (3,3) : ajout du suffixe à 3 lettres du préfixe à 3 lettres etc.
On en déduit la complexité de l'algorithme : O(n^3) (O(n^2) étapes de complexité O(n)).
L'un des problèmes, c'est qu'on ne sait pas exactement à quel endroit il faut ajouter la lettre suivante ; à chaque étape, on doit recommencer depuis la racine de l'arbre.
On peut le résoudre facilement en gardant en mémoire la liste des feuilles ajoutées à l'itération précédente. On peut représenter cette liste graphiquement comme des liens qui sautent d'une branche de l'arbre à l'autre.
A FAIRE : un dessin.
Ces flèches (entre un noeud représentant un mot de longueur n et le noeud représentant son suffixe de longueur n-1) s'appellent des "liens suffixes". L'algorithme que nous venons de présenter n'est pas trés intéressant (sa version naïve est beaucoup plus lente que l'algorithme initial, la version avec les liens suffixes a la même complexité), mais la notion de lien suffixe l'est : nous venons de voir qu'elle permettait d'accélérer l'algorithme. Elle jouera une place centrale dans l'algorithme d'Ukkonen -- il suffira "juste" de prendre des tries compacts et de regarder les modifications effectuées à chaque étape de l'algorithme, pour repérer celles qui sont inutiles.
Voici une implémentation possible :
#!/usr/bin/python
# (c) 2003 Vincent Zoonekynd <zoonek@math.jussieu.fr>
# Distributed freely -- use at your own risk
# Algorithme d'Ukkonen quadratique
# Il est quadratique pour deux raisons :
# - d'autre part, je n'utilise pas le "court-circuit" qui permet de le
# rendre lineaire (il est donc quadratique en temps)
# - d'une part, je construit un arbre non compresse (je le compresse apres)
# (il est donc quadratique en espace)
# L'algorithme normal, un peu plus complique, est lineaire en temps et
# en espace.
# Choix d'une chaine, au hasard
from random import *
chaine = ""
for i in range(10):
chaine = chaine + sample( ('a','a','a','b','b','c'), 1 )[0]
chaine = chaine + "$"
chaine="CCGTATATACG$"
print "Chaine : " + chaine
# Routine d'affichage
def affiche_arbre (a,n=0) :
for i in a.keys():
print n*" |" + "--" + i
affiche_arbre(a[i], n+1)
# Routine de compression
def compresse_arbre (a) :
done = 0
while not done:
done = 1
for i in a.keys():
if len(a[i].keys()) == 1:
j = a[i].keys()[0]
a[ i + j ] = a[i][j]
del(a[i])
done = 0
for i in a.keys():
compresse_arbre(a[i])
# Routine d'affichage d'un arbre compresse
import copy
def affiche_arbre_compresse (a):
b = copy.deepcopy(a)
compresse_arbre(b)
affiche_arbre(b)
# Construction de l'arbre des suffixes
racine = {}
liens = [racine]
for i in range(len(chaine)):
print "Inserting prefix " + chaine[0:(i+1)]
nouveaux_liens = []
for j in range(i+1):
print " Inserting factor " + j*" " + chaine[j:(i+1)]
if not liens[j].has_key(chaine[i]) :
liens[j][chaine[i]] = {}
nouveaux_liens.append( liens[j][chaine[i]] )
nouveaux_liens.append(racine)
liens = nouveaux_liens
print "\nIteration " + str(i)
affiche_arbre_compresse(racine)
# Compression de cet arbre
compresse_arbre(racine)
# Affichage
print "\nVoici l'arbre de suffixes de la chaine '" + chaine + "'\n"
affiche_arbre(racine)
# Exemple d'application
# Les sous-chaines repetees sont les (suffixes des) noeuds internes
print "\nLa (les) plus grande(s) sous-chaine(s) repetee(s) de '" + chaine + "' sont :"
resultat = []
def recherche (a, c=""):
for i in a.keys():
if( len(a[i].keys()) >0 ):
resultat.append(c+i)
recherche(a[i], c+i)
recherche(racine)
l = max(map(len,resultat))
for i in resultat:
if len(i) == l:
print " " + i
C'est une construction linéaire (en temps comme en espace) de l'arbre des suffixes basé sur la construction précédente. Notre algorithme a un énorme défaut : il est quadratique en temps (il est aussi quadratique en espace, mais c'est facile à corriger : il suffit de construire directement un arbre compact). En effet, si la chaine est de taille n, il y a n itérations et la k-ième itération insère k feuilles.
Lors de la construction de l'arbre que nous venons de présenter, on constate qu'à certaines étapes, on ne fait rien.
Pour bien comprendre, je suggère de faire l'exercice suivant.
Exercice : Construire l'arbre des suffixes de ACCGACCT comme précédemment (arbre des suffixes de A, puis de AC, puis de ACC, puis de ACCG, etc.). A chaque étape, indiquer quelle est l'opération effectuée : Ajouter une feuille, Etirer une feuille, Couper une arrête, Ne rien faire. Que remarque-t-on dans ce tableau ? Peut-on utiliser cette remarque pour accélérer l'algorithme ?
A A = Add a leaf EA E = Extend a leaf EEI S = Split an edge and Add a leaf EESA I = Implicit node / Do nothing EEEEI EEEESI EEEEEII EEEEEASA
(Il doit y avoir une erreur dans ce tableau, mais ça ressemble vaguement à ça : il y a plein de E...)
L'idée qui vient à l'esprit (c'est un morceau de l'algorithme d'Ukkonen), c'est de ne pas faire les E : ils resteront jusqu'à la fin (once a leaf node, always a leaf node). Au lieu de mettre (a,b) sur ces arrêtes, et de remplacer à chaque étape par (a,b+1), (a,b+2), etc., on met directement (a,infini) et on ne s'en préoccupe plus. (Mais attention : on sera quand même amené à couper l'arrête en question de temps à autre.)
Le noeud où on commence le travail (après tous les E) est le noeud actif.
L'autre remarque, c'est que les étapes où on ne fait rien sont juste à la fin. Sur cet exemple, il n'y en a pas beaucoup, mais avec des chaines plus longues, on peut en avoir bien plus. L'idée est donc de ne pas faire ces étapes : dès qu'on en rencontre une, on s'arrête et on passe à l'itération suivante. Le noeud correspondant à cette étape s'appelle le noeud final : c'est le noeud actif de l'itération suivante.
A chaque itération, il suffit donc d'ajouter des feuilles entre le noeud actif et le noeud final (i.e., tant qu'on peut en ajouter).
Dans la présentation naïve de l'algorithme d'Ukkonen, on savait exactement à quels noeuds il fallait ajouter une feuille, car on avait une liste les énumérant. Ici, nous allons retrouver ces noeuds à l'aide de liens suffixes.
Un lien suffixe dans un trie compressé est un lien (on pourrait dire "une arrête ajoutée à l'arbre" -- mais alors, ce n'est plus un arbre mais juste un graphe non orienté) d'un noeud interne vers son plus grand suffixe strict qui est aussi un noeud interne.
Par exemple, dans un trie contenant ACCTGC, pour construire le lien suffixe de ACCTGC, on regarderait si ses suffixes stricts, CCTGC, CTGC, TGC, GC, C ou "vide" sont dans le trie : le lien suffixe est un lien de ACCTGC vers le plus grand de ces suffixes qui est dans le trie. (Dans le pire des cas, il n'y en a aucun, le seul suffixe présent est "vide" et le lien suffixe pointe vers la racine du trie.)
Remarque : dans l'algorithme d'Ukkonen, tous les noeuds internes auront un lien suffixe (mais pas les feuilles, contrairement à ce qui se passait dans notre introduction).
Remarque : Les liens suffixes ne servent en fait que lors de la construction de l'arbre, plus après (mais on doit quand-même pouvoir imaginer des applications qui les utilisent quand-même).
Complexité : A chaque étape, l'arbre grossit, or sa taille est linéaire : il y a bien O(n) étapes. Grace aux liens suffixes, chacune de ces étapes s'effectue en temps constant.
Recherche d'une sous-chaine dans une base de données :
Exemple : on dispose d'un ensemble d'étiquettes (courtes séquences d'ADN) et de plein de gènes, on aimerait les associer : construire les arbres de suffixes des gènes, fusionner ces arbres et chercher les étiquettes dans le nouvel arbre. Exemple : complétion automatique dans un shell
Recherche approchée d'une chaine dans une base
Plus longue sous-chaine commune :
Application : détection de la contamination d'une séquence d'ADN ; on cherche une sous-chaine commune de longueur au moins N entre la séquence d'ADN et le texte obtenu en concaténénant les sources possibles de contamination (vecteur, amorces, génome de l'hôte, etc.)
Plus longue sous-chaine répétée
Plus longue sous-chaine palindromique
Alignement de deux séquences, très longues, en un temps linéaire :
construire l'arbre des suffixes de x ou y, en déduire les sous-chaines communes maximales de x et y (MUM, Maximum Unique Match), les ordonner, en déduire une "suite maximale" de MUMs (ça se fait en O(n log n)), combler les brèches.
Recherche des répétitions en tandem :
O(n)
MGA :
(Multiple Genome Aligner -- c'est une manière de calculer des alignements multiples avec des séquences très longues).
Détection du plagiat :
Recherche des chevauchements entre deux textes, analyse statistique de la pertinence du résultat. Utile pour les enseignants qui voudraient repérer ceux de leurs étudiants qui abusent d'Internet. http://www.plagiarism.org/
Recherche, non pas d'une chaine dans un texte, mais d'une PSSM :
On descend dans l'arbre tant que le score de l'alignement entre la PSMM et le texte ne descend pas en dessous d'un certains seuil. Remarque : Très souvent, les PSMM ont une partie centrale très conservée alors que les extrémités sont plus variables ; on doit pouvoir améliorer l'algorithme en commençant la recherche à partir du milieu de la chaine -- en remplaçant l'arbre des suffixes par l'arbre des affixes.
Compression : l'algorithme de Ziv-Lempel repose sur les arbres de suffixes (d'où l'intérêt de construire l'arbre des suffixe en lisant une chaine de gauche à droite...)
A FAIRE : URL de SuffixTrees.ppt
Application des liens suffixes (dans les discussions précédentes, on n'utilisait les liens suffixes que lors de la construction de l'arbre des suffixes, mais plus après) :
Pour trouver la plus grande sous-chaine d'une chaine donnée qui apparait dans un texte : construire l'arbre de suffixes du texte, avec les liens suffixes ; chercher la chaine dans l'arbre ; quand ça ne marche plus, suivre le lien suffixe.
C'est une structure plus simple, qui demande moins de mémoire (c'est toujours linéaire) mais la recherche d'une chaine dans le texte est un peu plus longue.
But : Quand on cherche une chaine de caractères dans un texte, on cherche en fait des préfixes de plus en plus grands. Par exemple, pour chercher "maison", on commence par chercher "m", puis "ma", puis "mai", etc. On parcourt donc la chaine à chercher de droite à gauche. On aimerait pouvoir la parcourrir dans les deux sens, par exemple en cherchant "s", "so", "iso", "aiso", "aison", "maison".
Idée : Construire l'arbre des suffixes du texte et du texte renversé ; Ajouter des arrêtes entre les sommets de l'un les sommets correspondants de l'autre (par exemple, une arrête entre le noeud "abc" de l'arbre du texte et le noeud "cba" du texte inversé).
Avec les arbres de suffixes non compacts, c'est facile, mais quadratique. On peut aussi faire ça avec les arbres de suffixes compacts, quitte à ajouter les noeuds qui manquent.
Si on compresse ces liens, on obtient un graphe (orienté) dont les arrêtes sont décorées par les lettres de l'alphabet et coloriées (une couleur pour l'arbre des suffixes du texte, une autre pour celui de son inverse).
La construction naïve est en O(n log n), mais il existe un algorithme linéaire.
Application : recherche de motifs palindromiques (p.ex., une épingle à cheveux)
stem = .4 hairpin = @stem CGCG reverse(@stem)
WOTD (Write-Only, Top-Down)
1. Calculer l'ensemble des suffixes
2. Les mettre dans l'ordre selon leur premier caractère
3. Pour chacun des groupes ainsi obtenus :
3a. Si c'est un singleton, créer une feuille
3b. Sinon, trouver le plus grans préfixe commun de ces suffixes,
créer le noeud interne correspondant, reprendre en 2 avec le reste
des suffixes.
Complexité : O(n^2) en pire cas, O(n log n) en moyenne, O(n) en espace.Lazy WOTD :
Idem, mais on ne construit pas l'arbre entier dès le début -- on ne construit au fur et à mesure, uniquement si on en a besoin.
Ils utilisent énormément de mémoire : c'est linéaire, certes, mais si le texte est énorme, c'est quand-même trop...
Leur construction utilise des morceaux de mémoire dans tous les sens : il faut que tout l'arbre soit en mémoire (Si on veut faire un arbre de suffixes avec une base de données de quelques Go...)
Il existe quelques variantes de (la construction de) ces arbres qui tient compte de la réalité des accès mémoire.
SST : Sparse Suffix Tree DST : Distributed Suffix Tree (aka Paged Suffix Tree)
Un vieil article de DDJ (TRÈS CLAIR) http://dogma.net/markn/articles/suffixt/suffixt.htm http://courses.cs.vt.edu/~algnbio/algnbio_2001/handouts/Ukkonen_Trees.ps Description rapide de l'algorithme, implémentation en JavaScript, Applications http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Tree/Suffix/ http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Tree/Suffix/ Recherche dans un texte Recherche de la sous-chaine répétée la plus longue Recherche de la sous-chaine commune la plus longue (soit marquer les noeuds pour savoir s'ils appartiennent à la première ou à la deuxième, soit construire l'arbre des suffixes de mot1$mot2#). Recherche du plus long palindrome (texte$reverse(texte)#) Une implémentation http://bioinformatics.rit.edu/~tex/publications/suffixtree/ Autre exemple d'application des arbres de suffixes : indexer des données qui ne sont pas textuelles. Ici, l'angle entre les AA d'une protéine, pour rechercher des similitudes de structure tertiaire. http://www14.informatik.tu-muenchen.de/personen/taeubig/ProSt/PAST.html Question 12: Describe which (i,j) iterations would be excuted or skipped in the Ukkonnen Algorithm to construct a suffix-tree for the following 8-base pair sequence : ACCGACCT Explain by showing the suffix tree after each i-iteration.
L'algorithme de McCreight est une autre construction linéaire de l'arbre des suffixes, qui utilise aussi des liens suffixes -- je trouve l'algorithme d'Ukkonen plus simple à comprendre et surtout plus utile...
0. On insère les suffixes dans l'arbre un par un, en commençant par les plus grands.
1. On choisit d'utiliser des arbres compacts, qui occuppent un espace linéaire (les labels des noeuds ne sont pas écrits explicitement, on se contente de mettre la position de début et de fin de ce label dans le texte).
2. Lors de la construction de l'arbre, il faut trouver l'endroit où ajouter la branche, puis ajouter la branche. Comme l'arbre est compact, on n'ajoute qu'une seule branche : cela s'effectue en temps constant.
3. Il reste à trouver l'endroit où il faut ajouter la branche. Si on procède naïvement, cela se fait en temps linéaire. On peut néanmoins améliorer les choses en ajoutant des arrêtes à l'arbre (qui ne sera plus un arbre) afin de préparer l'insertion des suffixes suivants. (Ces arrêtes ne resteront pas figées lors du déroulement de l'algorithme : elles pourront être modifi
Trouver l'endroit où greffer la branche suivante, revient à trouver la tête, i.e., le plus grand prefixe du mot à insérer qui est déjà présent dans l'arbre. On écrit cette tête sous la forme "a b" où a est une lettre et b un mot. On sait alors que la tête suivante commencera par b.
L'idée est alors la suivante : pour chaque noeud correspondant à un préfixe de la forme "a b", ajouter un "lien suffixe" vers le noeud "b". On ajoute ces liens suffixes tout en construisant l'arbre et on les utilise lorsqu"on passe de l'insertion d'un suffixe à l'insertion du suffixe suivant.
A chaque étape on ajoute au plus deux noeuds (soit une feuille, soit un noeud et une feuille), il n'y a qu'un nombre constant de liens suffixes à créer ou à modifier.
Une chaine de Markov est un processus (i.e., une suite de variables aléatoires) Xn telle que
P[ X(n+1)=a | X(n)=b ]
ne dépende pas de n (on dit souvent que le processus a une mémoire courte : si on voit X(n+1) comme une décision prise le jour n+1, elle ne dépend que de ce qui s'est passé la veille -- les évènements antérieurs sont complètement occultés).
Imaginons par exemple une araignée se déplaçant sur les arrètes d'un tétraèdre, allant d'un sommet à l'autre. Si on note 1, 2, 3, 4 les sommets du tétraèdre, X(n) la position de l'araignée après n unités de temps, si l'arraignée choisit au hasard le prochain sommet qu'elle va visiter, on a
P[ X(n+1) = a | X(n) = a ] = 0 P[ X(n+1) = a | X(n) = b ] = 1/3, si a <> b,
i.e., la probabilité que l'araignée reste sur le même sommet est nulle et la probabilité qu'elle visite l'un des trois autres sommets est 1/3 (pour chaque sommet).
On représente généralement une telle chaine de Markov par une matrice de transition (les colonnes correspondent à X(n), les lignes à X(n+1) et les entrées à la probabilité de la transition).
0 1/3 1/3 1/3 1/3 0 1/3 1/3 1/3 1/3 0 1/3 1/3 1/3 1/3 0
On a alors
P_n = M^n * P_0
où M est la matrice de transition et P_k le vecteur donnant la probabilité d'être dans un état donné au temps k.
On peut compliquer un peu ce modèle en mettant des mouches en certains sommets (leur odeur va attirer la mouche et l'inciter à rester), ou de la colle (l'araignée restera bloquée sur ce sommet).
0 1/4 1/6 0
1/4 0 1/6 0
(mouche) 1/2 1/2 1/2 0
(colle) 1/4 1/4 1/6 1On peut aussi représenter une chaine de Markov par un graphe orienté : les sommets correspondent aux états (dans l'exemple de l'araignée, ils sont numérotés par le numéro du sommet du tétraèdre) et les arrêtes aux transitions ; sur les arrêtes on indique la probabilité de la transition.
A FAIRE : dessin
Exercice : simuler une trajectoire de l'araignée, soit à l'aide d'un ordinateur, soit à l'aide de dés. (Lors d'un cours, on peut demander à chacun des étudiants de simuler une trajectoire -- il faudra peut-être prévoir des dés avec un nombre de faces peu standard ou changer les probabilités précédentes. On pourra ensuite comparer toutes ces trajectoires : par exemple, quelle est la proportion des trajectoires dans lesquelles l'araignée mange la mouche ? quelle est la durée de vie moyenne de la mouche ?)
Remarque : ça ressemble beaucoup aux automates finis. Un automate probabiliste, c'est comme une chaine de Markov, mais les labels sont sur les transitions, pas sur les états.
On constate que les descendants d'un fumeur ont 60% de "chances" de fumer, alors que les descendants d'un non fumeur ont 20% de chances de fumer. Sachant que la moitié de la population fume, peut-on s'attendre à une disparition des fumeurs ?
On peut modéliser cette situation par une chaine de Markov à deux états, "fumeur" et "non fumeur", dont les transitions correspondent au passage à la génération suivante. Le vecteur de probabilité P_k correspond alors à la proportion de fumeurs et de non-fumeurs dans la k-ième génération.
Exercice : écrire la matrice de transition (en détaillant les calculs, avec les probabilités conditionnelles). A l'aide d'un ordinateur, calculer P_100. Conclure.
Remarque : on peut rapprocher cet exemple de ceux décrivant l'évolution d'une population face à un prédateur à l'aide d'équations différentielles. Cela permet de voir les chaines de Markov comme un analogue discret de certains modèles évolutifs.
A FAIRE : mettre un lien ou un dessin...
On peut décrire un texte anglais par une suite de probabilités : la probabilité qu'un mot suive un autre, i.e., les probabilités du genre
P[ on lit "destruction" sachant qu'on vient de lire "mass" ] = 0.002.
ce qu'on écrirait
P[ X(n+1) = "destruction" | X(n) = "mass" ] = 0.002
Bien évidemment, il y en a beaucoup (si on se limite aux mots les plus courrants, on pourra en prendre 2000, et comme ces probabilités font intervenir deux mots, il y en aurait 4 millions).
Exercice : Ecrire un programme un programme qui calcule ses probabilités à l'aide d'un corpus de textes. Ecrire un autre programme qui effectue une simulation, i.e., qui écrive un texte à l'aide de cette chaine de Markov. Le texte est-il correct ? Modifier les programmes précédents pour qu'ils utilisent une chaine de Markov d'ordre deux, i.e., en utilisant les probabilités
P[ X(n+2) | X(n+1), X(n) ].
Cet exemple est détaillé (je crois) dans le livre "Programming Pearls".
Remarque : en prenant des lettres (ou des couples de lettres, on oubliant complètement les espaces) à la place des mots, on peut facilement reconnaître la langue d'un texte.
Remarque : On peut ainsi construire une chaine de Markov qui reconnait une (ou plusieurs) séquence(s) d'acides aminés.
# Construction d'une matrice de Markov
markov <- function (x) {
x <- strsplit(x,'')[[1]]
x <- factor(x)
aa <- strsplit("ACDEFGHIKLMNPSQRTVWY01",'')[[1]]
n <- length(aa)
m <- matrix(0, nr=n, nc=n)
colnames(m) <- aa
rownames(m) <- aa
m["1","1"] <- 1
m["0", x[1]] <- m["0", x[1]] +1
for (i in 1:(length(x)-1)) {
m[ x[i], x[i+1] ] <- m[ x[i], x[i+1] ] +1
}
m[ x[length(x)], "1" ] <- m[ x[length(x)], "1" ] +1
# C'est une matrice d'effectifs, on la transforme en une matrice de
# probabilités : la somme de chaque ligne doit être égale à 1
m <- m +.001 # Pas n
m <- m/apply(m,1,sum)
print(round(m, digits=2))
invisible(m)
}
x <-"MAKGVAVLNSSEGVTGTIFFTQEGDGVTTVSGTVSGLKPGLHGFHVHALGDTTNGCMSTGPHFNPDGKTHGAPEDANRHAGDLGNITVGDDGTATFTITDCQIPLTGPNSIVGRAVVVHADPDDLGKGGHELSLATGNAGGRVACGIIGLQG"
m <- markov(x)Exercice : utiliser cette chaine de Markov pour construire de nouvelles chaines...
On peut calculer les paramètres d'une chaine de Markov par simple comptage (i.e., calcul de proportion, ou calcul de probabilité), mais cela risque de produire des probabilités nulles : on peut les éviter en utilisant les "pseudocount" de Laplace
n_a + 1
P(a) = --------------
somme(n_i+1)ou des "m-estimateurs"
n_a + m
P(a) = ----------------
somme(n_i) + m(on dit qu'il y a m instances virtuelles)
(l'énoncé n'est pas de moi)
In the Land of Oz they never have two nice days in a row. If they have a nice day they are just as likely to have snow as rain in the next . If they have snow (or rain) they have an even chance of having the same in the next day. If there is a change from snow or rain, only half of the time is this change to a nice day. Vous avez prevu d'arriver au pays d'Oz un lundi mais pas dans n'importe quelles conditions : vous n'entreprenez le voyage qui s'il a plu au pays d'Oz le dimanche sinon vous repoussez votre voyage a la semaine suivante. 1. Quelle est la probabilite qu'il ne fasse pas beau le lundi de votre arrivee ? 2. Quelle est la probabilite qu'il fasse beau le samedi qui suit votre arrivee ? 3. Vous arrivez enfin et il neige. Quelle est la probabilite qu'il fasse beau le samedi ? 4. Le mardi vous dormez toute la journee et le mercredi il neige encore. Quelle est la probabilite qu'il ait fait beau pendant votre sommeil ? 5. Sur un temps assez long, peut-on prevoir la proportion de jours de neige, de pluie et de soleil au pays d'Oz ? Dans l'affirmative faut-il attendre tres longtemps pour que ces previsions soient bonnes ?
Tout cela pour dire que les chaines de Markov sont juste une autre présentation de la notion de probabilité conditionnelle, à la fois graphique et algorithmique.
On considere une population diploïde à deux allèles A et a, à generations non recouvrantes. Il y a trois génotypes possibles, b1=AA, b2=Aa et b3=aa. Il y a six types de couples possibles : E1 = AA × AA, E2 = AA × Aa, E3 = Aa × Aa, E4 = Aa × aa, E5 = AA × aa, E6 = aa × aa.
Modéliser cette situation par une chaîne de Markov en supposant que les accouplements se font entre frêres et soeurs (i.e., entre enfants issus d'un même couple).
(Indication : on pourra choisir les différents types de couples, E1, E2, ..., E6 comme états.)
Quelle est la fréquence limite des divers génotypes ?
Un modèle de Markov Caché (HMM, Hidden Markov Model), c'est une chaine de Markov dont les états sont labellés de manière probabiliste et dont on ne voit pas les états, mais juste les labels. En d'autres termes, c'est une chaine de Markov dont on ne voit pas les états mais dont les états émettent des symboles, de mani-re probabiliste.
Reprenons l'exemple de l'araignée se promenant sur son tétraèdre. Quand elle est coincée sur le sommet enduit de colle, elle pousse un cri de détresse avec une probabilité de 0.8 (et reste silencieuse, pétrifiée de peur, avec une probabilité de 0.2) ; quand elle est sur le sommet avec les mouches, elle pousse un crit de joie avec une probabilité de 0.8 (et reste silencieuse avec une probabilité de 0.2) ; quand elle est sur un autre sommet, elle pousse un cri de joie avec une probabilité de 0.1 (peut-être rêve-t-elle au sommet avec les mouches), un cri de détresse avec une probabilité de 0.1 (peut-être fait-elle un cauchemar) et reste silencieuse avec une probabilité de 0.8. On peut représenter ces probabilités sous la forme d'une matrice
silence .8 .8 .2 .2 cri de détresse .1 .1 0 .8 cri de joie .1 .1 .8 0
Imaginons qu'on ne ne voie pas l'araignée mais qu'on l'entende seulement : à l'aide de la nature de son cri à chaque instant, que peut-on dire de son trajet sur le tétraèdre ?
Le tableau suivant résume le vocabulaire des HMM.
tétraèdre matrice de transition comportement de matrice des probabilités d'émission l'araignée position de états l'araignée cris de l'araignée symboles, labels
Exemple d'application : reconnaissance de la parole. On construit un HMM pour chaque mot et, étant donné une suite de phonèmes, on essaye de trouver le HMM (ie, le mot) qui a le plus de chances de le produire.
Exemple d'application : reconnaissance de l'écriture manuscrite. C'est la même situation que précédemment : on essaye de construire un HMM pour chaque lettre et ensuite, étant donné un gribouilli, on essaye de sélectionne le HMM (ie, la lettre) le plus proche. (Ce problème se rencontre réellement, par exemple par la poste, pour le tri automatique du courrier (en particulier pour déchiffrer le code postal) ou par les PDAs, capables de reconnaître une écriture manuscrite simplifiée.)
Exemple d'application : analyse de texte (les états représentent la fonction des mots (sujet, complément, etc.), les symboles les mots).
Etant donné un corpus de textes d'un auteur (ou d'un genre donné, ou sur un sujet donné), construire un HMM le décrivant. Ensuite, étant donné un nouveau texte (par exemple un manuscrit récemment découvert, dont on cherche à identifier l'auteur), regarder s'il est raisonnable d'affirmer qu'il est du même auteur. Cela ne sert pas qu'à alimenter le débat "Ben Johnson est-il l'auteur des pièces de Shakespeare", mais aussi à classer automatiquement des documents électroniques, par exemple des pages Web, ou des articles de recherche.
Exemple d'application : appartenance d'une séquence d'ADN ou d'acides aminés à une famille.
Exemple d'application : prédiction des parties codantes et non codantes d'une séquence d'ADN : la chaine a deux états (codant, non codant) et les symboles sont des acides aminés.
Autre exemple : la segmentation. On a un seul HMM, mais les états sont typés (i.e., coloriés) et on cherche la suite de types (couleurs) la plus probable pour une séquence donnée. Exemple : trouver la classe grammaticale d'un mot.
Il y en a trois.
1. Trouver la probabilité d'une séquence de symboles : forward algorithm.
On peut transformer cette probabilité en un score d'"alignement" entre la chaine de Markov et la suite de symboles (par exemple, on peut construire une chaine de Markov qui décrive les parties codantes et une autre les parties non codantes et choisir celle qui décrit le mieux une séquence donnée ; idem pour les îlots CpG, etc.)
2. Trouver le chemin (la suite d'états) la plus probable, étant donnée une séquence de symboles : algorithme de Viterbi.
3. Trouver les paramètres d'une chaine de Markov (matrice de transition, matrice des probabilités d'émission -- on connait déjà la topologie de la chaine de Markov) : forward-backward, aka Baum-Welch, aka EM.
Le problème est le suivant : on dispose d'une chaine de Markov cachée et d'une suite de symboles produite (peut-être) par cette chaine. On cherche la probabilité d'observer cette suite avec cette chaine. Cela permet en particulier de choisir entre plusieurs chaines de Markov (par exemple, des chaines de Markov décrivant des domaines protéiques différents ou des familles de protéines) celle qui décrit le mieux la suite (par exemple, à quelle famille la protéine appartient). On parle parfois de calcul de score, voire d'alignement d'une suite de symboles sur une chaine de Markov.
Si on connaissait la suite d'états, ce serait très facile. Mais justement, on n'a que la suite de symboles, pas les états.
L'idée est de procéder "récursivement", en calculant la probabilité d'observer les t premiers symboles sachant qu'on arrive à l'état e.
Comme d'habitude, pour éviter les problèmes d'instabilité numérique, on remplace les probabilités par leur logarithme et les produits de probabilités par la somme de leur logarithme.
NON TESTÉ
# Matrice de transition : mt[i,j] est la probabilité de passer de i à j
mt <- matrix(c( 1/4,1/2,1/4,
1/8,3/4,1/8,
7/10,1/10,1/5),
nr=3, byrow=T)
rownames(mt) <- c('p1','p2','p3')
colnames(mt) <- rownames(mt)
mt
# Matrice de probabilités d'émission
pe <- matrix(c( 1/4,1/4,1/4,1/4,
1/3,1/6,1/6,1/3,
1/5,3/10,1/4,1/4),
nr=3, byrow=T)
rownames(pe) <- rownames(mt)
colnames(pe) <- c('A','C','G','T')
pe
# Suite de symboles observés
x <- "TTACGGATCGGGTATC"
x <- strsplit(x,'')[[1]]
# Algorithme forward : calcul du score de l'alignement d'un HMM et d'une séquence
forward <- function (mt,pe,x) {
if(length(x)==1) {
x <- strsplit(x,'')[[1]]
}
k <- length(x)
n <- dim(mt)[1]
m <- matrix(-Inf, nr=n, nc=k)
pe <- log(pe)
mt <- log(mt)
m[,1] <- pe[,x[1]]
for(i in 2:k) {
for (j in 1:n) {
m[j,i] <- sum( m[,i-1] + mt[,j] + pe[j,x[i]] )
}
}
m
}
forward(mt,pe,x)A FAIRE : score local (calcul du score avec l'algorithme forward, mais
Le problème des scores calculés précédemment, c'est qu'ils dépendent de la longueur de la séquence. Pour faciliter les comparaisons, on peut considérer des (logarithmes de) rapports de vraissemblance :
Probabilité d'observer cette séquence avec ce HMM
log -----------------------------------------------------------
Probabilité d'observer cette séquence avec un HMM trivialMais il n'y a pas de choix unique sur le HMM trivial : on peut le construire à l'aide des fréquences des acides aminés présents dans la séquence, à l'aide des fréquences des acides aminés présents dans le génome de l'espèce considérée, etc.
Très souvent, ce problème ne se pose pas, car on veut comparer deux modèles : il suffit de calculer le quotient de leur vraissemblance.
Le problème est le suivant : on connait la chaine de Markov et la suite de symboles qu'elle a produit et on cherche la suite d'états la plus probable.
L'idée est de procéder "récursivement" en cherchant, pour chaque instant t et chaque état e, quelle est la suite d'états la plus probable produisant les t premiers symboles et s'achevant à l'état e.
On s'apperçoit rapidement que le problème n'est pas récursif, mais qu'il relève de la programmation dynamique (si on trace le graphe des appels de fonctions, on constate qu'il ne forme pas un arbre) : on va stocker, dans un tableau, la probabilité que les t premiers symboles aient été produits par une suite d'états se terminant par e puis, une fois le tableau construit, le lire dans l'autre sens pour trouver le chemin le plus probable.
Pour éviter les problèmes d'instabilité numérique, on utilisera des logarithmes de probabilité, qu'on ajoutera, au lieu de multiplier des probabilités.
NON TESTÉ
# Matrice de transition : mt[i,j] est la probabilité de passer de i à j
mt <- matrix(c( 1/4,1/2,1/4,
1/8,3/4,1/8,
7/10,1/10,1/5),
nr=3, byrow=T)
rownames(mt) <- c('p1','p2','p3')
colnames(mt) <- rownames(mt)
mt
# Matrice de probabilités d'émission
pe <- matrix(c( 1/4,1/4,1/4,1/4,
1/3,1/6,1/6,1/3,
1/5,3/10,1/4,1/4),
nr=3, byrow=T)
rownames(pe) <- rownames(mt)
colnames(pe) <- c('A','C','G','T')
pe
# Suite de symboles observés
x <- "TTACGGATCGGGTATC"
x <- strsplit(x,'')[[1]]
viterbi <- function (mt,pe,x) {
if(length(x)==1) {
x <- strsplit(x,'')[[1]]
}
pe <- log(pe)
mt <- log(mt)
n <- dim(mt)[1] # Nombre d'états
k <- length(x) # longueur de la chaine de symboles émise
# Construction de la matrice de programmation dynamique
m <- matrix(-Inf, nr=n, nc=k)
m[,1] <- pe[,x[1]]
for (i in 2:k) {
for (j in 1:n) {
m[j,i] <- max(m[,i-1] + mt[,j] + pe[,x[i]])
}
}
rownames(m) <- rownames(mt)
colnames(m) <- x
print(m)
# Reconstruction du chemin
chemin <- rep(NA, k)
chemin[k] <- which( m[,k] == max(m[,k]) )[1]
for (i in k:2) {
j <- chemin[i]
chemin[i-1] <- which( m[j,i] == max(m[,i-1] + mt[,j] + pe[,x[i]]) )[1]
}
ch <- rownames(pe)[chemin]
print(ch)
image(t(m))
lines( seq(0,1,length=k), seq(0,1,length=n)[chemin] )
invisible(list(matrice=m, chemin=ch))
}
r <- viterbi(mt,pe,x)Voir par exemple :
http://www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/viterbi_algorithm/s1_pg1.html
Un exemple de topologie très populaire pour examiner des protéines : des états d'appariement, des états d'insertion, des états de délétion (qui n'émettent aucun symbole), un état initial et un état final.
A FAIRE : un dessin
A l'aide d'un alignement multiple (ou d'un PSSM), on peut facilement obtenir les probabilités de transition et d'émission : les chaines de Markov avec une telle topoloie sont simplement une autre manière de coder les PSSM.
A FAIRE : le même dessin que ci-dessus, avec des probabilités, et le PSSM correspondant.
Le fait que certaines probabilités du modèle obtenu soient exactement égales à zéro est problématique : il interdit absolument de reconnaître une séquence qui diffèrerait d'un caractère. Pour éviter cela, on impose que ces probabilités soient au moins égales à une certaine valeur. C'est la régularisation d'une chaine de Markov.
Pour construire un modèle de Markov caché, i.e., pour trouver quelles sont les propabilités (initiales, de transition et d'émission) qui permettent de décrire le mieux un phénomène, nous utiliserons l'algorithme EM (Expectation-Maximization). Comme cet algorithme est très général, nous allons le présenter dans un autre contexte, sans aucun rapport, mais plus simple.
Le problème est le suivant : on dispose d'un nuage de points et on essaye de les regrouper en plusieurs classes.
Par exemple, on a mesuré la longueur, l'envergure, le poids de quelques dizaines de papillons, ce qui nous donne quelques dizaines de points dans un espace de dimension trois ; on sait qu'il y a trois espèces de papillons différentes mais on n'a pas été capable de les distinguer lors des mesures, peut-on regrouper ensemble les papillons d'une même espèce ?
Autre exemple : on a mesuré le nombre d'ARNm de quelques milliers de gènes dans des cellules normales et des cellules cancéreuses, ce qui nous donne quelques milliers de points dans un espace de dimension deux ; on essaye de regrouper les gènes ayant le même profil d'expression.
Autre exemple : on a une photo d'un spot de biopuce et on aimerait classer les pixels en deux groupes, selon leur appartenance au spot (brillants) ou au fond (moins brillants).
Voici une manière de procéder : répartir les points dans trois groupes, au hasard ; calculer le centre de gravité (ça veut dire, le point moyen) de chacun des groupes ; changer les points de groupe en les mettant dans le groupe dont le centre est le plus proche ; itérer.
C'est l'algorithme des k-moyennes.
C'est le principe général de l'algorithme EM : on cherche la valeur de certains paramètres (ici, les coordonnées du point moyen de chaque groupe) ; on commence par prendre des valeurs au hasard ; on fait quelques calculs d'espérance (i.e., des calculs de moyenne : ici, on met les observations dans le groupe le plus proche et on calcule le point moyen de ce groupe) ; on en déduit de nouvelles estimations des paramètres ; on itère.
On peut reformuler cet algorithme ainsi :
0. Définir les classes au hasard, i.e., choisir des points m1,m2,...
et mettre chaque point x_i dans la classe G_j du m_j le plus proche.
1. E-step : Calculer la vraissemblance, i.e.
Produit P(x_i \in G_j)
où G_j est le groupe dans lequel on a mis x_i, qui contient les
points plus proche de m_j que des autres.
Si les X_i suivent une distribution normale (en dimension 1, pour
simplifier, mais ça ne change rien), on a
P(x_i \in G_j) = exp( - (x_i - m_j)^2 / 2 sigma )
2. M-step : calculer ces paramètres pour qu'ils maximisent la
vraissemblance. Ici, cela revient à minimiser la somme des
(x_i - m_j)^2 (les sigma disparaissent), i.e., à remplacer m_j
par la moyenne des x_i qui sont dans ce groupe.
3. Mettre chaque point dans sa classe la plus probable.
Itérer.Autre exemple d'application d'EM :
http://www.dcs.kcl.ac.uk/teaching/units/csmacmb/DOC/lecture16b.pdf
Remarque. Cet algorithme est en fait bien plus général : il permet de calculer des estimateurs dumaximum de vraissemblance.
Exercice : utiliser l'algorithme EM pour trouver les paramètres d'un mélange de gaussiennes (Indication : c'est une généralisation de l'algorithme des k-moyennes, avec deux différences : d'une part, les différents groupes n'ont pas nécessairement la même variance ; d'autre part, chaque groupe contient une proportion p_j des points). Exercice supplémentaire : et si la variance des groupes n'est pas une matrice scalaire ? Exercice supplémentaire : et si on ajoute
(Normalement, vous avez poussé un cri en lisant les deux questions précédentes, car elles font exploser le nombre de paramètres à estimer. On considère généralement qu'il faut au moins dix fois plus de données que de paramètres. De plus, un modèle trop compliqué risque de mieux décrire les données simplement par hasard : on parle de "surmodélisation" ("overfitting", en anglais).)
On peut aussi utiliser cet algorithme pour la reconstruction d'arbres phylogénétiques et en cartographie.
http://www.dcs.kcl.ac.uk/teaching/units/csmacmb/DOC/lecture15b.pdf
Jusqu'à présent, on supposait que le modèle de Markov était connu : nous allons maintenant voir comment on peut le construire. On se donne : un ensemble de suites, pour lesquelles on aimerait construire un HMM, et un ensemble d'états. Il faut trouver : les probabilités initiales, les probabilités de transition et les probabilités d'émission.
L'idée est la suivante : on prend des probabilités au hasard, on calcule les scores des différentes séquences pour ce HMM, on en déduit de meilleures estimations des probabilités et on itère.
Cet algorithme calcule un HMM qui maximise la vraissemblance des séquences qu'on lui a donn augmenter ses chances d'avoir un maximum global, on peut lancer l'algorithme plusieurs fois (avec des probabilités initiales différentes). Une autre id estimations à chaque étape (de moins en moins de bruit au fur et à mesure que l'on avance) -- c'est le recuit simulé.
Un autre problème, c'est que l'ensemble de séquences utilisées pour l'apprentissage de l'algorithme peuvent contenir des séquences très semblables, qui risquent produire un HMM biaisé. Pour corriger ce biais, on peut calculer un arbre phylogénétique (enraciné) sur ces séquences et leur accorder un poids : on part de la racine avec un poids de 1 et à chaque embranchement, on divise le poids entre les deux branches (c'est exactement comme les lois de Kirchoff, pour le calcul de l'intensité d'un courrant électrique).
Une autre manière de corriger ce biais consiste à affecter des poids aux différentes séquences d'autant plus importants que leur score est mauvais -- les poids changent donc à chaque itération. C'est très proche du "boosting", en statistique, pour ceux qui connaissent (ça consiste à "stabiliser" un estimateur en le calculant sur de nombreux échantillons de bootstrap -- mais, contrairement au bagging, ces échantillons de bootstrap ne sont pas choisis au hasard : les observations mal classées/prédites ont une plus grande chance d'être choisies).
Il y a encore d'autres manières de corriger ce biais, par exemple à l'aide de poids calculés pour chaque colonne, poids d'autant plus faibles que l'acide aminé en question est présent souvent.
Lors de l'apprentissage d'une chaine de Markov cachée, au lieu de considérer les acides aminés un par un, on peut les regarder en groupes : les acides aminés hydrophiles, les hydrophobes, etc. Ces classes peuvent se recouper : un même acide aminé peut être dans plusieurs classes. On ne cherche plus des probabilités de transition d'un acide aminé à un autre, mais d'une classe à une autre -- on peut en déduire les probabilités de transition d'un AA à un autre.
On peut aussi voir ça comme une forme de régularisation.
On peut calculer un alignement multiple comme suit : construire un HMM qui reconnaisse l'ensemble des séquences (en prenant la topologie "classique", présentée plus haut) et trouver le chemin le plus probable à l'aide de l'algorithme de Viterbi.
Je ne connais rien à ce qui suit.
Hybrids of HMMs and neural nets
dynamic Bayesian nets
factorial HMMs
Boltzmann trees
hidden Markov random fields <-- je connais
(reprendre l'exemple de la
reconnaissance de l'écriture
manuscrite)
http://www.cs.tufts.edu/comp/150B/lectures/RNA2DStructureLecture.ppt
Les bases des brins d'ARN ont tendance à s'apparier, donnant à la molécule une certaine structure tridimensionnelle.
On peut représenter cette structure par un dotplot de la séquence avec elle-même. On peut aussi écrire la séquence le long d'un cercle et tracer les cordes correspondant aux bases appariées.
On suppose généralement qu'il n'y a pas de pseudo-noeud, i.e., que les cordes de la représentation graphique ne se coupent pas.
On peut aussi représenter un pseudo-noeud ainsi :
--------------------
| ------------------ |
|| ---------------- ||
||| -------------- |||
|||| ------------ ||||
||||| ---------- |||||
|||||| -------- ||||||
||||||| ------ |||||||
|||||||| ||||||||
GCGCGCGCATATATGCGCGCGCATATAT
|||||| ||||||
||||| -------- |||||
|||| ---------- ||||
||| ------------ |||
|| -------------- ||
| ---------------- |
------------------
C'est toujours de la programmation dynamique -- je ne détaille pas.
Algorothme d'Ussinov-Jacobson Algorithme de Zuker
Un langage est un ensemble de mots (un mot, c'est une suite de lettres, une lettre c'est un élément d'un ensemble fixé une fois pour toute). Par exemple, "tous les mots de la langue anglaise" ou "tous les mots écrits avec les lettres a et b" ou "tous les mots commençant par 'aaabab'" ou "tous les mots contenant 'bbbababa'" ou "tous les mots contenant autant de a que de b".
Une grammaire, c'est un ensemble de règles permettant de décrire un langage, de la forme :
S --> aXaab S --> X X --> aXa X --> b aX --> Xb
Ici, les majuscules sont des "symboles non terminaux", les minuscules désignent des lettres : on commence par S (comme "Start") et on applique les règles, au choix, dans n'importe quel ordre, jusqu'à obtenir un mot sans symbole non terminal -- le langage décrit par la grammaire est l'ensemble des mots que l'on peut obtenir de cette manière.
On regroupe généralement les règles ayant la même prémice (la barre verticale se lit "ou") :
S --> aXaab | X X --> aXa | b aX --> Xb
Voici un exemple :
règle résultat ---------------------------------- S --> aXaab aXaab aX --> Xb Xbaab X --> aXa aXabaab X --> aXa aaXaabaab aX --> Xb aXbaabaab X --> b abbaabaab
Ainsi, le mot "abbaabaab" fait donc partie du langage décrit par la grammaire ci-dessus.
N. Chomsky (ça n'est pas un informaticien, mais un linguiste) distingue une hiérarchie de quatre classes de grammaires.
Une grammaire est régulière (on parle aussi d'expression régulière -- oui, ce sont les mêmes que sur votre éditeur de texte préféré) si ses règles sont de la forme
W --> a V
ou
W --> V a
où V et W sont des symboles non terminaux et a est une (ou plusieurs) lettre(s).
Les langages réguliers peuvent être décrits par des automates finis. En particulier, les expressions bien parenthésées ne peuvent pas être décrites par une grammaire régulière.
Une CFG (Context-Free Grammar) est une grammaire dont les règles sont de la forme
W --> S
où W est un symbole non terminal et S une suite quelconque de symboles non terminaux et de lettres. On les implémente à l'aide d'automates à pile.
Exemple : les palindromes (exemple : babab)
S --> a S a | b S b | aa | bb
Par contre, les copies (exemple : baba) ne sont pas décrites par une CFG.
On rencontre aussi des grammaires sensibles au contexte
a W b --> a S b
et même des grammaires quelconques
a W b --> S
On peut utiliser une grammaire (CFG) de manière aléatoire, afin d'obtenir des mots aléatoires du langage qu'elle décrit. On peut formaliser cela en ajoutant une probabilité à chacune des branches des alternatives, i.e., en choisissant une probabilité pour chacune des transformations.
Ainsi, on peut transformer la grammaire
S --> aXaab | X X --> aXa | b | Xb
en une grammaire stochastique :
S --> aXaab avec une probabilité 0.7 S --> X avec une probabilité 0.3 X --> aXa avec une probabilité 0.2 X --> b avec une probabilité 0.5 X --> Xb avec une probabilité 0.3
Les grammaires stochastiques généralisent les chaines de Markov cachées -- les chaines de Markov cachées correspondent au cas des langages réguliers.
On peut décrire une structure secondaire d'ARN par une grammaire :
S --> aS | cS | uS | gS (non appariement) S --> Sa | Sc | Su | Sg (non appariement) S --> aSu | uSa | cSg | gSc (appariement) S --> SS (bifurcation)
On procède alors comme pour les modèles de Markov cachés : on part d'un ensemble de séquences dont on connait la structure et on en déduit les probabilités du modèle ; après cette phase d'apprentissage, on peut soumettre une nouvelle séquence à la grammaire stochastique et lui demander le chemin le plus probable, i.e., la structure secondaire la plus probable.
L'apprentissage se fait toujours avec un algorithme EM.
L'alignement de la grammaire avec une séquence se fait grace à l'algorithme de Cocke-Younger-Kasami (CYK).
(Comme pour les chaines de Markov, on peut aussi s'intéresser à un troisième problème : calculer le score d'un alignement entre une séquence et une grammaire stochastique, afin de choisir, parmi plusieurs grammaires, la plus proche de la séquence étudiée.)
http://www.dcs.kcl.ac.uk/teaching/units/csmacmb/DOC/lecture18b.pdf http://www.imb-jena.de/RNA.html http://scor.lbl.gov/index.html http://www.rnabase.org/metaanalysis/
On peut décrire quelques déplacement élémentaires dans l'espace des structures secondaires (ajouter, retirer ou déplacer un appariement) et ainsi utiliser des algorithmes usuels d'optimisation (recuit simulé, algorithmes génétiques, etc.) pour trouver une structure optimale.
http://www.santafe.edu/%7Ewalter/Papers/kinfold.pdf
Généralement, les bases ne sont pas trop conservées, par contre la structure secondaire l'est. Il ne faut donc pas tant chercher les "séquences semblables" que les séquences dont la structure secondaire est semblable.
Présentation de la bioinformatique pour les non-informaticiens (pour ceux qui veulent utiliser la bioinformatique, pas en faire), qui se contenteront d'utiliser les outils sur le Web. Je crois me souvenir que le sous-titre est "Copy-paste for dummies".
Livre très gros et très cher (mais la bibliothèque en a heureusement un exemplaire) qui contient à peu près tout.
C'est visiblement LA référence pour tout ce qui relève des manipulations de chaines ou de textes, en particulier en biologie.
Je ne l'ai pas lu.
Je ne l'ai pas lu non plus, mais il parait que ce livre contient une présentation très claire des algorithmes d'alignement.
Notes assez succintes, mais souvent claires.
Exemple d'application du test du Chi2 pour trouver les cadres de lecture (les différents nucléotides n'apparaissent pas avec la même fréquence en première, deuxième ou troisième position d'un codon -- on peut rafiner et prendre les codons eux-mêmes).
Modélisation d'une séquence nucléique par une chaine de Markov d'ordre 0 (on se donne juste la fr 1, 2, 3, etc. et application au calcul de la fréquence théorique d'apparition d'un mot donn logiciels de recherche dans les bases.
Autre application des chaines de Markov : reconnaissance des cadres de lecture (on construit sept chaines de Markov : six qui décrivent des séquences codantes dans les six cadres de lecture, une septième qui décrit les séquences non codantes ; ces chaines sont construites à l'aide de gènes ou de non-gènes connus ; étant donnée une nouvelle chaine, on regarde quelle est la chaine de Markov qui la décrit le mieux).
Algorithmes de recherche de chaine : naïf, Hancart, KMP, Boyer-Moore.
Algorithmes d'alignement : naïfs, programmation dynamique, alignement global, semi-global ("bestfit"), local.
Z-score d'un alignement : calculer le score de l'alignement après avoir mélangé les séquences (plusieurs fois), calculer la moyenne et l'écart-type, comparer le score obtenu avec cette moyenne : s'il est à plus de 5 écarts-types, c'est bon.
Matrices de similarité. Inconvénients de la matrice PAM : hypothèses simplificatrices (indépendance des mutations (certaines parties de la séquence peuvent être changées sans problème, alors que d'autres o=sont vitales pour les propriétés de la protéine), symétrie des mutations) ; utilisation de séquences très proches (biais : les mutations les plus observées sont silencieuses...). Matrice BLOSUM : utilisation de séquences plus éloignées.
Le document explique comment les calculer -- mais tous les symboles mathématiques ont disparu dans la conversion en PDF, c'est incompréhensible...
Pour les alignements avec brèches, les différentes matrices donnent des résultats comparables -- par contre, les pénalités de bréches sont très importantes.
Présentation des matrices PAM et BLOSUM (pour PAM, ça me semble trop complexe, pour ne pas dire faux...)
Présentation du dotplot, de l'alignement global, local, semi-global.
Un cours d'algorithmique pour biologistes, AVEC DES EXERCICES (pour biologistes).
1. Recherche d'une sous-chaine dans une chaine circulaire, en adaptant un algorithme de recherche normal. Calcul de la complexité de l'algorithme naïf.
2. Calculs de complexité sur l'algorithme de Boyer-Moore.
3. Bornes supérieure et inférieure sur le nombre de comparaisons dans l'algorithme KMP.
Trouver une chaine de taille 10 et une famille de textes (de tailles de plus en plus grandes) telles que KMP requiert moins de comparaisons que BM (ou montrer qu'il n'en existe pas).
Idem dans l'autre sens (KMP>BM).
4. Bornes inférieure et supérieure sur la taille d'un arbre de suffixes (ou, plus précisément, du trie compact d'un ensemble de mots dont aucun n'est préfixe d'un autre).
Construction, à la mainn d'un arbre de suffixes.
5. Algorithme d'Ukkonen et calculs probabilistes sur l'arbre obtenu.
6. Algorithme d'Ukkonen
7. Programmation dynamique pour calculer l'alignement de deux séquences. Montrer qu'il peut y avoir plusieurs alignements optimaux. Montrer que le nombre d'alignement optimaux peut être exponentiel.
8. Recherche de répétitions inexactes optimales (c'est toujours de la programmation dynamique). (Procéder comme pour l'alignement local, en initialisant la table de programmation dynamique avec des zéros sur la première ligne, la première colonne et la diagonale).
9. Alignement de deux séquences protéiques (avec des beaux dessins de tableaux de programmation dynamique, dévorés avec des flèches colorées).
10. Alignements multiples.
Calcul d'un score (SP : score des paires) d'un alignement multiple à l'aide d'une fonction de score (qui permet de calculer le score d'un alignement de deux séquences).
Description de l'algorithme d'alignement d'une séquence et d'un profil PSSM (c'est comme pour aligner deux séquences : on fait un tableau).
11. Algorithme de l'étoile pour l'alignement multiple : dans l'ensemble de séquences à aligner, on cherche la séquence "la plus au milieu" et on aligne les autre séquences sur celle-là.
Ils parlent aussi d'alignement phylogénétique.
12. Après des rappels sur la construction des matrices PAM, on nous demande d'en déconstruire une, i.e., de retrouver les effectifs qui ont permis de la construire, puis d'en déduire une autre matrice PAM.
PAM20 = (PAM10)^2 (les log-odds se multiplient)
L'autre exercice nous demande de construire une chaine de Markov (avec un état par acide aminé, plus un état initial et un état final) qui reconnait une protéine, puis de calculer la vraissemblance correspondante. La dernière question nous demande de construire une chaine de Markov qui reconnaisse deux séquences.
# Construction d'une matrice de Markov
markov <- function (x) {
x <- strsplit(x,'')[[1]]
x <- factor(x)
aa <- strsplit("ACDEFGHIKLMNPSQRTVWY01",'')[[1]]
n <- length(aa)
m <- matrix(0, nr=n, nc=n)
colnames(m) <- aa
rownames(m) <- aa
m["1","1"] <- 1
m["0", x[1]] <- m["0", x[1]] +1
for (i in 1:(length(x)-1)) {
m[ x[i], x[i+1] ] <- m[ x[i], x[i+1] ] +1
}
m[ x[length(x)], "1" ] <- m[ x[length(x)], "1" ] +1
# C'est une matrice d'effectifs, on la transforme en une matrice de
# probabilités : la somme de chaque ligne doit être égale à 1
m <- m +.001
m <- m/apply(m,1,sum)
print(round(m, digits=2))
invisible(m)
}
x <-"MAKGVAVLNSSEGVTGTIFFTQEGDGVTTVSGTVSGLKPGLHGFHVHALGDTTNGCMSTGPHFNPDGKTHGAPEDANRHAGDLGNITVGDDGTATFTITDCQIPLTGPNSIVGRAVVVHADPDDLGKGGHELSLATGNAGGRVACGIIGLQG"
m <- markov(x)
# Calcul de la vraissemblance
logv <- function (m,x) {
x <- strsplit(x,'')[[1]]
x <- factor(x)
v <- 0
v <- v + log(m["0",x[1]])
for (i in 1:(length(x)-1)) {
v <- v + log(m[x[i],x[i+1]])
}
v <- v + log(m[x[length(x)],"1"])
v
}
logv(m,x)
y <- "MGKAVVVLSSSEGVSGTVYFAQEGDGPTTVTGNVSGLKPGLHGFHVHALGDTTNGCMSTGPHYNPNGKEHGAPEDDVRHAGDLGNITVGDDGTATFTIIDSQIPLSGPNSIVGRAVVVHAEPDDLGRGGHELSKTTGNAGGRVACGIIGLQG"
logv(m,y)
m <- (markov(x) + markov(y))/2
logv(m,x)
logv(m,y)
z <- paste( sample(strsplit("ACDEFGHIKLMNPSQRTVWY",'')[[1]], 150, replace=T),
collapse='')
logv(m,z)Exemples d'examen : 1. Trouver un algorithme linéaire qui recherche des sous-chaines répétées. 2. Optimiser l'algorithme de Boyer-Moore dans un cas particulier : AAACAAA (ou, plus généralement, A^kCA^k) 3. Trouver un algorithme qui recherche des sous-chaines semblables (gènes paralogues). 4. Expliquer comment construire un arbre d'Aho-Corasick (un arbre pour reconnaître un ensemble de mots, avec des liens d'échec).
Exercices des années antérieures :
1. Utilisation de sites internet pour récupérer des séquences.
Construction d'un arbre phylogénétique :
atgccttacaggcta atgccttatgggctg atgccctacaggcta atgccctacaggctg atgccttacaggctg
2. Calculs de probabilités avec des sites de restriction).
Réarranger ABCDEFGHIJ en DGCJEIBAFH à l'aide de renversements et de transpositions.
3. Construction de cartes de restriction à l'aide de la longueur des fragments de restriction (PQ-tree : c'est de la théorie des graphes). Cas d'un chromosome circulaire.
http://courses.cs.vt.edu/~algnbio/algnbio_2002/handouts/homeworks/hw3.php.html
4. Graphes d'intervalles
http://courses.cs.vt.edu/~algnbio/algnbio_2002/handouts/homeworks/hw4.php.html
5. Alignement local
6. Alignement local
7. Fasta
8. Blast
9. PSSM
Un exercice sur les graphes.
HMM :
An open reading frame (ORF) for a gene begins with a Start codon --- ATG --- and ends with one of three Stop codons --- TAA, TAG, or TGA. Of course, the length n of an ORF must be divisible by 3. Suppose that in Your Favorite Genome (YFG), it has been found that the frequency of the individual stop codons varies as follows: 0.25 is the frequency of TAA, 0.25 is the frequency of TAG, and 0.5 is the frequency of TGA. Further suppose that the frequency of each base in YFG is 0.25. Design a hidden Markov model (HMM) to recognize ORFs in YFG. You will need to make some assumptions. Give your assumptions and explain how they fit into your HMM.
10a. Appliquer l'algorithme de Viterbi : on a observé la chaine TTACGGATCGGGTATC à la sortie de la chaine de Markov cachée de matrice de transition
mt <- matrix(c( 1/4,1/2,1/4,
1/8,3/4,1/8,
7/10,1/10,1/5),
nr=3, byrow=T)
rownames(mt) <- c('p1','p2','p3')
colnames(mt) <- rownames(mt)
mtet de probabilités d'émission
pe <- matrix(c( 1/4,1/4,1/4,1/4,
1/3,1/6,1/6,1/3,
1/5,3/10,1/4,1/4),
nr=3, byrow=T)
rownames(pe) <- rownames(mt)
colnames(pe) <- c('A','C','G','T')
peQuelle est la suite d'états la plus probable ?
10b. On donne 5 séquences de dix caractères.
TCTTACATCT ACTTCCATTT ACTGCGATTT CCCTACATGT CACTGCAAGT
Construire un arbre phylogénétique à l'aide de l'algorithme UPGMA.
Construire un arbre phylogénétique à l'aide de l'algorithme NJ (Neighbor Joining).
Construire l'arbre phylogénétique de moindre parcimonie (le problème général est NP-complet, mais ici, on peut utiliser une recherche exhaustive ou effectuer un raisonnement ad hoc).
http://courses.cs.vt.edu/~algnbio/algnbio_2000/handouts/PearsonLipman_PNAS.pdf
Présentation de l'algorithme Fasta.
Autres articles sur Fasta :
http://courses.cs.vt.edu/~algnbio/FASTA.php.html
http://courses.cs.vt.edu/~algnbio/algnbio_2002/handouts/atga_letter.ps
http://www.student.math.uwaterloo.ca/~cs762/Notes/index.php
http://www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/main.html
Présentation très claire (avec quelques dessins) des modèles de Markov cachés et des algorithmes associés : algorithme forward pour le calcul du score de l'alignement d'un HMM et d'une séquence, algorithme de Viterbi pour calculer la suite d'états la plus probable, algorithme Forward-Backward pour l'apprentissage d'un HMM (pour ce dernier, ils se contentent de dire que c'est plus compliqué).
http://www.cse.ucsc.edu/research/compbio/ismb99.handouts/KK185FP.html
Présentation des PSMM et des HMM.
http://courses.cs.vt.edu/%257Ealgnbio/algnbio_2001/handouts/AltschulEtAlNAR_1997.pdf http://courses.cs.vt.edu/%257Ealgnbio/algnbio_2001/handouts/KarlinAltschulPNAS_1990.pdf http://courses.cs.vt.edu/%257Ealgnbio/algnbio_2001/handouts/KarlinAltschulPNAS_1993.pdf Interprétation des scores d'un blast : http://www.ncbi.nlm.nih.gov/BLAST/tutorial/Altschul-1.html
Vincent Zoonekynd
<zoonek@math.jussieu.fr>
latest modification on Wed Sep 8 23:06:37 CEST 2004