Découverte :
Métagénomique et Graphique de Recrutement : chercher l'aiguille dans la botte de foin !

Le massacre des Albigeois

La croisade contre les Albigeois, vers 1208. Manuscrit tiré des "Chronique de France ou de Saint Denis", 1325-1350. Source : Wikimedia Commons.

Un intérêt certain se développe pour la compréhension d'écosystèmes complexes comme les eaux, les sols ou encore les microbiomes (plus de 90% des cellules du corps humain sont en fait les bactéries qui peuplent son tube digestif). Plusieurs difficultés se posent :

  1. la plupart des organismes présents ne sont pas facilement cultivables et donc individualisables.
  2. tous les organismes ne sont pas présents en quantité égale et, par extension, d'importance égale dans le fonctionnement de l'écosystème, sans corrélation stricte entre abondance et importance.
  3. certains de ces écosystèmes sont extrêmement riches et des centaines d'espèces peuvent se cacher dans quelques grammes de terreau ou quelques millilitres d'océan.

Face à une tâche aussi ingrate que passionnante, les chercheurs ont trouvé une solution dans la célèbre phrase associée à la Croisade des Albigeois : «Tuez-les tous, Dieu reconnaîtra les siens», en l'adaptant légèrement :

«Séquencez-les tous, un doctorant fera le tri»

Ainsi naquît la Métagénomique et de grands projets comme le séquençage du rumen des vaches ou du microbiome humain, le projet METASOIL, le Global Ocean Survey et plus récemment le projet TARA.

Il convient de distinguer deux possibilités en métagénomique :

  1. séquencer un marqueur génétique capable d'identifier les organismes présents à l'échelle du genre voire de l'espèce comme le marqueur V9 (correspondant à l'ADNr 18S pour les eucaryotes, 16S pour les procaryotes)
  2. abuser des techniques modernes de séquençage très haut débit : extraire la totalité des acides nucléiques de l'échantillon, tout découper, tout séquencer. Pour l'analogie, imaginons une virée à la Bibliothèque Nationale, un découpage intensif de l'ensemble des pages de chacun des livres (dont certains sont -au grand bonheur des nombreux fans d'Harry Potter- en plusieurs copies), la réalisation d'un gros tas relativement stable, et pour finir se lancer, avec espoir, témérité et un peu de ruban adhésif, dans la reconstitution de tous ces chefs d’œuvre (si si, même ceux qu'on n'aime pas, comme 50 Shades of Grey ou la saga Twillight).

Fort heureusement, rares sont les chercheurs qui osent se lancer en quête d'une compréhension exhaustive des populations d'un écosystème complexe. Si la quête se révèle relativement réalisable lorsqu'on ne séquence qu'un marqueur universel et qu'il s'agit d'identifier "qui est là", tenter de reconstituer les génomes via une étape d'assemblage se révèle un défi. Bien souvent et c'est ce que nous allons aborder ici, le chercheur lambda ne s'intéressera qu'à détecter l'abondance de son picoeucaryote adoré ou de ses bactéries favorites. Pour cela un outil astucieux et très pratique existe : le Graphique de Recrutement (librement traduit du Recruitment Plot que vous trouverez dans la littérature)

L'idée étant simplement de recruter -par alignement- des reads de séquençage du métagénome sur un ou plusieurs génomes de référence. Et de distinguer ces reads par leur pourcentage d'identité avec la séquence cible. Lorsqu'une quantité significative de reads s'alignent sur le génome d'une souche de référence, et ce avec un pourcentage d'identité suffisamment important, on peut émettre l'hypothèse que cet organisme ou une variété d'organismes très proches sont présent dans l'échantillon considéré.

