Abstract

Secvențierea de mare randament, în special a exomilor, este un instrument de diagnosticare popular, dar este dificil de determinat care instrumente sunt cele mai bune pentru a analiza aceste date. În acest studiu, folosim rezultatele NIST Genome in a Bottle ca o resursă nouă pentru validarea pipeline-ului nostru de analiză a exomului. Utilizăm șase alinieri diferiți și cinci apelatori de variante diferiți pentru a determina care dintre cele 30 de conducte, din totalul de 30, funcționează cel mai bine pe un exom uman care a fost utilizat pentru a ajuta la generarea listei de variante detectate de Consorțiul Genome in a Bottle. Dintre aceste 30 de conducte, am constatat că Novoalign împreună cu GATK UnifiedGenotyper a prezentat cea mai mare sensibilitate, menținând în același timp un număr redus de falsuri pozitive pentru SNV-uri. Cu toate acestea, este evident că indelii sunt în continuare dificil de gestionat pentru orice pipeline, niciunul dintre instrumente nereușind să atingă o sensibilitate medie mai mare de 33% sau o valoare predictivă pozitivă (PPV) mai mare de 53%. În cele din urmă, așa cum era de așteptat, s-a constatat că aliniatorii pot juca un rol la fel de vital în detectarea variantelor ca și apelanții de variante în sine.

1. Context

În ultimii câțiva ani s-au făcut multe progrese în tehnologiile de secvențiere de mare capacitate. Datorită acestor progrese, este acum posibil să se detecteze un număr mare de variante potențial cauzatoare de boli , și, în câteva cazuri, datele de secvențiere de generație următoare (NGS) au fost chiar utilizate în scopuri de diagnosticare . Acest lucru se datorează parțial evoluției tehnologiilor de secvențiere din ultimii ani, dar și numărului de îmbunătățiri aduse diverselor instrumente bioinformatice utilizate pentru a analiza munții de date produse de instrumentele NGS .

Când se caută mutații la un pacient, un flux de lucru tipic este de a secvenția exomul acestuia cu un secvențiator Illumina, de a alinia datele brute la genomul uman de referință și apoi de a identifica variantele de nucleotide unice (SNV) sau inserțiile și delețiile scurte (indels) care ar putea provoca sau influența fenotipul de interes . În timp ce acest lucru este destul de simplu, nu este la fel de simplu să se decidă care sunt cele mai bune instrumente de utilizat în fiecare etapă a procesului de analiză. Există un număr mare de instrumente care sunt utilizate în diferite etape intermediare, dar cele două etape cele mai importante din întregul proces sunt alinierea citirilor brute la genom și apoi căutarea variantelor (de exemplu, SNV-uri și indeli) . În acest studiu, ne propunem să ajutăm bioinformaticianul de astăzi prin elucidarea combinației corecte a instrumentului de aliniere a citirilor scurte și a instrumentului de apelare a variantelor pentru procesarea datelor de secvențiere a exomului produse de instrumentele NGS.

În trecut au fost efectuate o serie de astfel de studii, dar toate au avut dezavantaje de o formă sau alta. În mod ideal, ar trebui să se dispună de o listă cu fiecare variantă cunoscută conținută într-o probă, astfel încât, atunci când se execută o conductă de instrumente de analiză, să o puteți testa pentru a ști cu certitudine că aceasta funcționează corect. Cu toate acestea, în trecut nu exista o astfel de listă, astfel încât validarea a trebuit să fie realizată prin metode mai puțin complete. În unele cazuri, validarea a fost efectuată prin generarea de date simulate, astfel încât să se creeze un set de adevărați pozitivi (TP) și adevărați negativi (TN) cunoscuți . Deși acest lucru oferă în mod convenabil o listă a fiecărui TP și TN din setul de date, nu reușește să reprezinte cu exactitate biologia. Alte metode de validare a conductelor de apelare a variantelor includ utilizarea matricelor de genotipare sau a secvențierii Sanger pentru a obține o listă de TP și de false pozitive (FP) . Aceste metode au avantajul de a furniza rezultate validate din punct de vedere biologic, dar au și dezavantajul de a nu fi exhaustive din cauza numărului limitat de puncte de pe rețelele de genotipare și a costului prohibitiv al validării Sanger atunci când se efectuează de mii de ori. În cele din urmă, niciunul dintre aceste studii nu a urmărit să analizeze efectul pe care alignerul de citire scurtă l-a avut asupra apelării variantelor. Prin urmare, efectul în amonte al performanței alignerului nu a putut fi evaluat independent.

