Approche génomique : étude de l'évolution de la cascade de régulation de la compétence chez les Streptocoques


Introduction

La transformation génétique naturelle est un mécanisme de transfert horizontal de gènes largement répandu chez les bactéries et qui nécessite la mise en place d'une machinerie d'entrée de l'ADN. Pour réaliser la transformation les bactéries doivent être dans un état physiologique particulier : la compétence. C'est un état génétiquement programmé et transitoire. Ce processus est étudié dans un certain nombre de batéries modèles, notamment chez  Streptococcus pneumoniae et Streptococcus mutans.  Chez Streptococcus pneumoniae, il a été démontré que l'état de compétence est régulé par un système à deux composants, ComDE, ComD étant l'histidine kinase (HK) et ComE le régulateur de réponse (RR). La détection d'un petit peptide (CSP, competence stimulating peptide) par le senseur ComD, conduit à son auto-phosphorylation et au transfert du groupement phosphoryl au régulateur de réponse ComE. ComE phosphorylé active directement l'expression des gènes précoces de la compétence, dont le facteur sigma alternatif ComX qui active la transcription des gènes de compétences dits tardifs. Cette réponse est très rapide car la transcription des gènes précoces atteind un maximum 5-7 minutes après l'adjonction de CSP dans le milieu. Les résulats obtenus chez  Streptococcus pneumoniae ont servi de modèle pour l'étude de ce processus chez S. mutans. Cependant dans cet organisme, un délai de 2h est observé entre l'adjonction de CSP et la transcription des gènes précoces.

Pour tenter de comprendre cette différence, une étude va être entreprise pour étudier l'évolution de cette régulation dans les génomes de Streptocoques et déterminer jusqu'à quel point le "modèle Streptococcus pneumoniae" peut être utilisé comme paradigme pour les autres espèces de Sterptocoques.

Dans un premier temps, nous allons analyser les 16 génomes complets de Streptocoques disponibles dans les bases de données. Pour simplifier l'analyse, nous n'avons retenu qu'un seul représentant de chaque espèce. Le gène codant pour le facteur sigma alternatif ComX est présent dans l'ensemble de ces génomes, ce qui suggère une conservation de la  régulation des gènes tardifs de la compétence. Afin d'étudier la régulation des gènes précoces, nous allons analyser la distribution des homologues au système ComDE dans ces génomes ainsi que leurs relations phylogénétiques.


Jeux de données

Le programme BlastP a été utilisé contre les protéines déduites des gènes identifiés sur les génomes complets avec comme séquences sondes les protéines ComE et ComD de S. pneumoniae. Nous avons utilisé une valeur seuil de e-value de 1e-05. Nous avons réalisé des alignements et des arbres afin d'éliminer les séquences partielles ou très divergentes des séquences ComE ou ComD. Les séquences retenues ont été extraites et sont collectées dans deux fichiers au format FASTA :  ComD et  ComE. Nous n'avons pas retenu de séquences dans les génomes de S. agalactiae et S. suis. Les informations associées à ces séquences sont compilées dans le tableau 1 (format Excel). Ce tableau est trié par groupes taxonomiques, par nom d'espèce et par position des gènes sur le chromosome. En l'absence de données expérimentales, nous prédisons que les partenaires d'un même système (HK et RR) sont codés par des gènes voisins sur le chromosome. Afin de faciliter l'étude de la coévolution des partenaires des systèmes, nous avons exclu  des  fichiers FASTA  les systèmes incomplets.
En plus du système, Com, les régulons Blp et Fas ont été étudiés expérimentalement.
  • Blp est impliqué dans la production de bactériocines. Un petit peptide inducteur (BlP) est codé par le gène blpC. Il est probablement maturé et exporté par le  transporteur ABC BlpAB. En réponse au BlP, l'histidine kinase BlpH s'auto-phosphoryle et transmet son phosphate au régulateur de réponse BlpR.
  • Fas (fibronectin/fibrinogen binding/haemolytic activity/streptokinase regulator) a été identifié chez Streptococcus pyogenes. Le système se compose d'un régulateur de réponse (FasA), de deux histidines kinases (FasB et FasC) et probablement d'un petit ARN non transcrit. Ce système ne serait pas impliqué dans la régulation par quorum sensing mais dans la régulation de l'expression de facteurs de virulences chez les Streptococcus pathogènes du groupe A (GAS).
Les séquences de Staphylococcus epidermidis et Staphylococcus aureus seront utilisées pour enraciner les arbres.

Analysez le tableau 1. Quelles observations pouvez-vous faire?

