Abstract

High-throughput sequencing, vooral van exomen, is een populair diagnostisch hulpmiddel, maar het is moeilijk om te bepalen welke hulpmiddelen het beste zijn in het analyseren van deze gegevens. In deze studie gebruiken we de NIST Genome in a Bottle resultaten als een nieuwe bron voor de validatie van onze exoom analyse pipeline. We gebruiken zes verschillende aligners en vijf verschillende variant callers om te bepalen welke pijplijn, van de in totaal 30, het beste presteert op een menselijk exoom dat werd gebruikt om de lijst van varianten te genereren die door het Genome in a Bottle Consortium werden gedetecteerd. Van deze 30 pijplijnen vonden we dat Novoalign in combinatie met GATK UnifiedGenoty de hoogste gevoeligheid vertoonde met behoud van een laag aantal fout-positieven voor SNV’s. Het is echter duidelijk dat indels voor geen enkele pijplijn nog steeds moeilijk te hanteren zijn, aangezien geen van de tools een gemiddelde gevoeligheid van meer dan 33% of een Positive Predictive Value (PPV) van meer dan 53% behaalde. Tenslotte werd, zoals verwacht, gevonden dat aligners een even vitale rol kunnen spelen in variant detectie als variant callers zelf.

1. Achtergrond

In de afgelopen jaren is er veel vooruitgang geboekt met high-throughput sequencing technologieën. Dankzij deze vooruitgang is het nu mogelijk een groot aantal potentiële ziekteveroorzakende varianten op te sporen en in enkele gevallen zijn sequencinggegevens van de volgende generatie (NGS) zelfs voor diagnostische doeleinden gebruikt. Dit is gedeeltelijk te danken aan de ontwikkelingen in sequencing-technologieën in de afgelopen jaren, maar ook aan het aantal verbeteringen dat is aangebracht in de verschillende bio-informaticahulpmiddelen die worden gebruikt om de bergen gegevens te analyseren die door NGS-instrumenten worden geproduceerd .

Bij het zoeken naar mutaties in een patiënt bestaat een typische workflow erin hun exoom te sequencen met een Illumina-sequencer, de ruwe gegevens uit te lijnen op het menselijke referentiegenoom, en vervolgens enkelvoudige nucleotide-varianten (SNV’s) of korte invoegingen en verwijderingen (indels) te identificeren die mogelijk het fenotype van interesse kunnen veroorzaken of beïnvloeden. Hoewel dit vrij eenvoudig is, is het niet eenvoudig om voor elke fase van de analysepijplijn te bepalen welke hulpmiddelen het best kunnen worden gebruikt. Er zijn een groot aantal tools die in verschillende tussenstappen worden gebruikt, maar de twee belangrijkste stappen in het hele proces zijn het uitlijnen van de ruwe reads op het genoom en vervolgens het zoeken naar varianten (d.w.z. SNVs en indels) . In deze studie willen we de bioinformaticus van vandaag helpen door de juiste combinatie van short read alignment tool en variant calling tool op te helderen voor het verwerken van exoom sequencing data geproduceerd door NGS instrumenten.

Er zijn in het verleden een aantal van deze studies uitgevoerd, maar ze hadden allemaal nadelen van de een of andere vorm. Idealiter zou men een lijst moeten hebben van elke bekende variant in een monster, zodat wanneer een pijplijn van analyse-instrumenten wordt uitgevoerd, men deze kan testen om met zekerheid te weten of zij correct werkt. In het verleden bestond zo’n lijst echter niet, zodat de validatie met minder volledige methoden moest worden uitgevoerd. In sommige gevallen werd de validatie uitgevoerd door gesimuleerde gegevens te genereren, zodat een reeks bekende true positives (TP) en true negatives (TN) ontstond. Hoewel dit gemakkelijk een lijst van elke TP en TN in de dataset oplevert, geeft het de biologie niet accuraat weer. Andere methoden om pijplijnen voor het aanroepen van varianten te valideren zijn het gebruik van genotyperingsarrays of Sanger-sequencing om een lijst van TP’s en vals-positieven (FP) te verkrijgen. Deze methoden hebben het voordeel dat ze biologisch gevalideerde resultaten opleveren, maar ze hebben ook het nadeel dat ze niet alomvattend zijn vanwege het beperkte aantal spots op genotyperingsarrays en de prohibitieve kosten van Sanger-validatie wanneer die duizenden keren wordt uitgevoerd. Tenslotte was geen van deze studies gericht op het effect van de short read aligner op het oproepen van varianten. Bijgevolg kon het stroomopwaartse effect van aligner prestaties niet onafhankelijk worden beoordeeld.