În acest studiu, avem avantajul unei liste de variante pentru o femeie anonimă din Utah (ID subiect: NA12878, secvențiată inițial pentru proiectul 1000 Genomes ) care a fost validată experimental de către consorțiul Genome in a Bottle (GiaB) condus de NIST. Această listă de variante a fost creată prin integrarea a 14 seturi de date diferite de la cinci secvențiatoare diferite și ne permite să validăm orice listă de variante generată de conductele noastre de analiză a exomului . Noutatea acestei lucrări constă în validarea combinației corecte de aliniere și de apelare a variantelor în raport cu un set de date de variante cuprinzător și determinat experimental: NIST-GiaB.

Pentru a efectua analiza noastră, vom utiliza unul dintre seturile de date exomice utilizate inițial pentru a crea lista NIST-GiaB. Am ales doar unul dintre exomii generați inițial de Illumina TruSeq pentru că am dorit să oferim un scenariu de utilizare standard pentru cineva care dorește să efectueze o analiză NGS și, în timp ce secvențierea genomului întreg continuă să scadă în preț, secvențierea exomului este încă o alternativă populară și viabilă . De asemenea, este important de remarcat faptul că, conform lui Bamshad et al., în prezent, numărul preconizat de SNV-uri per exom european-american este de 20 283 ± 523 . În ciuda acestui fapt, numărul total de SNV-uri găsite în lista NIST-GiaB cu potențialul de a exista în setul de date exome TruSeq a fost de 34 886, ceea ce este semnificativ mai mare decât se aștepta. Acest lucru se datorează probabil faptului că, în timp ce kitul exomic a fost utilizat pentru a genera datele NIST-GiaB, acesta a fost, de asemenea, completat de secvențierea întregului genom.

În cele din urmă, am luat în considerare un număr mare de aliniere și de apelatoare de variante, dar în cele din urmă am ales cele 11 instrumente pe baza prevalenței, popularității și relevanței pentru setul nostru de date (de exemplu, SNVMix, VarScan2 și MuTect nu au fost utilizate deoarece sunt destinate utilizării pe probe derivate din tumori). Analiza noastră propriu-zisă implică compararea a șase aliniere (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK și Novoalign) și a cinci apelatoare de variante (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup și SNPSVM ). În acest studiu încercăm, de asemenea, să determinăm cât de mult efect, dacă este cazul, are alignerul asupra apelării variantelor și care aliniatori au cele mai bune performanțe atunci când se utilizează o probă normală de exom Illumina. Din câte știm, acesta este primul raport care validează toate combinațiile posibile (în total 30 de pipe-line-uri) ale unei game largi de aliniere și de apelare a variantelor.

2. Methods

2.1. Methods

2.1. Seturi de date