Alignement multiples des séquences homologues a ComE de S. pneumoniae

Réaliser l'alignement multiple des séquences homologues à ComE avec le programme muscle au travers du logiciel seaview.

Ouvrir votre fichier séquences dans seaview (menu File).

Dans le menu Align sélectionner Align all (vérifier le programme sélectionné dans Alignment options), une nouvelle fenêtre s'ouvre avec la progression de l'alignement, cliquer sur OK pour accepter l'alignement dans l'éditeur.

Vérifier la qualité de l'alignement et effectuer des corrections manuellement si nécessaire. Puis sauvegarder l'alignement au format  Fasta. 

Attention au nom de fichier, pas d'espace, d'accentuation etc. Conseil : ComE_muscle.fst;

Construction des arbres en utilisant la méthode de distance BioNJ

Distance Poisson

Dans le menu Trees, sélectionner Distance Methods, choisir la méthode BioNJ et la distance Poisson. Vérifier que l'option ignore all gap sites est cochée. Vous pouvez demander une analyse Bootstrap avec 100 itérations.

L'arbre calculé s'affiche dans un éditeur d'arbre. Vous pouvez choisir entre plusieurs représentations et afficher ou non les longueurs de branches et les valeurs de bootstrap. 

  • Re-root permet d'enraciner l'arbre à l'aide d'un sous-arbre ou d'une séquence (selectionner le noeud correspondant).
  • Swap permet de faire tourner autour d'un noeud, les deux sous-arbres définis par ce noeud.

Pour les autres options se référer au menu Help.

Sauvegarder votre arbre avec le menu File, Save rooted tree. Vous pouvez également utiliser l'option Save to Trees menu pour conserver l'arbre avec le fichier alignement (il apparaîtra dans le menu Trees).

Format Newick

Ouvrir le fichier arbre avec un éditeur de texte. Vous découvrirez comment les arbres sont codés en informatique, le format Newick qui est utlisé par la plupart des programmes de phylogénie et interprété pour vous fournir une visualisation graphique d'arbre. Les groupes de séquences sont séparés par des parenthèses, les différents groupes par des virgules, le point virgule indique la fin de l'arbre. Les longueurs de branches sont indiquées par des valeurs précédées de deux points et les valeurs de bootstrap sont indiquées juste après la parenthèse fermante délimitant les groupes.

Distance Kimura

Calculer l'arbre avec la distance Kimura

Editer les arbres (Swap) de façon à ce que l'ordre des feuilles sur les topologies soient les plus proches possibles.

Comparez les deux topologies, en particulier les longueurs des branches, les valeurs de bootstraps et les incongruences?

Comparaison des bipartitions des arbres obtenues avec les distances de Poisson et Kimura

Nous allons rechercher les bipartitions communes à un ensemble d'arbre. Une bipartition est obtenue par la coupure d'une branche qui sépare les feuilles en deux sous-ensembles.

Lancer le logiciel Rstudio. Ensuite dans le menu File, choisir New File et R Markdown.

Ceci va vous permettre de conserver vos script dans un fichier que vous pourrez recharger dans R Studio la fois prochaine. La fenêtre de R Studio va être divisée en quatre. En haut, à gauche le fichier que vous venez d'ouvrir et dans lequel vous allez taper vos commandes R pour les conserver. En dessous, la console qui indiquera le déroulement de vos commandes et où vouspouvez aussi taper directement des commandes si vous ne voulez pas les sauvegarder. En haut, à droite, une partie environnement qui vous permet de voir les variables que vous avez créées dans votre espace R ainsi que leur valeur et en dessous, une partie où s'afficheront les graphiques.

Placez-vous dans le répertoire contenant les fichiers que vous avez sauvegardés. Deux façons de le faire. Soit, vous allez dans le menu Session puis dans Set Working Directory et vous parcourez votre répertoire jusqu'au dossier correspondant. Soit, dans le fichier R que vous venez de créer vous tapez la commande setwd(). Entre parenthèse doit figurer le chemin pour accéder au dossier contenant les fichiers sur lesquels vous allez travailler.  Pour connaître le chemin de votre dossier. Clic droit sur le dossier. Ouvrir Propriété. Vous trouverez l'emplacement. Copier/coller dans R. Attention modifier les \ en /. 

