Abstract

Höggenomströmningssekvensering, särskilt av exomer, är ett populärt diagnostiskt verktyg, men det är svårt att avgöra vilka verktyg som är bäst på att analysera dessa data. I den här studien använder vi NIST Genome in a Bottle-resultaten som en ny resurs för validering av vår pipeline för exomanalys. Vi använder sex olika aligners och fem olika variant callers för att avgöra vilken pipeline, av de totalt 30, som presterar bäst på ett mänskligt exom som användes för att hjälpa till att generera listan över varianter som upptäckts av Genome in a Bottle-konsortiet. Av dessa 30 pipelines fann vi att Novoalign i kombination med GATK UnifiedGenotyper uppvisade den högsta känsligheten samtidigt som vi bibehöll ett lågt antal falskt positiva SNV:er. Det är dock uppenbart att indels fortfarande är svåra att hantera för någon pipeline, eftersom inget av verktygen uppnådde en genomsnittlig känslighet högre än 33 % eller ett positivt prediktivt värde (PPV) högre än 53 %. Slutligen konstaterades det som väntat att aligners kan spela en lika viktig roll för upptäckt av varianter som variant callers själva.

1. Bakgrund

Under de senaste åren har det gjorts många framsteg inom tekniken för sekvensering med hög genomströmning. Tack vare dessa framsteg är det nu möjligt att upptäcka ett stort antal potentiella sjukdomsorsakande varianter , och i några fall har data från nästa generations sekvensering (NGS) till och med använts för diagnostiska ändamål . Detta beror delvis på utvecklingen av sekvenseringstekniken under de senaste åren, men också på det antal förbättringar som gjorts av de olika bioinformatiska verktyg som används för att analysera de berg av data som produceras av NGS-instrument .

När man söker efter mutationer hos en patient är ett typiskt arbetsflöde att sekvensera deras exom med en Illumina-sekvenser, anpassa rådata till det mänskliga referensgenomet och sedan identifiera enskilda nukleotidvarianter (SNV) eller korta insättningar och deletioner (indels) som eventuellt skulle kunna orsaka eller påverka den aktuella fenotypen . Detta är ganska okomplicerat, men att bestämma vilka verktyg som är bäst lämpade att använda i varje skede av analysen är det inte. Det finns ett stort antal verktyg som används i olika mellansteg, men de två viktigaste stegen i hela processen är att anpassa de obearbetade avläsningarna till genomet och sedan söka efter varianter (dvs. SNV:er och indels) . I den här studien syftar vi till att hjälpa dagens bioinformatiker genom att belysa den korrekta kombinationen av verktyg för anpassning av korta läsningar (short read alignment tool) och verktyg för att söka efter varianter (variant calling tool) för att bearbeta exomsekvenseringsdata som produceras av NGS-instrument.

Ett antal av dessa studier har utförts tidigare, men alla har haft nackdelar på ett eller annat sätt. Helst bör man ha en lista över varje känd variant som finns i ett prov, så att man när en pipeline av analysverktyg körs kan testa den för att med säkerhet veta att den fungerar korrekt. Tidigare fanns det dock ingen sådan förteckning, så valideringen fick utföras med mindre fullständiga metoder. I vissa fall utfördes valideringen genom att generera simulerade data för att skapa en uppsättning kända sant positiva (TP) och sant negativa (TN) resultat. Även om detta på ett bekvämt sätt ger en förteckning över alla TP och TN i datamängden, ger det en dålig bild av biologin på ett korrekt sätt. Andra metoder för att validera pipelines för variantframställning är att använda genotypningsarrayer eller Sanger-sekvensering för att få fram en lista över TP och falskt positiva (FP) . Dessa metoder har fördelen att de ger biologiskt validerade resultat, men de har också nackdelen att de inte är heltäckande på grund av det begränsade antalet fläckar på genotypningsarrayer och den oöverkomliga kostnaden för Sangervalidering när den utförs tusentals gånger. Slutligen var ingen av dessa studier inriktad på att undersöka vilken effekt short read aligner hade på variantbestämningen. Följaktligen kunde inte effekten uppströms av alignerns prestanda bedömas oberoende.