Genomul uman de referință hg19 a fost descărcat din browserul UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) și a fost utilizat pentru a efectua alinierile. Exomul uman, SRR098401, a fost descărcat de la Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Pentru adnotare și calibrare, au fost utilizate dbSNP137 fără site-uri după versiunea 129, HapMap 3.3, Human Omni 2.5 BeadChip și listele de seturi de indeluri Mills și 1000 G gold standard (toate de la ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. The Pipeline

Figura 1 prezintă fluxul de lucru utilizat în acest studiu, care este similar cu cel prezentat în ghidul Best Practices realizat de The Broad Institute . Acesta implică o serie de etape pentru a se asigura că fișierele de aliniere produse sunt de cea mai bună calitate, precum și alte câteva etape pentru a garanta că variantele sunt numite corect. În primul rând, citirile brute au fost aliniate la hg19, iar apoi duplicatele PCR au fost eliminate din aliniere. Apoi, pentru a ajuta la identificarea indelurilor mai târziu în cadrul procesului, s-a efectuat realinierea citirilor în jurul indelurilor. Ultimul pas al procesării alinierii a fost efectuarea unei etape de recalibrare a scorului de calitate a bazelor, care ajută la ameliorarea prejudecăților și inexactităților inerente ale scorurilor emise de secvențiatoare. Din păcate, în ciuda acestor etape, rata de aliniere a fiecărui aligner a fost semnificativ mai mică decât se aștepta, astfel încât, pentru a compensa acest lucru, a fost utilizat setul de instrumente fastx pentru a filtra citirile de calitate scăzută (tabelul 1). Citirile de calitate scăzută au fost definite ca fiind acele citiri care au avut cel puțin jumătate din scorurile de calitate sub 30. După procesarea alinierii, au fost efectuate apelarea variantelor și filtrarea variantelor.

Figura 1

Schema de analiză a datelor utilizate. Pentru a se asigura că sunt create alinieri de cea mai bună calitate, citirile sunt mai întâi filtrate și apoi aliniate la genomul uman de referință, hg19. Apoi, duplicatele PCR sunt eliminate, citirile sunt aliniate în jurul indelilor putative, iar scorurile de calitate a bazelor sunt recalibrate. În cele din urmă, variantele sunt numite și validate în raport cu lista de variante NIST-GiaB.

Cele șase instrumente utilizate pentru a genera alinieri au fost Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK și Novoalign, iar cele cinci instrumente utilizate pentru a genera variante au fost FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup și SNPSVM, după cum se poate observa în tabelul 2.

2.3. Filtrare

Datele brute au fost achiziționate de la SRA (SRR098401), împărțite cu fastq-dump și filtrate cu ajutorul setului de instrumente fastx. Mai exact, fastq-dump a folosit indicatoarele -split_files și -split_spot, iar fastq_quality_filter a fost rulat cu următoarele argumente: -Q 33 -q 30 -p 50. Apoi, citirile au fost împerecheate corespunzător cu un script personalizat.

2.4. Aliniere

Aliniatoarele au folosit argumente implicite, cu excepția cazului în care a fost folosit un argument pentru fire atunci când a fost disponibil. Comenzile utilizate sunt după cum urmează.

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

Pentru a asigura un calcul corect al adâncimii de acoperire, s-a utilizat modulul Picard Tools CalculateHsMetrics cu următoarele argumente:(1)java -jar CalculateHsMetrics.jar I=NA12878.ALN.BQSR.bam O=ALN.O.log R=genome/hg19.fa TI=genome/genome/truseq_exome.bed BI=genome/truseq_exome.bed VALIDATION_STRINGENCY=SILENT PER_TARGET_COVERAGE=ALN.ptc.bedEste important de reținut că fișierul de pat al exomului TruSeq trebuie să aibă antetul din fișierul de aliniere SAM atașat în prealabil pentru ca acest modul să funcționeze. Mai mult, coloana 6 trebuie mutată în coloana 4, iar coloana 5 trebuie să fie eliminată din fișierul de pat TruSeq.

2.6. Prelucrarea fișierelor de aliniere

Procesarea fișierelor de aliniere (fișiere SAM/BAM) a necesitat următorii pași pentru toate alinierele:(1)Conversia SAM în BAM cu SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Sortarea fișierelor BAM cu ajutorul modulului 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=coordinate(3)Eliminarea duplicatelor PCR utilizând modulul 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)Read Group adăugat la fișierele de aliniere utilizând modulul 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)Realinierea în jurul indelilor folosind modulele GATK RealignerTargetCreator și 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)Recalibrarea bazei folosind modulele GATK BaseRecalibrator și 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. Apelarea variantelor

