Une comparaison des pipelines d’appel de variants en utilisant Genome in a Bottle comme référence
Abstract
Le séquençage à haut débit, en particulier des exomes, est un outil de diagnostic populaire, mais il est difficile de déterminer quels outils sont les meilleurs pour analyser ces données. Dans cette étude, nous utilisons les résultats du NIST Genome in a Bottle comme une nouvelle ressource pour la validation de notre pipeline d’analyse des exomes. Nous utilisons six aligneurs différents et cinq appelants de variants différents pour déterminer quel pipeline, sur les 30 au total, est le plus performant sur un exome humain qui a été utilisé pour aider à générer la liste des variants détectés par le Consortium Genome in a Bottle. Parmi ces 30 pipelines, nous avons constaté que Novoalign, associé à GATK UnifiedGenotyper, présentait la plus grande sensibilité tout en maintenant un faible nombre de faux positifs pour les SNV. Cependant, il est évident que les indels sont encore difficiles à gérer pour n’importe quel pipeline, aucun des outils n’atteignant une sensibilité moyenne supérieure à 33 % ou une valeur prédictive positive (VPP) supérieure à 53 %. Enfin, comme prévu, il a été constaté que les aligneurs peuvent jouer un rôle aussi vital dans la détection des variants que les appelants de variants eux-mêmes.
1. Contexte
Au cours des dernières années, de nombreux progrès ont été réalisés dans les technologies de séquençage à haut débit. Grâce à ces progrès, il est maintenant possible de détecter un grand nombre de variants potentiels causant des maladies , et, dans quelques cas, les données de séquençage de nouvelle génération (NGS) ont même été utilisées à des fins de diagnostic . Cela est dû en partie aux développements des technologies de séquençage au cours des dernières années, mais aussi au nombre d’améliorations apportées aux différents outils bioinformatiques utilisés pour analyser les montagnes de données produites par les instruments NGS .
Lors de la recherche de mutations chez un patient, un flux de travail typique consiste à séquencer son exome avec un séquenceur Illumina, à aligner les données brutes sur le génome humain de référence, puis à identifier les variantes nucléotidiques uniques (SNV) ou les insertions et délétions courtes (indels) qui pourraient éventuellement causer ou influencer le phénotype d’intérêt . Si cela est assez simple, il n’en va pas de même pour le choix des meilleurs outils à utiliser à chaque étape du pipeline d’analyse. Il existe un grand nombre d’outils qui sont utilisés dans diverses étapes intermédiaires, mais les deux étapes les plus importantes de l’ensemble du processus sont l’alignement des lectures brutes sur le génome et la recherche de variants (c’est-à-dire les SNV et les indels). Dans cette étude, nous visons à aider le bioinformaticien d’aujourd’hui en élucidant la combinaison correcte de l’outil d’alignement des lectures courtes et de l’outil d’appel de variants pour traiter les données de séquençage de l’exome produites par les instruments NGS.
Un certain nombre de ces études ont été réalisées dans le passé, mais elles présentaient toutes des inconvénients sous une forme ou une autre. Idéalement, il faudrait disposer d’une liste de toutes les variantes connues contenues dans un échantillon, de sorte que lorsqu’un pipeline d’outils d’analyse est exécuté, on puisse le tester pour savoir avec certitude qu’il fonctionne correctement. Cependant, dans le passé, une telle liste n’existait pas, de sorte que la validation devait être effectuée par des méthodes moins complètes. Dans certains cas, la validation était effectuée en générant des données simulées afin de créer un ensemble de vrais positifs (TP) et de vrais négatifs (TN) connus. Bien que cette méthode permette d’obtenir une liste de tous les TP et TN de l’ensemble de données, elle ne permet pas de représenter fidèlement la biologie. D’autres méthodes de validation des pipelines d’appel de variants comprennent l’utilisation de matrices de génotypage ou le séquençage Sanger pour obtenir une liste de TP et de faux positifs (FP). Ces méthodes ont l’avantage de fournir des résultats biologiquement validés, mais elles ont aussi l’inconvénient de ne pas être complètes en raison du nombre limité de spots sur les matrices de génotypage et du coût prohibitif de la validation Sanger lorsqu’elle est effectuée des milliers de fois. Enfin, aucune de ces études ne visait à examiner l’effet de l’aligneur de lecture courte sur l’appel de variants. Par conséquent, l’effet en amont de la performance de l’aligneur n’a pas pu être évalué de manière indépendante.
Dans cette étude, nous avons l’avantage de disposer d’une liste de variants pour une femme anonyme de l’Utah (ID du sujet : NA12878, séquencé à l’origine pour le projet 1000 Génomes ) qui a été validée expérimentalement par le Consortium Genome in a Bottle (GiaB) dirigé par le NIST. Cette liste de variants a été créée en intégrant 14 ensembles de données différents provenant de cinq séquenceurs différents, et elle nous permet de valider toute liste de variants générée par nos pipelines d’analyse d’exome. La nouveauté de ce travail est de valider la bonne combinaison d’aligneurs et d’appelants de variants contre un ensemble de données de variants complet et déterminé expérimentalement : NIST-GiaB.
Pour effectuer notre analyse, nous utiliserons l’un des ensembles de données d’exome initialement utilisés pour créer la liste NIST-GiaB. Nous n’avons choisi qu’un seul des exomes générés à l’origine par Illumina TruSeq parce que nous voulions fournir un scénario de cas d’utilisation standard pour quelqu’un qui souhaite effectuer une analyse NGS, et alors que le séquençage du génome entier continue à baisser en prix, le séquençage de l’exome est toujours une alternative populaire et viable . Il est également important de noter que, selon Bamshad et al. le nombre attendu de SNV par exome européen-américain est actuellement de 20 283 ± 523 . Malgré cela, le nombre total de SNV trouvés dans la liste NIST-GiaB et susceptibles d’exister dans l’ensemble de données d’exome TruSeq était de 34 886, ce qui est considérablement plus élevé que prévu. Cela est probablement dû au fait que si le kit d’exome a été utilisé pour générer les données NIST-GiaB, il a également été complété par un séquençage du génome entier.
En dernier lieu, nous avons considéré un grand nombre d’aligneurs et d’appelants de variants, mais nous avons finalement choisi les 11 outils en fonction de leur prévalence, de leur popularité et de leur pertinence pour notre ensemble de données (par exemple, SNVMix, VarScan2 et MuTect n’ont pas été utilisés car ils sont destinés à être utilisés sur des échantillons dérivés de tumeurs). Notre analyse comprend la comparaison de six aligneurs (Bowtie2, BWA sampe, BWA mem, CUSHAW3, MOSAIK et Novoalign) et de cinq appelants de variants (FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup et SNPSVM). Dans cette étude, nous essayons également de déterminer l’effet, s’il y en a un, de l’aligneur sur l’appel de variante et quels aligneurs sont les plus performants lorsqu’on utilise un échantillon normal d’exome Illumina. À notre connaissance, il s’agit du premier rapport qui valide toutes les combinaisons possibles (total de 30 pipelines) d’un large éventail d’aligneurs et d’appelants de variants.
2. Méthodes
2.1. Jeux de données
Le génome de référence humain hg19 a été téléchargé à partir du navigateur UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) et a été utilisé pour effectuer les alignements. L’exome humain, SRR098401, a été téléchargé à partir de Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). À des fins d’annotation et d’étalonnage, les listes dbSNP137 sans sites après la version 129, HapMap 3.3, Human Omni 2.5 BeadChip et Mills et 1000 G gold standard indel set ont été utilisées (toutes sur ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).
2.2. Le pipeline
La figure 1 montre le flux de travail utilisé dans cette étude, qui est similaire à celui décrit dans le guide des meilleures pratiques produit par le Broad Institute . Cela implique un certain nombre d’étapes pour s’assurer que les fichiers d’alignement produits sont de la plus haute qualité, ainsi que plusieurs autres pour garantir que les variants sont appelés correctement. Tout d’abord, les lectures brutes ont été alignées sur hg19, puis les doublons PCR ont été supprimés de l’alignement. Ensuite, pour faciliter l’identification des indels plus tard dans le pipeline, un réalignement des lectures a été effectué autour des indels. La dernière étape du traitement de l’alignement a consisté à effectuer une étape de recalibrage du score de qualité des bases, ce qui permet d’améliorer le biais inhérent et les inexactitudes des scores émis par les séquenceurs. Malheureusement, malgré ces étapes, le taux d’alignement de chaque aligneur était nettement inférieur aux prévisions. Pour compenser, la boîte à outils fastx a été utilisée pour filtrer les lectures de faible qualité (tableau 1). Les lectures de faible qualité ont été définies comme étant celles dont au moins la moitié des scores de qualité étaient inférieurs à 30. Après le traitement de l’alignement, l’appel de variants et le filtrage des variants ont été effectués.
|
||||||||||||||||||||||||||||||||||||||||
Les six outils utilisés pour générer des alignements étaient Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK et Novoalign, et les cinq outils utilisés pour générer des variants étaient FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup et SNPSVM, comme on peut le voir dans le tableau 2.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2.3. Filtrage
Les données brutes ont été acquises à partir du SRA (SRR098401), divisées avec fastq-dump, et filtrées à l’aide de la boîte à outils fastx. Plus précisément, fastq-dump a utilisé les drapeaux -split_files et -split_spot, et fastq_quality_filter a été exécuté avec les arguments suivants : -Q 33 -q 30 -p 50. Ensuite, les lectures ont été correctement appariées avec un script personnalisé.
2.4. Alignement
Les aligneurs ont utilisé des arguments par défaut, sauf lorsqu’un argument filaire était utilisé lorsqu’il était disponible. Les commandes utilisées sont les suivantes.
2.4.1. Bowtie2
(1)bowtie2 -p 10 -x $INDEX -1 raw_data/read_1_filtered.fastq -2 raw_data/read_2_filtered.fastq -S alignments/NA12878.bt2.sam
2.4.2. BWA Sampe
(1)bwa aln -t 10 genome/hg19.fa raw_data/read_1_filtered.fastq > alignments/NA12878.R1.sai(2)bwa aln -t 10 genome/hg19.fa raw_data/read_2_filtered.fastq > alignments/NA12878.R2.sai(3)bwa sampe genome/hg19.fa alignements/NA12878.R1.sai alignements/NA12878.R2.sai raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > alignements/NA12878.bwa-sampe.sam
2.4.3. BWA Mem
(1)bwa mem -t 10 genome/hg19.fa raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > alignements/NA12878.bwa-mem.sam
2.4.4. CUSHAW3
(1)cushaw3 align -r $INDEX -t 10 -o alignements/NA12878.CUSHAW3.sam -q raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq
2.4.5. MOSAIK
(1)MosaikBuild -q raw_data/read_1_filtered.fastq -q2 raw_data/read_2_filtered.fastq -st illumina -outalignments/NA12878.MOSAIK.mkb(2)MosaikAligner -in alignments/NA12878.MOSAIK.mkb -out alignments/NA12878.MOSAIK -p 10 -ia genome/hg19.dat -j genome/hg19_15 -annpe tools/MOSAIK/src/networkFile/2.1.78.pe.ann -annse tools/MOSAIK/src/networkFile/2.1.78.se.ann
2.4.6. Novoalign
(1)novoalign -d $INDEX -f raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq -o SAM -c 10 > alignments/NA12878.novoalign.sam
2.5. Calcul de la profondeur de couverture des alignements
Pour assurer un calcul correct de la profondeur de couverture, le module CalculateHsMetrics de Picard Tools a été utilisé avec les arguments suivants :(1)java -jar CalculateHsMetrics.jar I=NA12878.ALN.BQSR.bam O=ALN.O.log R=genome/hg19.fa TI=genome/truseq_exome.bed BI=genome/truseq_exome.bed VALIDATION_STRINGENCY=SILENT PER_TARGET_COVERAGE=ALN.ptc.bedIl est important de noter que l’en-tête du fichier d’exome TruSeq doit être précédé du fichier d’alignement SAM pour que ce module fonctionne. De plus, la colonne 6 doit être déplacée vers la colonne 4, et la colonne 5 doit être retirée du fichier de lit TruSeq.
2.6. Traitement des fichiers d’alignement
Le traitement des fichiers d’alignement (fichiers SAM/BAM) a nécessité les étapes suivantes pour tous les aligneurs :(1)Conversion de SAM en BAM avec SAMtools view :(a)samtools view -bS alignments/NA12878.ALN.sam -o alignements/NA12878.ALN.bam(2)Tri des fichiers BAM à l’aide du module Picard Tools, SortSam :(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignements/NA12878.ALN.bam OUTPUT=alignements/NA12878.ALN.sorted.bam SORT_ORDER=coordinate(3)Suppression des doublons PCR à l’aide du module Picard Tools, MarkDuplicates :(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.sorted.bam O=alignments/NA12878.ALN.dups_removed.bam REMOVE_DUPLICATES=true M=alignements/métriques(4)Groupe de lecture ajouté aux fichiers d’alignement à l’aide du module Picard Tools, AddOrReplaceReadGroups :(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignements/NA12878.ALN.dups_removed.bam O=alignements/NA12878.ALN.RG.bam SO=coordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Réalignement autour des indels en utilisant les modules GATK RealignerTargetCreator et IndelRealigner :(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T RealignerTargetCreator -R genome/hg19.fa -I alignments/NA12878.ALN.RG.bam -known genome/mills.vcf -o tmp/ALN.intervals(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T IndelRealigner -R genome/hg19.fa -I alignments/NA12878.ALN.RG.bam -known genome/mills.vcf -o alignments/NA12878.ALN.indels.bam -maxReadsForRealignment 100000 -maxReadsInMemory 1000000 -targetIntervals tmp/ALN.intervals(6)Recalibrage de la base en utilisant les modules BaseRecalibrator et PrintReads de GATK :(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T BaseRecalibrator -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -knownSitesgenomedbsnp_137.hg19.excluding_sites_after_129.only_standard_chroms.vcf -o tmp/NA12878.ALN.grp(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T PrintReads -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -BQSR tmp/NA12878.ALN.grp -o alignments/NA12878.ALN.BQSR.bam
2.7. Appel de variants
Les arguments par défaut ont été utilisés pour chaque appel de variants, à moins qu’il ne contienne un drapeau « threads » ou « parallel », auquel cas cela a également été utilisé. De plus, les indels ont été appelés séparément des SNVs lorsque cela était possible. Plus précisément, les commandes utilisées sont les suivantes.
2.7.1. FreeBayes
(1)freebayes -f genome/hg19.fa -i -X -u -v vcf_files/NA12878.ALIGNER.freebayes.raw.snv.vcf alignments/NA12878.ALIGNER.BQSR.bam(2)freebayes -f genome/hg19.fa -I -X -u -v vcf_files/NA12878.ALIGNER.freebayes.raw.indel.vcf alignments/NA12878.ALIGNER.BQSR.bam
2.7.2. GATK HaplotypeCaller
(1)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R genome/hg19.fa -I alignments/NA12878.ALIGNER.BQSR.bam -dbsnp $DBSNP -o vcf_files/NA12878.ALIGNER.HC.raw.vcf -stand_call_conf 50
2.7.3. GATK UnifiedGenotyper
(1)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R genome/hg19.fa -nt 10 -I alignments/NA12878.ALIGNER.BQSR.bam -o vcf_files/NA12878.ALIGNER.UG.raw.snv.vcf -glm SNP -D $DBSNP(2)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R genome/hg19.fa -nt 10 -I alignments/NA12878.ALIGNER.BQSR.bam -o vcf_files/NA12878.ALIGNER.UG.raw.indel.vcf -glm INDEL -D $MILLS
2.7.4. SAMtools Mpileup
(1)samtools mpileup -uf genome/hg19.faalignments/NA12878.ALIGNER.BQSR.bam | bcftools view -bvcg – > vcf_files/NA12878.ALIGNER.mpileup.bcf && bcftools view vcf_files/NA12878.ALIGNER.mpileup.bcf > vcf_files/NA12878.ALIGNER.mpileup.raw.vcf
2.7.5. SNPSVM
(1)java -XX:ParallelGCThreads=10 -jar tools/SNPSVM/snpsvm.jar predict -R genome/hg19.fa -B alignments/NA12878.ALIGNER.BQSR.bam -M tools/SNPSVM/models/default.model -V vcf_files/NA12878.ALIGNER.SNPSVM.raw.vcfDu fait de l’inexistence des drapeaux CIGAR requis dans le fichier d’alignement, SNPSVM n’a pas réussi à appeler les variants pour CUSHAW3, et SAMtools mpileup n’a pas pu appeler les variants sur les alignements MOSAIK pour la même raison. De plus, en raison du fait que SNPSVM ne détecte que les SNV, aucun indel n’a été signalé pour ce programme.
2.8. Filtration des variants
La filtration variait en fonction de l’appelant de variants utilisé. Dans les cas de GATK HaplotypeCaller et de GATK UnifiedGenotyper, les modules GATK, VariantRecalibrator et ApplyRecalibrator, ont été utilisés pour filtrer les SNV en utilisant HapMap 3.3, l’Omni 2.5 SNP BeadChip et dbSNP 137 sans les données de 1000 Genome comme ensembles d’entraînement. Pour SNPSVM, des scores QUAL ≥ 4 et des valeurs DP ≥ 6 ont été utilisés. Pour FreeBayes et SAMtools, des scores QUAL ≥ 20 et des valeurs DP ≥ 6 ont été utilisés.
2.9. Comparaison des variants
Pour la comparaison des variants, USeq 8.8.1 a été utilisé pour comparer les SNV partagés entre tous les ensembles de données. Pour comparer les indels, l’outil vcflib vcfintersect a été utilisé. Le fichier TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, obtenu le 11 décembre 2013, a été utilisé pour restreindre les comparaisons aux emplacements qui pouvaient être capturés par le kit pull down exome utilisé pour le séquençage de SRR098401. Ce fichier peut être obtenu auprès d’Illumina ici : http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Pour s’assurer que les variants étaient représentés de manière identique entre les différents ensembles d’appels, l’outil vcflib vcfallelicprimitives a été utilisé pour prétraiter les fichiers vcf.
2.10. Calculs statistiques
Positif vrai (TP). C’est une mutation qui a été détectée par le pipeline testé et qui existe dans la liste NIST-GiaB.
Faux positif (FP). Il s’agit d’une mutation qui a été détectée par le pipeline testé mais qui n’existe pas dans la liste NIST-GiaB.
Vraiment négatif (TN). Il s’agit d’une mutation qui n’a pas été détectée par le pipeline testé et qui est une mutation qui n’existe pas dans la liste NIST-GiaB.
Faux négatif (FN). Il s’agit d’une mutation qui n’a pas été détectée par le pipeline testé mais qui existe dans la liste NIST-GiaB :
3. Résultats et discussion
3.1. Préfiltrage des variants
Lorsqu’on effectue une analyse de variants, l’un des nombreux pièges à prendre en considération est l’espace de séquence de l’exome (tel que défini par le kit de capture d’exome) et la façon dont il peut affecter les résultats de l’analyse. Dans ce cas, nous disposions d’un seul exome (SRR098401) qui a été extrait à l’aide du kit d’exome TruSeq d’Illumina et séquencé sur un HiSeq 2000. Dans cette optique, nous voulions nous assurer que nous mesurions la capacité des outils bioinformatiques à faire leur travail et non l’efficacité de la trousse de capture d’exome TruSeq d’Illumina. En d’autres termes, nous voulons seulement essayer d’appeler les variants qui sont censés être présents dans les exons tels que définis par le kit de capture d’exome. Pour cette raison, nous utilisons le fichier de lit fourni par Illumina, et non un fichier de lit d’annotation générique, par exemple, RefSeq pour hg19. Nous avons constaté que pour cet individu particulier, selon la liste NIST-GiaB, il devrait y avoir un total de 34 886 SNV et 1 473 indels dans les régions définies par le fichier de lit TruSeq.
Une fois que nous avons filtré les variants qui n’étaient pas situés dans les régions définies par le fichier de lit de l’exome TruSeq d’Illumina, nous sommes passés de centaines de milliers de variants putatifs (données non présentées) à, en moyenne, environ 23 000 variants (SNV et indels) par pipeline (tableau 3). Il s’agit d’une étape importante pour les chercheurs pour commencer, car elle réduit considérablement l’espace de recherche de variants potentiellement intéressants.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3.2. Résultats bruts des variants comparés à GiaB
Un aspect que nous voulions comprendre en faisant cette comparaison était l’importance du filtrage des variants détectés par ces outils. La raison en est qu’idéalement, on aimerait avoir un niveau de sensibilité aussi élevé que possible afin que les mutations d’intérêt ne soient pas perdues dans le processus de filtrage. Il nous incombe donc de déterminer si cette étape est nécessaire ou non et dans quelle mesure elle l’est, puisqu’il ressort clairement des résultats du NIST-GiaB et de l’examen de Bamshad et al. que la sensibilité pourrait être un problème.
Comme nous pouvons le voir dans le tableau 3, le filtrage est plus nécessaire pour certains appelants de variantes que pour d’autres lorsqu’il s’agit de SNV, et il est absolument nécessaire pour les indels. Dans la plupart des cas, le nombre de variants TP est proche ou supérieur à notre nombre attendu d’environ 20 000 , mais, d’autre part, dans certains cas, le nombre de FP est très élevé.
Il est clair qu’il y a beaucoup de variation dans les chiffres générés par chaque pipeline. Cependant, on peut trouver quelques points communs dans les nombres qui proviennent probablement des origines algorithmiques de chaque outil. FreeBayes produit à la fois le plus grand nombre de variantes non filtrées et le plus grand nombre de FP. Il est probable que nous ne voyons ce genre de performance de la part de cet outil que parce que, bien qu’il ne soit pas le seul appelant de variantes basé sur l’inférence bayésienne, il est unique dans son interprétation des alignements. C’est-à-dire qu’il s’agit d’un appelant basé sur les haplotypes qui identifie les variants sur la base de la séquence des lectures elles-mêmes au lieu de l’alignement, ce dernier étant la façon dont fonctionne le UnifiedGenotyper de GATK.
De plus, nous voyons les aligneurs basés sur Burrows-Wheeler avoir des performances très similaires les uns aux autres : Bowtie2, BWA mem, et BWA sampe obtiennent des résultats similaires sur toute la ligne. On pourrait supposer que cela est probablement dû au fait que tous ces outils utilisent des algorithmes similaires lorsqu’ils effectuent leur tâche désignée. Cette observation est soutenue par le fait que MOSAIK (alignements gappés utilisant l’algorithme de Smith-Waterman) et CUSHAW3 (une approche hybride de semis) ont tous deux des algorithmes sous-jacents très différents et produisent par la suite des résultats très différents.
Cette différence de résultats en corrélation avec différents algorithmes est mieux vue dans les résultats de SNPSVM. Parmi les appelants de variants, c’est le seul qui utilise des machines à vecteurs de support et la construction de modèles pour générer des appels SNV. Il semblerait que, même si elle présente l’inconvénient de ne pas être aussi sensible que d’autres méthodes, elle a l’avantage d’être extrêmement précise, quel que soit l’aligneur utilisé. Cela suggère que l’on est capable de sauter l’étape de filtrage complètement en utilisant cet appel de variante.
En ce qui concerne les indels, aucun aligneur ne semble se démarquer parmi les autres comme étant celui qui gère bien ce type de mutation. En fait, si l’on considère le nombre de FP, il est clair que c’est l’appelant de variante qui joue le plus grand rôle dans la précision de l’identification des indels. De plus, il n’y a pas de données pour les pipelines CUSHAW3 plus SNPSVM ou MOSAIK plus SAMtools mpileup car les fichiers d’alignement ne contiennent pas les chaînes CIGAR nécessaires pour que les appelants de variants fonctionnent en aval. Enfin, la raison pour laquelle il n’y a pas de données indel pour SNPSVM est que cet outil est uniquement utilisé pour l’identification des SNV.
3.3. Résultats filtrés par rapport à GiaB
Comme dans le tableau 4, les pratiques de filtrage standard parviennent à éliminer un grand nombre de SNV FP pour chaque pipeline ; il semble toutefois que ces filtres soient un peu trop agressifs dans la plupart des cas pour la détection des SNV, mais pas assez stricts pour les indels. Ceci est rendu évident en regardant les différences dans le nombre de FP rapportées dans chaque ensemble de données. Par exemple, Bowtie2 avec Freebayes voit une suppression de 72 570 FP SNVs (une réduction de 98%) mais seulement une suppression de 1 736 FP indels (une réduction de 70%). Il convient également de noter que les filtres utilisés dépendaient du pipeline et que, pour la plupart, chaque pipeline a produit des réductions similaires des FP SNV et indel. La seule exception ici est le nombre de variants identifiés à partir des alignements CUSHAW3 par rapport aux autres alignements : globalement le nombre de SNV TP est plus faible, le nombre de SNV FP est plus élevé, et c’est le seul aligneur qui produit plus de 1000 indels FP après filtrage.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Compte tenu du fait que le filtrage réduit considérablement le nombre de variants TP, il pourrait être judicieux, à l’exception des pipelines utilisant CUSHAW3 et FreeBayes, de sauter cette étape lors de la recherche de variants rares à fort impact. Au lieu de cela, on pourrait passer plus de temps sur un processus de filtrage qui est basé sur la biologie plutôt que sur les statistiques. Par exemple, il peut être plus judicieux d’investir du temps dans l’identification d’une petite liste de variants susceptibles d’avoir un impact élevé : mutations de sites d’épissage, indels qui provoquent des décalages de cadre, mutations de troncation, mutations de perte d’arrêt, ou mutations dans des gènes dont on sait qu’ils sont biologiquement pertinents pour le phénotype d’intérêt.
3.4. TPR et sensibilité moyenne
Comme on peut le voir dans le tableau 5, la valeur prédictive positive (VPP) de chaque outil, à l’exception de CUSHAW3, varie de 91% à 99,9% pour les SNV, mais la sensibilité moyenne est très faible (environ 50%). Cet écart peut être dû à un certain nombre de raisons, mais la plus probable est la profondeur variable de la couverture entre les exons. Nous pouvons voir qu’en plus de la faible sensibilité des SNV, la sensibilité des indels est faible (environ 30%) ; cependant, la VPP pour les indels est considérablement plus faible (35,86% à 52,95%). Cela pourrait être dû à l’une des raisons suivantes : les indels très courts sont difficiles à détecter par NGS conventionnel , la représentation des indels par différents appelants de variants peut amener les outils à prétendre de manière incorrecte que deux indels sont différents, ou les outils d’alignement produisent des représentations différentes du même indel .
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Peut-être que l’explication la plus probable pour les deux types de mutations est la question de la profondeur. Comme c’est le cas pour toute étude d’analyse de variants, une augmentation de la profondeur de couverture entraîne une augmentation de la sensibilité, mais il est impossible de garantir une bonne profondeur de couverture en raison de l’incapacité des kits de capture d’exome à tirer uniformément les exons . En outre, aucun kit de capture d’exome ne couvre tous les exons. En effet, il a été démontré que l’analyse des variantes du séquençage du génome entier à une profondeur moyenne inférieure à celle d’un exome donne de meilleurs résultats en raison de l’uniformité de cette profondeur. Ainsi, il est probable qu’un grand nombre de variants manquent du fait que la liste NIST-GiaB a été créée à partir d’une compilation de données de séquençage exomique et génomique. En fin de compte, pour obtenir une sensibilité appropriée, il faudra éventuellement effectuer un séquençage du génome entier, mais le coût en est actuellement prohibitif pour la plupart des laboratoires. Heureusement, ce coût continue de baisser, et nous verrons bientôt un passage progressif de l’analyse des exomes à l’analyse plus complète du génome entier.
3.5. Sensibilité en fonction de la profondeur
Parce que la sensibilité reflète l’une des mesures de performance les plus importantes d’un outil et que la plupart des outils luttent pour atteindre une sensibilité supérieure à 50%, nous aimerions explorer davantage comment la profondeur affectait la sensibilité d’appel de variants. Nous avons examiné un certain nombre de combinaisons différentes d’outils pour déterminer quels étaient les meilleurs pipelines, appelants de variantes et aligneurs. Pour la figure 2, nous avons pris les cinq meilleures combinaisons d’appelants de variantes et d’aligneurs, déterminées par leur sensibilité et leur taux de faux positifs (FPR). En d’autres termes, nous avons sélectionné celles qui avaient le plus grand nombre de SNV TP appelés en plus du plus petit nombre de SNV FP. Après inspection, la chose qui ressort immédiatement est que la sensibilité est plus faible que prévu. Tous les pipelines fonctionnent à peu près au même niveau : ils identifient la plupart de leurs variants lorsqu’une profondeur d’environ 150x a été atteinte, ce qui indique que cette profondeur est probablement suffisante et que le nombre de variants manquants est probablement dû au fait que certains exons ont une couverture inférieure à la moyenne. Notez que quatre des cinq pipelines les plus performants ont GATK UnifiedGenotyper comme appelant de variants, ce qui démontre sa performance supérieure indépendamment de l’aligneur utilisé, comme le montre la figure 3(b).
(a)
(b)
(a)
(b)
En plus d’examiner les cinq meilleurs pipelines, nous avons déterminé qu’il serait utile d’effectuer la même analyse sur le meilleur aligneur couplé à chaque appelant de variante (Figure 3(a)), ainsi que sur le meilleur appelant de variante couplé à chaque aligneur (Figure 3(b)). Comme pour les pipelines, le meilleur aligneur a été identifié comme celui qui a produit le plus grand nombre de SNV TP et le plus petit nombre de SNV FP – dans ce cas, BWA mem. Bien que nous ayons le meilleur alignement sur lequel travailler, nous constatons toujours une différence assez importante entre les appelants de variantes, ce qui est probablement attribuable aux différents algorithmes qu’ils utilisent (Figure 3(a)). Cependant, dans le cas de l’appelant de variante le plus performant, GATK UnifiedGenotyper, il semble y avoir moins de variation parmi les quatre meilleurs aligneurs, ce qui indique qu’il est assez performant dans la plupart des situations, les exceptions étant CUSAHW3 et MOSAIK.
3.6. Variantes partagées entre les meilleurs pipelines
En dernier lieu, nous avons voulu savoir à quel point les ensembles d’appels de variantes étaient uniques entre les différents pipelines. Pour ce faire, nous nous sommes à nouveau concentrés sur les cinq meilleurs pipelines d’appel de variants : Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper, et Novoalign plus UnifiedGenotyper. Comme on peut le voir dans la figure 4, il y a une grande quantité de chevauchement entre les cinq pipelines différents en question, avec 15 489 SNV (70 %) partagés sur un total de 22 324 variants distincts. Cependant, on pourrait également affirmer que cela est dû en grande partie au fait que quatre des cinq pipelines utilisent le UnifiedGenotyper comme appelant de variants. Cette notion est corroborée par le fait que le plus grand nombre de variantes uniques à un pipeline, 367, appartient à la combinaison BWA sampe plus HaplotypeCaller. Il est également intéressant de noter que le deuxième plus grand nombre de SNV uniques appartient également à l’aligneur BWA sampe, il est donc possible que le nombre élevé de SNV uniques soit mieux attribué à l’aligneur qu’à l’appelant de variants.
4. Conclusions
Nous avons constaté que parmi les trente pipelines différents testés, Novoalign plus GATK UnifiedGenotyper présentait la plus grande sensibilité tout en maintenant un faible nombre de FP pour les SNV. Parmi les aligneurs, BWA mem a toujours obtenu les meilleures performances, mais les résultats varient encore beaucoup en fonction de l’appelant de variante utilisé. Naturellement, il s’ensuit que le meilleur appelant de variante, GATK UnifiedGenotyper, a produit des résultats similaires indépendamment de l’aligneur utilisé. Cependant, il est évident que les indels sont toujours difficiles à gérer pour n’importe quel pipeline, aucun d’entre eux n’atteignant une sensibilité moyenne supérieure à 33 % ou une VPP supérieure à 53 %. En plus de la faible performance globale que nous constatons dans la détection des indels, la sensibilité, quel que soit le type de mutation, est un problème pour tous les pipelines présentés dans cet article. Le nombre attendu de SNV pour l’exome de NA12878 est de 34 886, mais même en utilisant l’union de tous les variants identifiés par les cinq meilleurs pipelines, le plus grand nombre identifié était très faible (22 324). Il semble que, bien que toujours très utile, l’analyse de l’exome a ses limites, même lorsqu’il s’agit de quelque chose d’aussi apparemment simple que la détection de SNV.
Disclosure
Adam Cornish est un étudiant diplômé dans le laboratoire de Chittibabu Guda avec une formation en informatique et en génomique. Chittibabu Guda (professeur associé) a une formation interdisciplinaire en biologie moléculaire et computationnelle. Il a publié un certain nombre de méthodes informatiques avec une variété d’applications dans la recherche biomédicale, depuis 2001.
Conflit d’intérêts
Les auteurs ne sont pas au courant d’intérêts concurrents.
Contribution des auteurs
Adam Cornish a conçu l’étude, a effectué toutes les analyses, a fait les figures et a écrit le papier. Chittibabu Guda a fourni des commentaires essentiels sur les améliorations à apporter au document et des contributions sur les analyses elles-mêmes, et a révisé en profondeur le document. Tous les auteurs ont lu et approuvé le document final.
Remerciements
Ce travail a été soutenu par les fonds de développement à Chittibabu Guda de l’Université du Nebraska Medical Center (UNMC). Les auteurs tiennent à remercier les créateurs de Novoalign pour avoir mis le logiciel à disposition sous forme de version d’essai, ainsi que le Bioinformatics and Systems Biology Core facility de l’UNMC pour le soutien de l’infrastructure facilitant cette recherche.