I den här studien har vi fördelen av en lista med varianter för en anonym kvinna från Utah (ämnes-ID: NA12878, ursprungligen sekvenserad för 1000 Genomes-projektet ) som har validerats experimentellt av det NIST-ledda konsortiet Genome in a Bottle (GiaB). Denna lista över varianter skapades genom att integrera 14 olika dataset från fem olika sekvenserare, och den gör det möjligt för oss att validera vilken lista över varianter som helst som genereras av våra pipelines för exomanalys . Det nya i detta arbete är att validera den rätta kombinationen av aligners och variant callers mot en omfattande och experimentellt bestämd datamängd med varianter: NIST-GiaB.

För att utföra vår analys kommer vi att använda ett av de exomdataset som ursprungligen användes för att skapa NIST-GiaB-listan. Vi valde endast ett av de ursprungliga Illumina TruSeq-genererade exomen eftersom vi ville tillhandahålla ett standardanvändningsscenario för någon som vill utföra NGS-analyser, och medan sekvensering av hela genomet fortsätter att sjunka i pris är exomsekvensering fortfarande ett populärt och livskraftigt alternativ . Det är också viktigt att notera att enligt Bamshad et al. är det förväntade antalet SNV:er per europeisk-amerikanskt exom för närvarande 20 283 ± 523 . Trots detta var det totala antalet SNV som hittades i NIST-GiaB-listan med potential att finnas i TruSeq-exomdatasetet 34 886, vilket är betydligt högre än förväntat. Detta beror sannolikt på att även om exomkitet användes för att generera NIST-GiaB-data kompletterades det också med helgenomsekvensering.

Slutningsvis övervägde vi ett stort antal aligners och variant callers men valde i slutändan de 11 verktygen baserat på prevalens, popularitet och relevans för vårt dataset (t.ex. användes inte SNVMix, VarScan2 och MuTect eftersom de är avsedda att användas på prover som härrör från tumörer). Själva analysen innebär att vi jämför sex aligners (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK och Novoalign) och fem variant callers (FreeBayes , GATK HaplotypeCaller , GATK UnifiedGenotyper , SAMtools mpileup och SNPSVM ). I den här studien försöker vi också fastställa hur stor effekt, om någon, aligneraren har på variant calling och vilka alignerare som presterar bäst när man använder ett normalt Illumina exomprov. Såvitt vi vet är detta den första rapporten som validerar alla möjliga kombinationer (totalt 30 pipelines) av ett brett utbud av aligners och variant callers.

2. Metoder

2.1. Datamängder

