La fonction de vraisemblance, définie comme la probabilité des données sachant le modèle, est l'outil central des procédures d'inférence statistique basées sur un modèle probabiliste, comme c'est le cas en phylogénie moléculaire - le modèle Markovien décrit le taux instantané auquel les différents types de changement de nucléotides/amino-acides se produisent, et l'on dérive aisément les probabilités de changements de long terme, puis les probabilités d'obtention de tel ou tel pattern de variation entre séquences (Galtier et al 2005, Bryant et al 2005). L'utilisation standard de la fonction de vraisemblance est sa maximisation sur l'espace des paramètres: on considère comme les plus vraisemblables les valeurs de paramètres qui confèrent une probabilité maximale aux données. En phylogénie moléculaire, la difficulté est d'estimer correctement le paramètre discret qu'est la topologie de l'arbre - l'espace des arbres permet de définir des distances, mais pas d'ordre, donc pas d'intervale, et pas d'encadrement de valeurs. L'algorithme PhyML (Guindon & Gascuel 2003 Syst Biol 52:696), développé au LIRMM, est un des standards actuels de l'optimisation de topologie au maximum de vraisemblance. Deux développements nouveaux ont vu le jour au cours du projet MODEL_PHYLO.
D'une part, la méthode PhyML a été étendue au cas de modèles non-homogènes et non-stationnaires, c'est-à-dire de modèles autorisant la composition en bases des séquences à évoluer (Galtier & Gouy 1998 Mol Biol Evol 15:871). Contrairement au cas stationnaire, la position de la racine de l'arbre influence la vraisemblance sous ce type de modèle. L'algorithme PhyML étant par nature non-enraciné, il a dû être adapté, notamment via la définition des vraisemblances conditionnelles supérieures aux noeuds. Le modèle de Galtier & Gouy, qui prévoyait un taux de GC d'équilibre indépendant par branche de l'arbre, est par ailleurs modifié: le biais GC affecté à chaque branche est tiré d'un mélange fini, ce qui accélère encore les calculs. Ce travail permet donc la reconstruction rapide d'arbres phylogénétiques même lorsque le pourcentage de A, C G et T varie entre séquences, et l'estimation des composition en bases ancestrales (Boussau & Gouy 2006).
D'autre part, l'algorithme PhyML lui-même a été amélioré grace à l'introduction de mouvements topologiques de plus grande envergure. A la base, PhyML explore l'espace des topologies au travers de mouvements de type NNI (pour Nearest Neighbor Interchange), par lesquels l'arbre ((A,B)(C,D)) devient ((A,C)(B,D)) ou ((A,D)(B,C)), A, B , C, et D étant des groupes de taxons (feuilles) inchangés. Un NNI consiste donc à intervertir les (groupes de) taxons autour d'une branche interne donnée. L'intérêt algorithmique est qu'une bonne partie du calcul de la vraisemblance de l'arbre initial peut être utilisée pour le nouvel arbre, et notamment tout ce qui est interne aux groupes A, B, C et D. La pratique a toutefois montré que PhyML_NNI était parfois piégé dans des maxima locaux. Pour remédier à ce problème, une nouvelle classe de mouvements, les SPR (pour Subtree Pruning Regrafting) ont été introduits. Il s'agit de mouvements dans lesquels un sous-arbre est coupé, puis collé à une nouvelle position, éventuellement distante de sa position initiale - un NNI est un SPR de courte distance. L'introduction des SPR ne peut se faire que de manière raisonnée puisque (i) leur nombre est quadratique en le nombre de taxa (linéaire pour les NNI), et (ii) un SPR demande plus de re-calcul qu'un NNI. Un filtre basé sur les méthodes des distances (Gascuel & Steel 2006) a donc été mis en place pour sélectionner les SPR prometteurs. Le nouvel algorithme améliore significativement les performances de PhyML à un coût acceptable en termes de temps de calcul (Hordijk & Gascuel 2005).
L'approche Bayésienne, et l'intégration numérique permise par les algorithmes de Monte-Carlo (par exemple, Metropolis-Hastings), permettent l'implémentation de modèles complexes que l'on ne peut pas traiter au maximum de vraisemblance en raison d'un trop grand nombre de paramètres, ou d'un calcul analytique impossible de la fonction de vraisemblance. L'approche bayésienne nécessite toutefois la définition d'hypothèses a priori sur la valeur des paramètres, sous forme de lois de probabilités (ce qui peut être vécu comme un atout ou comme un prix à payer, selon les circonstances). Les algorithmes MCMC explorent l'espace des paramètres à la manière du recuit-simulé, les mouvements vers un état moins probable étant autorisés.
L'approche bayésienne a été introduite en phylogénie moléculaire à la fin des années 90 (par ex. le programme MrBayes), et a sous-tendu un important effort de modélisation de l'évolution des molécules. Notamment, Lartillot & Philippe (2004) ont proposé un modèle prenant en compte les contraintes site-spécifiques. Dans ce modèle (appelé CAT), la fréquence stationnaire des nucléotides ou amino-acides varie entre sites - un peu comme si chaque position de la molécule possédait un ensemble d'états "autorisés" - alors que les modèles classiques affectent la même matrice de substitution à tous les sites. Sur le plan de la modélisation de l'évolution des séquences, deux voies majeures ont été explorées dans le cadre de MODEL_PHYLO.
Nous avons tout d'abord mis en place un modèle non-stationnaire sous lequel la composition en bases ou en amino-acides des séquences peut varier, au même titre que dans Galtier & Gouy 1998 et Boussau & Gouy 2006 (cf ci-dessus). Le cadre bayésien autorise une modélisation des variations de fréquences stationnaires statistiquement plus satisfaisante (car plus plausible biologiquement, plus économe en paramètres, et plus robuste à des variations dans l'échantillonnage des taxons). Plutôt que d'affecter un spectre de fréquence stationnaire par branche de l'arbre (Galtier & Gouy 1998), ce qui amène vite à une sur-paramétrisation lorsque la taille du jeu de données augmente, nous proposons ici qu'un nombre fini (aléatoire) de changements de spectre (points de cassure) se produisent dans la phylogénie, le taux de changements étant un paramètre libre estimable. Chaque point de cassure définit un changement discret des fréquences d'équilibre du processus de substitution (Blanquart & Lartillot 2006). De cette manière, on obtient un modèle dans lequel les processus de spéciation et de variations de composition nucléotidiques se trouvent totalement découplés a priori.
Plus récemment, nous avons réussi à combiner ce processus de points de cassure (qui décrit les variations globales de composition en bases/amino-acides, typiquement dues à des effets mutationnels) avec le modèle CAT (qui modélise la variation site-spécifique, reflétant typiquement les contraintes sélectives). Le modèle CAT-BP (Blanquart & Lartillot, soumis) prend ainsi en compte, de manière relativement économique en nombre de paramètres, les sources de variations majeures affectant la composition en bases/amino-acides, à la fois le long des séquences et au cours du temps - deux des biais les plus critiques des modèles classiques de phylogénie moléculaire.
Les modèles décrits dans le paragraphe précédent sont de type empirique et non-paramétrique : à savoir, ils sont capables d'accommoder les variations des processus d'évolution observées, en particulier, le long des séquences, mais sans toutefois chercher à les expliquer. Ceci est très utile pour reconstruire les phylogénies, mais ne permet pas toujours d'accéder à une compréhension approfondie des modalités et des causes de l'évolution moléculaire.
Pour ces raisons, en collaboration avec un étudiant de Montréal (Nicolas Rodrigue), nous avons construit des modèles d'évolution des protéines explicitement basés sur la contrainte conformationnelle. L'idée est de comprendre quelle proportion des substitutions observées dans les séquences protéiques peuvent s'expliquer suivant un mode quasi-neutre, c'est-à-dire, comme étant le résultat d'une sélection stabilisante dominée par la contrainte conformationnelle. De tels modèles posent une double difficulté : (1) être capable de conduire des calculs de vraisemblance en prenant en compte les interdépendances entre positions induites par les contacts entre résidus dans l'espace tri-dimensionnel ; (2) correctement modéliser la thermodynamique conformationnelle des protéines.
L'essentiel du travail de Nicolas Rodrigue a consisté en la construction de méthodes de Monte Carlo à même de résoudre le premier point (Rodrigue et al, 2005, 2007). Par ailleurs, en vue d'affronter le second point avec les outils adéquats, une série de méthodes statistiques de comparaisons de modèles, basées sur l'évaluation du facteur de Bayes, ont également été implémentées (Rodrigue et al, 2006). En collaboration avec deux autres étudiants, de Montpellier (Cécile Bonnard) et de Montréal (Claudia Kleinman), des développements méthodologiques ont également été menés, afin de construire de meilleures potentiels statistiques, mimant l'énergie libre conformationnelle des protéines (Kleinman et al, 2006). Tous ces développements permettent maintenant d'effectuer un test systématique de différentes formes de potentiels thermodynamiques, dans le cadre Monte Carlo écrit par Nicolas Rodrigue, et d'évaluer ainsi leur pertinence dans la description de l'évolution des protéines réelles.
L'approche Bayésienne en phylogénie est attractive notamment grâce à sa flexibilité, et donc sa capacité à implémenter des modèles réalistes, donc complexes, de l'évolution des séquences, comme illustré ci-dessus. Cet intérêt d'ordre fondamental se heurte toutefois à des difficultés d'ordre technique. Toute chaîne de Markov bien formée de type Métropolis-Hastings converge sur le long terme, mais on n'a typiquement aucune garantie sur la vitesse de convergence. Par ailleurs, la comparaison de modèles (donc d'hypothèses biologiques) est moins naturelle en bayésien qu'en maximum de vraisemblance. Des efforts algorithmiques et informatiques importants ont dû être mis en ¿uvre pour traduire les modèles (équations) nouvellement proposées en méthodes (programmes) utilisables. Ces efforts théoriques sont réutilisables pour d'autres projets, et ont été valorisés par des publications (Lartillot 2006, Lartillot & Philippe 2006, Rodrigue et al 2007), ainsi que par des implémentations mises à disposition de la communauté des phlogénéticiens (cf infra).
Le calcul de vraisemblance ou bayésien en phylogénie moléculaire intègre sur l'ensemble des états pris à chaque site aux noeuds internes de l'arbre (états ancestraux). De ce fait, cette approche ne permet pas (directement) d'estimer les séquences ancestrales, ou la nature et la position dans l'arbre des changements de nucléotides (d'amino-acides). Reconstruire l'histoire des substitutions est pourtant un objectif biologique d'intérêt, en particulier en vue de mieux comprendre les forces évolutives qui sous-tendent l'évolution des macro-molécules. Ziheng Yang, notamment, a montré comment reconstruire les séquences ancestrales dans le cadre de l'approche statistique en phylogénie (Yang et al. 1995 Genetics 14:1641). La probabilité postérieure d'un état X donné à un noeud interne N donné pour un site S donné s'écrit comme le rapport de la probabilité de S conditionnellement à N=X par la proabilité inconditionnelle de S - on parle d'estimation "bayésienne empirique".
Dans le cadre de MODEL_PHYLO, nous avons étendu ces aspects théoriques au cas de la cartographie de substitutions: nous fournissons un estimateur du nombre (et de la nature) des substitutions de nucléotides ou d'amino-acides qui sont intervenues dans chacune des branches de l'arbre (donné) et pour chacun des sites. Nous montrons que dans ce contexte, et contrairement au cas des reconstructions de séquences ancestrales, les reconstructions marginales (où chaque branche de l'arbre est considérée séparément) et conjointes (où toutes les branches sont considérées simultanément) sont identiques, et nous fournissons une procédurede calcul rapide et flexible prenant en compte les spécificités du processus évolutif de la molécule analysée, et l'incertitude sur les états ancestraux.
Cet outil a été appliqué au problème de la détection de groupes de sites évoluant de manière corrélée (Galtier & Dutheil 2007), et à la caractérisation des contraintes biochimiques site-spécifiques. Pour la première de ces deux applications, nous définissons la proximité de deux sites comme le coefficient de corrélation entre leurs deux vecteurs de substitution - vecteurs récapitulant le nombre estimé de changements s'étant produits dans chaque branche de l'arbre, pour un site donné (Dutheil et al 2005). Puis les sites sont groupés sur la base de leur proximité par une méthode de clustering, et la significativité de ces groupes est évaluée par une procédure se simulation tenant compte de l'effet tests multiples (Dutheil & Galtier soumis). Le deuxième objectif est atteint en exhibant, grâce à la cartographie de substitution, des sites qui conservent significativement, au cours de leur évolution, une propriété biochimique donnée (par exemple la charge, ou le volume, ou la polarité). Ces sites sont détectés par un résidu élevé après régression de la vitesse d'évolution totale (toutes substitutions) sur le vitesse d'évolution biochimique (substitutions non-conservatives, Dutheil, soumis).
Le groupe dirigé par M. Gouy est depuis longtemps impliqué dans le développement de banques de données de séquences, et plus particulièrement de banques intégrant des informations phylogénétiques (alignements et arbres). Ces banques (HOVERGEN, HOGENOM et HOMOLENS) constituent une source de donnée de choix pour de nombreuses analyses évolutives. Cependant, un des principaux problèmes posés pour leur maintenance est lié à l'énorme quantité de séquences disponibles, quantité entraînant des temps de calcul extrêmement importants lors leur mise à jour. Par ailleurs, des outils de classification automatique sont également requis, en particulier dans le cadre de projets de séquençage de génomes complets.
Dans ce contexte, nous avons développé HoSeqI (Homologous Sequence Identification), un système permettant d'automatiser l'identification de séquences dans de grandes banques de familles de gènes homologues (Arigon et al. 2006). HoSeqI propose une interface qui permet à un utilisateur d'identifier rapidement une nouvelle séquence dans une banque de familles de gènes homologues choisie, c'est-à-dire de trouver la famille à laquelle appartient la séquence analysée, puis, de visualiser l'alignement et l'arbre phylogénétique obtenus. Il est ainsi facile de localiser la nouvelle séquence dans l'arbre de la famille de gènes homologues.
À partir des algorithmes développés pour HoSeqI, une autre application, MultiHoSeqI, a été implémentée afin d'ajouter rapidement les séquences d'un ou de plusieurs génomes aux familles présentes dans les banques de familles homologues. Le test en vraie grandeur de cette application a été effectué pour deux génomes de bactéries du genre Frankia, dont tous les gènes ont été rajoutés à la banque HOGENOM de cette façon. Enfin, l'application ChiSeqI, a également été développée sur la base des algorithmes utilisés avec HoSeqI (Arigon et al. 2007). Ce programme est dédié à la détection automatique de séquences chimères et à l'identification taxonomique basée sur l'utilisation de séquences d'ARNr 16S. Cet outil permet d'identifier l'espèce la plus proche de l'espèce de provenance d'une séquence d'ARNr 16S, ceci après avoir vérifié qu'il ne s'agissait pas d'une chimère.
Toutes ces méthodes ont pour objectif de reconstruire l'histoire des molécules, et donc potentiellement celle des génomes/espèces qui les hébergent, mais aussi de mieux comprendre les mécanismes de l'évolution moléculaire. Ces méthodes ont donc été appliquées à nombre de jeux de données, tout d'abord dans un souci de validation, puis pour répondre aux questions biologiques spécifiques posées dans MODEL_PHYLO: phylogénies anciennes et histoire de la thermophilie.
Les méthodes développées ont été confrontées à des jeux de données réel pour validation. Ainsi, lors d'une analyse de génomes mitochondriaux d'arthropodes, le modèle CAT-BP permet de replacer l'abeille, à la composition en bases (et en amino-acides) extrême, au sein des insectes, alors que les modèles classiques la positionne faussement comme groupe frère des tiques (qui présentent le même biais compositionnel, Blanquart & Lartillot, soumis).
Les méthodes de détection de la coévolution entre sites ont été appliquées à plusieurs alignements d'ARN et protéines. L'analyse de 80 ARN risomiques bactériens à permis de valider la méthode au sens où elle détecte une bonne proportion des paires Watson-Crick (bien documentées) qui stabilisent les hélices (Dutheil et al 2005). L'analyse de jeux de données protéiques a révélé un pattern de coévolution très variable selon les molécules. L'hémoglobine des vertébrés, par exemple, fournit très peu de signal, alors que la protéine bactérienne MAP, ou encore le gène SRK des brassicaceae (plantes), amènent à la détection d'un nombre importants de groupes de sites coévoluant, de taille et propriétés biochimiques variable. La majorité des groupes détectés font sens sur le plan structural, comme l'a révélé l'examen de leur position sur la structure tridimensionnelle des protéines. Une sophistication de la méthode permet de distinguer la coévolution liée à de la compensation vraie (un changement petit®gros compénsé par un gros®petit à un site voisin, par exemple) d'autres modes de coévolution moins intuitifs, telle que la corélation des vitesses d'évolution (Dutheil & Galtier, soumis).
Le modèle CAT a été utilisé pour une analyse phylogénomique à l'échelle des métazoaires (Brinkamnn et al 2005, Lartillot et al 2007). Une comparaison rigoureuses des différents modèles disponibles a permis de montrer que CAT s'ajuste aux données de manière optimale, par rapport au modèle à une matrice WAG, et qu'il diminue significativement les artefacts méthodologiques. Le modèle CAT reconnaît qu'un site peut osciller fréquemment entre peu d'états (par exemple: E->D->E->N->E->D->N), alors que le modèle WAG considère ce scénario comme très improbable, tous les sites étant censés suivre le processus "moyen". De ce fait le modèle WAG soutiendra une topologie où les réversions de ce type sont rares, alors que CAT réalise qu'il s'agit de sites à évolution rapide et à pouvoir résolutif faible, et leur donne correctement un poids implicite moindre. Cette innovation méthodologique a permis de démontrer la monophylie des Ecdysozoaires (animaux à mue: arthtropodes et nématodes) et des Lophotrochozoaires (Annélides, Mollusques, Plathelmintes), alors que les méthodes classiques positionnent typiquement les groupes à évolution rapide (Nématodes, Plathelmintes) comme groupes frères, à la base des Protostomiens. Il s'agit d'un pas important vers la résolution de la phylogénie des animaux - un des derniers groupes d'intérêt qui résiste encore aux données moléculaires.
Quelques dizaines d'espèces procaryotes actuelles connues vivent à très haute température (>80°, et jusqu'à 120°). Il a été proposé que ce caractère thermophile pourrait être ancestral, la vie étant apparue à très haute température (Woese 1987). La compréhension de l'évolution de la thermophilie est donc une question importante dans le contexte de la reconstruction des premières étapes de l'histoire de la vie. Par chance, le caractère thermophile est bien corrélé avec le taux de G+C dans les tiges de ARN ribosomiques - un taux élevé est nécessaire pour la stabilité à haute temprétaure. Grâce au programme rapide de reconstruction de phylogénies sous modèle non-stationnaire NHPhyML, Boussau & Gouy (2006) ont pu confirmer l'étude de Galtier et al (1999 Science 283:220), et montrer que l'ancêtre commun universel, quelle que soit la position de la racine, n'était probablement pas un organisme hyperthermophile, puisque son taux de GC ribosomique est modéré. Calteau et al (2005) ont exploré l'importance des transferts latéraux de gènes dans l'évolution de la thermophilie, et confirmé que des échanges génétiques avaient eu lieu entre archaea et bacteria (les deux grands groupes procaryotes) - il est plausible d'imaginer que l'hyperthermophilie ait été "inventée" par les archaea, et transmise secondairement vers certaines bacteria par transfert de gènes. Enfin, Boussau et al (soumis) montrent que les deux phylums hyperthermophiles bactériens Aquificales et Thermotogales sont probablement groupes frères dans la phylogénie des bactéries, et auraient donc possiblement une origine thermophile unique.