M1MABS BGPG Projets
From silico.biotoul.fr
Contents |
Projet 1 : Détermination du génome coeur
Nous avons vu en cours que la taille et le nombre de gènes entre génomes de différentes souches d'une même espèce bactérienne pouvaient varier. Le génome coeur (ou core genome) a été défini comme correspondant à l'ensemble des gènes communs à toutes les souches d'une même espèce (ensemble des gènes orthologues).
Comment pratiquement et rapidement peut-on identifier les gènes orthologues entre deux espèces ?
La technique la plus couramment utilisée consiste à analyser les résultats d'un programme appelé Blast qui permet de comparer rapidement une séquence donnée (séquence requête ou query) à un ensemble de séquences constituant une banque. Le résultat du Blast est un alignement local indiquant les régions présentant des similarités entre les deux séquences. Cette similarité est quantifiée par le pourcentage d'identité trouvé entre les deux séquences alignées, le score de l'alignement (appelé bit score) et la Evalue qui correspond au nombre attendu d’alignements qui par chance aurait un score ≥ au score obtenu avec nos deux "vraies" séquences. Ce nombre est évalué en regardant les scores des alignements générés par comparaison de séquences aléatoires ayant même longueur et même composition que la séquence requête. Plus la Evalue est faible plus l'alignement obtenu est significatif.
Pour rechercher les gènes orthologues, la comparaison se fait souvent au niveau de leurs produits protéiques. Nous utiliserons donc ici, par abus de langage, les termes d'orthologues et de paralogues pour des séquences protéiques, ceci pour simplifier l'écriture.
Dans les fichiers joints, vous trouverez donc le résultat de la comparaison de chaque protéine d'un génome (chaque séquence sera à son tour séquence requête) avec l'ensemble des protéines d'un autre génome constituant donc la banque de séquences. Le programme utilisé pour la comparaison est BlastP (version du blast pour comparer une séquence protéique à une banque protéique).
Nomenclature utilisée et exemple
Code à 5 lettres se référant au nom de l'espèce suivi d'une cinquième lettre assignée à la souche : EcolA: Ecol pour Escherichia coli et la dernière lettre se référant à la souche (A pour K12 MG1655, B pour O157:H7, D pour CFT073, L pour K12 DH10B, M pour SECEC SMS-3-5, R pour IAII, U pour UMN026 et X pour BL21)
Pour plus d'informations sur les souches : https://www-abcdb.biotoul.fr/#/entry/displaytaxonomytree/Taxon_Name/EcolA (remplacer EcolA dans l'URL par le génome qui vous intéresse).
EcolB_EcolD.blt : chaque séquence du génome de la souche EcolB01 sera une séquence query qui sera comparée à l'ensemble des séquences protéiques de la souche EcolD01. Comme résultats on ne gardera que le premier hit (donc la séquence de EcolD01 qui possède la plus petite Evalue). On parlera alors de Best Hit.
Si on prend l'exemple de la séquence EcolB01.ECS002 son best hit dans EcolD01 est la séquence EcolD01.THRA. On notera BH(EcolB01.ECS002, EcolD01.THRA) la Evalue est de 0 et le score de 1590.
EcolD_EcolB.blt : chaque séquence du génome de la souche EcolD01 sera une séquence query qui sera comparée à l'ensemble des séquences protéiques de la souche EcolB0101. Comme résultats on ne gardera que le premier hit.
Si on prend l'exemple de la séquence EcolD01.THRA, son best hit est la séquence EcolB01.ECS002. On a donc BH(EcolD01.THRA,EcolB01.ECS002)
Dans cet exemple on a donc BH(EcolB01.ECS002, EcolD01.THRA) = BH(EcolD01.THRA,EcolB01.ECS002), les deux séquences sont réciproquement best hit l'une de l'autre. Dans ce cas on parle de Best Bidirectional Hits appelé dans la littérature BBH.
La première condition pour que deux protéines soient orthologues est qu'elles soient BBH. Ensuite, pour que cette relation soit 1:1, il ne faut pas qu'une des deux séquences possède une séquence paralogue dont le score d'alignement serait supérieur à celui du BBH.
Nous devons donc vérifier qu'il n'y a pas eu de duplication du gène dans un des deux génomes après leur divergence (qui serait donc plus récente que la divergence) et qui résulterait dans la présence d'une séquence paralogue (dans le même génome) ayant un score de d'alignement de Blast supérieur à celui obtenu entre les séquences BBH de nos deux génomes. Pour cela, pour chaque génome, chacune de ces séquences protéiques a été comparée à l'ensemble des autres séquences protéiques et le premier et donc meilleur hit obtenu est appelé Best Paralog, noté BP.
Les résultats de ces comparaisons sont trouvés dans les fichiers portant deux fois le nom de la souche.
Si nous reprenons l'exemple précédant :
EcolB_EcolB.blt : chaque séquence protéique de la souche EcolB01 est comparée à toutes les autres protéines de cette souche.
Le best hit de la séquence EcolB01.ECS002 est EcolB01.ECS4869, après le premier hit trivial qui est celui de la séquence requête avec elle-même. Ces deux séquences sont donc best paralog et le score de l'alignement obtenu est 341.
EcolD_EcolD.blt : chaque séquence protéique de la souche EcolD01 est comparée à toutes les autres protéines de cette souche. Le best hit de la séquence EcolD01.THRA est EcolD01.METL. Ces deux séquences sont donc best paralog et le score de l'alignement obtenu est 344.
Les deux scores BP(EcolB01.ECS002,EcolB01.ECS4869) = 341 et BP(EcolD01.THRA,EcolD01.METL) = 344 < BBH((EcolB01.ECS002, EcolD01.THRA) = 1590
Il n'existe donc pas de séquences paralogues dans aucun des deux génomes possédant un score d'alignement supérieur au BBH. Donc nous pouvons conclure que EcolB01.ECS002 et EcolD01.THRA sont deux séquences orthologues 1:1.
Données fournies
8 souches ont été analysées : EcolA, EcolB, EcolD, EcolL,EcolM, EcolR, EcolU et EcolX.
Pour chacune d'entre elle, il y a 8 fichiers : Comparaison de chaque protéine de la souche avec les autres protéines de la même souche (fichier ou le nom de la souche apparaît deux fois), comparaison de chaque protéine de la souche avec l'ensemble des protéines des 7 autres souches (dans le nom du fichier, le nom de la première souche est celle dont les séquences sont utilisées comme query et le deuxième est le nom de la souche dont les séquences ont servi à constituer la banque).
Dans chaque cas seul le best hit est reporté.
Le fichier contient le résultat du BlastP sous forme de tableau avec les colonnes correspondant aux données suivantes:
- nom séquence query ,
- nom séquence banque,
- % d'identité,
- longueur de l'alignement,
- nombre de substitutions,
- nombre de gaps (indels),
- position début alignement séquence query,
- position fin alignement séquence query,
- position début alignement séquence banque,
- position fin alignement séquence banque,
- Evalue,
- bit score.
Pour trouver le génome coeur il va donc falloir comparer ces 8 fichiers et trouver les protéines orthologues dans les 8 génomes. Dans l'étape précédente, on a trouvé que les deux séquences EcolB01.ECS002 et EcolD01.THRA étaient orthologues. Maintenant, il va falloir rechercher le BBH de chacune d'entre-elles dans les autres génomes, vérifier que la séquence identifiée comme BBH dans le nouveau génome que l'on parte de EcolB01.ECS002 ou de EcolD01.THRA est bien la même et qu'il n'y a pas de paralogues dans ce nouveau génome de plus fort score Blast.
Exemple si on regarde la souche EcolL.
On a BBH(EcolB01.ECS002,EcolL01.THRA) = 1595 et on a BBH(EcolD01.THRA,EcolL01.THRA). Le best paralog est BP(EcolL01.THRA, EcolL01.METL) = 344. Donc dans les trois souches pas de paralogues de score plus élevé. Donc ici les 3 protéines sont orthologues.
On fait de même avec les cinq génomes restants. Si on a trouvé une protéine orthologue dans chacun des génomes, alors la protéine appartient au génome coeur.
Travail demandé
Faire le programme qui permet d'extraire le génome coeur en appliquant les principes énoncés ci-dessus et donner le résultat obtenu sur les 7 génomes analysés, c'est à dire le nombre et la liste des protéines constituant ce génome coeur.
Rendre un rapport décrivant votre démarche et intégrant l'algorithme du programme. Le code du programme (commenté) sera envoyé séparément.
Projet 2 : Création et analyse des profils phylogénétiques
Pour rappel : un profil phylogénétique décrit la présence ou l'absence d'un gène (d'une protéine) dans un ensemble de génomes comparés à un génome de référence. Ces gènes (ou protéines) doivent présenter une relation d'orthologie avec le gène (la protéine) du génome de référence.
Comment pratiquement et rapidement peut-on identifier les gènes orthologues entre deux espèces ?
La technique la plus couramment utilisée consiste à analyser les résultats d'un programme appelé Blast qui permet de comparer rapidement une séquence donnée (séquence requête ou query) à un ensemble de séquences constituant une banque. Le résultat du Blast est un alignement local indiquant les régions présentant des similarités entre les deux séquences. Cette similarité est quantifiée par le pourcentage d'identité trouvé entre les deux séquences alignées, le score de l'alignement (appelé bit score) et la Evalue qui correspond au nombre attendu d’alignements qui par chance aurait un score ≥ au score obtenu avec nos deux "vraies" séquences. Ce nombre est évalué en regardant les scores des alignements générés par comparaison de séquences aléatoires ayant même longueur et même composition que la séquence requête. Plus la Evalue est faible plus l'alignement obtenu est significatif.
Pour rechercher les gènes orthologues, la comparaison se fait souvent au niveau de leurs produits protéiques. Nous utiliserons donc ici, par abus de langage, les termes d'orthologues et de paralogues pour des séquences protéiques, ceci pour simplifier l'écriture.
Dans les fichiers joints, vous trouverez donc le résultat de la comparaison de chaque protéine d'un génome (chaque séquence sera à son tour séquence requête) avec l'ensemble des protéines d'un autre génome constituant donc la banque de séquences. Le programme utilisé pour la comparaison est BlastP (version du blast pour comparer une séquence protéique à une banque protéique).
Nomenclature utilisée et exemple
Code à 5 lettres se référant au nom de l'espèce suivi d'une cinquième lettre assignée à la souche : EcolA: Ecol pour Escherichia coli et la dernière lettre se référant à la souche (A pour K12 MG1655, B pour O157:H7, D pour CFT073, L pour K12 DH10B, M pour SECEC SMS-3-5, R pour IAII, U pour UMN026 et X pour BL21)
EcolB_EcolD.blt : chaque séquence du génome de la souche EcolB01 sera une séquence query qui sera comparée à l'ensemble des séquences protéiques de la souche EcolD01. Comme résultats on ne gardera que le premier hit (donc la séquence de EcolD01 qui possède la plus petite Evalue). On parlera alors de Best Hit.
Si on prend l'exemple de la séquence EcolB01.ECS002 son best hit dans EcolD01 est la séquence EcolD01.THRA. On notera BH(EcolB01.ECS002, EcolD01.THRA) la Evalue est de 0 et le score de 1590.
EcolD_EcolB.blt : chaque séquence du génome de la souche EcolD01 sera une séquence query qui sera comparée à l'ensemble des séquences protéiques de la souche EcolB0101. Comme résultats on ne gardera que le premier hit.
Si on prend l'exemple de la séquence EcolD01.THRA, son best hit est la séquence EcolB01.ECS002. On a donc BH(EcolD01.THRA,EcolB01.ECS002)
Dans cet exemple on a donc BH(EcolB01.ECS002, EcolD01.THRA) = BH(EcolD01.THRA,EcolB01.ECS002), les deux séquences sont réciproquement best hit l'une de l'autre. Dans ce cas on parle de Best Bidirectional Hits appelé dans la littérature BBH.
La première condition pour que deux protéines soient orthologues est qu'elles soient BBH. Ensuite, pour que cette relation soit 1:1, il ne faut pas qu'une des deux séquences possède une séquence paralogue dont le score d'alignement serait supérieur à celui du BBH.
Nous devons donc vérifier qu'il n'y a pas eu de duplication du gène dans un des deux génomes après leur divergence (qui serait donc plus récente que la divergence) et qui résulterait dans la présence d'une séquence paralogue (dans le même génome) ayant un score de d'alignement de Blast supérieur à celui obtenu entre les séquences BBH de nos deux génomes. Pour cela, pour chaque génome, chacune de ces séquences protéiques a été comparée à l'ensemble des autres séquences protéiques et le premier et donc meilleur hit obtenu est appelé Best Paralog, noté BP.
Les résultats de ces comparaisons sont trouvés dans les fichiers portant deux fois le nom de la souche.
Si nous reprenons l'exemple précédant :
EcolB_EcolB.blt : chaque séquence protéique de la souche EcolB01 est comparée à toutes les autres protéines de cette souche.
Le best hit de la séquence EcolB01.ECS002 est EcolB01.ECS4869, après le premier hit trivial qui est celui de la séquence requête avec elle-même. Ces deux séquences sont donc best paralog et le score de l'alignement obtenu est 341.
EcolD_EcolD.blt : chaque séquence protéique de la souche EcolD01 est comparée à toutes les autres protéines de cette souche. Le best hit de la séquence EcolD01.THRA est EcolD01.METL. Ces deux séquences sont donc best paralog et le score de l'alignement obtenu est 344.
Les deux scores BP(EcolB01.ECS002,EcolB01.ECS4869) = 341 et BP(EcolD01.THRA,EcolD01.METL) = 344 < BBH((EcolB01.ECS002, EcolD01.THRA) = 1590
Il n'existe donc pas de séquences paralogues dans aucun des deux génomes possédant un score d'alignement supérieur au BBH. Donc nous pouvons conclure que EcolB01.ECS002 et EcolD01.THRA sont deux séquences orthologues 1:1.
Données fournies
271 génomes ont été sélectionnés pour constituer les profils.
Pour chacun d'entre eux, il y a 3 fichiers :
- comparaison de chaque protéine du génome avec les autres protéines du même génome (fichier ou le nom de génome apparaît deux fois),
- comparaison de chaque protéine du génome avec l'ensemble des protéines d'EcolA (dans le nom du fichier, le nom de la première souche est celle dont les séquences sont utilisées comme query et le deuxième est le nom de la souche dont les séquences ont servi à constituer la banque)
- comparaison de chaque protéine d'EcolA avec l'ensemble des protéines codés par le génome.
Le fichier contient le résultat du BlastP sous forme de tableau avec les colonnes correspondant aux données suivantes:
- nom séquence query,
- nom séquence banque (subject),
- % d'identité,
- longueur de l'alignement,
- nombre de substitutions,
- nombre de gaps (indels),
- position début alignement séquence query,
- position fin alignement séquence query,
- position début alignement séquence banque,
- position fin alignement séquence banque,
- Evalue,
- bit score.
Travail demandé
Par comparaison de chacun des autres génomes avec le génome de référence (EcolA), identifier les protéines orthologues 1:1. A partir de ces résultats construire les profils phylogénétiques en codant 1 la présence de la protéine orthologue et 0 son absence. A partir de ce tableau, construire une matrice de pseudo-distances entre les profils phylogénétiques en utilisant l'indice de Jaccard. On aura pris soin avant de construire la matrice de pseudo-distances d'éliminer les résultats triviaux (protéines orthologues présentes dans tous les génomes). Appliquer ensuite une méthode de classification automatique pour trouver les clusters de profils similaires. Analyser deux clusters pris au choix au niveau fonctionnel (les protéines appartiennent-elles a une même processus biologique ?).
Pour cela, il faudra programmer l'identification des protéines orthologues et la construction des profils phylogénétiques. Vous utiliserez ensuite R pour calculer la matrice de pseudo-distances. La classification automatique pourra être réalisée avec la méthode et les outils de votre choix. Pour l'analyse fonctionnelle des clusters sélectionnés vous pourrez vous aider de la recherche d'annotations sur-représentées (programme sook de mobyle disponible à l'adresse http://silico.biotoul.fr/cgi-bin/mobyle/portal.py). Pour le trouver aller dans annotation puis enrichment. Les annotations Gene Ontology ainsi que les mots-clés UniProt-SwissProt sont disponibles pour la souche EcolA (K-12/MG1655).
Rendre un rapport décrivant votre démarche et intégrant l'algorithme du(des) programme(s) ainsi que la description de vos résultats. Le code du(des) programme(s) sera envoyé séparément.
