Abstract

A nagy áteresztőképességű szekvenálás, különösen az exom szekvenálás népszerű diagnosztikai eszköz, de nehéz meghatározni, hogy mely eszközök a legjobbak ezen adatok elemzésében. Ebben a tanulmányban a NIST Genome in a Bottle eredményeit használjuk új forrásként az exom-elemző csővezetékünk validálásához. Hat különböző igazítót és öt különböző variánshívót használunk annak meghatározására, hogy az összesen 30 csővezeték közül melyik teljesít a legjobban egy olyan emberi exomon, amelyet a Genome in a Bottle konzorcium által kimutatott variánsok listájának létrehozásához használtak. E 30 pipeline közül azt találtuk, hogy a Novoalign a GATK UnifiedGenotyperrel együtt mutatta a legnagyobb érzékenységet, miközben alacsony volt a hamis pozitív SNV-k száma. Nyilvánvaló azonban, hogy az indeleket még mindig nehéz kezelni, mivel egyik eszköz sem ért el 33%-nál magasabb átlagos érzékenységet vagy 53%-nál magasabb pozitív prediktív értéket (PPV). Végül, ahogyan az várható volt, kiderült, hogy az illesztőprogramok ugyanolyan fontos szerepet játszhatnak a variánsok felderítésében, mint maguk a variánshívók.

1. Háttér

Az elmúlt néhány évben számos előrelépés történt a nagy áteresztőképességű szekvenálási technológiák terén. Ezeknek az előrelépéseknek köszönhetően ma már nagyszámú potenciális betegséget okozó variáns detektálására van lehetőség , és néhány esetben az újgenerációs szekvenálás (NGS) adatait még diagnosztikai célokra is felhasználták . Ez részben a szekvenálási technológiák elmúlt évekbeli fejlődésének köszönhető, de annak is, hogy számos fejlesztés történt a különböző bioinformatikai eszközökben, amelyeket az NGS-eszközök által termelt adathegyek elemzésére használnak.

Mutációk keresésekor egy betegnél a tipikus munkafolyamat az exom szekvenálása Illumina szekvenálóval, a nyers adatok összehangolása a humán referencia genommal, majd az egynukleotid variánsok (SNV-k) vagy rövid beillesztések és deléciók (indelek) azonosítása, amelyek esetleg a kívánt fenotípust okozhatják vagy befolyásolhatják . Míg ez meglehetősen egyszerű, az elemzési folyamat egyes szakaszaiban használandó legjobb eszközök kiválasztása már nem az. Számos eszközt használnak a különböző köztes lépésekben, de az egész folyamat két legfontosabb lépése a nyers leolvasások összehangolása a genommal, majd a variánsok (azaz SNV-k és indelek) keresése. Ebben a tanulmányban az a célunk, hogy segítsük a mai bioinformatikusokat azáltal, hogy megvilágítjuk a rövid olvasatok összehangolására szolgáló eszköz és a variánshívó eszköz helyes kombinációját az NGS eszközök által előállított exom szekvenálási adatok feldolgozásához.

A múltban számos ilyen vizsgálatot végeztek, de mindegyiknek valamilyen formában hátrányai voltak. Ideális esetben rendelkezni kellene egy listával a mintában található minden ismert variánsról, hogy az elemzőeszközök csővezetékének futtatásakor tesztelni lehessen azt, hogy biztosan lehessen tudni, hogy az megfelelően működik. A múltban azonban nem létezett ilyen lista, így a validálást kevésbé teljes körű módszerekkel kellett elvégezni. Egyes esetekben a validálást szimulált adatok generálásával végezték, hogy létrehozzák az ismert igaz pozitív (TP) és igaz negatív (TN) eredmények készletét. Ez ugyan kényelmesen biztosítja az adathalmazban lévő összes TP és TN listáját, de a biológia pontos reprezentálására nem alkalmas. A variánshívó pipelinek validálásának egyéb módszerei közé tartozik a genotipizáló tömbök vagy a Sanger-szekvenálás használata a TP-k és a hamis pozitívumok (FP) listájának összeállításához. Ezek előnye, hogy biológiailag validált eredményeket adnak, de hátrányuk, hogy a genotipizáló tömbökön található foltok korlátozott száma és a Sanger-vizsgálat több ezer alkalommal történő elvégzésének megfizethetetlen költségei miatt nem átfogóak. Végül pedig e tanulmányok egyikének sem volt célja, hogy megvizsgálja, milyen hatással van a rövid leolvasások összehangolója a variánsok megnevezésére. Következésképpen az illesztő teljesítményének upstream hatását nem lehetett függetlenül értékelni.

