Abstrakt

Vysokokapacitní sekvenování, zejména exomů, je oblíbeným diagnostickým nástrojem, ale je obtížné určit, které nástroje jsou při analýze těchto dat nejlepší. V této studii používáme výsledky projektu NIST Genome in a Bottle jako nový zdroj pro ověření naší pipeline pro analýzu exomů. Používáme šest různých zarovnávačů a pět různých voličů variant, abychom určili, která pipeline z celkového počtu 30 funguje nejlépe na lidském exomu, který byl použit k vytvoření seznamu variant detekovaných konsorciem Genome in a Bottle. Z těchto 30 pipeline jsme zjistili, že Novoalign ve spojení s GATK UnifiedGenotyper vykazuje nejvyšší citlivost při zachování nízkého počtu falešně pozitivních SNV. Je však zřejmé, že indely jsou stále obtížně zpracovatelné pro jakoukoli pipeline, přičemž žádný z nástrojů nedosáhl průměrné citlivosti vyšší než 33 % nebo pozitivní prediktivní hodnoty (PPV) vyšší než 53 %. Nakonec bylo podle očekávání zjištěno, že srovnávače mohou hrát při detekci variant stejně zásadní roli jako samotné volače variant

1. Pozadí

V posledních několika letech došlo k mnoha pokrokům v technologiích vysokokapacitního sekvenování. Díky těmto pokrokům je nyní možné detekovat velké množství potenciálních variant způsobujících onemocnění , a v několika případech byla data sekvenování nové generace (NGS) dokonce použita pro diagnostické účely. To je částečně způsobeno vývojem sekvenačních technologií v posledních několika letech, ale také řadou vylepšení různých bioinformatických nástrojů používaných k analýze hor dat produkovaných přístroji NGS .

Při hledání mutací u pacienta je typickým pracovním postupem sekvenování jeho exomu pomocí sekvenátoru Illumina, zarovnání nezpracovaných dat s lidským referenčním genomem a následná identifikace jednonukleotidových variant (SNV) nebo krátkých inzercí a delecí (indelů), které by případně mohly způsobit nebo ovlivnit fenotyp, o který se zajímáme . Zatímco toto je poměrně jednoduché, rozhodování o nejlepších nástrojích, které se mají použít v každé fázi analytického potrubí, už tak jednoduché není. Existuje velké množství nástrojů, které se používají v různých mezikrocích, ale dva nejdůležitější kroky v celém procesu jsou zarovnání surových čtení ke genomu a následné vyhledávání variant (tj. SNV a indelů) . V této studii si klademe za cíl pomoci dnešnímu bioinformatikovi objasněním správné kombinace nástroje pro zarovnání krátkých čtení a nástroje pro volání variant pro zpracování dat exomového sekvenování produkovaných přístroji NGS.

V minulosti byla provedena řada těchto studií, ale všechny měly v té či oné podobě nedostatky. V ideálním případě je třeba mít k dispozici seznam všech známých variant obsažených ve vzorku, aby bylo možné při spuštění pipeline analytických nástrojů otestovat a s jistotou vědět, že fungují správně. V minulosti však žádný takový seznam neexistoval, takže se validace musela provádět méně úplnými metodami. V některých případech se validace prováděla generováním simulovaných dat, aby se vytvořil soubor známých pravdivě pozitivních (TP) a pravdivě negativních (TN) . Tento postup sice pohodlně poskytuje seznam všech TP a TN v souboru dat, ale špatně reprezentuje biologii. Jiné metody validace pipeline pro volání variant zahrnují použití genotypovacích polí nebo Sangerova sekvenování k získání seznamu TP a falešně pozitivních (FP) . Tyto metody mají tu výhodu, že poskytují biologicky ověřené výsledky, ale mají také tu nevýhodu, že nejsou komplexní kvůli omezenému počtu míst na genotypovacích maticích a neúměrně vysokým nákladům na Sangerovu validaci, pokud se provádí tisíckrát. A konečně, žádná z těchto studií se nezaměřila na zkoumání vlivu zarovnávače krátkých čtení na volání variant. V důsledku toho nebylo možné nezávisle posoudit vliv výkonu aligátoru na vyšší úroveň.