Argumentele implicite au fost utilizate pentru fiecare apelator de variante, cu excepția cazului în care acesta conținea un indicator „threads” sau „parallel”, caz în care a fost utilizat și acesta. În plus, indels au fost apelate separat de SNV-uri atunci când a fost posibil. În mod specific, comenzile utilizate sunt următoarele.

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.vcfDin cauza inexistenței stegulețelor CIGAR necesare în fișierul de aliniere, SNPSVM nu a reușit să apeleze variantele pentru CUSHAW3, iar SAMtools mpileup nu a putut apela variantele pe alinierile MOSAIK din același motiv. De asemenea, datorită faptului că SNPSVM detectează numai SNV-uri, nu au fost raportate indels pentru acest program.

2.8. Filtrarea variantelor

Filtrarea a variat în funcție de apelatorul de variante utilizat. În cazul GATK HaplotypeCaller și GATK UnifiedGenotyper, modulele GATK, VariantRecalibrator și ApplyRecalibration, au fost utilizate pentru a filtra SNV-urile folosind HapMap 3.3, Omni 2.5 SNP BeadChip și dbSNP 137 fără datele 1000 Genome ca seturi de antrenament. Pentru SNPSVM, au fost utilizate scoruri QUAL ≥ 4 și valori DP ≥ 6. Pentru FreeBayes și SAMtools, s-au utilizat scoruri QUAL ≥ 20 și valori DP ≥ 6.

2.9. Compararea variantelor

Pentru compararea variantelor, s-a utilizat USeq 8.8.1 pentru a compara SNV-urile partajate între toate seturile de date. Pentru a compara indels, a fost utilizat instrumentul vcflib vcfintersect. Fișierul TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, obținut în 11 decembrie 2013, a fost utilizat pentru a restricționa comparațiile la locațiile care au putut fi capturate de kitul de extragere a exomului utilizat în secvențierea lui SRR098401. Acest fișier poate fi obținut de la Illumina aici: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Pentru a se asigura că variantele au fost reprezentate în mod identic între diferite seturi de apeluri, instrumentul vcflib vcfallelicprimitives a fost utilizat pentru a preprocesa fișierele vcf.

2.10. Calcule statistice

True Positive (TP). Este o mutație care a fost detectată de pipeline-ul testat și este una care există în lista NIST-GiaB.

False Positive (FP). Este o mutație care a fost detectată de către conducta testată, dar care nu există în lista NIST-GiaB.

True Negative (TN). Este o mutație care nu a fost detectată de pipeline-ul testat și care nu există în lista NIST-GiaB.

False Negative (FN). Este o mutație care nu a fost detectată de pipeline-ul testat, dar este una care există în lista NIST-GiaB:

3. Rezultate și discuții

3.1. Prefiltrarea variantelor

Când se efectuează analiza variantelor, una dintre numeroasele capcane care trebuie luate în considerare este spațiul secvențial al exomului (așa cum este definit de kitul de captare a exomului) și modul în care acesta poate afecta rezultatele analizei. În acest caz, am avut un singur exom (SRR098401) care a fost extras cu ajutorul kitului Illumina TruSeq exome și secvențiat pe un HiSeq 2000. Având în vedere acest lucru, am vrut să ne asigurăm că măsurăm capacitatea instrumentelor bioinformatice de a-și face treaba și nu cât de bine a funcționat kitul de captare a exomului Illumina TruSeq. Altfel spus, dorim să încercăm să apelăm doar variantele care ar trebui să fie prezente în exoni, așa cum sunt definite de kitul de extragere. Din acest motiv, folosim fișierul bed furnizat de Illumina, nu un fișier bed de adnotare generic, de exemplu, RefSeq pentru hg19. Am constatat că pentru acest individ anume, în conformitate cu lista NIST-GiaB, ar trebui să existe un total de 34 886 SNV și 1 473 indels în regiunile definite de fișierul de pat TruSeq.