Dans la partie fichier, les commandes R doivent être tapez entre ```{r } et ```. Ceci va définir un Chunk qui pourra être évaluer en allant dans  Run. Run Current Chunk évaluera les commandes du Chunk dans lequel le curseur est placé. Cela permet d'éxécuter les commandes du Chunk dont le résultat de l'exécution apparaît dans la partie Console. S'il y a des erreurs dans votre commande, elles vous serontt indiquées. Il faudra alors corriger la commande dans la partie fichier et réévaluer celle-ci en faisant de nouveau Run Current Chunk.
Il est conseillé de découper le script en plusieurs Chunk pour les évaluer au fur et à mesure.
A la fn du TP, n'oublier pas d'enregistrer votre fichier (File, Save).

```{r }
library(ape)
setwd("C:/Users/fichant/Documents/Phylogenie/Phylogenomic/cours-TD/TD_phylo/2012")
```
Pour travailler en phylogénie, il va falloir charger une librairie spécifique, la librairie ape (cf ci-dessus)

La commande read.tree() permet de lire un arbre dans le format Newick et plot.phylo() permet de déssiner cet arbre. 

Lire l'arbre obtenu avec la distance de Poisson :

tP <- read.tree(file = "ComE_BioNJ_Poisson_edited.ph")
plot.phylo(tP, show.node.label = TRUE, edge.color = "blue", edge.width = 1, no.margin = FALSE, cex=0.8)
title('BioNJ Poisson')
axisPhylo(side = 1)

La variable tP contiendra les informations concernant votre arbre, avec le nom des feuilles, les longeurs de branches et les labels associés aux noeuds. Vous pouvez voir son contenu en allant dans la partie environnement et en cliquant sur la variable ou en tapant son nom dans la partie console. L'arbre apparaît dans la partie graphique.

Lire l'arbre obtenu avec la distance de Kimura :
tK <- read.tree(file = "ComE_BioNJ_Kimura_edited.ph")
plot.phylo(tK, show.node.label = TRUE, edge.color = "blue", edge.width = 1, no.margin = FALSE,
cex=0.8)
title(
'BioNJ Kimura')
axisPhylo(side = 1)

Calculer l'arbre consensus strict (les partitions communes à tous les arbres)
tC = consensus(tP, tK, p=1);
plot.phylo(tC,
cex=0.8, main='Consensus strict between Poisson and Kimura');

Congruence des arbres obtenues avec les distances de Poisson et Kimura

Nous allons utiliser une méthode graphique qui permet de dessiner les deux arbres en vis-à-vis et de relier leur feuilles.

ntips = length(tP$tip.label) # nombre de feuilles dans l'arbre

Creer un tableau associant les feuilles à relier entre les deux arbres (ici elles sont toutes identiques)

association <- matrix(ncol=2, nrow=ntips)
association[,1]<-  tP$tip.label
association[,1] <- sort(association[,1])
association[,2] <- tK$tip.label
association[,2] <- sort(association[,2])
cophyloplot(tP, tK, assoc=association, show.tip.label=FALSE, col = "red", use.edge.length=TRUE)

Pour améliorer le rendu, nous allons utiliser une version modifiée des fonctions mycophyloplot.R et myplotCophylo2.R. Sauvegardez ces fonctions dans votre répertoire courant et les lire dans R avec les commandes source("mycophyloplot.R") et source("myplotCophylo2.R").
Puis utilisez la commande :
mycophyloplot(tP, tK, assoc=association, show.tip.label=FALSE, col = "red", use.edge.length=TRUE)
layout(matrix(1:2, 1, 2))
#partitionne la fenêtre graphique (ici en 2  1 ligne et 2 colonnes)
plot.phylo(tC, cex=0.6)
title ('consensus strict')
mycophyloplot(tP, tK, assoc=association, show.tip.label=FALSE, use.edge.length=TRUE)
title('Poisson versus Kimura')
layout(1)

Vous pouvez sauvegader votre graphique. Faire Export et choisir pdf. On peut avoir une prévisualisation.

Relation entre les distances d'arbres Poisson/ Kimura

La méthode NJ permet de reconstruire un arbre à partir d'une matrice de distances. Afin de comparer les distances inférées sur les arbres, nous allons extraire et  les comparer pour les deux méthodes de distances.

mDisP <- cophenetic(tP)
mDisK <- cophenetic(tK)

Les matrices sont obtenues dans l'ordre des feuilles des arbres, or celui peut varier entre les deux arbres. Comme nous voulons calculer une corrélation entre les paires de distances calculées avec les deux méthodes, nous devons tout d'abord trier les matrices de distances mDisP et mDisK sur les lignes et les colonnes par rapport à l'ordre alphabétique du nom des séquences. Par cette procédure, seules les valeurs des matrices ont été réordonnées mais pas les noms des lignes et des colonnes.

