Didacticiel :
Introduction à l'analyse des SNPs

Introduction

Un SNP (Single Nucleotid Polymorphism) est la mutation d'un seul nucléotide à une position donnée, entre individus d'une même population. Ces substitutions ponctuelles permettent d'identifier des sous-populations.

Vous avez sûrement étudié en cours l'exemple de la drépanocytose, une maladie causée par la mutation ponctuelle d'un gène. L'origine de cette maladie génétique peut être observée en analysant les SNPs dans la population humaine. Analyser ces derniers peut ainsi permettre de trouver des prédispositions à certaines maladies, mais peut également être utile entre autres, en génétique des populations ou en cancérologie. L'arrivée des séquenceurs de deuxième génération (les fameux NGS) a facilité leur étude.

Pour la suite, nous partirons du principe que vous êtes familiers avec les concepts liés à l'analyse de données NGS (reads, mapping, format SAM/BAM, fastq, profondeur de séquençage etc.) sinon la suite risque d'être un peu ardue.

L'analyse

Pour analyser vos SNPs, je vous propose le pipeline ci-dessous. Ce dernier a pour but d'identifier les SNPs présents dans vos échantillons puis de les filtrer. À la fin de l'article nous verrons comment les annoter par rapport à une base de données.

Pipeline d'analyse de SNPs

Pipeline d'analyse de SNPs

 

L'analyse commence après alignement de vos reads contre votre génome de référence via votre mapper préféré. À la suite de cette étape vous obtenez classiquement un fichier au format BAM. La première étape propre à l'analyse des SNPs consiste à utiliser la commande mpileup, fournie avec l'outil samtools. Elle s'utilise de la manière suivante :

D'autres options peuvent bien sûr être rajoutées (cf. la doc).

Si vous possédez plusieurs fichiers BAM (correspondant à plusieurs échantillons mappés contre le même génome), vous pouvez utiliser cette commande sur tous vos fichier d'un coup, il est d'ailleurs fortement conseillé de procéder ainsi. D'une part vous n'aurez qu'un fichier à gérer par la suite et d'autre part, l'algorithme de  mpileup utilise les reads de tous vos échantillons pour ses calculs. La commande  mpileup produit un fichier au format pileup (plus trop utilisé actuellement). Plus généralement, l'option -u est utilisée afin de produire un fichier binaire (format BCF). Ce dernier est le format compressé de VCF : Variant Call Format (en gros BCF est à VCF ce que BAM est à SAM).

 

Exemple d'un fichier au format vcf (source : https://en.wikipedia.org/wiki/Variant_Call_Format )

Exemple d'un fichier au format vcf (source : https://en.wikipedia.org/wiki/Variant_Call_Format )

 

Ci-dessus un exemple de fichier VCF. Ce format se compose de 3 parties :

  1.   Les premières lignes, commençant par "##", contiennent la description des champs. Dans l'exemple ci-dessus on peut lire "##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">" ce qui indique que le champ DP ne doit contenir qu'un seul chiffre ("Number=1") entier ("Type=Integer") et ce chiffre correspond à la profondeur du read ("Description="Reads Depth"").
  2. L'en-tête (header) : Le nom des colonnes, commençant par "#"
  3. Les données elles-mêmes contenant différents champs séparés par une tabulation.

 

Plus d'infos sur ce format sur la page wikipedia d'où provient la capture d'écran. La documentation officielle du format VCF est disponible ici.

Revenons à notre pipelinempileup ne fait qu' "entasser" les reads alignés sur votre génome (du verbe anglais to pileup). Afin de supprimer les faux positifs dus à des erreurs de séquençage (c'est-à-dire un mismatch entre un read et le génome de référence sans pour autant que cela soit un SNP), il est possible d'utiliser l'outil bcftools. Celui-ci vous permettra également d'obtenir un fichier au format VCF qui sera utile pour la suite.

Par exemple, pour convertir le fichier output.bcf  obtenu précédemment, il suffit d'écrire :

L'option -v sert à afficher uniquement les potentiels variants.

D'une manière générale  bcftools permet de manipuler les fichier VCF et BCF ( bcftools  est aux formats VCF et BCF ce que  samtools est aux formats SAM et BAM).

À noter :

La commande  view est faite pour fonctionner avec  samtools ≤ 0.1.19 (la version de  samtools disponible dans les paquets Ubuntu à l'heure actuelle). Si vous avez une version de samtools et  bcftools plus récente il faudra passer par la commande call (l'étape de détection des SNPs se nomme SNP calling).  Ainsi la commande sera :

Pourquoi -m ? Eh ben d'après le site officiel de bcftools :

 Users are now required to choose between the old samtools calling model (-c/--consensus-caller) and the new multiallelic calling model (-m/--multiallelic-caller). The multiallelic calling model is recommended for most tasks.