Ingrédients de la recette du jour : Le Graphique de Recrutement

  1. Quelques dizaines de millions de "reads" d'une taille honnête (pas moins de 100 paires de bases), au mieux une banque Illumina "Pair-End" de 2 x 100 pb chevauchants voire de longs reads 454 ou IonTorrent.
  2.  Quelques génomes de référence, représentatifs de la diversité du groupe dont vous souhaitez détecter la présence
  3.  Un gros, gros ordinateur (ou bien mieux un cluster de calcul)
  4. Quelques compétences en "lançage de BLAST" et dans un langage de programmation "User Friendly" (type Bash, Python ou Perl) ainsi qu'un peu de R pour faire de beaux graphiques, mais Excel reste utilisable pour les frileux
  5. Beaucoup d'enthousiasme !

Difficulté : 3 / 5     Rien d'extravagant ici, mais nécessite une certaine souplesse devant une ligne de commande

Méthode en 3 étapes :

1/ Aligner l'ensemble des reads du métagénome sur l'ensemble des génomes "de référence" à l'aide de BLAST, en autorisant des "hits" aussi faibles que 55% d'identité

L'idée étant de trouver les paramètres appropriés :

blastall -p blastn -i myReads.fasta -d myDatabase -o outfile.tab -m 8 -G 8 -E 6 -r 5 -q -4 -W 8 -e 1 -F "m L" -U T

Avec -W 8 : la taille du "mot", -G 8 -E 6 -r 5 -q 4 : les coûts d'ouverture et d'élongation de gap, la récompense pour un match et la pénalité pour un mismatch, -e 1 -F "m L" -U T : quelques filtres pour alléger les fichiers de sortie, et non moins important -m 8 pour une sortie tabulée facile à manipuler. Ces paramètres ne sont qu'une suggestion, il conviendra de les adapter à ses propres données.