V této studii máme výhodu seznamu variant pro anonymní ženu z Utahu (ID subjektu: NA12878, původně sekvenovaná pro projekt 1000 Genomes ), který byl experimentálně ověřen konsorciem Genome in a Bottle (GiaB) pod vedením NIST. Tento seznam variant byl vytvořen integrací 14 různých datových sad z pěti různých sekvenátorů a umožňuje nám validovat jakýkoli seznam variant vytvořený našimi pipeline pro analýzu exomu . Novinkou této práce je ověření správné kombinace srovnávačů a voličů variant na rozsáhlém a experimentálně určeném souboru dat variant: NIST-GiaB.

K provedení naší analýzy budeme používat jednu ze sad exomových dat původně použitých k vytvoření seznamu NIST-GiaB. Vybrali jsme pouze jeden z původních exomů vytvořených technologií Illumina TruSeq, protože jsme chtěli poskytnout standardní scénář použití pro někoho, kdo chce provést analýzu NGS, a zatímco cena sekvenování celého genomu stále klesá, sekvenování exomů je stále populární a životaschopnou alternativou . Je také důležité poznamenat, že podle Bamshada et al. je v současné době očekávaný počet SNV na jeden evropsko-americký exom 20 283 ± 523 . Navzdory tomu byl celkový počet SNV nalezených v seznamu NIST-GiaB s potenciálem existence v souboru dat exomu TruSeq 34 886, což je výrazně více, než se očekávalo. To je pravděpodobně způsobeno skutečností, že ačkoli byla k vytvoření dat NIST-GiaB použita sada exomů, byla také doplněna sekvenováním celého genomu.

Nakonec jsme zvažovali velké množství aligátorů a voličů variant, ale nakonec jsme vybrali 11 nástrojů na základě rozšířenosti, popularity a relevance pro náš soubor dat (např. nástroje SNVMix, VarScan2 a MuTect nebyly použity, protože jsou určeny pro použití na vzorcích získaných z nádorů). Naše vlastní analýza zahrnuje srovnání šesti srovnávačů (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK a Novoalign) a pěti voličů variant (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup a SNPSVM ). V této studii se také snažíme určit, jak velký vliv, pokud vůbec nějaký, má zarovnávač na volání variant a které zarovnávače mají nejlepší výsledky při použití normálního vzorku exomu Illumina. Pokud je nám známo, jedná se o první zprávu, která ověřuje všechny možné kombinace (celkem 30 pipelines) široké škály alignerů a variant callerů.

2. Metody

2.1. Způsoby analýzy a vyhodnocování variant

2.2. Metody

2.3. Způsoby analýzy a vyhodnocování variant Datové sady