Ebben a tanulmányban az az előnyünk, hogy egy névtelen utahi nő esetében (alany azonosítója: NA12878, eredetileg az 1000 Genomes projekthez szekvenálták ) egy olyan variánslistára támaszkodhatunk, amelyet a NIST által vezetett Genome in a Bottle (GiaB) konzorcium kísérletileg validált. Ez a variánslista öt különböző szekvenálótól származó 14 különböző adatkészlet integrálásával jött létre, és lehetővé teszi számunkra, hogy az exómaelemző pipeline-aink által generált bármely variánslistát validáljuk . E munka újdonsága az, hogy az igazítók és a variánshívók megfelelő kombinációját egy átfogó és kísérletileg meghatározott variánsadatkészleten validáljuk: NIST-GiaB.

Az elemzés elvégzéséhez a NIST-GiaB-lista létrehozásához eredetileg használt exom-adatkészletek egyikét fogjuk használni. Azért választottuk csak az egyik eredeti Illumina TruSeq által generált exomot, mert egy standard felhasználási forgatókönyvet akartunk adni annak, aki NGS-elemzést szeretne végezni, és bár a teljes genom szekvenálás ára folyamatosan csökken, az exom szekvenálás még mindig népszerű és életképes alternatíva . Azt is fontos megjegyezni, hogy Bamshad et al. szerint jelenleg az SNV-k várható száma európai-amerikai exomonként 20 283 ± 523 . Ennek ellenére a NIST-GiaB listán talált, a TruSeq exom-adatkészletben potenciálisan létező SNV-k száma összesen 34 886 volt, ami jelentősen magasabb a vártnál. Ez valószínűleg annak köszönhető, hogy bár a NIST-GiaB adatok előállításához az exome kitet használtuk, azt egész genom szekvenálással is kiegészítettük.

Végül számos alignert és variánshívót vettünk figyelembe, de végül a 11 eszközt az elterjedtség, a népszerűség és az adatállományunk szempontjából való relevancia alapján választottuk ki (pl. az SNVMix, VarScan2 és MuTect nem került felhasználásra, mivel azokat tumorból származó mintákra szántuk). Maga az elemzésünk hat alignert (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK és Novoalign) és öt variánshívót (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup és SNPSVM ) hasonlít össze. Ebben a tanulmányban azt is megpróbáljuk meghatározni, hogy az igazítónak mekkora hatása van – ha van egyáltalán – a variánshívásra, és hogy mely igazítók teljesítenek a legjobban egy normál Illumina exome minta használata esetén. Tudomásunk szerint ez az első olyan jelentés, amely az alignerek és variánshívók széles skálájának minden lehetséges kombinációját (összesen 30 pipeline) validálja.

2. Módszerek

2.1. Vizsgálati módszerek

2.1. Vizsgálati módszerek

2. Módszerek

Adathalmazok

A hg19 emberi referencia genomot a UCSC böngészőjéből töltöttük le (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/), és azt használtuk az illesztések elvégzéséhez. A humán exomot, az SRR098401-et a Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) adatbázisból töltöttük le. Annotációs és kalibrációs célokra a 129-es verzió utáni helyek nélküli dbSNP137, a HapMap 3.3, a Human Omni 2.5 BeadChip, valamint a Mills és az 1000 G arany standard indel készletlistákat használtuk (mind a ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz oldalról).

2.2. Annotáció és kalibráció. A csővezeték