Det mänskliga referensgenomet hg19 hämtades från UCSC:s webbläsare (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) och användes för att utföra anpassningarna. Det mänskliga exomet, SRR098401, hämtades från Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). För annotering och kalibrering användes dbSNP137 utan platser efter version 129, HapMap 3.3, Human Omni 2.5 BeadChip samt Mills och 1000 G gold standard indel set lists (alla från ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. Pipeline

Figur 1 visar det arbetsflöde som användes i den här studien, som liknar det som beskrivs i den vägledning om bästa praxis som tagits fram av The Broad Institute . Detta inbegriper ett antal steg för att säkerställa att de alinjeringsfiler som produceras är av högsta kvalitet samt ytterligare flera steg för att garantera att varianterna kallas korrekt. Först anpassades råavläsningarna till hg19 och sedan avlägsnades PCR-duplikat från anpassningen. För att underlätta identifieringen av indels senare i pipelinen utfördes sedan en omjustering av avläsningar runt indels. Det sista steget i anpassningsbehandlingen var att utföra en omkalibrering av baskvalitetspoängen, vilket bidrar till att förbättra de inneboende biaserna och felaktigheterna i de poäng som utfärdas av sekvenserare. Trots dessa steg var tyvärr anpassningsgraden för varje alignerare betydligt lägre än förväntat, så för att kompensera detta användes verktygslådan fastx för att filtrera bort läsningar av låg kvalitet (tabell 1). Lågkvalitativa läsningar definierades som de läsningar som hade minst hälften av sina kvalitetspoäng under 30. Efter justeringsprocessen utfördes variantkallning och variantfiltrering.

Figur 1

Skematisk bild av den dataanalysledning som använts. För att se till att det skapas anpassningar av högsta kvalitet filtreras först läsningarna och anpassas sedan till det mänskliga referensgenomet hg19. Därefter tas PCR-duplikat bort, läsningar anpassas runt förmodade indels och baskvalitetsvärden omkalibreras. Slutligen kallas varianter och valideras mot NIST-GiaB-listan över varianter.

De sex verktyg som användes för att generera anpassningar var Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK och Novoalign, och de fem verktyg som användes för att generera varianter var FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup och SNPSVM, vilket framgår av tabell 2.

Tool Type Version Reference Bowtie2 Aligner 2.1.0 BWA sampe Aligner 0.7.5a BWA mem Aligner 0.7.5a CUSHAW3 Aligner 3.0.3 MOSAIK Aligner 2.2.3 Novoalign Aligner 3.02.07 N/A FreeBayes Genotyper v9.9.2-19-g011561f GATK HaplotypeCaller Genotyper 2.7-2 N/A GATK UnifiedGenotyper Genotyper 2.7-2 SAMtools mpileup Genotyper 0.1.19 SNPSVM Genotyper 0.01

Tabell 2

Detta är de elva olika verktyg som användes och som utgjorde de 30 olika pipelines (sex aligners fem variant callers). Programvaruversioner finns också med för att säkerställa reproducerbarhet.

2.3. Filtrering

Rå data förvärvades från SRA (SRR098401), delades med fastq-dump och filtrerades med hjälp av verktygslådan fastx. Specifikt använde fastq-dump flaggorna -split_files och -split_spot, och fastq_quality_filter kördes med följande argument: -Q 33 -q 30 -p 50. Därefter parades läsningarna korrekt med ett anpassat skript.

2.4. Aligning

Alignerna använde standardargument utom när ett argument för trådar användes när det fanns tillgängligt. De kommandon som används är följande.

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 alignments/NA12878.R1.sai alignments/NA12878.R2.sai raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > alignments/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.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 raw_data/read_2_filtered.fastq -o SAM -c 10 > alignments/NA12878.novoalign.sam

2.5. Beräkning av täckningsdjup för anpassningar

För att säkerställa en korrekt beräkning av täckningsdjupet användes Picard Tools-modulen CalculateHsMetrics med följande argument:(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.bedDet är viktigt att notera att TruSeq-exombäddfilen måste ha rubriken från SAM-justeringsfilen i förväg för att den här modulen ska fungera. Vidare måste kolumn 6 flyttas till kolumn 4 och kolumn 5 måste tas bort från TruSeq-bäddfilen.

2.6. Behandling av anpassningsfiler

Behandlingen av anpassningsfilerna (SAM/BAM-filer) krävde följande steg för alla anpassare:(1)SAM till BAM-konvertering med SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2) Sortering av BAM-filer med hjälp av Picard Tools-modulen 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)Borttagning av PCR-dubbletter med hjälp av Picard Tools-modulen 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)Läsgrupp som läggs till i anpassningsfiler med hjälp av Picard Tools-modulen AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=koordinat RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Justering runt indels med hjälp av GATK-modulerna RealignerTargetCreator och 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)Basåterkalibrering med hjälp av GATK-modulerna BaseRecalibrator och 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. Variantkallning

Standardargument användes för varje variantkallare såvida den inte innehöll en ”threads”- eller ”parallell”-flagga, i vilket fall detta också användes. Dessutom kallades indels separat från SNVs när det var möjligt. Specifikt används följande kommandon:

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 $D $DBSNP(2)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar.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.vcfPå grund av att det inte fanns nödvändiga CIGAR-flaggor i anpassningsfilen misslyckades SNPSVM med att kalla varianter för CUSHAW3, och SAMtools mpileup kunde inte kalla varianter för MOSAIK-anpassningar av samma anledning. På grund av att SNPSVM endast upptäcker SNV:er rapporterades inga indels för detta program.

2.8. Variantfiltrering