L'option -m est donc la nouvelle recommandation des auteurs. L'option -c est plutôt utilisée pour rechercher des sites bialléliques.

La différence entre  mpileup et  bcftools peut vous sembler floue. Pour résumer,  mpileup rassemble tous les reads à la même position sur le génome et bcftools détecte les SNPs de manière plus fine en supprimant les faux positifs. Dans la pratique, ces deux commandes sont utilisées conjointement via le symbole pipe ( |) sous Unix. C'est d'ailleurs cette façon de faire qui est indiquée dans la doc de mpileup.

Ainsi si vous ne deviez retenir qu'une seule commande, retenez celle-là :

 

Filtres

Maintenant que vous avez identifié vos SNPs, vous pouvez les compter, les regarder, mais vous pouvez aussi vous demander que faire (de plus) avec tous ces SNPs ? Comment ne garder que les plus "pertinents" ? (La difficulté ici est de définir "pertinent").

Vous pouvez filtrer vos SNPs selon :

  • Le nombre de fois qu'ils sont lus;
  • La qualité de l'alignement;
  • Le taux d’allèles à une position donnée

 

C'est par exemple ce qui est effectué dans Seo et. al (2012)  sur des SNPs humains (cf. paragraphe "Somatic single nucleotide and short-indel discovery").

Afin de vous aider à choisir vos filtres, vous pouvez vous poser les deux questions suivantes :

  • Avez-vous des contrôles positifs ? Si vous êtes sûrs que vous devez retrouver tel ou tel SNPs, est-ce que vous les retrouvez après avoir appliqué vos filtres ?;
  • Quel est le contexte de votre étude ? Est-ce que vous vous attendez à avoir "beaucoup" ou "peu" de SNPs ?

 

Pour filtrer vos SNPs, vous pouvez par exemple utiliser snpsift.

Utilisation :

Par exemple, dans la publication citée ci-dessus, les auteurs préconisent de ne garder que les SNPs lus avec une profondeur 3. Ce qui se traduit par la commande suivante :

 

Plus d'info sur http://snpeff.sourceforge.net/SnpSift.html#filter

Sachez que certains outils, tels que GATK, offrent des filtres plus fin pour analyser les SNPs. Par contre ce dernier est plus compliqué à utiliser : l'installation n'est pas simple et avant de commencer à analyser vos SNPs, prévoyez de nombreuses heures pour lire la doc.

L'annotation

Le pipeline présenté tout au long de cet article s'arrête après l'étape de filtrage. Ceci dit, une fois vos SNPs obtenus, vous pouvez également les comparer à des bases de données afin de les annoter.

SnpSift permet de le faire grâce à la commande annotate. L'idée est de regarder si vos SNPs sont présents dans la base de données choisie et si oui de leur rajouter certaines informations. Prenons par exemple la base de données dbSNP qui est une base de données de SNPs humains hébergée par le NCBI. Mettons que si l'un de vos SNPs est présent chez dbSNP alors vous voulez rajouter l'information à propos de la fréquence de l'allèle basée sur l'étude 1 000 genomes (cette information est contenue dans le champ CAF d'après le glossaire). Pour réaliser ceci vous pouvez utiliser la commande suivante :

Dans la commande ci-dessus 00-All.vcf correspond à la base de données.

Conclusion

Pour résumer, une analyse de SNPs commence toujours par un mapping, suivi d'un SNP calling (avec samtools et bcftools) afin d'obtenir un fichier VCF. Ce dernier peut éventuellement être filtré dans le but d'éliminer les allèles moins fréquents par exemple. Enfin, la partie annotation des SNPs est la plus compliquée car elle est dépendante des données et plus particulièrement du génome et du contexte de l'étude.

Comme évoqué dans l'introduction, les NGS ont grandement facilité l'étude des SNPs. La probabilité de détecter un SNP rare dépend du nombre d'individus séquencés. Ainsi, grâce au projet 100,000 genomes (séquençage de 100,000 génomes humains pour étudier des maladies rares) on peut s'attendre à découvrir de nouveaux SNPs chez l'humain. Pour finir, sachez que les SNPs sont le plus souvent étudiés dans le cadre des maladies génétiques mais rien ne vous empêche de les analyser dans d'autres domaines comme la phylogénie par exemple.

 

Merci à Julien Fouret pour son aide ainsi qu'à LelouarBérénice, et NiGoPol pour les commentaires et discussions lors de l’édition de cet article.

  • À propos de
  • Bordelais d'origine, je pratique le pierre-feuille-ciseaux en club.

Catégorie: Didacticiel | Tags: , , ,