In deze studie hebben we het voordeel van een lijst van varianten voor een anonieme vrouw uit Utah (subject ID: NA12878, oorspronkelijk gesequenced voor het 1000 Genomes project ) die experimenteel werd gevalideerd door het NIST-geleide Genome in a Bottle (GiaB) Consortium. Deze lijst van varianten werd gecreëerd door het integreren van 14 verschillende datasets van vijf verschillende sequencers, en het stelt ons in staat om elke lijst van varianten gegenereerd door onze exoom analyse pijplijnen te valideren. De nieuwigheid van dit werk is het valideren van de juiste combinatie van aligners en variant callers tegen een uitgebreide en experimenteel bepaalde variant dataset: NIST-GiaB.

Om onze analyse uit te voeren zullen we gebruik maken van een van de exoom datasets die oorspronkelijk gebruikt zijn om de NIST-GiaB lijst te maken. We hebben slechts één van de oorspronkelijke Illumina TruSeq-gegenereerde exomen gekozen omdat we een standaard use-casescenario wilden bieden voor iemand die NGS-analyse wil uitvoeren, en terwijl whole genome sequencing in prijs blijft dalen, is exome sequencing nog steeds een populair en levensvatbaar alternatief . Het is ook belangrijk op te merken dat, volgens Bamshad et al., het verwachte aantal SNV’s per Europees-Amerikaans exoom momenteel 20 283 ± 523 is. Desondanks bedroeg het totale aantal SNV’s dat in de NIST-GiaB-lijst werd gevonden en potentieel aanwezig is in de TruSeq-exoomdataset 34 886, wat aanzienlijk hoger is dan verwacht. Dit is waarschijnlijk te wijten aan het feit dat, terwijl de exoomkit werd gebruikt om NIST-GiaB-gegevens te genereren, deze ook werd aangevuld met whole genome sequencing.

Ten slotte hebben we een groot aantal aligners en variant callers overwogen, maar uiteindelijk hebben we de 11 tools gekozen op basis van prevalentie, populariteit en relevantie voor onze dataset (zo werden SNVMix, VarScan2 en MuTect niet gebruikt omdat ze zijn bedoeld voor gebruik op tumor-afgeleide monsters). Onze analyse zelf omvat een vergelijking van zes aligners (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK , en Novoalign) en vijf variant callers (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup , en SNPSVM ). In deze studie proberen we ook te bepalen hoeveel van een effect, indien aanwezig, de aligner heeft op variant calling en welke aligners het beste presteren bij gebruik van een normale Illumina exoom monster. Voor zover wij weten, is dit het eerste rapport dat alle mogelijke combinaties (totaal van 30 pijplijnen) van een breed scala van aligners en variant callers valideert.

2. Methoden

2.1. Datasets