Filtreringen varierade beroende på vilken variantkälla som användes. I fallen GATK HaplotypeCaller och GATK UnifiedGenotyper användes GATK-modulerna VariantRecalibrator och ApplyRecalibration för att filtrera SNV:er med HapMap 3.3, Omni 2.5 SNP BeadChip och dbSNP 137 utan 1000 Genome-data som träningsuppsättningar. För SNPSVM användes QUAL-poäng ≥ 4 och DP-värden ≥ 6. För FreeBayes och SAMtools användes QUAL-poäng ≥ 20 och DP-värden ≥ 6.

2.9. Variantjämförelse

För variantjämförelse användes USeq 8.8.1 för att jämföra SNV:er som delas mellan alla dataset. För att jämföra indels användes vcflib-verktyget vcfintersect. TruSeq hg19 exome bed-filen truseq_exome_targeted_regions.hg19.bed.chr, som erhölls den 11 december 2013, användes för att begränsa jämförelserna till platser som kunde fångas upp av det exome pull down-kit som användes vid sekvenseringen av SRR098401. Denna fil kan erhållas från Illumina här: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. För att säkerställa att varianter representerades identiskt mellan olika call sets användes verktyget vcflib vcfallelicprimitives för att förbehandla vcf-filer.

2.10. Statistiska beräkningar

True Positive (TP). Det är en mutation som upptäcktes av den pipeline som testas och som finns i NIST-GiaB-listan.

Falskt positivt (FP). Det är en mutation som upptäcktes av den testade pipelinen men som inte finns i NIST-GiaB-listan.

True Negative (TN). Det är en mutation som inte upptäcktes av den testade pipelinen och som inte finns i NIST-GiaB-listan.

Falskt negativt (FN). Det är en mutation som inte upptäcktes av den testade pipelinen men som finns i NIST-GiaB-listan:

3. Resultat och diskussion

3.1. Förfiltrering av varianter

När man utför variantanalys är en av de många fallgroparna som måste beaktas exomsekvensutrymmet (enligt definitionen i exomfångstkitet) och hur det kan påverka analysresultaten. I det här fallet hade vi ett enda exom (SRR098401) som extraherades med hjälp av Illumina TruSeq exomkit och sekvenserades på en HiSeq 2000. Med detta i åtanke ville vi försäkra oss om att vi mätte bioinformatikverktygens förmåga att utföra sitt arbete och inte hur väl Illumina TruSeq exome capture kit fungerade. Det vill säga, vi vill bara försöka ringa in varianter som ska finnas i exonerna enligt definitionen i pull down-kitet. Därför använder vi den bäddfil som tillhandahålls av Illumina, inte en generisk bäddfil för annotation, till exempel RefSeq för hg19. Vi fann att det för denna specifika individ, enligt NIST-GiaB-listan, borde finnas totalt 34 886 SNV:er och 1 473 indels inom de regioner som definieras av TruSeq-bäddfilen.

När vi väl hade filtrerat bort varianter som inte fanns i de regioner som definierades av Illumina TruSeq exom-bäddfilen, gick vi från hundratusentals förmodade varianter (data visas inte) till i genomsnitt cirka 23 000 varianter (SNV:er och indels) per pipeline (tabell 3). Detta är ett viktigt steg för forskare att börja med, eftersom det avsevärt minskar sökutrymmet för potentiellt intressanta varianter.

3.2. Råvariantresultat jämfört med GiaB

En aspekt som vi ville förstå när vi gjorde denna jämförelse var vikten av att filtrera varianter som upptäckts av dessa verktyg. Anledningen till detta är att man idealt sett vill ha en så hög känslighetsnivå som möjligt så att de mutationer som är av intresse inte går förlorade i filtreringsprocessen. Det åligger oss därför att avgöra om detta steg är nödvändigt eller inte och i vilken grad det är nödvändigt, eftersom det framgår tydligt av NIST-GiaB-resultaten och Bamshad et al:s granskning att känsligheten kan vara ett problem.

Som vi kan se i tabell 3 behövs filtrering mer för vissa variant callers än för andra när det gäller SNV:er, och det är absolut nödvändigt för indels. I de flesta fall är antalet TP-varianter nära eller högre än vårt förväntade antal på cirka 20 000 , men å andra sidan är antalet FPs i vissa fall mycket högt.