17 commentaires sur “Introduction à l'analyse des SNPs

  1. Bonjour,

    Merci pour votre article, fort intéressant.

    Je tiens juste à indiquer que le mapping des reads contre un génome de référence n\'est pas la seul solution pour détecter des SNP. On peut même réaliser ces détections sans génome de référence comme avec DiscoSNP https://gatb.inria.fr/software/discosnp/ pour ne cité que celui que je connais le mieux.

  2. Salut,
    Décidément on en apprend tous les jours 🙂
    Moi ça me paraitrait plus rassurant de mapper mes reads contre mon génome de référence pour être sur que la majorité s\'aligne dessus, mais si le génome de référence n\'est pas dispo ça peut être un bon outil.

    • On peut aussi regarder du coté de MindTheGap https://github.com/GATB/MindTheGap (de la même équipe qui a crée DiscoSNP)qui lui utilise un génome de référence et des reads NGS, pour détecte les SNP, les insertion et les délétions.

  3. Bonjour,

    Existe t-il une méthode pour déterminer les SNPs sans génome de référence en mappant des génomes non des reads ? j\'ai trouvé un seul mais avec référence : https://harvest.readthedocs.io/en/latest/

  4. Mapper des génomes ?
    Sans génome de référence il y a disco SNP cf le 1er commentaire 😉

    • DiscoSNP mappe des reads, je voulais savoir s\'il existe un programme utilisant input = des génomes assemblés en fasta?

  5. À ma connaissance non. Ce que tu peux faire à la rigueur c\'est de l\'alignement multiple et repérer les SNPs \"à la main\" http://bioinfo-fr.net/alignements-multiples-quels-logiciels-choisir

    • Merci, mais je ne pense pas que ces outils d\'alignement seront bien adapté à de grandes données je vais peut être teste mafft.

      • Bonjour,

        Je serai toi je partirai d\'abord sur Clustal Omega qui est plus récent est mieux adapté aux gros jeux de données.
        Si ça ne marche pas, alors j\'envisagerai ensuite MAFFT ou/et MUSCLE.
        Bon courage.

  6. Salut à tous !

    Première fois que je poste sur le site 😛 Or ce n\'est pas faute de m\'intéresser à tout ce que vous faites (Merci ^_^ !) J\'attends impatiemment d\'avoir quelques vraies journées devant moi pour vous lire entièrement (ou en grande partie du moins !)

    Je suis étudiante en thèse de génétique des populations et je vais être amenée très prochainement à analyser des données brutes de RAD-Seq (et plus précisément de Pool-seq, aoutch 😉 !) Pour le moment j\'ai surtout travaillé avec Stacks, que je m\'étonne de ne pas voir listé ici (et sauf erreur, il n\'a pas non plus été évoqué dans un autre article du blog ?), y a t\'il une raison particulière à cela ?

    Sachant que de mon côté mon petit challenge personnel sera de filtrer mes SNP en sachant que je n\'ai qu\'un pool (donc un seul barcode) pour chacune de mes 18 espèces, et zéro génome de référence ! Je pensais partir sur une pipeline type Velvetoptimiser, BWA/Bowtie2, puis implémenter tout ça dans Popoolation2, un logiciel qui permet justement de faire du SNP calling sur du pool. Peut-être que certains ici ont fait quelque chose de similaire 🙂 ?

    En tout cas bravo encore, j\'essaierai de vous voir aux JOBIM puisque vous y serez 🙂 !

    Chrys

    • Salut Chrys,

      Merci pour ton commentaire, il fait chaud au coeur !
      Je ne vais pas pouvoir t\'aider concernant ta requête, ce n\'est absolument pas ma spécialité. Mais si tu n\'obtiens pas de réponse d\'ici JOBIM? je t\'encourage fortement à venir nous rencontrer, il y aura au moins une experte (Isabelle) capable de te répondre.

      Et puis tu l\'as dis : nous n\'avons pas encore d\'article sur ton sujet. Qui sait, peut-être que ça pourrait être toi qui écrira le premier sur ce sujet ici ? 🙂

      Pour info : on donnera une brève présentation lors du workshop JeBiF le lundi 27 juin et on aura ensuite un poster bioinfo-fr.net tout au long de la conf 🙂

      Au plaisir de discuter avec toi !

    • Salut,
      Pour ce qui est de l\'article ici l\'idée c\'était de faire une introduction aux analyses de SNPs, au format VCF et compagnie, pas du tout d\'être exhaustif et de dérouler la liste de tous les logiciels existants.
      Après comme l\'a dit Yoann, se focaliser sur quelques logiciels qui ont des particularités intéressantes, ça pourrait faire le sujet d\'un autre article 😉

      T\'as 18 espèces différentes et tous tes reads ont le même index ? C\'est...bizarre comme façon de faire, bon courage !

    • Salut!

      Le travail que tu vas effectuer est relativement similaire à ce que je fais actuellement en stage de Bioinfo à l\'INRA. J\'vais du coup apporter ma petite pierre à cet article en espérant que ça puisse (t\')aider quelques personnes.

      De mon côté je bosse également sur des données brutes de RAD-Seq et concrètement je dois tester plusieurs méthodes de détection de SNP, les comparer, voir les points positifs et négatifs de chacune et déterminer s\'il y en a pas une meilleure dans tout ça. Bien sûr ce que je vais citer s\'applique plus particulièrement à mon contexte, c-a-d l\'étude d\'une population d\'arbres (sans génome de référence / organisme modèle). Soit un nombre assez conséquent de séquences...

      J\'ai débuté mon travail à partir de 8 fichiers .fastq (8 individus séquencés). Ce que j\'ai testé:

      - Stacks. Effectivement, un bon outil lorsqu\'on a un paquet de données. L\'avantage de ce programme est la rapidité de traitement (1 nuit pour 20Go de données), même sans génome de référence. J\'ai pas encore vérifié les SNPs qu\'il trouve mais j\'ai amplifié par PCR un lot de loci qu\'il a déterminé et la plupart sont existants.

      - Pyrad. Très similaire à Stacks seulement celui-ci permet de travailler à partir de séquences de longueurs variables. Un avantage comme un inconvénient car il demande énormément de RAM... et de temps de calcul (1 semaine de traitement sur cluster :-|). Les résultats sont comparables à Stacks, plutôt positifs dans l\'ensemble.

      - CLC / (BWA-SAMTOOLS) / Reads2SNP. Avec ces 4 outils j\'ai testé une méthode différente des deux autres, c-a-d construire un génome de ref à partir de toutes les séquences (CLC) puis mapper ces séquences sur cette ref (BWA). L\'idée peut paraître séduisante mais les résultats le sont moins. La plupart des loci possédant des SNPs (détectés par R2S) de ne sont pas amplifiés lorsque j\'ai essayé de vérifier par PCR (quelque chose à voir avec la ref de CLC je suppose).

      - DiscoSNP. Probablement pas adapté aux données sur lesquelles je travaille. Je n\'ai pas pu exploiter les .vcf obtenus.

      - GATK. Pas adapté aux données RAD-Seq. Mais un must test.

      Les trois premières méthodes m\'ont donné +500 000 SNPs sur +50 000 loci (à partir de +80 millions de reads). A partir de là j\'ai écris pas mal de scripts pour filtrer un peu tout ça, essayer de ressortir les SNPs les plus plausibles.

      Voilà pour moi 🙂

  7. Bonjour
    C\' est vraiment un bon article.
    Ca me fait toujours plaisir de vous lire.
    Continuez comme ca... vous aidez beaucoup de gens, vu que nous pouvons avoir une idee du sujet au moins en francais. Car la plupart de ces sujets sont en anglais et ce n\'est pas facile pour les anglophobes 😀

  8. ClemBuntu :

    Désolée du retard de ma réponse, ça faisait un moment que je n\'avais pas vu venir me balader ici.

    Non non, j\'ai un (et même deux, pour avoir un réplicat) barcode par espèce bien sûr 🙂 mais effectivement du coup on n\'a que deux \"samples\" par espère car en fait pour nous le but du RAD est uniquement de trouver des SNP pertinents. Ils seront génotypés individuellement par spectrométrie de masse ensuite. On a fait ça pour des questions de coût, mais ça reste prometteur, on aura ou pas confirmation lors du retour de génotypage ;).

    Quentin J. :

    Effectivement j\'avais un peu plongé le nez dans les soft dont tu parles et étais arrivée aux même conclusions que toi ^^. Au final Stacks est un bon SNP caller même si dans mon cas il ne peut se passer d\'une grosse phase de vérification avec un aligneur (je réaligne, avec BWA, l\'ensemble des reads bruts contre mes séquences consensus extraites et compilées en fasta)

    Aurais-tu par hasard fait ton stage à l\'INRA de Bordeaux 😉 ?

    Bonne journée à tous !

  9. Salut,
    je suis heureux de decouvrir ce site il est très edifiant. je suis nouveau en bioinfo, je travaille sur ubuntu mais je ne parviens pas à deployer et à utiliser cette commande ainsi que ses differentes options : java -jar SnpSift.jar filter [options] \'expression‘

    besoin d\'aide !

    bonne journée à tous

    C-drick

    • Salut,
      Pour t\'aider j\'aurais besoin de plus d\'info ; ta version du logiciel, de Java et surtout ton message d\'erreur 😉

Laisser un commentaire