Introduction à l'analyse des SNPs

Introduction

Un SNP (Single Nucleo­tid Poly­mor­phism) est la muta­tion d'un seul nucléo­tide à une posi­tion don­née, entre indi­vi­dus d'une même popu­la­tion. Ces sub­sti­tu­tions ponc­tuelles per­mettent d'identifier des sous-popu­la­tions.

Vous avez sûre­ment étu­dié en cours l'exemple de la dré­pa­no­cy­tose, une mala­die cau­sée par la muta­tion ponc­tuelle d'un gène. L'origine de cette mala­die géné­tique peut être obser­vée en ana­ly­sant les SNPs dans la popu­la­tion humaine. Ana­ly­ser ces der­niers peut ain­si per­mettre de trou­ver des pré­dis­po­si­tions à cer­taines mala­dies, mais peut éga­le­ment être utile entre autres, en géné­tique des popu­la­tions ou en can­cé­ro­lo­gie. L'arrivée des séquen­ceurs de deuxième géné­ra­tion (les fameux NGS) a faci­li­té leur étude.

Pour la suite, nous par­ti­rons du prin­cipe que vous êtes fami­liers avec les concepts liés à l'analyse de don­nées NGS (reads, map­ping, for­mat SAM/​BAM, fastq, pro­fon­deur de séquen­çage etc.) sinon la suite risque d'être un peu ardue.

L'analyse

Pour ana­ly­ser vos SNPs, je vous pro­pose le pipe­line ci-des­sous. Ce der­nier a pour but d'identifier les SNPs pré­sents dans vos échan­tillons puis de les fil­trer. À la fin de l'article nous ver­rons com­ment les anno­ter par rap­port à une base de don­nées.

Pipeline d'analyse de SNPs
Pipe­line d'analyse de SNPs

 

L'analyse com­mence après ali­gne­ment de vos reads contre votre génome de réfé­rence via votre map­per pré­fé­ré. À la suite de cette étape vous obte­nez clas­si­que­ment un fichier au for­mat BAM. La pre­mière étape propre à l'analyse des SNPs consiste à uti­li­ser la com­mande

mpileup

, fournie avec l'outil

samtools

. Elle s'utilise de la manière suivante :

samtools mpileup -uf ref.fasta fichier.bam > output.bcf

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 pipeline

mpileup

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 :

bcftools view -v output.bcf > output.vcf

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 :

bcftools call -v -m output.bcf > output.vcf

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à :

samtools mpileup -uf ref.fasta fichier.bam  | bcftools call -v -m > output.vcf

 

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 :

java -jar SnpSift.jar filter [options] 'expression‘

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 :

java -jar SnpSift.jar filter "(DP >= 3)" -f Outputfile.vcf > Outputfile_filter.vcf

 

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 :

java -jar SnpSift.jar annotate 00-All.vcf myFile.vcf -info 'CAF' > annotatedSample.vcf

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.



Pour continuer la lecture :


Commentaires