Det är tydligt att det finns en stor variation i de siffror som genereras av varje pipeline. Man kan dock hitta vissa gemensamma drag i siffrorna som sannolikt härrör från varje verktygs algoritmiska ursprung. FreeBayes producerar både det största antalet ofiltrerade varianter och det högsta antalet FPs. Det är troligt att vi endast ser denna typ av prestanda från detta verktyg på grund av det faktum att det inte är det enda verktyg för att hitta varianter som bygger på Bayesiansk inferens, men att det är unikt i sin tolkning av anpassningar. Det vill säga, det är en haplotypsbaserad sökare som identifierar varianter baserat på sekvensen av själva läsningarna i stället för på anpassningen, vilket är hur GATK:s UnifiedGenotyper fungerar.

Därutöver ser vi att de Burrows-Wheeler-baserade anpassningsverktygen presterar mycket likt varandra: Bowtie2, BWA mem och BWA sampe uppnår liknande resultat över hela linjen. Man kan anta att detta sannolikt beror på att alla dessa verktyg använder liknande algoritmer när de utför sin uppgift. Denna observation stöds av det faktum att MOSAIK (gapped alignments med hjälp av Smith-Waterman-algoritmen) och CUSHAW3 (en hybridmetod för sådd) båda har mycket olika underliggande algoritmer och därefter ger mycket olika resultat.

Denna skillnad i resultat som korrelerar med olika algoritmer syns bäst i SNPSVM-resultaten. Av variantkallarna är det den enda som använder stödvektormaskiner och modellbygge för att generera SNV-kallarna. Det verkar som om den visserligen har nackdelen att den inte är lika känslig som andra metoder, men den gynnas av att den är extremt noggrann oavsett vilken alignerare som används. Detta tyder på att man kan hoppa över filtreringssteget helt och hållet när man använder denna variant caller.

Med avseende på indels tycks ingen alignerare sticka ut bland de övriga som en som hanterar denna typ av mutation väl. Faktum är att när man tittar på antalet FPs är det tydligt att det är variant caller som spelar störst roll för noggrannheten i indelidentifieringen. Dessutom finns det uppgifter för varken CUSHAW3 plus SNPSVM eller MOSAIK plus SAMtools mpileup pipelines på grund av att anpassningsfilerna inte innehåller de CIGAR-strängar som krävs för att variantkallarna ska fungera nedströms. Slutligen är anledningen till att det inte finns några indel-data för SNPSVM att detta verktyg endast används för identifiering av SNV:er.

3.3. Filtrerade resultat jämfört med GiaB

Som framgår av tabell 4 lyckas standardfiltreringsmetoderna ta bort ett stort antal FP-SNV för varje pipeline; det verkar dock som om dessa filter är lite för aggressiva i de flesta fall för SNV-detektering, men inte tillräckligt strikta för indels. Detta blir tydligt när man tittar på skillnaderna i antalet FPs som rapporteras i varje dataset. I Bowtie2 med Freebayes kan man t.ex. se en borttagning av 72 570 FP SNVs (en minskning med 98 %) men endast en borttagning av 1 736 FP indels (en minskning med 70 %). Det bör också noteras att de filter som användes var pipelineberoende och att de för det mesta inom varje pipeline gav liknande minskningar av SNV och indel FPs. Det enda undantaget här är antalet varianter som identifierades från CUSHAW3-anpassningarna jämfört med andra anpassningar: totalt sett är antalet TP SNV lägre, antalet FP SNV högre, och det är den enda alignerare som producerar mer än 1 000 FP indels efter filtrering.

Med tanke på att filtrering avsevärt minskar antalet TP-varianter kan det vara klokt att, med undantag för pipelines som använder CUSHAW3 och FreeBayes, hoppa över detta steg när man söker efter sällsynta varianter med stor inverkan. Istället kan man lägga mer tid på en filtreringsprocess som bygger på biologi snarare än statistik. Det kan till exempel vara vettigare att investera tid i att identifiera en liten lista över varianter som sannolikt kommer att ha stor inverkan: mutationer på skarvställen, indels som orsakar frameshifts, trunkeringsmutationer, stop-loss-mutationer eller mutationer i gener som är kända för att vara biologiskt relevanta för den aktuella fenotypen.

3.4. Genomsnittlig TPR och känslighet