După ce am filtrat variantele care nu se aflau în regiunile definite de fișierul de pat al exomului Illumina TruSeq, am trecut de la sute de mii de variante presupuse (datele nu sunt prezentate) la, în medie, aproximativ 23 000 de variante (SNV și indels) pe conductă (tabelul 3). Acesta este un pas important pentru cercetători pentru început, deoarece reduce semnificativ spațiul de căutare a variantelor potențial interesante.

3.2. Rezultatele privind variantele brute în comparație cu GiaB

Un aspect pe care am dorit să îl înțelegem atunci când am făcut această comparație a fost importanța filtrării variantelor detectate de aceste instrumente. Motivul este că, în mod ideal, s-ar dori un nivel de sensibilitate cât mai ridicat posibil, astfel încât mutațiile de interes să nu se piardă în procesul de filtrare. Prin urmare, este de datoria noastră să determinăm dacă acest pas este sau nu necesar și în ce măsură este necesar, deoarece din rezultatele NIST-GiaB și din analiza Bamshad et al. reiese clar că sensibilitatea ar putea fi o problemă.

După cum se poate vedea în tabelul 3, filtrarea este necesară mai mult pentru unele dispozitive de apelare a variantelor decât pentru altele atunci când vine vorba de SNV-uri, și este absolut necesară pentru indeli. În cele mai multe cazuri, numărul de variante TP este apropiat sau mai mare decât numărul nostru așteptat de aproximativ 20 000 , dar, pe de altă parte, în unele cazuri, numărul de FP este foarte mare.

Este clar că există o mare variație în numerele generate de fiecare pipeline. Cu toate acestea, se pot găsi unele puncte comune în aceste numere, care provin probabil din originile algoritmice ale fiecărui instrument. FreeBayes produce atât cel mai mare număr de variante nefiltrate, cât și cel mai mare număr de FP-uri. Este probabil ca acest tip de performanță să fie obținut doar de acest instrument datorită faptului că, deși nu este singurul instrument de apelare a variantelor bazat pe inferența Bayesiană, este unic în interpretarea sa a alinierilor. Adică, este un apelator bazat pe haplotipuri care identifică variantele pe baza secvenței de citire în sine în loc de aliniere, acesta din urmă fiind modul în care funcționează UnifiedGenotyper al GATK.

În plus, observăm că aliniatoarele bazate pe Burrows-Wheeler au performanțe foarte asemănătoare între ele: Bowtie2, BWA mem și BWA sampe obțin rezultate similare pe toată linia. Se poate presupune că acest lucru se datorează probabil faptului că toate aceste instrumente utilizează algoritmi similari atunci când își îndeplinesc sarcina desemnată. Această observație este susținută de faptul că MOSAIK (alinieri lacunare care utilizează algoritmul Smith-Waterman) și CUSHAW3 (o abordare hibridă de însămânțare) au ambii algoritmi de bază foarte diferiți și, ulterior, produc rezultate foarte diferite.

Această diferență în rezultatele corelate cu diferiți algoritmi se observă cel mai bine în rezultatele SNPSVM. Dintre dispozitivele de apelare a variantelor, acesta este singurul care utilizează mașini vectoriale de suport și construirea de modele pentru a genera apeluri SNV. S-ar părea că, deși are dezavantajul de a nu fi la fel de sensibil ca alte metode, beneficiază de faptul că este extrem de precis, indiferent de alignerul utilizat. Acest lucru sugerează că se poate sări peste etapa de filtrare cu totul atunci când se utilizează acest apelant de variante.

În ceea ce privește indelurile, niciun aligner nu pare să se evidențieze printre celelalte ca fiind unul care gestionează bine acest tip de mutație. De fapt, atunci când se analizează numărul de FP, este clar că apelantul de variante este cel care joacă cel mai mare rol în acuratețea identificării indelurilor. În plus, nu există date pentru conductele CUSHAW3 plus SNPSVM și MOSAIK plus SAMtools mpileup din cauza faptului că fișierele de aliniere nu conțin șirurile CIGAR necesare pentru ca apelatorii de variante să funcționeze în aval. În cele din urmă, motivul pentru care nu există date indel pentru SNPSVM este că acest instrument este utilizat exclusiv pentru identificarea SNV-urilor.