Lidský referenční genom hg19 byl stažen z prohlížeče UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) a byl použit k provedení zarovnání. Lidský exom, SRR098401, byl stažen z Archivu čtení sekvencí (SRA) (http://www.ncbi.nlm.nih.gov/sra). Pro účely anotace a kalibrace byly použity dbSNP137 bez míst po verzi 129, HapMap 3.3, Human Omni 2.5 BeadChip a seznamy sad zlatých standardů indelů Mills a 1000 G (vše z ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2 Analýza a kalibrace. Pipeline

Obrázek 1 ukazuje pracovní postup použitý v této studii, který je podobný postupu uvedenému v příručce osvědčených postupů vypracované institutem The Broad Institute . Zahrnuje řadu kroků, které mají zajistit, aby vytvořené soubory zarovnání byly co nejkvalitnější, a také několik dalších kroků, které mají zaručit správné pojmenování variant. Nejprve byla surová čtení zarovnána do hg19 a poté byly ze zarovnání odstraněny PCR duplikáty. Dále bylo provedeno zarovnání čtení kolem indelů, aby se usnadnila pozdější identifikace indelů. Posledním krokem zpracování zarovnání bylo provedení kroku rekalibrace skóre kvality báze, který pomáhá zmírnit přirozené zkreslení a nepřesnosti skóre vydávaného sekvenátory. Bohužel i přes tyto kroky byla míra zarovnání jednotlivých zarovnávačů výrazně nižší, než se očekávalo, takže k vyrovnání této skutečnosti byla použita sada nástrojů fastx k odfiltrování čtení s nízkou kvalitou (tabulka 1). Jako čtení s nízkou kvalitou byla definována ta čtení, která měla alespoň polovinu skóre kvality nižší než 30. Po zpracování zarovnání bylo provedeno volání variant a filtrování variant.

Obrázek 1

Schéma použité pipeline analýzy dat. Aby bylo zajištěno vytvoření zarovnání nejvyšší kvality, jsou čtení nejprve filtrována a poté zarovnána k lidskému referenčnímu genomu hg19. Poté jsou odstraněny PCR duplikáty, čtení jsou zarovnána kolem domnělých indelů a skóre kvality bází je rekalibrováno. Nakonec jsou varianty pojmenovány a ověřeny podle seznamu variant NIST-GiaB.

Šest nástrojů použitých ke generování zarovnání bylo Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK a Novoalign a pět nástrojů použitých ke generování variant bylo FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup a SNPSVM, jak je vidět v tabulce 2.

2.3. Filtrování

Surová data byla získána ze SRA (SRR098401), rozdělena pomocí fastq-dump a filtrována pomocí sady nástrojů fastx. Konkrétně fastq-dump používal příznaky -split_files a -split_spot a fastq_quality_filter byl spuštěn s následujícími argumenty: -Q 33 -q 30 -p 50. Poté byly čtení správně spárovány pomocí vlastního skriptu.

2.4. Zarovnávání

Zarovnávače používaly výchozí argumenty s výjimkou případů, kdy byl použit argument vlákna, pokud byl k dispozici. Byly použity následující příkazy.

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. 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. Výpočet hloubky pokrytí zarovnání

Pro zajištění správného výpočtu hloubky pokrytí byl použit modul Picard Tools CalculateHsMetrics s následujícími argumenty:(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.bedDůležité je poznamenat, že soubor TruSeq exome bed musí mít předřazenou hlavičku ze souboru SAM alignment, aby tento modul fungoval. Dále musí být sloupec 6 přesunut do sloupce 4 a sloupec 5 musí být ze souboru TruSeq bed odstraněn.

2.6. Zpracování zarovnávacích souborů

Zpracování zarovnávacích souborů (soubory SAM/BAM) vyžadovalo u všech zarovnávačů následující kroky:(1)Převod SAM na BAM pomocí SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Třídění souborů BAM pomocí modulu Picard Tools, SortSam:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam SORT_ORDER=souřadnice(3)Odstranění duplicit PCR pomocí modulu 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)Skupina čtení přidaná do souborů zarovnání pomocí modulu 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)Zarovnání kolem indelů pomocí modulů GATK RealignerTargetCreator a 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)Rekalibrace báze pomocí modulů GATK BaseRecalibrator a 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. Volání variant

Pro každý volací program variant byly použity výchozí argumenty, pokud neobsahoval příznak „threads“ nebo „parallel“, v takovém případě byl použit i ten. Kromě toho byly indely volány odděleně od SNV, pokud to bylo možné. Konkrétně byly použity následující příkazy.

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.vcfVzhledem k neexistenci požadovaných příznaků CIGAR v souboru zarovnání se SNPSVM nepodařilo vyvolat varianty pro CUSHAW3 a SAMtools mpileup nemohl ze stejného důvodu vyvolat varianty na zarovnáních MOSAIK. Také vzhledem k tomu, že SNPSVM detekuje pouze SNV, nebyly pro tento program hlášeny žádné indely.

2.8. Filtrace variant

Filtrace se lišila v závislosti na použitém programu pro volání variant. V případě programů GATK HaplotypeCaller a GATK UnifiedGenotyper byly moduly GATK, VariantRecalibrator a ApplyRecalibration, použity k filtrování SNV s použitím HapMap 3.3, Omni 2.5 SNP BeadChip a dbSNP 137 bez dat 1000 Genome jako trénovacích sad. Pro SNPSVM bylo použito skóre QUAL ≥ 4 a hodnoty DP ≥ 6. Pro FreeBayes a SAMtools bylo použito skóre QUAL ≥ 20 a hodnoty DP ≥ 6.

2.9. Porovnání variant

Pro porovnání variant byl použit USeq 8.8.1 pro porovnání SNV sdílených mezi všemi soubory dat. Pro porovnání indelů byl použit nástroj vcflib vcfintersect. Soubor TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, získaný 11. prosince 2013, byl použit k omezení porovnání na místa, která mohla být zachycena soupravou exome pull down použitou při sekvenování SRR098401. Tento soubor lze získat od společnosti Illumina zde: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Aby bylo zajištěno identické zastoupení variant mezi různými soubory volání, byl k předzpracování souborů vcf použit nástroj vcflib vcfallelicprimitives