mDisProw <-  mDisP[order(rownames(mDisP)),]
mDisProwcol <- mDisProw[,order(colnames(mDisProw))]
mDisKrow <-  mDisK[order(rownames(mDisK)),]
mDisKrowcol <- mDisKrow[,order(colnames(mDisKrow))]

Comme ces matrices sont symétriques par rapport à la diagonale et que les valeurs de la diagonale sont nulles, nous allons utiliser la fonction as.dist() pour extraire la moitié inferieure gauche de la matrice de distance et éviter ainsi de prendre en compte deux fois les mêmes paires de valeurs ainsi que les valeurs contenues sur la diagonale.

plot(as.dist(mDisProwcol),as.dist( mDisKrowcol), pch=20, xlab='Poisson', ylab='Kimura'); 
abline(0,1, col='red');

Construction des arbres en utilisant une méthode du maximum de vraisemblance

Utiliser la méthode PhyML implémentée dans seaview avec les paramètres par défaut.

Comparer le nouvel arbre avec ceux obtenus précédemment avec la méthode de distance. Voyez-vous des différences? Si oui, cela est-il étonnant ?

Relation entre les distances d'arbres NJ obtenu avec le modèle Kimura et PhyML obtenu avec le modèle LG

tK <- read.tree(file = "ComE_BioNJ_Kimura_edited.ph") #uniquement si vous n'avez plus la variable tK dans votre environnement
tLG <- read.tree(file = "ComE_PhyML_LG_edited.ph")

Retapez les 3 commandes suivantes si vous n'avez plus la matrice mDisK dans votre environnement. Sinon fait juste au dessus. Donc juste faire le traitement pour les distances entre les séquences dans l'arbre obtenu avec PhyML (tLG)

mDisK <- cophenetic(tK)
mDisKrow <-  mDisK[order(rownames(mDisK)),]
mDisKrowcol <- mDisKrow[,order(colnames(mDisKrow))]

mDisLG <- cophenetic(tLG)
mDisLGrow <-  mDisLG[order(rownames(mDisLG)),]
mDisLGrowcol <- mDisLGrow[,order(colnames(mDisLGrow))]
plot(as.dist(mDisKrowcol),as.dist( mDisLGrowcol), pch=20, xlab='Kim', ylab='LG', main='cophenetic distances')
abline(0,1,col="red")

Présentez les deux arbres en vis-à-vis

layout(matrix(1:2, 1, 2))
plot.phylo(tK, cex=0.6)
title('NJ Kim')

plot.phylo(tLG, cex=0.6, direction='leftwards')
title('PhyML LG')
layout(1)

Comparaison des bipartitions des trois arbres

tP <- read.tree(file = "ComE_BioNJ_Poisson_edited.ph")
tK <- read.tree(file = "ComE_BioNJ_Kimura_edited.ph")
tLG <- read.tree(file = "ComE_PhyML_LG_edited.ph")
tC <- consensus(tP, tK, tLG, p=1)
plot.phylo(tC, cex=0.8)
title('Consensus strict, Poisson, Kimura, PAM, and LG')

Coloriage de l'arbre en fonction de la classe des protéines Lecture du tableau en format csv

tableau <- read.table("Gene_description_CleanUp-1.csv", header=TRUE, sep=';')
Code couleur pour chaque classe:
codes <- matrix(c('COME', 'BLPR', 'FASA', 'blue', 'red', 'green'), 3,2)
Association de la classe de la protéine (colonne Classification) à une couleur en gardant l'ordre des feuilles dans l'arbre tC.
col <- codes[match(tableau[match(tC$tip, tableau$Gene_Name), 12], codes[,1]), 2]
plot.phylo(tC, cex=0.8, tip.color=col);
legend("bottomleft", legend=codes[,1], col=codes[,2], pch=20)

Coloriage de chacun des 3 arbres obtenus
layout(matrix(1:3,1,3))
col_tP <- codes[match(tableau[match(tP$tip, tableau$Gene_Name), 12], codes[,1]), 2]
plot.phylo(tP,cex=0.6, tip.color=col_tP)
title('BioNJ Poisson')
col_tK <- codes[match(tableau[match(tK$tip, tableau$Gene_Name), 12], codes[,1]), 2]
plot.phylo(tK,cex=0.6, tip.color=col_tK)
title('BioNJ Kimura')
col_tLG <- codes[match(tableau[match(tLG$tip, tableau$Gene_Name), 12], codes[,1]), 2]
plot.phylo(tLG,cex=0.6, tip.color=col_tLG)
title('PhyML LG')
layout(1)