3.3. Rezultate filtrate în comparație cu GiaB

Ca în tabelul 4, practicile standard de filtrare reușesc să elimine un număr mare de SNV-uri FP pentru fiecare conductă; cu toate acestea, se pare că aceste filtre sunt un pic prea agresive în majoritatea cazurilor pentru detectarea SNV-urilor, dar nu sunt suficient de stricte pentru indeli. Acest lucru devine evident atunci când se analizează diferențele în ceea ce privește numărul de FP-uri raportate în fiecare set de date. De exemplu, Bowtie2 cu Freebayes constată o eliminare a 72 570 de SNV-uri FP (o reducere de 98%), dar numai o eliminare a 1 736 de indeli FP (o reducere de 70%). De asemenea, trebuie remarcat faptul că filtrele utilizate au fost dependente de pipeline și, în cea mai mare parte, în cadrul fiecărui pipeline au produs reduceri similare ale SNV și indel FP. Singura excepție aici este numărul de variante identificate din alinierile CUSHAW3 în comparație cu alte alinieri: în general, numărul de SNV TP este mai mic, numărul de SNV FP este mai mare și este singurul aligner care produce mai mult de 1.000 de indeli FP după filtrare.

Datorită faptului că filtrarea reduce semnificativ numărul de variante TP, ar putea fi înțelept ca, cu excepția conductelor care utilizează CUSHAW3 și FreeBayes, să se sară peste acest pas atunci când se caută variante rare, cu impact ridicat. În schimb, s-ar putea dedica mai mult timp unui proces de filtrare care se bazează mai degrabă pe biologie decât pe statistică. De exemplu, ar putea avea mai mult sens să se investească timp pentru a identifica o listă mică de variante care sunt susceptibile de a fi de mare impact: mutații ale situsului de îmbinare, indels care cauzează schimbări de cadre, mutații de trunchiere, mutații de stop-loss sau mutații în gene despre care se știe că sunt relevante din punct de vedere biologic pentru fenotipul de interes.

3.4. Valoarea predictivă pozitivă (PPV) medie și sensibilitatea

După cum se poate observa în tabelul 5, valoarea predictivă pozitivă (PPV) pentru fiecare instrument, cu excepția CUSHAW3, variază între 91 % și 99,9 % pentru SNV-uri, dar sensibilitatea medie este foarte scăzută (aproximativ 50 %). Această discrepanță s-ar putea datora mai multor motive, dar cel mai probabil este adâncimea variabilă a acoperirii între exoni. Putem observa că, în plus față de sensibilitatea scăzută a SNV, sensibilitatea indel este scăzută (aproximativ 30 %); cu toate acestea, PPV pentru indel este considerabil mai mică (35,86 % până la 52,95 %). Acest lucru s-ar putea datora oricăruia dintre următoarele motive: indelurile foarte scurte sunt greu de detectat prin NGS convențional , reprezentarea indelurilor de către diferite apelatoare de variante poate face ca instrumentele să afirme în mod incorect că două indeluri sunt diferite sau instrumentele de aliniere produc reprezentări diferite ale aceluiași indel .

Poate că cea mai probabilă explicație pentru ambele tipuri de mutații este problema profunzimii. La fel ca în cazul oricărui studiu de analiză a variantelor, o creștere a profunzimii de acoperire duce la o creștere a sensibilității, dar este imposibil să se garanteze o bună profunzime de acoperire din cauza incapacității kiturilor de captare a exomului de a extrage în mod uniform exoni . În plus, niciun kit de captare a exomului nu acoperă fiecare exon. Într-adevăr, s-a demonstrat că analiza variantelor din secvențierea întregului genom la o adâncime medie mai mică decât un exom are rezultate mai bune datorită uniformității acestei adâncimi. Astfel, este probabil ca un număr mare de variante să lipsească din cauza faptului că lista NIST-GiaB a fost creată dintr-o compilație de date de secvențiere exomică și genomică. În cele din urmă, pentru a obține o sensibilitate adecvată, va trebui în cele din urmă să se efectueze secvențierea întregului genom, dar acest lucru este în prezent prohibitiv din punct de vedere al costurilor pentru majoritatea laboratoarelor. Din fericire, acest cost continuă să scadă și, în curând, vom asista la o trecere treptată de la analiza exomului la analiza mai completă a genomului întreg.