Humaan referentie genoom hg19 werd gedownload van de UCSC browser (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) en werd gebruikt om de uitlijningen uit te voeren. Het menselijk exoom, SRR098401, werd gedownload van de Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Voor annotatie- en kalibratiedoeleinden werden dbSNP137 zonder sites na versie 129, HapMap 3.3, Human Omni 2.5 BeadChip, en Mills en 1000 G gold standard indel set lists gebruikt (alle van ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. De pijplijn

Figuur 1 toont de workflow gebruikt in deze studie, die vergelijkbaar is met die geschetst in de Best Practices gids geproduceerd door The Broad Institute . Dit omvat een aantal stappen om ervoor te zorgen dat de uitlijning bestanden geproduceerd zijn van de hoogste kwaliteit, alsmede een aantal meer om te garanderen dat de varianten correct worden genoemd. Eerst werden ruwe leest uitgelijnd met hg19, en vervolgens PCR duplicaten werden verwijderd uit de uitlijning. Om later in de pijplijn te helpen bij de identificatie van indels, werd de uitlijning rond de indels uitgevoerd. De laatste stap in de uitlijning was het herijken van de score van de basiskwaliteit, wat helpt om de inherente bias en onnauwkeurigheden van de scores van sequencers te verbeteren. Helaas, ondanks deze stappen, de uitlijning tarief van elke aligner was aanzienlijk lager dan verwacht, dus om dit te compenseren, werd de fastx toolkit gebruikt om te filteren uit lage kwaliteit leest (tabel 1). Lezingen van lage kwaliteit werden gedefinieerd als die lezingen die ten minste de helft van hun kwaliteit scores onder de 30 hadden. Na alignment verwerking, variant calling en variant filtering werden uitgevoerd.

Figuur 1

Schema van de data-analyse pijplijn gebruikt. Om ervoor te zorgen dat de hoogste kwaliteit alignments worden gemaakt, worden gelezen eerst gefilterd en vervolgens uitgelijnd met het menselijk referentiegenoom, hg19. Vervolgens PCR duplicaten worden verwijderd, gelezen worden uitgelijnd rond vermoedelijke indels, en base kwaliteit scores worden opnieuw gekalibreerd. Ten slotte worden varianten genoemd en gevalideerd tegen de NIST-GiaB lijst van varianten.

De zes hulpmiddelen die zijn gebruikt om uitlijningen te genereren waren Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK, en Novoalign, en de vijf hulpmiddelen die zijn gebruikt om varianten te genereren waren FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup, en SNPSVM, zoals te zien is in tabel 2.

2.3. Filtering

Ruwe gegevens werden verkregen van de SRA (SRR098401), gesplitst met fastq-dump, en gefilterd met behulp van de fastx toolkit. Meer specifiek werden bij fastq-dump de -split_files en -split_spot flags gebruikt, en fastq_quality_filter werd uitgevoerd met de volgende argumenten -Q 33 -q 30 -p 50. Daarna werden de gelezen correct gepaird met een aangepast script.

2.4. Aligning

Aligners gebruikten standaard argumenten behalve wanneer een threads argument werd gebruikt waar beschikbaar. De gebruikte commando’s zijn als volgt.

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 genoom/hg19.fa uitlijningen/NA12878.R1.sai uitlijningen/NA12878.R2.sai raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > uitlijningen/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 > alignments/NA12878.bwa-mem.sam

2.4.4. CUSHAW3

(1)cushaw3 align -r $INDEX -t 10 -o alignments/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. Alignment Depth of Coverage Calculation

Om een correcte berekening van de dekkingsdiepte te verzekeren, werd de Picard Tools module CalculateHsMetrics gebruikt met de volgende argumenten:(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.bedHet is belangrijk op te merken dat het TruSeq exoombedbestand de header van het SAM-uitlijningsbestand moet hebben om deze module te laten werken. Verder moet kolom 6 worden verplaatst naar kolom 4, en kolom 5 moet worden verwijderd uit de TruSeq bed file.

2.6. Alignment File Processing

Verwerking van de alignment files (SAM/BAM files) vereist de volgende stappen voor alle aligners:(1)SAM naar BAM conversie met SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)BAM-bestand sorteren met behulp van de Picard Tools module, SortSam:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam SORT_ORDER=coordinate(3)PCR duplicaatverwijdering met de 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=alignments/metrics(4)Leesgroep toegevoegd aan alignmentbestanden met de module Picard Tools, AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=coordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Uitlijning rond indels met behulp van de GATK-modules RealignerTargetCreator en 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)Base recalibratie met behulp van de GATK modules BaseRecalibrator en PrintReads:(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. Varianten aanroepen

De standaard argumenten werden gebruikt voor elke variant aanroeper, tenzij deze een “threads” of “parallel” vlag bevatte, in welk geval dat ook werd gebruikt. Bovendien werden indels waar mogelijk apart van SNV’s aangeroepen. De gebruikte commando’s zijn als volgt:

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.vcfDoor het ontbreken van de vereiste CIGAR vlaggen in het uitlijningsbestand, slaagde SNPSVM er niet in varianten aan te roepen voor CUSHAW3, en SAMtools mpileup kon om dezelfde reden geen varianten aanroepen op MOSAIK uitlijningen. Omdat SNPSVM alleen SNV’s detecteert, werden er ook geen indels gerapporteerd voor dit programma.

2.8. Variant Filtratie

Filtratie varieerde afhankelijk van de variant caller die werd gebruikt. In het geval van GATK HaplotypeCaller en GATK UnifiedGenotyper, werden de GATK modules, VariantRecalibrator en ApplyRecalibration, gebruikt om SNV’s te filteren met behulp van HapMap 3.3, de Omni 2.5 SNP BeadChip, en dbSNP 137 zonder 1000 Genome data als training sets. Voor SNPSVM werden QUAL scores ≥ 4 en DP waarden ≥ 6 gebruikt. Voor FreeBayes en SAMtools werden QUAL-scores ≥ 20 en DP-waarden ≥ 6 gebruikt.

2.9. Variantenvergelijking

Voor variantenvergelijking werd USeq 8.8.1 gebruikt om SNV’s te vergelijken die tussen alle datasets werden gedeeld. Voor de vergelijking van indels werd het vcflib-hulpprogramma vcfintersect gebruikt. Het TruSeq hg19 exoombedbestand truseq_exome_targeted_regions.hg19.bed.chr, verkregen in 11 december 2013, werd gebruikt om de vergelijkingen te beperken tot locaties die konden worden vastgelegd door de exoom pull down kit gebruikt bij de sequencing van SRR098401. Dit bestand kan hier van Illumina worden verkregen: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Om ervoor te zorgen dat varianten identiek vertegenwoordigd waren tussen verschillende call sets, werd de vcflib tool vcfallelicprimitives gebruikt om vcf files.preprocess

2.10. Statistische berekeningen

True Positive (TP). Dit is een mutatie die door de geteste pijplijn is gedetecteerd en die voorkomt in de NIST-GiaB-lijst.

False Positive (FP). Dit is een mutatie die door de geteste pijplijn is gedetecteerd, maar die niet voorkomt in de NIST-GiaB-lijst.

True Negative (TN). Dit is een mutatie die niet werd gedetecteerd door de geteste pijplijn en die niet voorkomt in de NIST-GiaB-lijst.

False Negative (FN). Dit is een mutatie die niet door de geteste pijplijn werd gedetecteerd, maar die wel in de NIST-GiaB-lijst voorkomt:

3. Resultaten en discussie

3.1. Bij het uitvoeren van variantenanalyses is een van de vele valkuilen waarmee rekening moet worden gehouden de exoomsequentieruimte (zoals gedefinieerd door de exoom capture kit) en hoe die de analyseresultaten kan beïnvloeden. In dit geval hadden we één exoom (SRR098401) dat met de Illumina TruSeq exoomkit was geëxtraheerd en op een HiSeq 2000 was gesequeneerd. Met dit in gedachten wilden we ervoor zorgen dat we het vermogen van de bio-informatische hulpmiddelen om hun werk te doen en niet hoe goed de Illumina TruSeq exoom capture kit werkte te meten. Dat wil zeggen, we willen alleen proberen om varianten die worden verondersteld aanwezig te zijn in de exonen zoals gedefinieerd door de pull down kit noemen. Om deze reden gebruiken we het bed file geleverd door Illumina, niet een generieke annotatie bed file, bijvoorbeeld, RefSeq voor hg19. We vonden dat voor dit specifieke individu, volgens de NIST-GiaB lijst, moet er een totaal van 34.886 SNV’s en 1.473 indels binnen de regio’s gedefinieerd door de TruSeq bed file.

Zodra we uitgefilterd varianten die niet waren gelegen in de regio’s gedefinieerd door de Illumina TruSeq exoom bed file, gingen we van honderdduizenden vermeende varianten (gegevens niet weergegeven) tot, gemiddeld, ongeveer 23.000 varianten (SNV’s en indels) per pijplijn (tabel 3). Dit is een belangrijke stap voor onderzoekers om mee te beginnen, omdat het de zoekruimte voor potentieel interessante varianten aanzienlijk vermindert.

3.2. Een aspect dat we wilden begrijpen bij het maken van deze vergelijking was het belang van het filteren van varianten die door deze tools werden gedetecteerd. De reden hiervoor is dat men idealiter een zo hoog mogelijk gevoeligheidsniveau zou willen hebben, zodat de mutaties van belang niet verloren gaan in het filterproces. Daarom moet worden nagegaan of en in hoeverre deze stap nodig is, aangezien uit de NIST-GiaB-resultaten en de review van Bamshad et al. duidelijk blijkt dat gevoeligheid een probleem kan zijn.

Zoals uit tabel 3 blijkt, is filtering voor sommige variant callers meer nodig dan voor andere als het gaat om SNV’s, en voor indels is filtering absoluut noodzakelijk. In de meeste gevallen ligt het aantal TP-varianten dicht bij of hoger dan ons verwachte aantal van ongeveer 20.000, maar anderzijds is in sommige gevallen het aantal FP’s zeer hoog.

Het is duidelijk dat er veel variatie is in de aantallen die door elke pijplijn worden gegenereerd. Toch zijn er enkele overeenkomsten te vinden in de getallen die waarschijnlijk voortkomen uit de algoritmische oorsprong van elk gereedschap. FreeBayes produceert zowel het grootste aantal ongefilterde varianten als het grootste aantal FPs. Waarschijnlijk zien we dit soort prestaties alleen bij dit gereedschap omdat het, hoewel het niet het enige varianten-oproepprogramma is dat gebaseerd is op Bayesiaanse inferentie, uniek is in zijn interpretatie van alignments. Dat wil zeggen, het is een haplotype-gebaseerde caller die varianten identificeert op basis van de volgorde van de gelezen data zelf in plaats van op basis van de uitlijning, wat de manier is waarop GATK’s UnifiedGenotyper werkt.

Daarnaast zien we dat de Burrows-Wheeler gebaseerde aligners zeer vergelijkbaar met elkaar presteren: Bowtie2, BWA mem, en BWA sampe behalen over de hele linie vergelijkbare resultaten. Men zou kunnen veronderstellen dat dit waarschijnlijk te wijten is aan het feit dat al deze tools gelijkaardige algoritmes gebruiken bij het uitvoeren van hun taak. Deze observatie wordt ondersteund door het feit dat MOSAIK (gapped alignments met behulp van het Smith-Waterman algoritme) en CUSHAW3 (een hybride seeding benadering) beide zeer verschillende onderliggende algoritmen hebben en daardoor zeer verschillende resultaten produceren.

Dit verschil in resultaten die correleren met verschillende algoritmen is het best te zien in de SNPSVM resultaten. Van de variant-callers is het de enige die gebruik maakt van support vector machines en modelbouw om SNV-calls te genereren. Het lijkt erop dat deze methode, hoewel ze het nadeel heeft niet zo gevoelig te zijn als andere methoden, het voordeel heeft dat ze uiterst accuraat is, ongeacht de aligner die wordt gebruikt. Dit suggereert dat men in staat is om de filterstap helemaal over te slaan bij gebruik van deze variant caller.

Met betrekking tot indels lijkt geen enkele aligner zich te onderscheiden van de rest als een die goed met dit type mutatie omgaat. Als we kijken naar het aantal FP’s, is het zelfs duidelijk dat de variantaanroeper de grootste rol speelt bij de nauwkeurigheid van de indel-identificatie. Bovendien zijn er noch voor CUSHAW3 plus SNPSVM noch voor MOSAIK plus SAMtools mpileup-pijplijnen gegevens beschikbaar, omdat de alignmentbestanden niet de noodzakelijke CIGAR-strings bevatten om de variant-callers downstream te laten functioneren. Ten slotte is de reden dat er geen indel-gegevens voor SNPSVM zijn, omdat dit hulpmiddel uitsluitend wordt gebruikt voor de identificatie van SNV’s.

3.3. Gefilterde resultaten vergeleken met GiaB

Zoals in tabel 4, slagen standaard filterpraktijken erin om een groot aantal FP SNV’s voor elke pijplijn te verwijderen; het lijkt er echter op dat deze filters in de meeste gevallen een beetje te agressief zijn voor SNV-detectie, maar niet streng genoeg voor indels. Dit wordt duidelijk als we kijken naar de verschillen in het aantal FPs dat in elke dataset wordt gerapporteerd. Bij Bowtie2 met Freebayes worden bijvoorbeeld 72.570 FP SNV’s verwijderd (een reductie van 98%), maar slechts 1.736 FP indels (een reductie van 70%). Ook moet worden opgemerkt dat de gebruikte filters pijplijn-afhankelijk waren en, voor het grootste deel, binnen elke pijplijn vergelijkbare reducties in SNV en indel FPs opleverden. De enige uitzondering hierop is het aantal varianten dat werd geïdentificeerd in de CUSHAW3-uitlijningen in vergelijking met andere uitlijningen: Over het geheel genomen is het aantal TP SNV’s lager, het aantal FP SNV’s hoger, en het is de enige aligner die na filtering meer dan 1.000 FP indels oplevert.

Gezien het feit dat filteren het aantal TP-varianten aanzienlijk vermindert, zou het verstandig kunnen zijn om, met uitzondering van pijplijnen die gebruik maken van CUSHAW3 en FreeBayes, deze stap over te slaan bij het zoeken naar zeldzame, high-impact varianten. In plaats daarvan zou men meer tijd kunnen besteden aan een filterproces dat meer gebaseerd is op biologie dan op statistiek. Het kan bijvoorbeeld zinvoller zijn om tijd te investeren in het identificeren van een kleine lijst van varianten die waarschijnlijk high-impact zijn: splice site mutaties, indels die frameshifts veroorzaken, truncatiemutaties, stop-loss mutaties, of mutaties in genen waarvan bekend is dat ze biologisch relevant zijn voor het fenotype van interesse.

3.4. Gemiddelde TPR en gevoeligheid

Zoals blijkt uit tabel 5 varieert de Positive Predictive Value (PPV) voor elk instrument, met uitzondering van CUSHAW3, van 91% tot 99,9% voor SNV’s, maar is de gemiddelde gevoeligheid zeer laag (rond 50%). Deze discrepantie kan verschillende oorzaken hebben, maar de meest waarschijnlijke is een variabele dekkingsdiepte over de exonen. We zien dat, naast de lage SNV-gevoeligheid, ook de indel-gevoeligheid laag is (rond 30%); de PPV voor indels is echter aanzienlijk lager (35,86% tot 52,95%). Dit kan aan een van de volgende redenen te wijten zijn: zeer korte indels zijn moeilijk te detecteren met conventionele NGS, de weergave van indels door verschillende variant callers kan ertoe leiden dat tools ten onrechte beweren dat twee indels verschillend zijn, of alignment tools produceren verschillende weergaven van dezelfde indel .

Misschien is de meest waarschijnlijke verklaring voor beide soorten mutaties de kwestie van de diepte. Zoals het geval is met elke variant analyse studie, een toename van de diepte van de dekking leidt tot een toename van de gevoeligheid, maar het is onmogelijk om een goede diepte van de dekking te garanderen als gevolg van het onvermogen van exoom capture kits om uniform pull down exons . Bovendien dekt geen enkele exome capture kit elk exon. Er is inderdaad aangetoond dat variantenanalyse van sequencing van het volledige genoom bij een gemiddelde diepte die kleiner is dan een exoom, beter presteert omdat die diepte uniform is. Het is dus waarschijnlijk dat een groot aantal varianten ontbreekt doordat de NIST-GiaB-lijst is samengesteld uit een compilatie van exoom- en genoomsequencinggegevens. Om de juiste gevoeligheid te bereiken zal men uiteindelijk whole genome sequencing moeten uitvoeren, maar dat is momenteel voor de meeste laboratoria een te dure aangelegenheid. Gelukkig dalen die kosten steeds verder en zullen we binnenkort een geleidelijke verschuiving zien van exoomanalyse naar de completere analyse van het volledige genoom.

3.5. Gevoeligheid als functie van diepte

Omdat gevoeligheid een van de belangrijkste prestatiemaatstaven van een tool is en de meeste tools moeite hebben om een gevoeligheid van meer dan 50% te bereiken, wilden we verder onderzoeken hoe diepte de gevoeligheid voor het oproepen van varianten beïnvloedde. We hebben gekeken naar een aantal verschillende combinaties van tools om te bepalen wat de beste pipelines, variant callers en aligners waren. Voor Figuur 2 hebben we de vijf beste combinaties van variant callers en aligners genomen, zoals bepaald door hun gevoeligheid en fout-positieve rate (FPR). Dat wil zeggen, we hebben die combinaties geselecteerd met het hoogste aantal TP SNV’s dat werd genoemd en het laagste aantal FP SNV’s. Wat bij inspectie onmiddellijk opvalt, is dat de gevoeligheid lager is dan verwacht. Alle pijplijnen presteren op ongeveer hetzelfde niveau: ze identificeren de meeste varianten tegen de tijd dat een diepte van ongeveer 150x bereikt is, wat erop wijst dat deze diepte waarschijnlijk voldoende is en dat het aantal ontbrekende varianten waarschijnlijk te wijten is aan het feit dat bepaalde exonen een lagere dan gemiddelde dekking hebben. Merk op dat vier van de vijf best presterende pijplijnen hebben GATK UnifiedGenotyper als hun variant caller, waaruit blijkt dat de superieure prestaties, ongeacht de aligner gebruikt zoals getoond in figuur 3 (b).

Figuur 2

Gevoeligheid als functie van de diepte voor de top vijf pijplijnen. De top vijf pijplijnen worden hier getoond met de diepte van elke SNV uitgezet tegen de gevoeligheid.

(a)

(b)

(a)
(b)

Naast het bekijken van de top vijf pijplijnen, besloten we dat het nuttig zou zijn om dezelfde analyse uit te voeren op de beste aligner gekoppeld aan elke variant caller (figuur 3(a)), evenals de beste variant caller gekoppeld aan elke aligner (figuur 3(b)). Net als bij de pijplijnen werd de beste aligner geïdentificeerd als die welke het hoogste aantal TP SNV’s en het laagste aantal FP SNV’s opleverde – in dit geval BWA mem. Ondanks de beste uitlijning om mee te werken, zien we nog steeds een vrij groot verschil tussen de variant-callers, wat waarschijnlijk toe te schrijven is aan de verschillende algoritmes die ze gebruiken (Figuur 3(a)). Echter, in het geval van de best presterende variant caller, GATK UnifiedGenotyper, lijkt er minder variatie te zijn tussen de top vier aligners, wat aangeeft dat het redelijk goed presteert in de meeste situaties, met als uitzonderingen CUSAHW3 en MOSAIK.

3.6. Gedeelde varianten tussen de top pijplijnen

Ten slotte wilden we weten hoe uniek de sets variante oproepen waren tussen de verschillende pijplijnen. Om dit te doen, richtten we ons opnieuw op de top vijf van pijplijnen voor het oproepen van varianten: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper, en Novoalign plus UnifiedGenotyper. Zoals kan worden gezien in Figuur 4, is er een grote mate van overlap tussen de vijf verschillende pijplijnen in kwestie, met 15.489 SNVs (70%) gedeeld op een totaal van 22.324 verschillende varianten. Men zou echter ook kunnen argumenteren dat dit grotendeels te wijten is aan het feit dat vier van de vijf pijplijnen de UnifiedGenotyper gebruiken als hun variant caller. Dit idee wordt bevestigd door het feit dat het grootste aantal varianten uniek voor een pijplijn, 367, behoort tot de BWA sampe plus HaplotypeCaller combinatie. Het is ook vermeldenswaard dat het op een na hoogste aantal unieke SNV’s ook behoort tot de BWA sampe aligner, dus het is mogelijk dat het hoge aantal unieke SNV’s beter kan worden toegeschreven aan de aligner dan aan de variant caller.

Figuur 4

Het snijpunt van de SNV’s geïdentificeerd door de top vijf pijplijnen.

4. Conclusies

Wij vonden dat van de dertig verschillende geteste pijplijnen Novoalign plus GATK UnifiedGenotyper de hoogste gevoeligheid vertoonden met behoud van een laag aantal FP voor SNV’s. Van de aligners presteerde BWA mem consequent het beste, maar de resultaten varieerden nog steeds sterk, afhankelijk van de gebruikte variant caller. Uiteraard volgt hieruit dat de beste variant caller, GATK UnifiedGenotyper, meestal vergelijkbare resultaten opleverde, ongeacht de gebruikte aligner. Het is echter duidelijk dat indels nog steeds moeilijk te verwerken zijn door welke pijplijn dan ook, aangezien geen van de pijplijnen een gemiddelde gevoeligheid van meer dan 33% of een PPV van meer dan 53% behaalde. Naast de lage algemene prestaties die we zien bij het detecteren van indels, is de gevoeligheid, ongeacht het mutatietype, een probleem voor elke pijplijn die in dit artikel wordt beschreven. Het verwachte aantal SNV’s voor het exoom van NA12878 is 34.886, maar zelfs bij gebruik van de unie van alle varianten geïdentificeerd door de top vijf pijplijnen, was het grootste aantal geïdentificeerd zeer laag (22.324). Het lijkt erop dat hoewel nog steeds zeer nuttig exoomanalyse zijn beperkingen heeft, zelfs als het gaat om zoiets schijnbaar eenvoudigs als SNV-detectie.

Disclosure

Adam Cornish is een afgestudeerde student in het lab van Chittibabu Guda met een opleiding in computerwetenschappen en genomica. Chittibabu Guda (universitair hoofddocent) heeft een interdisciplinaire achtergrond in moleculaire en computationele biologie. Hij heeft een aantal computationele methoden gepubliceerd met een verscheidenheid aan toepassingen in biomedisch onderzoek, sinds 2001.

Conflict of Interests

De auteurs zijn zich niet bewust van concurrerende belangen.

Authors’ Contribution

Adam Cornish ontwierp de studie, voerde alle analyses uit, maakte de figuren, en schreef de paper. Chittibabu Guda gaf essentiële feedback over verbeteringen aan het artikel en input over de analyses zelf en redigeerde het artikel grondig. Alle auteurs lazen en keurden de definitieve paper goed.

Acknowledgments

Dit werk werd ondersteund door de ontwikkelingsfondsen aan Chittibabu Guda van de Universiteit van Nebraska Medical Center (UNMC). De auteurs willen de makers van Novoalign bedanken voor het beschikbaar stellen van de software als proefversie en de Bioinformatics and Systems Biology Core faciliteit van UNMC voor de infrastructurele ondersteuning die dit onderzoek mogelijk maakte.