Comparer et discuter les trois arbres.

Recherche du modèle évolutif le plus adapté à nos données

Nous allons utiliser le programme prottest3 qui permet de selectionner le meilleur modéle d'évolution pour un alignement de séquences d'acides aminés.  Le programme est disponible comme ressource WEB (prottest2_server) ou comme application que vous pouvez lancer sur votre PC. Nous allons utiliser la version locale que l'on peut obtenir en allant sur le site drive.google.com pour télécharger l'archive (dernière version, cliquez sur la flèche vers le bas dans le bandeau).

Décompresser l'archive à l'aide du programme 7-Zip (clic droit sur le nom du fichier, extraire ici: la première opération vous décompresse le fichier; on obtient le fichier avec l'extension .tar. recommencer la même opération sur ce nouveau fichier. Cela vous extrait vos fichier de l'archive et vous construit le répertoire prottest3). Aller dans le répertoire prottest3 et lancer la commande runXProtTestHPC.bat.

  • menu File selectionner l'alignement à analyser (format FASTA).
  • menu Analysis, choisir Compute likelihood scores,  sélectionner les matrices JTT (Jones et al., 1992), LG (Le and Gascuel, 2008), WAG (Whelan and Goldman, 2001) et Blosum62 (Henikoff and Henikoff, 1992) et l'estimation des vitesses de variations (+G) avec 4 catégories et des fréquences des acides aminés (+F). Puis lancer la recherche...
  • menu Selection, choisir Results. Vous avez un tableau avec les différents modèles classés par ordre décroissant par rapport à leur adéquation avec l'alignement multiple utilisé.
  • Utiliser Export to main console pour copier les résultats dans la console et les enregistrer dans un fichier.

Dans le tableau Results vous avez plusieurs calculs réalisés suivant différents critères, notamment l'AIC qui signifie Akaike Information Criterion (AIC).

Lorsque l'on estime un modèle statistique, il est possible d'augmenter la vraisemblance du modèle en ajoutant un paramètre. Le critère d'information d'Akaike permet de pénaliser les modèles en fonction du nombre de paramètres afin de satisfaire le critère de parcimonie. On choisit alors le modèle avec le critère d'information d'Akaike le plus faible.

Dans notre contexte, c'est un estimateur qui correspond à la minimisation de la distance attendue entre un modèle vrai et son estimation. Les modèles correspondant aux valeurs minimales de l'AIC sont considérés comme les plus appropriés pour la reconstruction. Une même topologie de référence doit être utilisée pour tester les différents modèles.

Le critère d'information d'Akaike s'écrit comme la différence entre 2 fois le nombre de paramètres (k) deux fois la log-vraisemblance du modèle estimé.


formule AIC
                                                                                     k = nombre de paramètres libres du modèle

Alignement multiples des séquences homologues a ComD de S. pneumoniae

Réaliser l'alignements multiples des séquences homologues à ComD avec le programme muscle au travers du logiciel seaview comme vous l'avez fait pour ComE.

Vérifier la qualité de l'alignements et effectuer des corrections manuellement si nécessaire. Puis sauvegarder l'alignement.

Recherche du modèle évolutif le plus adapté à nos données

Utiliser prottest3 pour trouver le modèle le mieux adapté à ces séquences.

Construction de l'arbre en utilisant une méthode du maximum de vraisemblance

Utiliser la méthode PhyML implémentée dans seaview en modifiant si nécesaire les paramètres correspondant à la matrice et à la correction gamma. Réorganiser votre arbre si nécessaire (groupe externe, swap des feuilles par rapport à l'arbre obtenu avec ComE). Sauvegarder le en format Newick.

Comparaison des deux arbres obtenus sur ComE et ComD

tComE <- read.tree(file = "ComE_PhyML_LG_edited_ComD.ph")
tComD <- read.tree(file = "ComD_PhyML_LG_edited_ComE.ph")

codes2 <- matrix(c('COMD', 'BLPH', 'FASB', 'FASC', 'blue', 'red', 'green', 'lightgreen'), 4,2)
col2 <- codes2[match(tableau[match(tComD$tip, tableau$Gene_Name), 12], codes2[,1]), 2]

layout(matrix(1:2, 1,2))
opar <- par()
par(xpd=TRUE, mar=c(8,4,4,3))
plot.phylo(tComE, cex=0.5, tip.color=col_tLG)
legend(-0.5,-0.5, legend=codes[,1], col=codes[,2], pch=20, cex=0.55)
title('ComE')

plot.phylo(tComD, cex=0.5, direction='leftwards', tip.color=col2)
legend(1.5,-0.5, legend=codes2[,1], col=codes2[,2], pch=20, cex=0.55)
title('ComD')
layout(1)
par(opar)

Analyser les arbres. Conclusions.

Annotation des arbres avec iTOL (Interactive Tree Of Life)

Interactive Tree Of Life : http://itol.embl.de/

Dans la barre de menu en haut de la page, cliquez sur Upload.

Vous pouvez copier/coller votre arbre ou charger le fichier contenant l'arbre en format Newick (d'autres formats sont disponibles).