3.5. Sensibilitatea în funcție de profunzime

Pentru că sensibilitatea reflectă unul dintre cei mai importanți parametri de performanță ai unui instrument și majoritatea instrumentelor se străduiesc să atingă o sensibilitate mai mare de 50%, am dori să explorăm în continuare modul în care profunzimea a afectat sensibilitatea de apelare a variantelor. Am analizat o serie de combinații diferite de instrumente pentru a determina care sunt cele mai bune conducte, apelatoare de variante și aliniatoare. Pentru figura 2, am luat cele mai bune cinci combinații de apelanți de variante și de aliniere, așa cum au fost determinate de sensibilitatea și de rata falsurilor pozitive (FPR). Altfel spus, le-am selectat pe cele care au avut cel mai mare număr de SNV TP apelate în plus față de cel mai mic număr de SNV FP. La inspecție, lucrul care iese imediat în evidență este că sensibilitatea este mai mică decât se aștepta. Toate conductele performează aproximativ la același nivel: acestea identifică majoritatea variantelor în momentul în care se atinge o adâncime de aproximativ 150x, ceea ce indică faptul că această adâncime este probabil suficientă și că numărul de variante lipsă se datorează probabil faptului că anumiți exoni au o acoperire mai mică decât media. Rețineți că patru din cele cinci conducte cu cele mai bune performanțe au GATK UnifiedGenotyper ca apelant de variante, demonstrând performanța sa superioară indiferent de alignerul utilizat, după cum se arată în figura 3(b).

Figura 2

Sensibilitatea în funcție de adâncime pentru cele cinci conducte de top. Primele cinci conducte sunt prezentate aici cu adâncimea fiecărui SNV reprezentată grafic în funcție de sensibilitate.

(a)

(b)

(a)
(b)

Figura 3

Sensibilitatea în funcție de adâncime pentru aliniatorul de top și apelatorul de variante de top. (a) Rezultatele pentru adâncimea fiecărui SNV reprezentate grafic în funcție de sensibilitate pentru alignerul de top, BWA mem, asociat cu fiecare apelant de variante. (b) Rezultatele pentru adâncimea fiecărui SNV reprezentate în funcție de sensibilitate pentru cel mai bun apelant de variante, GATK UnifiedGenotyper, asociat cu fiecare aligner.

În plus față de examinarea celor mai bune cinci conducte, am stabilit că ar fi util să efectuăm aceeași analiză pe cel mai bun aligner asociat cu fiecare apelant de variante [Figura 3(a)], precum și pe cel mai bun apelant de variante asociat cu fiecare aligner [Figura 3(b)]. Ca și în cazul conductelor, cel mai bun aligner a fost identificat ca fiind cel care a produs cel mai mare număr de SNV TP și cel mai mic număr de SNV FP – în acest caz, BWA mem. În ciuda faptului că avem cea mai bună aliniere cu care să lucrăm, observăm totuși o diferență destul de mare între apelatorii de variante, care se poate atribui probabil algoritmilor diferiți pe care îi folosesc [Figura 3(a)]. Cu toate acestea, în cazul celui mai performant dispozitiv de apelare a variantelor, GATK UnifiedGenotyper, se pare că există o variație mai mică între primii patru aliniați, ceea ce indică faptul că acesta se comportă destul de bine în majoritatea situațiilor, excepțiile fiind CUSAHW3 și MOSAIK.

3.6. Variante partajate între cele mai bune conducte