Som framgår av tabell 5 varierar det positiva prediktiva värdet (PPV) för varje verktyg, med undantag för CUSHAW3, från 91 % till 99,9 % för SNV:er, men den genomsnittliga känsligheten är mycket låg (omkring 50 %). Denna diskrepans kan bero på ett antal orsaker, men den mest sannolika är varierande täckningsdjup mellan exoner. Vi kan se att förutom den låga SNV-känsligheten är indelkänsligheten låg (cirka 30 %), men PPV för indels är betydligt lägre (35,86 % till 52,95 %). Detta kan bero på någon av följande orsaker: mycket korta indels är svåra att upptäcka med konventionell NGS, representationen av indels av olika variant callers kan leda till att verktygen felaktigt påstår att två indels är olika, eller att anpassningsverktygen ger olika representationer av samma indel .

Den kanske mest sannolika förklaringen till båda typerna av mutationer är frågan om djup. Som i alla variantanalysstudier leder en ökning av täckningsdjupet till en ökning av känsligheten, men det är omöjligt att garantera ett bra täckningsdjup på grund av exomfångstkiteternas oförmåga att dra ner exoner på ett enhetligt sätt . Dessutom finns det inte ett enda kit för exomfångst som täcker alla exoner. Det har faktiskt visats att variantanalys av sekvensering av hela genomet med ett genomsnittligt djup som är lägre än ett exom ger bättre resultat på grund av det enhetliga djupet. Det är därför troligt att ett stort antal varianter saknas på grund av att NIST-GiaB-listan skapades från en sammanställning av exomiska och genomiska sekvenseringsdata. I slutändan måste man för att uppnå en korrekt känslighet så småningom utföra sekvensering av hela genomet, men detta är för närvarande kostnadsdrivande för de flesta laboratorier. Lyckligtvis fortsätter denna kostnad att sjunka, och vi kommer snart att se en gradvis övergång från exomanalys till den mer kompletta helgenomanalysen.

3.5. Känslighet som funktion av djup

Med tanke på att känslighet återspeglar ett av de viktigaste prestandamåtten för ett verktyg och att de flesta verktygen kämpar för att uppnå en känslighet som är högre än 50 %, skulle vi vilja utforska ytterligare hur djupet påverkade känsligheten för variantavrop. Vi tittade på ett antal olika kombinationer av verktyg för att fastställa vilka de bästa pipelines, variant callers och aligners var. I figur 2 tog vi de fem bästa kombinationerna av variantkallar och aligners enligt deras känslighet och falskt positiva frekvens (FPR). Det vill säga, vi valde de som hade det högsta antalet TP SNV:er som kallades utöver det lägsta antalet FP SNV:er. Vid en inspektion är det som genast sticker ut att känsligheten är lägre än förväntat. Alla pipelines presterar på ungefär samma nivå: de identifierar de flesta av sina varianter när ett djup på cirka 150x har uppnåtts, vilket tyder på att detta djup sannolikt är tillräckligt och att antalet saknade varianter troligen beror på att vissa exoner har lägre täckning än genomsnittet. Observera att fyra av de fem bäst presterande pipelines har GATK UnifiedGenotyper som variant caller, vilket visar på dess överlägsna prestanda oavsett vilken alignerare som används, vilket visas i figur 3(b).

Figur 2

Känslighet som en funktion av djup för de fem bästa pipelines. De fem bästa pipelines visas här med djupet för varje SNV plottat mot känsligheten.

(a)

(b)

(a)
(b)

Figur 3

Känslighet som en funktion av djupet för den bästa aligneraren och den bästa variantsökaren. (a) Resultat för djupet för varje SNV plottat mot känsligheten för den bästa aligneraren, BWA mem, parat med varje variant caller. (b) Resultat för djupet av varje SNV plottat mot känsligheten för den bästa variantutredaren, GATK UnifiedGenotyper, parat med varje alignerare.