L’arbre est tracé dans la fenêtre principale. Vous pouvez le déplacer en cliquant dessus et en déplaçant votre souris. A gauche, vous avez des icônes pour zoomer en avant et en arrière.  Opération également possible avec la roulette de la souris.

Modification de l’arbre

Les fonctions principales d’iTOL sont rassemblées dans le panneau de contrôle. Elles sont réparties en quatre catégories (Basic, Advanced, Datasets et Export).

Par défaut et pour des raisons esthétiques, iTOL place les branches avec le plus faible nombre de feuille en haut de l’arbre (mode rectangulaire). Il est possible d’annuler cette option avec Advanced->Leaf sorting.

Vous pouvez modifier l’aspect de l’arbre avec Basic->Display mode. Vous pouvez ajouter les valeurs de bootstraps avec Advanced->Bootstrap.

Export vous permet de générer une image de votre arbre et de son annotation.

Un clique gauche sur les branches de l’arbre, vous donne accès à des fonctions s’appliquant aux nœuds et aux branches. Il est par exemple possible de ré-enraciner l’arbre à la branche sélectionnée (Tree structure->Reroot the tree here).

Annotation des arbres

iTOL offre de nombreuses possibilités d’annotations des arbres. Chaque opération est décrite par un fichier texte, dont le nom doit avoir le suffixe ‘.txt’. L’annotation se fait par glissement (drag-dropping) du fichier dans la fenêtre principale d’iTOL ou par le chargement du fichier à l’aide du bouton + localisé en bas et à droite de la fenêtre.

Le format des fichiers d’annotations est le suivant :

Une première ligne avec un mot clé correspondant au type d’annotation :
TREE_COLORS

La définition du type de séparateur utilisé pour délimiter les informations (TAB, SPACE ou COMMA).
SEPARATOR TAB

Un ensemble de paramètres globaux qui vont définir l’aspect global de l’annotation.
Le mot clé DATA indique le début des annotations. 

Annotation de ComE

 Exemple 1 (RR_color_class_range_cover.txt) : coloriage des sous arbres correspondant aux familles de gènes

Chaque sous arbre est définit par le dernier ancêtre commun des feuilles correspondantes,  à l’aide de deux labels de feuilles séparés par ‘|’. Plusieurs combinaisons de pairs de feuilles donnent le même résultat. Le mot range définit le type de coloriage à réaliser (l’arrière-plan du sous arbre).

TREE_COLORS
SEPARATOR TAB
DATA
SaurA01.AGRA|SepiA01.AAO05237.1   range       #80f31f    AGRA
SgorA01.COME|SpneA01.COME           range       #3541fe   COME
SubeA01.CAR41219.1|SsalA01.RR09     range       #ca039c   BLPR
SsalA01.CCB96178.1|SpyoA01.AAK33322.1         range       #078eec   FASA

Après chargement du fichier, il est possible de modifier le type de coloriage avec Cover dans la fenêtre associée à cette annotation (label, clade, full, ou off). En cliquant sur la roue crantée, vous pouvez modifier la couleur de chaque clade (Datasets dans la fenêtre de contrôle).

Exemple 2 (RR_replace_tips_by_Code.txt) : remplacement du nom des feuilles.
A chaque feuille est associé un nouveau nom.

LABELS
SEPARATOR TAB
DATA
SaurA01.AGRA      SaurRR1
SepiA01.AAO05237.1             SepiRR1
SequA01.ACG62860.1            SequRR2

Exemple3 (RR_protein_domains.txt) : annotation des domaines Pfam.
La longueur de la protéine est indiquée et à la suite, chaque domaine est décrit par une forme (HH), un début et une fin dans la séquence, une couleur et un label.