Az 1. ábra mutatja az ebben a vizsgálatban használt munkafolyamatot, amely hasonló a Broad Institute által készített Best Practices útmutatóban leírtakhoz. Ez számos lépést tartalmaz annak biztosítására, hogy az előállított illesztési fájlok a lehető legjobb minőségűek legyenek, valamint számos további lépést a változatok helyes megnevezésének garantálására. Először a nyers leolvasásokat igazítottuk a hg19-hez, majd a PCR-duplikátumokat eltávolítottuk az igazításból. Ezután az indelek azonosításának megkönnyítése érdekében a csővezeték későbbi szakaszában az indelek körül újrarendezést végeztünk. Az összehangolás feldolgozásának utolsó lépése a bázisminőségi pontszámok újrakalibrálása volt, amely segít a szekvencerek által kiadott pontszámok eredendő torzításainak és pontatlanságainak enyhítésében. Sajnos e lépések ellenére az egyes igazítók igazítási aránya jelentősen alacsonyabb volt a vártnál, ezért ennek ellensúlyozására a fastx eszközkészletet használtuk a rossz minőségű olvasatok kiszűrésére (1. táblázat). Az alacsony minőségű olvasatokat úgy definiáltuk, hogy azok az olvasatok, amelyek minőségi pontszámainak legalább a fele 30 alatt volt. Az illesztés feldolgozása után a variánshívást és a variánsszűrést végeztük el.

Az igazítások létrehozásához használt hat eszköz a Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK és Novoalign, a variánsok létrehozásához használt öt eszköz pedig a FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup és SNPSVM volt, amint az a 2. táblázatban látható.

2.3. A szoftververziók. Szűrés

A nyers adatokat az SRA-ból (SRR098401) szereztük be, fastq-dump segítségével osztottuk fel, és a fastx eszközkészlet segítségével szűrtük. Konkrétan a fastq-dump a -split_files és -split_spot jelzőket használta, a fastq_quality_filter pedig a következő argumentumokkal futott: -Q 33 -q 30 -p 50. Ezután a leolvasásokat egy egyéni szkript segítségével megfelelően párosítottuk.

2.4. Aligning

Alignerek alapértelmezett argumentumokat használtak, kivéve, ha egy szál argumentumot használtak, ahol ez rendelkezésre állt. A használt parancsok a következők.

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 -o SAM -c 10 > alignments/NA12878.novoalign.sam

2.5. Alignment Depth of Coverage Calculation