2.10. Statistické výpočty

Pravdivý výsledek (TP). Jedná se o mutaci, která byla detekována testovanou pipeline a která existuje v seznamu NIST-GiaB.

False Positive (FP). Jedná se o mutaci, která byla zjištěna testovaným postupem, ale která neexistuje v seznamu NIST-GiaB.

Pravdivě negativní (TN). Jedná se o mutaci, která nebyla zjištěna testovaným postupem a která neexistuje v seznamu NIST-GiaB.

Falešně negativní (FN). Jedná se o mutaci, která nebyla testovanou soupravou detekována, ale která existuje v seznamu NIST-GiaB:

3. Výsledky a diskuse

3.1. Výsledky a diskuse Předfiltrování variant

Při provádění analýzy variant je jedním z mnoha úskalí, která je třeba vzít v úvahu, prostor sekvence exomu (definovaný soupravou pro zachycení exomu) a to, jak může ovlivnit výsledky analýzy. V tomto případě jsme měli k dispozici jeden exom (SRR098401), který byl získán pomocí sady Illumina TruSeq exome kit a sekvenován na přístroji HiSeq 2000. S ohledem na tuto skutečnost jsme se chtěli ujistit, že měříme schopnost bioinformatických nástrojů vykonávat svou práci, a ne to, jak dobře pracuje souprava pro zachycení exomu Illumina TruSeq. To znamená, že jsme se chtěli pokusit vyvolat pouze varianty, které by se měly vyskytovat v exonech definovaných soupravou pro stahování. Z tohoto důvodu používáme lůžkový soubor poskytnutý společností Illumina, nikoli obecný anotační lůžkový soubor, například RefSeq pro hg19. Zjistili jsme, že pro tohoto konkrétního jedince by se podle seznamu NIST-GiaB mělo v oblastech definovaných ložem souboru TruSeq nacházet celkem 34 886 SNV a 1 473 indelů.

Po odfiltrování variant, které se nenacházely v oblastech definovaných ložem souboru Illumina TruSeq exome, jsme se dostali ze stovek tisíc předpokládaných variant (údaje nejsou uvedeny) na průměrně asi 23 000 variant (SNV a indelů) na pipeline (tabulka 3). To je pro výzkumné pracovníky pro začátek důležitý krok, protože se tím výrazně zmenšuje prostor pro vyhledávání potenciálně zajímavých variant.

3.2. Surové výsledky variant ve srovnání s GiaB

Jedním z aspektů, které jsme chtěli při tomto srovnání pochopit, byl význam filtrování variant zjištěných těmito nástroji. Důvodem je to, že v ideálním případě by člověk chtěl mít co nejvyšší úroveň citlivosti, aby se mutace, které ho zajímají, neztratily v procesu filtrování. Proto je vhodné zjistit, zda je tento krok nutný a do jaké míry je nutný, protože z výsledků NIST-GiaB a přehledu Bamshad et al. je zřejmé, že citlivost může být problémem.

Jak vidíme v tabulce 3, filtrování je u některých variant callerů nutné více než u jiných, pokud jde o SNV, a u indelů je naprosto nezbytné. Ve většině případů je počet TP variant blízký nebo vyšší než námi očekávaný počet přibližně 20 000 , ale na druhou stranu je v některých případech počet FP velmi vysoký.

Je zřejmé, že v počtech generovaných jednotlivými pipeline jsou velké rozdíly. V číslech však lze najít některé společné rysy, které pravděpodobně vyplývají z algoritmického původu jednotlivých nástrojů. FreeBayes produkuje jak největší počet nefiltrovaných variant, tak nejvyšší počet FP. Je pravděpodobné, že takový výkon vidíme pouze u tohoto nástroje díky tomu, že ačkoli není jediným nástrojem pro volání variant založeným na bayesovské inferenci, je jedinečný ve své interpretaci zarovnání. To znamená, že se jedná o volající nástroj založený na haplotypech, který identifikuje varianty na základě sekvence samotných čtení namísto zarovnání, což je druhý způsob, jakým pracuje UnifiedGenotyper společnosti GATK.