Inom att titta på de fem bästa pipelinerna bestämde vi oss för att det skulle vara användbart att utföra samma analys på den bästa aligneraren kopplad till varje variantutredare (figur 3 a)), samt på den bästa variantutredaren kopplad till varje alignerare (figur 3 b)). Liksom för pipelines identifierades den bästa aligneraren som den som gav det högsta antalet TP SNV och det lägsta antalet FP SNV – i detta fall BWA mem. Trots att vi har den bästa anpassningen att arbeta med ser vi fortfarande en ganska stor skillnad mellan variantkallarna, vilket sannolikt beror på de olika algoritmer som de använder (figur 3(a)). När det gäller den bäst presterande variantkalleraren, GATK UnifiedGenotyper, verkar det dock finnas mindre variation mellan de fyra bästa aligners, vilket tyder på att den presterar ganska bra i de flesta situationer, med CUSAHW3 och MOSAIK som undantag.

3.6. Delade varianter bland de bästa pipelines

Slutningsvis ville vi veta hur unika uppsättningarna av variantutrop var mellan de olika pipelines. För att göra detta fokuserade vi återigen på de fem bästa pipelines för variantuppringning: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper och Novoalign plus UnifiedGenotyper. Som framgår av figur 4 finns det en stor mängd överlappning mellan de fem olika pipelines i fråga, med 15 489 SNV (70 %) som delas av totalt 22 324 distinkta varianter. Man kan dock också hävda att detta till stor del beror på att fyra av de fem pipelines använder UnifiedGenotyper som variant caller. Denna uppfattning bekräftas av det faktum att det största antalet varianter som är unika för en pipeline, 367, tillhör kombinationen BWA sampe plus HaplotypeCaller. Det är också värt att notera att det näst högsta antalet unika SNV:er också tillhör BWA sampe aligner, så det är möjligt att det höga antalet unika SNV:er bättre tillskrivs aligern än variant caller.

Figur 4

Skärningspunkten mellan SNV:erna som identifierats av de fem bästa pipelines.

4. Slutsatser

Vi fann att bland de trettio olika pipelines som testades uppvisade Novoalign plus GATK UnifiedGenotyper den högsta känsligheten samtidigt som antalet FP för SNVs var lågt. Av aligners presterade BWA mem konsekvent bäst, men resultaten varierade fortfarande mycket beroende på vilken variant caller som användes. Naturligtvis följer därav att den bästa variantkällaren, GATK UnifiedGenotyper, för det mesta gav liknande resultat oavsett vilken alignerare som användes. Det är dock uppenbart att indels fortfarande är svåra att hantera för någon pipeline, eftersom ingen av pipelinerna uppnådde en genomsnittlig känslighet på mer än 33 % eller en PPV på mer än 53 %. Förutom den låga allmänna prestandan när det gäller att upptäcka indels är känsligheten, oavsett mutationstyp, ett problem för varje pipeline som beskrivs i detta dokument. Det förväntade antalet SNV:er för NA12878:s exom är 34 886, men även när man använder föreningen av alla varianter som identifierats av de fem bästa pipelines var det största antalet identifierade mycket lågt (22 324). Det verkar som om exomanalysen, även om den fortfarande är mycket användbar, har sina begränsningar även när det gäller något så till synes enkelt som SNV-detektering.

Offentliggörande

Adam Cornish är doktorand i Chittibabu Gudas labb med utbildning inom datavetenskap och genomik. Chittibabu Guda (biträdande professor) har en tvärvetenskaplig bakgrund inom molekylär- och beräkningsbiologi. Han har publicerat ett antal beräkningsmetoder med olika tillämpningar inom biomedicinsk forskning sedan 2001.

Intressekonflikter

Författarna känner inte till några konkurrerande intressen.

Författarnas bidrag

Adam Cornish har utformat studien, utfört alla analyser, gjort figurerna och skrivit artikeln. Chittibabu Guda gav viktig återkoppling om förbättringar av dokumentet och bidrag till själva analyserna och redigerade dokumentet grundligt. Alla författare läste och godkände den slutliga uppsatsen.

Acknowledgments

Detta arbete stöddes av utvecklingsmedel till Chittibabu Guda från University of Nebraska Medical Center (UNMC). Författarna vill tacka skaparna av Novoalign för att de gjorde programvaran tillgänglig som testversion och Bioinformatics and Systems Biology Core facility vid UNMC för infrastrukturstödet som underlättade denna forskning.