A megfelelő lefedettségi mélység kiszámításához a Picard Tools CalculateHsMetrics modulját használtuk a következő argumentumokkal:(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.bed Fontos megjegyezni, hogy a modul működéséhez a TruSeq exome bed fájlnak rendelkeznie kell a SAM igazítási fájl fejlécével. Továbbá a 6. oszlopot át kell helyezni a 4. oszlopba, az 5. oszlopot pedig el kell távolítani a TruSeq ágyfájlból.

2.6. Igazítási fájlok feldolgozása

Az igazítási fájlok (SAM/BAM fájlok) feldolgozása a következő lépéseket igényelte minden igazító esetében:(1)SAM BAM-ba konvertálás SAMtools view segítségével:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)BAM fájlok rendezése a Picard Tools moduljával, SortSam:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam SORT_ORDER=koordinate(3)PCR duplikátumok eltávolítása a Picard Tools modul, MarkDuplicates használatával:(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)Read Group hozzáadása az igazítási fájlokhoz a Picard Tools modul, AddOrReplaceReadGroups segítségével:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=koordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Realignment az indelek körül a GATK modulok RealignerTargetCreator és IndelRealigner használatával:(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 rekalibrálás a GATK modulok BaseRecalibrator és PrintReads használatával:(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. Variánshívás

Az egyes variánshívókhoz az alapértelmezett argumentumokat használtuk, kivéve, ha tartalmazott “threads” vagy “parallel” jelzőt, amely esetben azt is használtuk. Ezenkívül az indeleket lehetőség szerint az SNV-ktől elkülönítve hívtuk meg. Konkrétan a használt parancsok a következők.

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.vcfAz igazítási fájlban a szükséges CIGAR-jelzők hiánya miatt az SNPSVM nem tudott változatokat hívni a CUSHAW3 esetében, és a SAMtools mpileup ugyanezen okból nem tudott változatokat hívni a MOSAIK igazításokon. Továbbá, mivel az SNPSVM csak az SNV-ket észleli, a program nem jelentett indeleket.

2.8. Variánsszűrés

A szűrés a használt variánshívótól függően változott. A GATK HaplotypeCaller és a GATK UnifiedGenotyper esetében a GATK VariantRecalibrator és ApplyRecalibration modulokat használtuk az SNV-k szűrésére a HapMap 3.3, az Omni 2.5 SNP BeadChip és a dbSNP 137 segítségével, 1000 Genome adatok nélkül, mint képzési készletek. Az SNPSVM esetében ≥ 4 QUAL pontszámokat és ≥ 6 DP értékeket használtunk. A FreeBayes és a SAMtools esetében QUAL pontszámokat ≥ 20 és DP értékeket ≥ 6 használtunk.

2.9. Variáns-összehasonlítás

A variánsok összehasonlításához az USeq 8.8.1 programot használtuk az összes adatkészlet között közös SNV-k összehasonlítására. Az indelek összehasonlításához a vcflib vcfintersect eszközt használtuk. A 2013. december 11-én kapott TruSeq hg19 exome bed file truseq_exome_targeted_regions.hg19.bed.chr fájlt használtuk az SRR098401 szekvenálásához használt exome pull down kit által megragadható helyekre való korlátozáshoz. Ez a fájl az Illuminától itt szerezhető be: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Annak biztosítása érdekében, hogy a különböző híváskészletek között a variánsok azonos módon legyenek reprezentálva, a vcflib vcfallelicprimitives eszközt használtuk a vcf-fájlok előfeldolgozására.

2.10. Statisztikai számítások

Igaz pozitív (TP). Olyan mutáció, amelyet a vizsgált csővezeték észlelt, és amely szerepel a NIST-GiaB listán.

False Positive (FP). Olyan mutáció, amelyet a tesztelt csővezeték észlelt, de amely nem létezik a NIST-GiaB listán.

True Negative (TN). Olyan mutáció, amelyet a vizsgált csővezeték nem észlelt, és amely nem szerepel a NIST-GiaB listán.

Hamis negatív (FN). Olyan mutáció, amelyet a vizsgált csővezeték nem észlelt, de amely létezik a NIST-GiaB listán:

3. Eredmények és megbeszélés

3.1. Eredmények és megbeszélés

3.1. Eredmények és megbeszélés

. A variánsok előszűrése

A variánselemzés elvégzésekor a sok buktató közül az egyik, amit figyelembe kell venni, az exom szekvencia tér (ahogyan azt az exom capture kit meghatározza) és az, hogy ez hogyan befolyásolhatja az elemzés eredményeit. Ebben az esetben egyetlen exom (SRR098401) állt rendelkezésünkre, amelyet az Illumina TruSeq exome kit segítségével nyertünk ki és szekvenáltunk egy HiSeq 2000-en. Ezt szem előtt tartva meg akartunk győződni arról, hogy a bioinformatikai eszközök munkaképességét mérjük, és nem azt, hogy mennyire jól működik az Illumina TruSeq exome capture kit. Vagyis csak olyan variánsokat akarunk megpróbálni hívni, amelyeknek a pull down kit által meghatározott exonokban kellene jelen lenniük. Ezért az Illumina által biztosított ágyfájlt használjuk, nem pedig egy általános annotációs ágyfájlt, például a RefSeq-et a hg19 esetében. Azt találtuk, hogy az adott egyén esetében a NIST-GiaB lista szerint összesen 34 886 SNV-nek és 1473 indelnek kellene lennie a TruSeq ágyfájl által meghatározott régiókban.

Amikor kiszűrtük azokat a variánsokat, amelyek nem az Illumina TruSeq exome ágyfájl által meghatározott régiókban helyezkedtek el, a több százezer feltételezett variánsról (adatok nem láthatóak) átlagosan körülbelül 23 000 variánsra (SNV-k és indelek) jutottunk csővezetékenként (3. táblázat). Ez fontos lépés a kutatók számára, mivel jelentősen csökkenti a potenciálisan érdekes variánsok keresési terét.

3.2. Nyers variánseredmények a GiaB-vel összehasonlítva

Az egyik szempont, amelyet az összehasonlítás során meg akartunk érteni, az ezen eszközök által észlelt variánsok szűrésének fontossága volt. Ennek oka, hogy ideális esetben a lehető legnagyobb érzékenységet szeretnénk elérni, hogy az érdeklődésre számot tartó mutációk ne vesszenek el a szűrési folyamat során. Ezért kötelességünk meghatározni, hogy szükséges-e ez a lépés, és milyen mértékben szükséges, mivel a NIST-GiaB eredményei és a Bamshad et al. áttekintése alapján egyértelmű, hogy az érzékenység kérdéses lehet.

Amint a 3. táblázatban láthatjuk, egyes variánshívók esetében a szűrésre nagyobb szükség van, mint másoknál, ha SNV-kről van szó, az indelek esetében pedig feltétlenül szükséges. A legtöbb esetben a TP variánsok száma megközelíti vagy meghaladja az általunk várt, körülbelül 20 000-es számot , másrészt viszont néhány esetben az FP-k száma nagyon magas.

Az egyes pipeline-ok által generált számok között nyilvánvalóan nagy a szórás. A számokban azonban találhatunk néhány közös vonást, amelyek valószínűleg az egyes eszközök algoritmikus eredetéből erednek. A FreeBayes produkálja a legtöbb szűretlen variánst és a legtöbb FP-t is. Valószínűleg csak azért látunk ilyen teljesítményt ettől az eszköztől, mert bár nem ez az egyetlen Bayes-féle következtetésen alapuló variánshívó, az igazítások értelmezésében egyedülálló. Vagyis ez egy haplotípus-alapú hívó, amely a variánsokat maguknak a leolvasásoknak a szekvenciája alapján azonosítja, nem pedig az igazítás alapján, ez utóbbi a GATK UnifiedGenotyper működésének módja.

Egyébként azt látjuk, hogy a Burrows-Wheeler-alapú igazítók nagyon hasonlóan teljesítenek egymáshoz: A Bowtie2, a BWA mem és a BWA sampe összességében hasonló eredményeket érnek el. Feltételezhetjük, hogy ez valószínűleg annak köszönhető, hogy ezek az eszközök mind hasonló algoritmusokat használnak a kijelölt feladatuk elvégzése során. Ezt a megfigyelést támasztja alá az a tény, hogy a MOSAIK (Smith-Waterman algoritmust használó gapped alignment) és a CUSHAW3 (hibrid vetítési megközelítés) mögöttes algoritmusai nagyon különbözőek, és ezt követően nagyon különböző eredményeket produkálnak.

A különböző algoritmusokkal korreláló eredményeknek ez a különbsége a SNPSVM eredményeiben látszik a legjobban. A variánshívók közül ez az egyetlen, amely támogató vektorgépeket és modellépítést használ az SNV-hívások generálásához. Úgy tűnik, hogy bár hátránya, hogy nem olyan érzékeny, mint más módszerek, előnye, hogy a használt alignertől függetlenül rendkívül pontos. Ez azt sugallja, hogy a szűrési lépést teljesen ki lehet hagyni, ha ezt a variánshívót használjuk.

Az indeleket tekintve úgy tűnik, hogy egyetlen illesztő sem emelkedik ki a többi közül, amelyik jól kezeli ezt a fajta mutációt. Valójában, ha az FP-k számát nézzük, egyértelmű, hogy a variánshívó az, amely a legnagyobb szerepet játszik az indelek azonosításának pontosságában. Ezenkívül sem a CUSHAW3 plus SNPSVM, sem a MOSAIK plus SAMtools mpileup pipelines esetében nincsenek adatok, mivel az illesztési fájlok nem tartalmazzák a variánshívók downstream működéséhez szükséges CIGAR karakterláncokat. Végül, az SNPSVM esetében azért nincsenek indel-adatok, mert ezt az eszközt kizárólag SNV-k azonosítására használják.

3.3. Szűrt eredmények a GiaB-vel összehasonlítva

A 4. táblázatban látható, hogy a standard szűrési gyakorlatoknak sikerül eltávolítaniuk az FP SNV-k nagy részét minden egyes csővezeték esetében; úgy tűnik azonban, hogy ezek a szűrők a legtöbb esetben kissé túl agresszívek az SNV-k felismeréséhez, de nem elég szigorúak az indelek esetében. Ez nyilvánvalóvá válik, ha megnézzük az egyes adatkészletekben jelentett FP-k számának különbségeit. Például a Bowtie2 Freebayes-szel 72 570 FP SNV eltávolítását látja (98%-os csökkenés), de csak 1736 FP indel eltávolítását (70%-os csökkenés). Azt is meg kell jegyezni, hogy az alkalmazott szűrők a csővezetéktől függtek, és az egyes csővezetékeken belül többnyire hasonló SNV és indel FP csökkentést eredményeztek. Az egyetlen kivétel itt a CUSHAW3 igazításokból azonosított variánsok száma a többi igazításhoz képest: összességében a TP SNV-k száma alacsonyabb, az FP SNV-k száma magasabb, és ez az egyetlen olyan illesztő, amely a szűrés után több mint 1000 FP indelt eredményez.

Mivel a szűrés jelentősen csökkenti a TP-variánsok számát, a ritka, nagy hatású variánsok keresésekor – a CUSHAW3 és FreeBayes pipeline-ok kivételével – érdemes lehet kihagyni ezt a lépést. Ehelyett érdemes több időt fordítani egy olyan szűrési folyamatra, amely inkább a biológián, mint a statisztikán alapul. Például több értelme lehet időt szánni a nagy valószínűséggel nagy hatású variánsok egy kis listájának azonosítására: splice site mutációk, frameshiftset okozó indelek, csonka mutációk, stop-loss mutációk, vagy olyan génekben található mutációk, amelyekről ismert, hogy biológiailag relevánsak az adott fenotípus szempontjából. Átlagos TPR és érzékenység

Amint az 5. táblázatban látható, a pozitív prediktív érték (PPV) minden eszköz esetében – a CUSHAW3 kivételével – 91% és 99,9% között mozog az SNV-k esetében, de az átlagos érzékenység nagyon alacsony (50% körüli). Ennek az eltérésnek számos oka lehet, de a legvalószínűbb az exonok közötti lefedettség változó mélysége. Láthatjuk, hogy az alacsony SNV-érzékenység mellett az indel-érzékenység is alacsony (30% körüli); az indelekre vonatkozó PPV azonban jóval alacsonyabb (35,86% és 52,95% között). Ennek a következő okok bármelyike lehet az oka: a nagyon rövid indeleket a hagyományos NGS-sel nehéz kimutatni , az indelek különböző variánshívók általi ábrázolása miatt az eszközök tévesen állíthatják, hogy két indel különböző, vagy az összehangoló eszközök ugyanazt az indelt különbözőképpen ábrázolják.

Talán a legvalószínűbb magyarázat mindkét típusú mutációra a mélység kérdése. Mint minden variánselemző vizsgálat esetében, a lefedettség mélységének növelése az érzékenység növekedéséhez vezet, de lehetetlen garantálni a jó lefedettségi mélységet, mivel az exome capture kitek nem képesek egységesen lehúzni az exonokat . Ezenkívül egyetlen exome capture kit sem fed le minden exont. Valójában kimutatták, hogy az exomnál alacsonyabb átlagos mélységű teljes genomszekvenálás variánselemzése az említett mélység egyenletessége miatt jobban teljesít. Így valószínű, hogy nagyszámú variáns hiányzik, mivel a NIST-GiaB listát exomi és genomi szekvenálási adatok összeállításából állították össze. Végső soron a megfelelő érzékenység eléréséhez előbb-utóbb teljes genomszekvenálást kell végezni, de ez jelenleg a legtöbb laboratórium számára megfizethetetlen. Szerencsére ez a költség folyamatosan csökken, és hamarosan fokozatos elmozdulást fogunk látni az exom-elemzésről a teljesebb teljes genom-elemzésre.

3.5. Érzékenység a mélység függvényében

Mivel az érzékenység egy eszköz egyik legfontosabb teljesítménymetrikáját tükrözi, és a legtöbb eszköz küzd az 50%-nál nagyobb érzékenység elérésével, szeretnénk tovább vizsgálni, hogy a mélység hogyan befolyásolja a variánshívás érzékenységét. Számos különböző eszközkombinációt vizsgáltunk, hogy meghatározzuk, melyek a legjobb pipelinek, variánshívók és illesztőprogramok. A 2. ábrához a variánshívók és illesztőprogramok öt legjobb kombinációját vettük alapul az érzékenységük és a hamis pozitív arányuk (FPR) alapján. Vagyis azokat választottuk ki, amelyeknél a legkevesebb FP SNV mellett a legtöbb TP SNV-t hívták ki. A vizsgálat során azonnal feltűnik, hogy az érzékenység alacsonyabb a vártnál. Az összes pipeline nagyjából azonos szinten teljesít: a legtöbb variánst a körülbelül 150x-es mélység eléréséig azonosítják, ami azt jelzi, hogy ez a mélység valószínűleg elegendő, és a hiányzó variánsok száma valószínűleg annak köszönhető, hogy bizonyos exonok az átlagosnál alacsonyabb lefedettséggel rendelkeznek. Megjegyzendő, hogy az öt legjobban teljesítő pipeline közül négyben a GATK UnifiedGenotyper a variánshívó, ami a 3. b) ábrán látható módon a használt alignertől függetlenül bizonyítja annak kiváló teljesítményét.

2. ábra

Szenzitivitás a mélység függvényében a legjobb öt pipeline esetében. Az öt legjobb csővezeték látható itt, ahol az egyes SNV-k mélységét ábrázoljuk az érzékenység függvényében.

(a)

(b)

(a)
(b)

3. ábra

Az érzékenység a mélység függvényében a top aligner és a top variant caller esetében. (a) Minden SNV mélységének eredményei az érzékenység függvényében ábrázolva a top aligner, a BWA mem és az összes variánshívó párosításában. (b) Minden SNV mélységének eredményei a legjobb variánshívó, a GATK UnifiedGenotyper, és az összes igazító párosítás érzékenységének függvényében.

Az öt legjobb pipeline vizsgálatán túlmenően úgy döntöttünk, hogy hasznos lenne ugyanezt az elemzést elvégezni a legjobb igazító és az összes variánshívó párosításával (3. ábra (a)), valamint a legjobb variánshívó és az összes igazító párosításával (3. ábra (b)). A csővezetékekhez hasonlóan a legjobb igazítót úgy azonosítottuk, mint amelyik a legtöbb TP SNV-t és a legkevesebb FP SNV-t eredményezte – ebben az esetben a BWA mem. Annak ellenére, hogy a legjobb igazítással dolgozhatunk, még mindig elég nagy különbséget látunk a variánshívók között, ami valószínűleg az általuk alkalmazott különböző algoritmusoknak tulajdonítható (3(a) ábra). A legjobban teljesítő variánshívó, a GATK UnifiedGenotyper esetében azonban úgy tűnik, hogy kevesebb eltérés van a legjobb négy igazító között, ami azt jelzi, hogy a legtöbb helyzetben meglehetősen jól teljesít, a CUSAHW3 és a MOSAIK.

3.6 kivételével. Közös variánsok a legjobb csővezetékek között

Végül arra voltunk kíváncsiak, hogy a különböző csővezetékek között mennyire egyediek a variánshívási készletek. Ehhez ismét a legjobb öt variánshívó pipeline-ra összpontosítottunk: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper és Novoalign plus UnifiedGenotyper. Amint az a 4. ábrán látható, a szóban forgó öt különböző pipeline között nagy az átfedés, az összesen 22 324 különböző variánsból 15 489 SNV (70%) közös. Ugyanakkor azzal is érvelhetünk, hogy ez nagyrészt annak köszönhető, hogy az öt pipeline közül négy az UnifiedGenotyper-t használja variánshívóként. Ezt az elképzelést alátámasztja az a tény, hogy a legnagyobb számú, 367 egyedi variáns a BWA sampe és a HaplotypeCaller kombinációjához tartozik. Azt is érdemes megjegyezni, hogy a második legtöbb egyedi SNV szintén a BWA sampe alignerhez tartozik, így lehetséges, hogy az egyedi SNV-k magas száma inkább az alignernek, mint a variánshívónak tulajdonítható.

4. ábra

Az első öt pipeline által azonosított SNV-k metszete.

4. Következtetések

Megállapítottuk, hogy a tesztelt harminc különböző pipeline közül a Novoalign plus GATK UnifiedGenotyper mutatta a legnagyobb érzékenységet, miközben az SNV-k alacsony FP-számát tartotta fenn. Az igazítók közül a BWA mem következetesen a legjobban teljesített, de az eredmények még mindig nagymértékben változtak a használt variánshívótól függően. Ebből természetesen következik, hogy a legjobb variánshívó, a GATK UnifiedGenotyper többnyire hasonló eredményeket produkált, függetlenül az alkalmazott illesztőprogramtól. Az azonban jól látható, hogy az indeleket még mindig nehéz kezelni, mivel egyik csővezeték sem ért el 33%-nál magasabb átlagos érzékenységet vagy 53%-nál magasabb PPV-t. Az indelek felismerésében tapasztalt alacsony általános teljesítmény mellett az érzékenység – a mutáció típusától függetlenül – az ebben a tanulmányban bemutatott összes csővezeték esetében problémát jelent. Az NA12878 exomjának várható SNV-száma 34 886, de még a legjobb öt pipeline által azonosított összes variáns egyesítése esetén is nagyon alacsony (22 324) volt a legnagyobb azonosított szám. Úgy tűnik, hogy bár még mindig nagyon hasznos az exómaelemzés, vannak korlátai, még akkor is, ha olyan látszólag egyszerű dologról van szó, mint az SNV-k felderítése.

Felvilágosítás

Adam Cornish Chittibabu Guda laboratóriumában végzett hallgató, informatikai és genomikai képzéssel. Chittibabu Guda (docens) interdiszciplináris háttérrel rendelkezik a molekuláris és számítógépes biológia területén. Számos számítási módszert publikált 2001 óta, amelyeknek számos alkalmazása van az orvosbiológiai kutatásban.

Érdekütközések ütközése

A szerzőknek nincs tudomása konkurens érdekekről.

A szerzők hozzájárulása

Adam Cornish tervezte a tanulmányt, elvégezte az összes elemzést, elkészítette az ábrákat, és megírta a cikket. Chittibabu Guda alapvető visszajelzéseket adott a dolgozat javításához, valamint inputot adott magukhoz az elemzésekhez, és alaposan megszerkesztette a dolgozatot. Minden szerző elolvasta és jóváhagyta a végleges papírt.

Köszönet

Ezt a munkát a University of Nebraska Medical Center (UNMC) Chittibabu Guda számára nyújtott fejlesztési támogatással támogatta. A szerzők köszönetet mondanak a Novoalign készítőinek a szoftver kipróbálható változatának rendelkezésre bocsátásáért, valamint az UNMC Bioinformatikai és Rendszerbiológiai magjának a kutatáshoz nyújtott infrastrukturális támogatásért.