Dále vidíme, že zarovnávače založené na Burrows-Wheelerovi si vedou velmi podobně: Bowtie2, BWA mem a BWA sampe dosahují podobných výsledků ve všech oblastech. Lze se domnívat, že je to pravděpodobně způsobeno tím, že všechny tyto nástroje využívají při plnění určeného úkolu podobné algoritmy. Toto pozorování je podpořeno skutečností, že MOSAIK (zarovnání s mezerami pomocí Smith-Watermanova algoritmu) a CUSHAW3 (hybridní přístup seeding) mají oba velmi odlišné základní algoritmy a následně poskytují velmi odlišné výsledky.

Tento rozdíl ve výsledcích korelujících s různými algoritmy je nejlépe vidět na výsledcích SNPSVM. Z voličů variant je to jediný, který ke generování volání SNV využívá stroje podpůrných vektorů a sestavování modelů. Zdá se, že ačkoli má nevýhodu v tom, že není tak citlivý jako ostatní metody, těží z toho, že je mimořádně přesný bez ohledu na použitý aligátor. To naznačuje, že při použití tohoto volače variant je možné zcela vynechat krok filtrování.

Co se týče indelů, zdá se, že žádný aligátor nevyniká mezi ostatními jako ten, který dobře zvládá tento typ mutací. Při pohledu na počet FP je totiž zřejmé, že právě variant caller hraje největší roli v přesnosti identifikace indelů. Navíc nejsou k dispozici údaje pro pipeline CUSHAW3 plus SNPSVM ani MOSAIK plus SAMtools mpileup, a to z důvodu, že zarovnávací soubory neobsahují potřebné řetězce CIGAR, aby mohly variant callery fungovat následně. A konečně, důvodem, proč nejsou k dispozici údaje o indelech pro SNPSVM, je skutečnost, že tento nástroj se používá výhradně pro identifikaci SNV.

3.3. Výsledky filtrování ve srovnání s GiaB

Jak je uvedeno v tabulce 4, standardními postupy filtrování se podařilo odstranit velké množství FP SNV pro každou pipeline; zdá se však, že tyto filtry jsou ve většině případů příliš agresivní pro detekci SNV, ale nejsou dostatečně přísné pro indely. To je zřejmé při pohledu na rozdíly v počtu FP uváděných v jednotlivých souborech dat. Například u Bowtie2 s Freebayesem bylo odstraněno 72 570 FP SNV (snížení o 98 %), ale pouze 1 736 FP indelů (snížení o 70 %). Je třeba také poznamenat, že použité filtry byly závislé na pipeline a většinou v rámci každé pipeline vedly k podobnému snížení SNV a indelů FP. Jedinou výjimkou je zde počet variant identifikovaných ze zarovnání CUSHAW3 ve srovnání s ostatními zarovnáními: celkově je počet TP SNV nižší, počet FP SNV je vyšší a je to jediný alignner, který po filtrování produkuje více než 1 000 FP indelů.

Vzhledem k tomu, že filtrování výrazně snižuje počet variant TP, mohlo by být rozumné, s výjimkou pipeline používajících CUSHAW3 a FreeBayes, tento krok při vyhledávání vzácných variant s vysokým dopadem vynechat. Místo toho je možné věnovat více času procesu filtrování, který je založen spíše na biologii než na statistice. Například může mít větší smysl investovat čas do identifikace malého seznamu variant, které budou mít pravděpodobně vysoký dopad: mutace v místě sestřihu, indely, které způsobují posun rámce, zkrácené mutace, stop-loss mutace nebo mutace v genech, o nichž je známo, že jsou biologicky relevantní pro fenotyp, který nás zajímá.

3.4. Varianty, o nichž je známo, že jsou biologicky relevantní pro fenotyp, který nás zajímá. Průměrná TPR a citlivost