În cele din urmă, am vrut să știm cât de unice au fost seturile de apelare a variantelor între diferitele conducte. Pentru a face acest lucru, ne-am concentrat din nou pe primele cinci conducte de apelare a variantelor: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper și Novoalign plus UnifiedGenotyper. După cum se poate observa în figura 4, există o cantitate mare de suprapunere între cele cinci pipe-line-uri diferite în cauză, cu 15 489 de SNV-uri (70%) partajate dintr-un total de 22 324 de variante distincte. Cu toate acestea, s-ar putea argumenta, de asemenea, că acest lucru se datorează în mare parte faptului că patru dintre cele cinci conducte utilizează UnifiedGenotyper ca apelant de variante. Această noțiune este coroborată de faptul că cel mai mare număr de variante unice pentru o conductă, 367, aparține combinației BWA sampe plus HaplotypeCaller. De asemenea, este demn de remarcat faptul că al doilea cel mai mare număr de SNV-uri unice aparține, de asemenea, alignerului BWA sampe, astfel încât este posibil ca numărul mare de SNV-uri unice să fie mai bine atribuit alignerului decât apelatorului de variante.

Figura 4

Intersecția SNV-urilor identificate de primele cinci pipe-line-uri.

4. Concluzii

Am constatat că, dintre cele treizeci de pipe-line-uri diferite testate, Novoalign plus GATK UnifiedGenotyper a prezentat cea mai mare sensibilitate, menținând în același timp un număr redus de FP pentru SNV-uri. Dintre aliniatoare, BWA mem a avut în mod constant cea mai bună performanță, dar rezultatele au variat totuși foarte mult în funcție de apelatorul de variante utilizat. În mod firesc, rezultă că cel mai bun apelant de variante, GATK UnifiedGenotyper, a produs în cea mai mare parte rezultate similare indiferent de alignerul utilizat. Cu toate acestea, este ușor de observat că indelii sunt în continuare dificil de gestionat pentru orice tip de pipe-line, niciunul dintre aceste pipe-line-uri nereușind să atingă o sensibilitate medie mai mare de 33% sau o VPP mai mare de 53%. În plus față de performanța generală scăzută pe care o observăm în detectarea indelilor, sensibilitatea, indiferent de tipul de mutație, reprezintă o problemă pentru fiecare conductă descrisă în această lucrare. Numărul preconizat de SNV-uri pentru exomul lui NA12878 este de 34 886, dar chiar și atunci când se utilizează uniunea tuturor variantelor identificate de primele cinci pipeline-uri, cel mai mare număr identificat a fost foarte mic (22 324). Se pare că, deși este încă foarte utilă, analiza exomului are limitările sale chiar și atunci când vine vorba de ceva atât de aparent simplu precum detectarea SNV.

Dezvăluiri

Adam Cornish este student absolvent în laboratorul lui Chittibabu Guda, cu pregătire în informatică și genomică. Chittibabu Guda (profesor asociat) are o pregătire interdisciplinară în domeniul biologiei moleculare și computaționale. El a publicat o serie de metode computaționale cu o varietate de aplicații în cercetarea biomedicală, începând din 2001.

Conflict de interese

Autorii nu au cunoștință de interese concurente.

Contribuția autorilor

Adam Cornish a proiectat studiul, a efectuat toate analizele, a realizat figurile și a scris articolul. Chittibabu Guda a oferit un feedback esențial privind îmbunătățirile aduse lucrării și contribuții la analizele în sine și a editat amănunțit lucrarea. Toți autorii au citit și aprobat lucrarea finală.

Recunoștințe

Această lucrare a fost susținută de fondurile de dezvoltare pentru Chittibabu Guda de la University of Nebraska Medical Center (UNMC). Autorii ar dori să mulțumească creatorilor Novoalign pentru că au pus la dispoziție software-ul ca versiune de încercare și facilitatea Bioinformatics and Systems Biology Core de la UNMC pentru suportul de infrastructură care a facilitat această cercetare.

.