DATASET_DOMAINS
SEPARATOR
TAB
DATASET_LABEL Domains_Pfam
COLOR   #0000ff
DATASET_SCALE 100-amino_acid_100-#0000ff  500-line at aa500-#ff0000        1000-aa1000-#00ff00
DATA
SaurA01.AGRA      238          HH|3|121|#7300ff|pfam00072  RE|149|237|#0000ff|pfam04397
SepiA01.AAO05237.1             238          HH|3|121|#7300ff|pfam00072              RE|149|238|#0000ff|pfam04397
StheA01.ABJ66794.1           245          HH|34|120|#7300ff|pfam00072 RE|145|237|#0000ff|pfam04397
SubeA01.CAR41219.1            246          HH|3|127|#7300ff|pfam00072               RE|146|233|#0000ff|pfam04397
SubeA01.FASA       245          HH|3|114|#7300ff|pfam00072  RE|152|242|#0000ff|pfam04397

Annotation de ComD

Exemple 4 (HK_color_class_strips.txt) : ajout d'une colonne avec les familles des protéines

Exemple 5 (HK_replace_tips_by_Code.txt) : remplacement du nom des feuilles par un code.

Exemple 6 (HK_protein_domains.txt) : annotation des domaines Pfam.

Construction d'un arbre espèces en utilisant un ensemble de protéines

Une des approches utilisées pour reconstruire l'arbre phylogénétique d'un ensemble d'espèces consiste à concaténer les alignements obtenus sur un ensemble de gènes ou de protéines conservés dans les génomes de ces espèces (méthode supermatrice). Nous allons illustrer cette approche avec les protéines des petites et grandes sous-unités des ribosomes. Nous avons utilisé les profils COGs de la base de données CD pour annoter les génomes complets de streptocoques à l'aide du programme RPS-BLAST. Nous avons retenu un échantillon de 43 familles de protéines (tableau 2) qui étaient annotées correctement dans la grande majorité des  17 génomes complets de streptocoques. Les séquences du génome de Lactococcus lactis seront utilisées pour enraciner l'arbre.
Pour chaque famille, les séquences des protéines ont été extraites (Ribo_proteins).

Illustration

A titre d'illustration, nous allons utiliser seaview pour aligner avec muscle les 5 premiers fichiers et nous allons calculer les arbres phylogénétiques sur ces fichiers.
Pour pouvoir réaliser facilement l'étape de concaténation, pensez à ouvrir vos fichiers à partir de la première fenêtre seaview. Dans tous les cas sauvegarder vos alignements.

Que pensez-vous des alignements obtenus?

Que pensez-vous des arbres obtenus?

Recherchons les bipartitions communes à tous ces arbres

t1 <- read.tree(file = "Tree1.ph")
t2 <- read.tree(file = "Tree2.ph")
t3 <- read.tree(file = "Tree3.ph")
t4 <- read.tree(file = "Tree4.ph")
t5 <- read.tree(file = "Tree5.ph")
tC = consensus(t1, t2, t3, t4,t5, p=1)
Qu'avez vous obtenus?

Supprimer l'arbre qui pose des problèmes et recommencer l'opération.

plot.phylo(tC, cex=0.8, main = 'Consensus strict')

Nous allons maintenant concaténer les 5 premiers alignements. Pour cela, placez-vous dans la fenêtre seaview contenant le premier alignement et dans le menu File choisissez Concatenate.

On selectionne successivement les 4 autres alignements. Pensez à cocher la case add gaps  pour que les alignements avec un nombre de séquences plus petit soient correctement entrés.

Que pensez-vous de l'alignement obtenu?

Nous allons utiliser le pgrogramme Gblocks pour éliminer automatiquement les positions avec une faible ou aucune information phylogénétique. Dans le menu Sites choisir  Create set, puis Gblocks. Dans le menu Gblocks cocher Allow gap positions within final blocks.
Relancer le calcul de l'arbre phylogénétique.


Alignement des 43 familles de protéines

Dans le fichier (super matrix) vous trouverez la concaténation des alignements des 43 familles de protéines.

Calculer un arbre en utilisant la méthode BioNJ et la matrice de distance Kimura (100 bootstraps).

Annotation de l'arbre des espèces avec iTOL

Exemple 7 (AE_replace_tips_by_species.txt) : remplacement du nom des feuilles par l'espèce

Exemple 8 (AE_color_taxonomy_range.txt) : coloriage des sous arbres correspondant rang taxonomique.

Exemple 9 (AE_family_distribution.txt) : distribution phylogénétique des familles.

Construction de la phylogénie des Streptocoques à partir des séquences de l'ARN 16S