Jak je vidět v tabulce 5, pozitivní prediktivní hodnota (PPV) se u každého nástroje, s výjimkou nástroje CUSHAW3, pohybuje od 91 % do 99,9 % pro SNV, ale průměrná citlivost je velmi nízká (kolem 50 %). Tento rozdíl může být způsoben řadou důvodů, ale nejpravděpodobnějším z nich je různá hloubka pokrytí napříč exony. Vidíme, že kromě nízké citlivosti SNV je nízká i citlivost indelů (kolem 30 %); PPV pro indely je však podstatně nižší (35,86 % až 52,95 %). To může být způsobeno některým z následujících důvodů: velmi krátké indely jsou běžným NGS těžko detekovatelné , reprezentace indelů různými variant callery může způsobit, že nástroje nesprávně tvrdí, že dva indely jsou různé, nebo nástroje pro zarovnávání vytvářejí různé reprezentace stejného indelu.

Možná nejpravděpodobnějším vysvětlením pro oba typy mutací je otázka hloubky. Stejně jako v případě každé studie analýzy variant vede zvýšení hloubky pokrytí ke zvýšení citlivosti, ale není možné zaručit dobrou hloubku pokrytí kvůli neschopnosti souprav pro zachycení exomu rovnoměrně stáhnout exony . Kromě toho žádná souprava pro zachycení exomu nepokrývá všechny exony. Bylo totiž prokázáno, že analýza variant sekvenování celého genomu s průměrnou hloubkou nižší než exom dosahuje lepších výsledků díky rovnoměrnosti uvedené hloubky. Je tedy pravděpodobné, že velký počet variant chybí vzhledem k tomu, že seznam NIST-GiaB byl vytvořen na základě kompilace dat z exomového a genomového sekvenování. Pro dosažení správné citlivosti bude nakonec nutné provést sekvenování celého genomu, což je však v současné době pro většinu laboratoří finančně nedostupné. Naštěstí tyto náklady stále klesají a brzy budeme svědky postupného přechodu od exomové analýzy k úplnější analýze celého genomu.

3.5. Analýza celého genomu. Citlivost jako funkce hloubky

Protože citlivost odráží jeden z nejdůležitějších ukazatelů výkonnosti nástroje a většina nástrojů se snaží dosáhnout citlivosti vyšší než 50 %, chtěli bychom dále prozkoumat, jak hloubka ovlivnila citlivost volání variant. Zkoumali jsme řadu různých kombinací nástrojů, abychom zjistili, jaké jsou nejlepší pipeline, variant callery a alignery. Pro obrázek 2 jsme vzali pět nejlepších kombinací variant callerů a alignerů podle jejich citlivosti a míry falešně pozitivních výsledků (FPR). To znamená, že jsme vybrali ty, které měly kromě nejvyššího počtu vyvolaných TP SNV také nejnižší počet FP SNV. Při prohlídce okamžitě vynikne, že citlivost je nižší, než se očekávalo. Všechny pipeline dosahují zhruba stejné úrovně: většinu variant identifikují v době, kdy je dosaženo hloubky přibližně 150x, což naznačuje, že tato hloubka je pravděpodobně dostatečná a že počet chybějících variant je pravděpodobně způsoben tím, že některé exony mají nižší než průměrné pokrytí. Všimněte si, že čtyři z pěti nejvýkonnějších pipeline mají jako volič variant nástroj GATK UnifiedGenotyper, což ukazuje jeho lepší výkonnost bez ohledu na použitý aligátor, jak ukazuje obrázek 3b).

Obrázek 2

Citlivost v závislosti na hloubce pro pět nejvýkonnějších pipeline. Je zde zobrazeno pět nejlepších pipeline s hloubkou každého SNV vynesenou v závislosti na citlivosti.

(a)

(b)

(a)
(b)

Obrázek 3

Citlivost v závislosti na hloubce pro špičkový srovnávač a špičkový volič variant. (a) Výsledky pro hloubku každého SNV vynesené v závislosti na citlivosti pro top aligner, BWA mem, spárovaný s každým variant callerem. (b) Výsledky pro hloubku každého SNV vynesené v závislosti na citlivosti pro nejlepší variant caller, GATK UnifiedGenotyper, spárovaný s každým aligátorem.