17 réponses à “Introduction à l'analyse des SNPs”

  1. Avatar de Natir

    Bon­jour,

    Mer­ci pour votre article, fort inté­res­sant.

    Je tiens juste à indi­quer que le map­ping des reads contre un génome de réfé­rence n'est pas la seul solu­tion pour détec­ter des SNP. On peut même réa­li­ser ces détec­tions sans génome de réfé­rence comme avec Dis­coSNP https://​gatb​.inria​.fr/​s​o​f​t​w​a​r​e​/​d​i​s​c​o​s​np/ pour ne cité que celui que je connais le mieux.

  2. Avatar de ClemBuntu
    ClemBuntu

    Salut,
    Déci­dé­ment on en apprend tous les jours 🙂
    Moi ça me parai­trait plus ras­su­rant de map­per mes reads contre mon génome de réfé­rence pour être sur que la majo­ri­té s'aligne des­sus, mais si le génome de réfé­rence n'est pas dis­po ça peut être un bon outil.

    1. Avatar de Natir

      On peut aus­si regar­der du coté de Mind­The­Gap https://​github​.com/​G​A​T​B​/​M​i​n​d​T​h​e​Gap (de la même équipe qui a crée DiscoSNP)qui lui uti­lise un génome de réfé­rence et des reads NGS, pour détecte les SNP, les inser­tion et les délé­tions.

  3. Avatar de Sassi

    Bon­jour,

    Existe t‑il une méthode pour déter­mi­ner les SNPs sans génome de réfé­rence en map­pant des génomes non des reads ? j'ai trou­vé un seul mais avec réfé­rence : https://​har​vest​.read​the​docs​.io/​e​n​/​l​a​t​e​st/

  4. Avatar de ClemBuntu
    ClemBuntu

    Map­per des génomes ?
    Sans génome de réfé­rence il y a dis­co SNP cf le 1er com­men­taire 😉

    1. Avatar de Sassi

      Dis­coSNP mappe des reads, je vou­lais savoir s'il existe un pro­gramme uti­li­sant input = des génomes assem­blés en fas­ta ?

    1. Avatar de Sassi

      Mer­ci, mais je ne pense pas que ces outils d'alignement seront bien adap­té à de grandes don­nées je vais peut être teste mafft.

      1. Avatar de Yoann M.
        Yoann M.

        Bon­jour,

        Je serai toi je par­ti­rai d'abord sur Clus­tal Ome­ga qui est plus récent est mieux adap­té aux gros jeux de don­nées.
        Si ça ne marche pas, alors j'envisagerai ensuite MAFFT ou/​et MUSCLE.
        Bon cou­rage.

  5. Avatar de Chrys

    Salut à tous !

    Pre­mière fois que je poste sur le site 😛 Or ce n'est pas faute de m'intéresser à tout ce que vous faites (Mer­ci ^_​^ !) J'attends impa­tiem­ment d'avoir quelques vraies jour­nées devant moi pour vous lire entiè­re­ment (ou en grande par­tie du moins !)

    Je suis étu­diante en thèse de géné­tique des popu­la­tions et je vais être ame­née très pro­chai­ne­ment à ana­ly­ser des don­nées brutes de RAD-Seq (et plus pré­ci­sé­ment de Pool-seq, aoutch 😉 !) Pour le moment j'ai sur­tout tra­vaillé avec Stacks, que je m'étonne de ne pas voir lis­té ici (et sauf erreur, il n'a pas non plus été évo­qué dans un autre article du blog ?), y a t'il une rai­son par­ti­cu­lière à cela ?

    Sachant que de mon côté mon petit chal­lenge per­son­nel sera de fil­trer mes SNP en sachant que je n'ai qu'un pool (donc un seul bar­code) pour cha­cune de mes 18 espèces, et zéro génome de réfé­rence ! Je pen­sais par­tir sur une pipe­line type Vel­ve­top­ti­mi­ser, BWA/​Bowtie2, puis implé­men­ter tout ça dans Popoolation2, un logi­ciel qui per­met jus­te­ment de faire du SNP cal­ling sur du pool. Peut-être que cer­tains ici ont fait quelque chose de simi­laire 🙂 ?

    En tout cas bra­vo encore, j'essaierai de vous voir aux JOBIM puisque vous y serez 🙂 !

    Chrys

    1. Avatar de Yoann M.
      Yoann M.

      Salut Chrys,

      Mer­ci pour ton com­men­taire, il fait chaud au coeur !
      Je ne vais pas pou­voir t'aider concer­nant ta requête, ce n'est abso­lu­ment pas ma spé­cia­li­té. Mais si tu n'obtiens pas de réponse d'ici JOBIM ? je t'encourage for­te­ment à venir nous ren­con­trer, il y aura au moins une experte (Isa­belle) 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 pour­rait être toi qui écri­ra le pre­mier sur ce sujet ici ? 🙂

      Pour info : on don­ne­ra une brève pré­sen­ta­tion lors du work­shop JeBiF le lun­di 27 juin et on aura ensuite un pos­ter bioin​fo​-fr​.net tout au long de la conf 🙂

      Au plai­sir de dis­cu­ter avec toi !

    2. Avatar de ClemBuntu
      ClemBuntu

      Salut,
      Pour ce qui est de l'article ici l'idée c'était de faire une intro­duc­tion aux ana­lyses de SNPs, au for­mat VCF et com­pa­gnie, pas du tout d'être exhaus­tif et de dérou­ler la liste de tous les logi­ciels exis­tants.
      Après comme l'a dit Yoann, se foca­li­ser sur quelques logi­ciels qui ont des par­ti­cu­la­ri­tés inté­res­santes, ça pour­rait faire le sujet d'un autre article 😉

      T'as 18 espèces dif­fé­rentes et tous tes reads ont le même index ? C'est…bizarre comme façon de faire, bon cou­rage !

    3. Avatar de Quentin J.
      Quentin J.

      Salut !

      Le tra­vail que tu vas effec­tuer est rela­ti­ve­ment simi­laire à ce que je fais actuel­le­ment en stage de Bioin­fo à l'INRA. J'vais du coup appor­ter ma petite pierre à cet article en espé­rant que ça puisse (t')aider quelques per­sonnes.

      De mon côté je bosse éga­le­ment sur des don­nées brutes de RAD-Seq et concrè­te­ment je dois tes­ter plu­sieurs méthodes de détec­tion de SNP, les com­pa­rer, voir les points posi­tifs et néga­tifs de cha­cune et déter­mi­ner s'il y en a pas une meilleure dans tout ça. Bien sûr ce que je vais citer s'applique plus par­ti­cu­liè­re­ment à mon contexte, c‑a-d l'étude d'une popu­la­tion d'arbres (sans génome de réfé­rence /​ orga­nisme modèle). Soit un nombre assez consé­quent de séquences…

      J'ai débu­té mon tra­vail à par­tir de 8 fichiers .fastq (8 indi­vi­dus séquen­cés). Ce que j'ai tes­té :

      - Stacks. Effec­ti­ve­ment, un bon outil lorsqu'on a un paquet de don­nées. L'avantage de ce pro­gramme est la rapi­di­té de trai­te­ment (1 nuit pour 20Go de don­nées), même sans génome de réfé­rence. J'ai pas encore véri­fié les SNPs qu'il trouve mais j'ai ampli­fié par PCR un lot de loci qu'il a déter­mi­né et la plu­part sont exis­tants.

      - Pyrad. Très simi­laire à Stacks seule­ment celui-ci per­met de tra­vailler à par­tir de séquences de lon­gueurs variables. Un avan­tage comme un incon­vé­nient car il demande énor­mé­ment de RAM… et de temps de cal­cul (1 semaine de trai­te­ment sur clus­ter :-|). Les résul­tats sont com­pa­rables à Stacks, plu­tôt posi­tifs dans l'ensemble.

      - CLC /​ (BWA-SAMTOOLS) /​ Reads2SNP. Avec ces 4 outils j'ai tes­té une méthode dif­fé­rente des deux autres, c‑a-d construire un génome de ref à par­tir de toutes les séquences (CLC) puis map­per ces séquences sur cette ref (BWA). L'idée peut paraître sédui­sante mais les résul­tats le sont moins. La plu­part des loci pos­sé­dant des SNPs (détec­tés par R2S) de ne sont pas ampli­fiés lorsque j'ai essayé de véri­fier par PCR (quelque chose à voir avec la ref de CLC je sup­pose).

      - Dis­coSNP. Pro­ba­ble­ment pas adap­té aux don­nées sur les­quelles je tra­vaille. Je n'ai pas pu exploi­ter les .vcf obte­nus.

      - GATK. Pas adap­té aux don­nées RAD-Seq. Mais un must test.

      Les trois pre­mières méthodes m'ont don­né +500 000 SNPs sur +50 000 loci (à par­tir de +80 mil­lions de reads). A par­tir de là j'ai écris pas mal de scripts pour fil­trer un peu tout ça, essayer de res­sor­tir les SNPs les plus plau­sibles.

      Voi­là pour moi 🙂

  6. Avatar de Lansana
    Lansana

    Bon­jour
    C' est vrai­ment un bon article.
    Ca me fait tou­jours plai­sir de vous lire.
    Conti­nuez comme ca… vous aidez beau­coup de gens, vu que nous pou­vons avoir une idee du sujet au moins en fran­cais. Car la plu­part de ces sujets sont en anglais et ce n'est pas facile pour les anglo­phobes 😀

  7. Clem­Bun­tu :

    Déso­lée du retard de ma réponse, ça fai­sait un moment que je n'avais pas vu venir me bala­der ici.

    Non non, j'ai un (et même deux, pour avoir un répli­cat) bar­code par espèce bien sûr 🙂 mais effec­ti­ve­ment du coup on n'a que deux "samples" par espère car en fait pour nous le but du RAD est uni­que­ment de trou­ver des SNP per­ti­nents. Ils seront géno­ty­pés indi­vi­duel­le­ment par spec­tro­mé­trie de masse ensuite. On a fait ça pour des ques­tions de coût, mais ça reste pro­met­teur, on aura ou pas confir­ma­tion lors du retour de géno­ty­page ;).

    Quen­tin J. :

    Effec­ti­ve­ment j'avais un peu plon­gé le nez dans les soft dont tu parles et étais arri­vée aux même conclu­sions que toi ^^. Au final Stacks est un bon SNP cal­ler même si dans mon cas il ne peut se pas­ser d'une grosse phase de véri­fi­ca­tion avec un ali­gneur (je réaligne, avec BWA, l'ensemble des reads bruts contre mes séquences consen­sus extraites et com­pi­lées en fas­ta)

    Aurais-tu par hasard fait ton stage à l'INRA de Bor­deaux 😉 ?

    Bonne jour­née à tous !

  8. Avatar de C-drick

    Salut,
    je suis heu­reux de decou­vrir ce site il est très edi­fiant. je suis nou­veau en bioin­fo, je tra­vaille sur ubun­tu mais je ne par­viens pas à deployer et à uti­li­ser cette com­mande ain­si que ses dif­fe­rentes options : java ‑jar SnpSift.jar fil­ter [options] 'expres­sion‘

    besoin d'aide !

    bonne jour­née à tous

    C‑drick

    1. Avatar de clembuntu
      clembuntu

      Salut,
      Pour t'aider j'aurais besoin de plus d'info ; ta ver­sion du logi­ciel, de Java et sur­tout ton mes­sage d'erreur 😉

Laisser un commentaire