2/ Nettoyer les résultats (à l'aide de Bash, Python ou Perl, ou autre...)

Il s'agit de conserver les hits BLAST qui vous semblent pertinents, non pas en termes de pourcentage d'identité, mais en terme de longueur de l'alignement. J'essaye par exemple d'avoir au moins 80% de la séquence de mon read qui soit effectivement alignée. À ce stade, vous voulez savoir si un organisme de référence donné est présent dans le métagénome. Il ne s'agit pas de faire de compétition entre vos organismes de référence, je vous suggère donc de conserver également l'information de reads qui s'alignent sur plusieurs de vos génomes (en l'occurence, considérant que vos génomes sont proches, cela devrait se produire régulièrement). Par contre, n'hésitez pas à ne conserver que le meilleur hit de chaque read pour chacune de  vos souches de référence.

Il s'agit ensuite d'extraire trois informations clefs : le génome de référence où le read s'est aligné, le pourcentage d'identité du hit et la localisation de ce hit sur le génome de référence (pas de pitié, faites la moyenne du Stop et du Start).

3/ Construction du Graphique de Recrutement

Il s'agit enfin, pour un organisme de référence donné de générer le graphique suivant : en abscisse la position sur le chromosome de votre organisme (s'il y en a plusieurs, il faudra faire plusieurs graphes ou bien les concaténer), et en ordonnée, pour chacun des reads son pourcentage d'identité. Petit conseil : il n'est pas utile d'afficher 2 millions de points sur un graphiques, sélectionnez tout au plus 100 000 points à afficher.

Résultats :

Graphique de Recrutement pour la Souche A (G.Farrant, Tous droits réservés)

Graphique de Recrutement pour la Souche B (G.Farrant, Tous droits réservés)

Dans ces exemples, j'ai fusionné les profils de recrutement des métagénomes de trois stations de prélèvement sur les génomes de deux souches bactériennes A et B. Cette combinaison de 3 conditions, associées à des couleurs différentes facilite la comparaison des stations. L'histogramme à droite du graphique représente l'abondance de points sur des tranches de pourcentage d'identité.

On observe que quasiment aucun read ne s’aligne sur la souche A avec une identité supérieure à 90% alors que sur la souche B, une quantité significative de reads s’aligne avec plus de 90% d’identité, ce qui permet de conclure que des organismes génétiquement proches de la souche B sont présents dans l’échantillon étudié. Les points recrutés entre 55 et 80% d'identité correspondraient à du recrutement non-spécifique.

On observe également quelques artefacts potentiellement intéressants tels des sites où le signal est très riche, l'intuition amène à penser qu'il s'agit de régions contenant des gènes de ménage largement conservés entre organismes (type ARN ribosomiques) ainsi que des régions "vides", dont une très marquée pour la Souche 2, qui doivent être interprétés (îlots génomiques ?).

Le Graphique de Recrutement est un premier outil d'analyse de la présence d'organismes de référence dans un métagénome. Il permet entre autres de définir un seuil de pourcentage d'identité prompt à discriminer les familles d'organismes. Il permet également de révéler des particularités riches d'interprétations. Pour plus d'infos et du subtilités, je vous suggère la publication du GlobalOceanSurvey (GOS) qui s'octroie la paternité de la méthode.

Je tiens à remercier Malicia, Hautbit, Yoann M., Nico M., Mica et Guillaume Collet pour leurs commentaires avisés.

  • À propos de
  • Post-doc pour l'institut Matís de Reykjavík (Islande) sur la biodiversité des lacs sous-glaciaires et des sources géothermales (projet MarieCurie AstroLakes). Doctorat de l'UPMC pour une thèse en écologie des picocyanobactéries marine à la Station Biologique de Roscoff, Un profil pltôt axé bioanalyse/bioinformatique, sur des données de génomique et métagénomique (assemblage, biodiversité, miTAGs...) Sinon... joueur de go (à mes heures perdues), choriste, amateur de bouquins et de musique, de randonnées islandaises... je milite pour un nouveau monde (sur une base extrêmiste modérée).

7 commentaires sur “Métagénomique et Graphique de Recrutement : chercher l'aiguille dans la botte de foin !

  1. Pour minimiser le nettoyage, j\'aurai effectuer un lançage de blast avec une valeur de -e bien inférieure à 1 😛

    Sinon, les sorties graphiques des graphes de recrutement sont de toute beauté. Possibilité de partage du code?

    • Oui, je suis assez d\'accord sur le principe, mais en vérité, je ne voulais pas prendre le risque de louper quelque chose et n\'étant pas limité par l\'espace disque, j\'ai préféré faire du nettoyage post-BLAST.

      En ce qui concerne les sorties graphiques, je ne suis pas en mesure de fournir le code, mais il n\'y a clairement rien d\'ésotérique vu que j\'ai pu le faire tout seul.

  2. Cette représentation est vraiment très intéressante, elle pourrait permettre d\'obtenir un visuel comparatif plus complet des génomes connus dans ce type de jeux de données qui sont quand même assez conséquents et \"fouillis\" à appréhender de prime abord.

    Quelques questions (je n\'ai pas encore lu la publi donc pardonne-moi si ça y est mentionné) :
    - Lorsqu\'on observe les 3 profils de 3 sites différents sur un même génome et un même graphe, inclues-tu une étape de normalisation de la quantité des reads ?
    - 90% d\'identité pour valider une souche me parait peu, comment as-tu fixé ton seuil ? (je serais plus encline à ne valider que l\'espèce avec un tel pourcentage)

    Je sens que je pourrais m\'inspirer de tout ceci pour une représentation complémentaire à certaines de mes analyses, merci en tout cas de l\'avoir mise en avant, il va falloir que je creuse !

    PS : Petite note de pinailleuse, il serait judicieux de mettre à jour la ligne de commande blast avec la nouvelle nomenclature des dernières versions (blastn -query, -db, -out, -word_size, -evalue, ...) (je sais, passer de l\'un à l\'autre est une gymnastique bien pénible, je me plante encore souvent...)

    • Tout d\'abord, merci pour ton intérêt.

      Alors non, il n\'y pas d\'étape de normalisation pour la conception des Graphiques de Recrutement car je me contente d\'afficher les 100 000 premiers \'hits\' de chaque stations. J\'ai pris le temps de voir les résultat en affichant la totalité de mes points, c\'est passablement illisible, et en échantillonnant les 100 000 points de façon aléatoire, le graphique est relativement stable à l’échantillonnage.
      La comparaison fine des stations est difficile sur le Graphique de recrutement, personnellement, je m\'en sers plus comme d\'une base pour détecter les sus-nommés artefacts et comparer leur présence dans les différentes stations.
      Comparer les stations de façon quantitative est difficile car on ne connaît pas systématiquement le volume ou la masse de l\'échantillon de métagénome. Par contre, on peut -et on doit- faire du semi-quantitatif et comparer les profils de distribution entre les stations. Typiquement, on peut faire du recrutement sur les génomes de 36 souches appartenant à 17 clades différents, il est possible, en dénombrant les meilleurs hits BLAST pour chacun des clades d\'en estimer la proportion dans le métagénome, et ce profil semi-quantitatif est quand à lui comparable entre les stations.

      Concernant la question du seuil, il n\'y pas de bonne réponse. Une stratégie possible consiste à établir arbitrairement un seuil, puis extraire les séquences correspondantes, les \'blaster\' sur nr et voir le ratio de séquences qui ont un \'best hit\' sur l\'organisme de référence, adapter le seuil, et recommencer, mais nous sommes d\'accord -je l\'espère- pour penser que LA souche de référence n\'existe que dans la littérature et que le milieu étudié contient au mieux un spectre relativement hétérogène d\'organismes proches. Finalement on en revient à la question de définition d\'espèce.

      Ravi que ça t\'ai inspiré. (Pour le pinaillage, je dirais \"You\'re raising an interesting concern\" ^^)

  3. Bonjour,

    Tout d\'abord merci beaucoup pour cette article ainsi que pour l\'ensemble de votre site qui est très intéressant.

    Je me permets de solliciter votre aide car j\'aimerais effectuer une manipulation similaire sur mes échantillons, mais hélas je débute en métagénomique.
    J\'aurais voulu savoir par quel moyen vous décidiez des valeurs pour les différents paramètres du blast. Ensuite auriez vous plus de renseignements sur la manière d\'effectuer le nettoyage sur les résultats du blast ? J\'ai essayé de trouver de la documentation ou des tutoriels, mais je ne trouve rien d\'accessible ou de très précis.

    Par avance merci beaucoup pour votre aide !

    • Bonjour,
      Désolé de n\'avoir pris le temps de répondre plus tôt.
      Ravi que cela vous plaise =)

      Bien, alors concernant BLAST, je pense qu\'il faudrait éditer un erratum. En fait, à l\'époque du début de ma thèse, j\'utilisais BLASTALL (une version plus ancienne, mais pas obsolète de BLAST. Avec les paramètres par défaut, BLASTALL et la dernière version du logiciel (appelée \"BLAST+\") se comporte comme MEGABLAST et génère des alignements génétiquement proches compris globalement entre 80 et 100% d\'identité, et donc inutilisables pour la construction desdits \"recruitment plots\", ces derniers se basent sur des alignements entre 50 et 100%.
      Les paramètres \"sensibles\" de blast présentés dans cet article permettent de forcer MEGABLAST à générer des alignements distants jusqu\'à ~50% d\'identité.

      En fait, j\'ignorais à l\'époque de la rédaction l\'existence d\'une option de BLAST+: -task blastn qui lui permet de se comporter comme un BLASTN normal. L\'utilisation des options \"sensibles\" de BLAST permettent de générer plus d\'alignements.

      J\'essaierai de faire un article plus détaillé sur la question.

      • Merci pour ces informations complémentaires !

Laisser un commentaire