Nous allons réaliser cette nouvelle phylogénie à l'aide des séquences des ARNr 16S.  Nous allons utiliser comme ressource la Ribosomal Database Project II (RDP).  Plusieurs outils d'analyses vous sont proposés. Utiliser le  Hierarchy browser (lien Browsers)  pour sélectionner les séquences d'ARNr 16S des streptococcus. Sélectionner seulement les séquences types.

Combien de séquences obtenez_vous?

Cocher toutes les séquences pour les sélectionner (cliquer dans le +).  Ensuite faire download pour qu'elles soient mises dans votre seqCart. On peut choisir en plus de sauvegarder ces séquences sur son ordinateur. Pour cela, choisir le format Fasta et garder l'option Remove common gaps sinon nous obtenons un alignement très long plein de gaps. Ouvrir le fichier avec Seaview (format Fasta). Vous allez voir que vos séquences sont identifiées suivant leur accession number dans la RDP et que vous avez donc perdu le nom de l'espèce. Il va donc falloir les renommer. Comme cela est fastidieux, vous trouverez ici le fichier à utiliser. 
  • Construction de l'arbre.
Vous pouvez soit construire l'arbre directement dans Seaview avec la distance HKY (la plus proche de Tajima&Nei), ou ouvrir ce dernier dans Seaview et  sauvegarder en format mase comme à l'exercice 1et réaliser l'arbre phylogénétique avec phylo_win en utilisant la distance de Tajima&Nei et la Neigbor Joining method. Les calculs étant assez rapides vous pouvez tenter 50 bootstrap. Votre topologie est-elle bien résolue?
Nous allons suppimer les séquences d'ARNr 16S qui sont redondantes pour certains organismes. Refaire l'arbre. Etes-vous satisfait de votre topologie.

  • Recherche du modèle évolutif le plus adapté à nos données.
Etant donné le grand nombre de modèles évolutifs disponibles pour le traitement des séquences d'acides nucléiques, des méthodes ont été développées pour permettre de choisir le modèle le plus adapté aux données. Elles utilisent PhyML comme méthode de reconstruction d'arbre. Deux tests seront utilisés :

      • Le rapport de vraisemblance ou LRT (Likelihood Ratio Test)
Deux modèles sont comparés (modèle correponsdant à l'hypothèse nulle et le modèle alternatifs) et le rapport des vraisemblances est calculé.

Rapport de vraisemblance

Le rapport de vraisemblance est d'autant plus grand que la vraisemblance du modèle alternatif est grande. Le modèle nul sera rejeté si R est supérieur au seuil fixé par l'utilisateur. Pour que le LRT est un sens il faut que les modèles soient imbriqués car parcours d'un arbre de décision et testé sur une même topologie de référence.

Exemple d'un arbre de décision (extrait du manuel de JModeltest).

 exemple arbre de décision
      • L'Akaike Information Criterion (AIC)
C'est une alernative au LRT qui présente l'avantage de pouvoir être appliquée à des hypothèses non imbriquées. C'est un estimateur qui correspond à la minimisation de la distance attendue entre un modèle vrai et son estimation. Les modèles correspondant aux valeurs minimales de l'AIC sont considérés comme les plus appropriés pour la reconstruction. Une même topologie de référence doit être utilisée pour tester les différents modèles.

formule AIC
k = nombre de paramètres libres du modèle

 Nous allons donc rechercher le modèle le plus approprié correspondant à nos données sur ComD pour réaliser notre arbre. Aller sur le site de Phylemon2. Créer un compte ou se connecter en tant qu'utilisateur anonyme. Rechercher l'application JMoldelTest. Copier/coller votre alignement. Pour un nombre de type de substitution égale à 3, utiliser les deux méthodes AIC et LRT. Pour cette dernière on effectuera l'option " Do dynamical likelihood ratio tests".
Quel modèle est le plus adéquate?
A quoi correspond ce modèle ? (voir le manuel de JModelTest dans le Help)
Est-il implémenté dans PhyML?
Quel est le second modèle le plus adapté?
Noter l'impact de la correction par la gamma distribution sur la vraisemblance du modèle.

Reconstruire l'arbre en choisissant le modèle proposé dans PhyML et en choisissant Best of NNI & SPR pour la recherche de l'arbre.
Comparer avec le précédent. La topologie est-elle bien résolue?

Nous allons donc essayer d'améliorer celle-ci en recherchant de façon exhaustive le modèle le plus approprié. Pour cela sur Phylemon2 rechercher l'option Phyml Best AIC Tree (v. 1.02b). Sauveragder l'arbre obtenu dans un fichier et importer le ensuite dans seaview.
Persiste-t-il encore des problèmes de résolution des branches?