Kromě zkoumání pěti nejlepších pipeline jsme zjistili, že by bylo užitečné provést stejnou analýzu pro nejlepší aligátor spojený s každým variant callerem (obrázek 3(a)) a také pro nejlepší variant caller spojený s každým aligátorem (obrázek 3(b)). Stejně jako u pipeline byl nejlepší aligner identifikován jako ten, který produkoval nejvyšší počet TP SNV a nejnižší počet FP SNV – v tomto případě BWA mem. Přestože máme k dispozici nejlepší zarovnání, stále vidíme poměrně velké rozdíly mezi volači variant, což lze pravděpodobně přičíst rozdílným algoritmům, které používají (obrázek 3(a)). V případě nejvýkonnějšího volače variant, GATK UnifiedGenotyper, se však zdá, že mezi čtyřmi nejlepšími zarovnávači jsou menší rozdíly, což naznačuje, že si ve většině situací vede poměrně dobře, přičemž výjimku tvoří CUSAHW3 a MOSAIK.

3.6. Sdílené varianty mezi nejlepšími pipeline

Nakonec jsme chtěli zjistit, jak jedinečné jsou sady volání variant mezi jednotlivými pipeline. Za tímto účelem jsme se opět zaměřili na pět nejlepších pipeline pro volání variant: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper a Novoalign plus UnifiedGenotyper. Jak je vidět na obrázku 4, existuje velké množství překryvů mezi pěti různými dotyčnými pipeline, přičemž 15 489 SNV (70 %) je společných z celkového počtu 22 324 odlišných variant. Lze však také tvrdit, že je to z velké části způsobeno tím, že čtyři z pěti pipeline používají jako volací program variant UnifiedGenotyper. Tuto domněnku potvrzuje skutečnost, že největší počet variant jedinečných pro jednu pipeline, 367, patří kombinaci BWA sampe plus HaplotypeCaller. Stojí také za povšimnutí, že druhý nejvyšší počet unikátních SNV patří rovněž aligátoru BWA sampe, takže je možné, že vysoký počet unikátních SNV lze lépe připsat aligátoru než variant calleru.

Obrázek 4

Průsečík SNV identifikovaných pěti nejlepšími pipeline.

4. Závěry

Zjistili jsme, že z třiceti různých testovaných pipeline Novoalign plus GATK UnifiedGenotyper vykazuje nejvyšší citlivost při zachování nízkého počtu FP pro SNV. Ze srovnávačů si trvale nejlépe vedl BWA mem, ale výsledky se přesto značně lišily v závislosti na použitém variant calleru. Z toho samozřejmě vyplývá, že nejlepší variant caller, GATK UnifiedGenotyper, většinou poskytoval podobné výsledky bez ohledu na použitý aligner. Je však snadno patrné, že indely jsou pro jakoukoli pipeline stále obtížně zpracovatelné, přičemž žádná z pipeline nedosáhla průměrné citlivosti vyšší než 33 % nebo PPV vyšší než 53 %. Kromě nízké celkové výkonnosti, kterou pozorujeme při detekci indelů, je citlivost, bez ohledu na typ mutace, problémem pro každou pipeline popsanou v tomto článku. Očekávaný počet SNV pro exom NA12878 je 34 886, ale i při použití spojení všech variant identifikovaných pěti nejlepšími pipeline byl největší počet identifikovaných variant velmi nízký (22 324). Zdá se, že i když je analýza exomu stále velmi užitečná, má svá omezení, a to i pokud jde o něco tak zdánlivě jednoduchého, jako je detekce SNV.

Zveřejnění

Adam Cornish je postgraduální student v laboratoři Chittibabu Guda se vzděláním v oblasti informatiky a genomiky. Chittibabu Guda (docent) má interdisciplinární vzdělání v oblasti molekulární a výpočetní biologie. Od roku 2001 publikoval řadu výpočetních metod s různými aplikacemi v biomedicínském výzkumu.

Konflikt zájmů

Autoři si nejsou vědomi žádných konkurenčních zájmů.

Přínos autorů

Adam Cornish navrhl studii, provedl všechny analýzy, vytvořil obrázky a napsal článek. Chittibabu Guda poskytl zásadní zpětnou vazbu ohledně vylepšení článku a vstupy k samotným analýzám a článek důkladně upravil. Všichni autoři si přečetli a schválili konečnou verzi článku.

Poděkování

Tato práce byla podpořena rozvojovými fondy pro Chittibabu Gudu z University of Nebraska Medical Center (UNMC). Autoři by rádi poděkovali tvůrcům programu Novoalign za poskytnutí zkušební verze softwaru a Bioinformatics and Systems Biology Core facility při UNMC za infrastrukturní podporu umožňující tento výzkum.