Abstract

Il sequenziamento ad alta velocità, specialmente degli esomi, è uno strumento diagnostico popolare, ma è difficile determinare quali strumenti sono i migliori per analizzare questi dati. In questo studio, usiamo i risultati del NIST Genome in a Bottle come una nuova risorsa per la convalida della nostra pipeline di analisi degli esomi. Usiamo sei diversi allineatori e cinque diversi chiamatori di varianti per determinare quale pipeline, delle 30 totali, esegue il meglio su un esoma umano che è stato utilizzato per aiutare a generare l’elenco delle varianti rilevate dal Genome in a Bottle Consortium. Di queste 30 pipeline, abbiamo trovato che Novoalign in combinazione con GATK UnifiedGenotyper ha esibito la più alta sensibilità mantenendo un basso numero di falsi positivi per SNVs. Tuttavia, è evidente che gli indel sono ancora difficili da gestire per qualsiasi pipeline con nessuno degli strumenti che raggiunge una sensibilità media superiore al 33% o un valore predittivo positivo (PPV) superiore al 53%. Infine, come ci si aspettava, si è scoperto che gli allineatori possono giocare un ruolo vitale nell’individuazione delle varianti tanto quanto i chiamanti delle varianti stesse.

1. Sfondo

Negli ultimi anni sono stati fatti molti progressi nelle tecnologie di sequenziamento ad alta velocità. Grazie a questi progressi, è ora possibile rilevare un gran numero di potenziali varianti che causano malattie e, in alcuni casi, i dati del sequenziamento di prossima generazione (NGS) sono stati addirittura utilizzati per scopi diagnostici. Ciò è in parte dovuto agli sviluppi delle tecnologie di sequenziamento negli ultimi anni, ma anche al numero di miglioramenti apportati ai vari strumenti bioinformatici utilizzati per analizzare le montagne di dati prodotti dagli strumenti NGS.

Quando si cercano mutazioni in un paziente, un flusso di lavoro tipico è quello di sequenziare il loro esoma con un sequenziatore Illumina, allineare i dati grezzi al genoma umano di riferimento, e quindi identificare le varianti a singolo nucleotide (SNV) o brevi inserzioni e delezioni (indel) che potrebbero eventualmente causare o influenzare il fenotipo di interesse. Mentre questo è abbastanza semplice, decidere i migliori strumenti da utilizzare in ogni fase della pipeline di analisi non lo è. Ci sono un gran numero di strumenti che vengono utilizzati in varie fasi intermedie, ma le due fasi più importanti dell’intero processo sono l’allineamento delle letture grezze al genoma e poi la ricerca di varianti (cioè, SNVs e indel). In questo studio, ci proponiamo di aiutare il bioinformatico di oggi chiarendo la corretta combinazione di strumento di allineamento di lettura breve e strumento di chiamata delle varianti per l’elaborazione dei dati di sequenziamento dell’esoma prodotti da strumenti NGS.

Un certo numero di questi studi è stato eseguito in passato, ma tutti hanno avuto inconvenienti di una forma o un’altra. Idealmente si dovrebbe avere un elenco di ogni variante nota contenuta in un campione in modo che quando una pipeline di strumenti di analisi viene eseguita, è possibile testarla per sapere con certezza che sta eseguendo correttamente. Tuttavia, in passato non esisteva un tale elenco, quindi la convalida doveva essere eseguita con metodi meno completi. In alcuni casi, la convalida è stata eseguita generando dati simulati in modo da creare un insieme di veri positivi (TP) e veri negativi (TN) noti. Mentre questo fornisce convenientemente un elenco di ogni TP e TN nel set di dati, fa un lavoro povero di rappresentare accuratamente la biologia. Altri metodi di convalida delle pipeline di chiamata delle varianti includono l’utilizzo di array di genotipizzazione o il sequenziamento Sanger per ottenere un elenco di TP e falsi positivi (FP). Questi hanno il vantaggio di fornire risultati biologicamente validati, ma hanno anche il rovescio della medaglia di non essere completi a causa del numero limitato di spot sugli array di genotipizzazione e il costo proibitivo della convalida Sanger quando viene eseguita migliaia di volte. Infine, nessuno di questi studi mirava a esaminare l’effetto dell’allineatore di letture brevi sulla chiamata delle varianti. Di conseguenza, l’effetto a monte delle prestazioni dell’allineatore non poteva essere valutato in modo indipendente.

In questo studio, abbiamo il vantaggio di un elenco di varianti per una femmina anonima dello Utah (ID soggetto: NA12878, originariamente sequenziato per il progetto 1000 Genomi) che è stato convalidato sperimentalmente dal Consorzio Genome in a Bottle (GiaB) guidato dal NIST. Questa lista di varianti è stata creata integrando 14 diversi set di dati da cinque diversi sequenziatori, e ci permette di convalidare qualsiasi lista di varianti generate dalle nostre pipeline di analisi dell’esoma. La novità di questo lavoro è di convalidare la giusta combinazione di allineatori e chiamanti di varianti contro un set di dati di varianti completo e determinato sperimentalmente: NIST-GiaB.

Per eseguire la nostra analisi useremo uno dei dataset di esomi originariamente utilizzati per creare la lista NIST-GiaB. Abbiamo scelto solo uno degli esomi originali generati da Illumina TruSeq perché volevamo fornire uno scenario di caso d’uso standard per qualcuno che desidera eseguire l’analisi NGS, e mentre il sequenziamento del genoma intero continua a scendere di prezzo, il sequenziamento dell’esoma è ancora un’alternativa popolare e valida. È anche importante notare che, per Bamshad et al., attualmente il numero previsto di SNV per esoma europeo-americano è 20.283 ± 523 . Nonostante questo, il numero totale di SNVs trovati nella lista NIST-GiaB con il potenziale di esistere nel dataset di esomi TruSeq era 34.886, che è significativamente più alto del previsto. Questo è probabilmente dovuto al fatto che mentre il kit esoma è stato utilizzato per generare i dati NIST-GiaB è stato anche integrato dal sequenziamento dell’intero genoma.

Infine, abbiamo considerato un gran numero di allineatori e variant caller ma alla fine abbiamo scelto gli 11 strumenti basati sulla prevalenza, la popolarità e la rilevanza per il nostro set di dati (ad esempio, SNVMix, VarScan2 e MuTect non sono stati utilizzati perché sono destinati all’uso su campioni derivati dal tumore). La nostra analisi prevede il confronto di sei allineatori (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK e Novoalign) e cinque chiamatori di varianti (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup e SNPSVM ). In questo studio cerchiamo anche di determinare quanto di un effetto, se del caso, l’allineatore ha sulla chiamata variante e quali allineatori eseguire meglio quando si utilizza un campione esoma Illumina normale. A nostra conoscenza, questo è il primo rapporto che convalida tutte le possibili combinazioni (totale di 30 pipeline) di una vasta gamma di allineatori e chiamanti di varianti.

2. Metodi

2.1. Datasets

Il genoma umano di riferimento hg19 è stato scaricato dal browser UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) ed è stato utilizzato per eseguire gli allineamenti. L’esoma umano, SRR098401, è stato scaricato dal Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Per scopi di annotazione e calibrazione, dbSNP137 senza siti dopo la versione 129, HapMap 3.3, Human Omni 2.5 BeadChip, e Mills e 1000 G gold standard indel set list sono stati utilizzati (tutti da ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. La pipeline

La figura 1 mostra il flusso di lavoro utilizzato in questo studio, che è simile a quello delineato nella guida Best Practices prodotta da The Broad Institute . Ciò comporta una serie di passaggi per garantire che i file di allineamento prodotti siano della massima qualità, nonché diversi altri per garantire che le varianti siano chiamate correttamente. In primo luogo, le letture grezze sono state allineate a hg19, e poi i duplicati PCR sono stati rimossi dall’allineamento. Successivamente, per aiutare l’identificazione degli indelebili nella pipeline, è stato eseguito il riallineamento delle letture intorno agli indelebili. L’ultimo passo dell’elaborazione dell’allineamento è stato quello di eseguire un passo di ricalibrazione del punteggio di qualità della base, che aiuta a migliorare la distorsione intrinseca e le imprecisioni dei punteggi emessi dai sequenziatori. Sfortunatamente, nonostante questi passaggi, il tasso di allineamento di ogni allineatore era significativamente più basso del previsto, quindi per compensare questo, il toolkit fastx è stato utilizzato per filtrare le letture di bassa qualità (Tabella 1). Le letture di bassa qualità sono state definite come quelle che avevano almeno la metà dei loro punteggi di qualità inferiori a 30. Dopo l’elaborazione dell’allineamento, sono state eseguite la chiamata delle varianti e il filtraggio delle varianti.

Figura 1

Schema della pipeline di analisi dati utilizzata. Per garantire che vengano creati allineamenti di altissima qualità, le letture vengono prima filtrate e poi allineate al genoma umano di riferimento, hg19. Successivamente, i duplicati PCR vengono rimossi, le letture vengono allineate intorno agli indel putativi e i punteggi di qualità delle basi vengono ricalibrati. Infine, le varianti sono chiamate e convalidate rispetto alla lista di varianti NIST-GiaB.

I sei strumenti utilizzati per generare gli allineamenti erano Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK e Novoalign, e i cinque strumenti utilizzati per generare le varianti erano FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup e SNPSVM, come si può vedere nella tabella 2.

2.3. Filtraggio

I dati grezzi sono stati acquisiti dal SRA (SRR098401), suddivisi con fastq-dump e filtrati utilizzando il toolkit fastx. In particolare, fastq-dump ha utilizzato i flag -split_files e -split_spot, e fastq_quality_filter è stato eseguito con i seguenti argomenti: -Q 33 -q 30 -p 50. Poi le letture sono state correttamente accoppiate con uno script personalizzato.

2.4. Allineamento

Gli allineatori hanno usato argomenti predefiniti tranne quando è stato usato un argomento threads dove disponibile. I comandi utilizzati sono i seguenti.

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. Calcolo della profondità di copertura dell’allineamento

Per assicurare un corretto calcolo della profondità di copertura, è stato utilizzato il modulo Picard Tools CalculateHsMetrics con i seguenti argomenti:(1)java -jar CalculateHsMetrics.jar I=NA12878.ALN.BQSR.bam O=ALN.O.log R=genoma/hg19.fa TI=genome/truseq_exome.bed BI=genome/truseq_exome.bed VALIDATION_STRINGENCY=SILENT PER_TARGET_COVERAGE=ALN.ptc.bedE’ importante notare che il file TruSeq exome bed deve avere l’intestazione del file di allineamento SAM aggiunta ad esso affinché questo modulo funzioni. Inoltre, la colonna 6 deve essere spostata nella colonna 4 e la colonna 5 deve essere rimossa dal file letto TruSeq.

2.6. Elaborazione dei file di allineamento

L’elaborazione dei file di allineamento (file SAM/BAM) richiede i seguenti passi per tutti gli allineatori:(1)Conversione da SAM a BAM con SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Ordinamento dei file BAM utilizzando il modulo 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)Rimozione duplicati PCR usando il modulo 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)Gruppo di lettura aggiunto ai file di allineamento usando il modulo 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)Riallineamento intorno agli indel usando i moduli GATK RealignerTargetCreator e IndelRealigner:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T RealignerTargetCreator -R genoma/hg19.fa -I allineamenti/NA12878.ALN.RG.bam – genoma noto/mills.vcf -o tmp/ALN.intervals(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T IndelRealigner -R genoma/hg19.fa -I allineamenti/NA12878.ALN.RG.bam -conosciuto genoma/mills.vcf -o allineamenti/NA12878.ALN.indels.bam -maxReadsForRealignment 100000 -maxReadsInMemory 1000000 -targetIntervals tmp/ALN.intervals(6)Ricalibrazione della base utilizzando i moduli GATK BaseRecalibrator e 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. Chiamata delle varianti

Gli argomenti predefiniti sono stati usati per ogni chiamante di varianti a meno che non contenesse un flag “threads” o “parallel”, nel qual caso è stato usato anche quello. Inoltre, gli indel sono stati chiamati separatamente dagli SNV quando possibile. Nello specifico, i comandi utilizzati sono i seguenti.

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 allineamenti/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 allineamenti/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 vista -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.vcfA causa dell’inesistenza dei flag CIGAR richiesti nel file di allineamento, SNPSVM non è riuscito a chiamare le varianti per CUSHAW3, e SAMtools mpileup non ha potuto chiamare le varianti sugli allineamenti MOSAIK per lo stesso motivo. Inoltre, dato che SNPSVM rileva solo gli SNV, non sono stati riportati indel per questo programma.

2.8. Filtrazione delle varianti

La filtrazione varia a seconda del variant caller utilizzato. Nel caso di GATK HaplotypeCaller e GATK UnifiedGenotyper, i moduli GATK, VariantRecalibrator e ApplyRecalibration, sono stati utilizzati per filtrare gli SNV utilizzando HapMap 3.3, Omni 2.5 SNP BeadChip e dbSNP 137 senza i dati di 1000 Genome come set di allenamento. Per SNPSVM, sono stati utilizzati punteggi QUAL ≥ 4 e valori DP ≥ 6. Per FreeBayes e SAMtools, sono stati utilizzati punteggi QUAL ≥ 20 e valori DP ≥ 6.

2.9. Confronto delle varianti

Per il confronto delle varianti, è stato utilizzato USeq 8.8.1 per confrontare gli SNV condivisi tra tutti i set di dati. Per confrontare gli indel, è stato utilizzato lo strumento vcflib vcfintersect. Il file TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, ottenuto nel dicembre 11, 2013, è stato utilizzato per limitare i confronti alle posizioni che potrebbero essere catturati dal kit exome pull down utilizzato nel sequenziamento di SRR098401. Questo file può essere ottenuto da Illumina qui: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Per garantire che le varianti fossero rappresentate in modo identico tra diversi set di chiamate, lo strumento vcflib vcfallelicprimitives è stato utilizzato per preprocessare i file vcf.

2.10. Calcoli statistici

True Positive (TP). Si tratta di una mutazione che è stata rilevata dalla pipeline in esame e che esiste nella lista NIST-GiaB.

Falso positivo (FP). Si tratta di una mutazione che è stata rilevata dalla pipeline in fase di test, ma che non esiste nella lista NIST-GiaB.

Vero Negativo (TN). Si tratta di una mutazione che non è stata rilevata dalla pipeline in esame e che non esiste nella lista NIST-GiaB.

Falso negativo (FN). Si tratta di una mutazione che non è stata rilevata dalla pipeline in esame ma che esiste nella lista NIST-GiaB:

3. Risultati e Discussione

3.1. Prefiltraggio delle varianti

Quando si esegue l’analisi delle varianti, una delle molte insidie che devono essere prese in considerazione è lo spazio di sequenza dell’esoma (come definito dal kit di cattura dell’esoma) e come può influenzare i risultati dell’analisi. In questo caso, avevamo un singolo esoma (SRR098401) che è stato estratto utilizzando il kit Illumina TruSeq exome e sequenziato su un HiSeq 2000. Con questo in mente, volevamo essere sicuri che stavamo misurando la capacità degli strumenti bioinformatici di fare il loro lavoro e non quanto bene ha funzionato il kit Illumina TruSeq di acquisizione dell’esoma. Cioè, vogliamo solo provare a chiamare le varianti che si suppone siano presenti negli esoni come definito dal kit di pull down. Per questo motivo, usiamo il file letto fornito da Illumina, non un generico file letto di annotazione, per esempio, RefSeq per hg19. Abbiamo trovato che per questo particolare individuo, secondo la lista NIST-GiaB, ci dovrebbe essere un totale di 34.886 SNVs e 1.473 indel all’interno delle regioni definite dal file TruSeq bed.

Una volta che abbiamo filtrato le varianti che non si trovavano nelle regioni definite dal file Illumina TruSeq exome bed, siamo passati da centinaia di migliaia di varianti putative (dati non mostrati) a, in media, circa 23.000 varianti (SNVs e indel) per pipeline (Tabella 3). Questo è un passo importante per i ricercatori per iniziare, in quanto riduce significativamente lo spazio di ricerca per le varianti potenzialmente interessanti.

3.2. Raw Variant Results Compared to GiaB

Un aspetto che volevamo capire quando abbiamo fatto questo confronto era l’importanza di filtrare le varianti rilevate da questi strumenti. Il motivo è che idealmente si vorrebbe avere un livello di sensibilità il più alto possibile in modo che le mutazioni di interesse non vadano perse nel processo di filtraggio. È quindi opportuno determinare se questo passo è necessario o meno e in che misura è necessario, poiché è chiaro dai risultati di NIST-GiaB e dalla revisione di Bamshad et al. che la sensibilità potrebbe essere un problema.

Come possiamo vedere nella tabella 3, il filtraggio è necessario più per alcuni variant caller che per altri quando si tratta di SNV, ed è assolutamente necessario per gli indel. Nella maggior parte dei casi, il numero di varianti TP è vicino o superiore al nostro numero previsto di circa 20.000, ma, d’altra parte, in alcuni casi il numero di FP è molto alto.

E’ chiaro che c’è molta variazione nei numeri generati da ogni pipeline. Tuttavia, si possono trovare alcuni punti in comune nei numeri che probabilmente derivano dalle origini algoritmiche di ogni strumento. FreeBayes produce sia il maggior numero di varianti non filtrate che il maggior numero di FP. È probabile che vediamo questo tipo di prestazioni solo da questo strumento a causa del fatto che, pur non essendo l’unico variant caller basato sull’inferenza bayesiana, è unico nella sua interpretazione degli allineamenti. Cioè, è un chiamatore basato sugli aplotipi che identifica le varianti basate sulla sequenza delle letture stesse invece che sull’allineamento, che è il modo in cui opera UnifiedGenotyper di GATK.

Inoltre vediamo che gli allineatori basati su Burrows-Wheeler hanno prestazioni molto simili tra loro: Bowtie2, BWA mem e BWA sampe ottengono risultati simili su tutta la linea. Si potrebbe supporre che questo è probabilmente dovuto al fatto che tutti questi strumenti utilizzano algoritmi simili quando eseguono il loro compito designato. Questa osservazione è supportata dal fatto che MOSAIK (allineamenti gapped utilizzando l’algoritmo Smith-Waterman) e CUSHAW3 (un approccio ibrido di semina) hanno entrambi algoritmi di base molto diversi e di conseguenza producono risultati molto diversi.

Questa differenza di risultati correlati a diversi algoritmi si vede meglio nei risultati SNPSVM. Tra i variant caller, è l’unico che utilizza macchine vettoriali di supporto e la costruzione di modelli per generare chiamate SNV. Sembrerebbe che, pur avendo lo svantaggio di non essere così sensibile come altri metodi, benefici di essere estremamente accurato indipendentemente dall’allineatore utilizzato. Questo suggerisce che si è in grado di saltare del tutto la fase di filtraggio quando si usa questo variant caller.

Per quanto riguarda gli indel, nessun allineatore sembra distinguersi tra gli altri come uno che gestisce bene questo tipo di mutazione. Infatti, quando si guarda al numero di FP, è chiaro che è il variant caller che gioca il ruolo maggiore nell’accuratezza dell’identificazione degli indel. Inoltre, non ci sono dati né per CUSHAW3 più SNPSVM né per MOSAIK più SAMtools mpileup pipeline a causa dei file di allineamento che non contengono le stringhe CIGAR necessarie ai variant caller per funzionare a valle. Infine, il motivo per cui non ci sono dati indel per SNPSVM è che questo strumento è utilizzato esclusivamente per l’identificazione di SNVs.

3.3. Risultati filtrati rispetto a GiaB

Come nella tabella 4, le pratiche di filtraggio standard riescono a rimuovere un gran numero di FP SNVs per ogni pipeline; tuttavia sembra che questi filtri siano un po’ troppo aggressivi nella maggior parte dei casi per la rilevazione degli SNV, ma non abbastanza severi per gli indel. Questo è reso ovvio quando si osservano le differenze nel numero di FP riportati in ogni set di dati. Per esempio, Bowtie2 con Freebayes vede una rimozione di 72.570 FP SNVs (una riduzione del 98%) ma solo una rimozione di 1.736 FP indel (una riduzione del 70%). Va anche notato che i filtri utilizzati erano dipendenti dalla pipeline e, per la maggior parte, all’interno di ogni pipeline hanno prodotto riduzioni simili in SNV e indel FP. L’unica eccezione è il numero di varianti identificate dagli allineamenti CUSHAW3 rispetto agli altri allineamenti: Nel complesso il numero di TP SNVs è inferiore, il numero di FP SNVs è superiore, ed è l’unico allineatore che produce più di 1.000 FP indel dopo il filtraggio.

Dato che il filtraggio riduce significativamente il numero di varianti TP, potrebbe essere saggio, ad eccezione delle pipeline che utilizzano CUSHAW3 e FreeBayes, saltare questo passaggio quando si cercano varianti rare e ad alto impatto. Invece, si potrebbe spendere più tempo in un processo di filtraggio che si basa sulla biologia piuttosto che sulla statistica. Per esempio, potrebbe avere più senso investire del tempo nell’identificazione di un piccolo elenco di varianti che sono probabilmente ad alto impatto: mutazioni del sito di giunzione, indel che causano frameshift, mutazioni di troncamento, mutazioni di stop-loss, o mutazioni in geni che sono noti per essere biologicamente rilevanti per il fenotipo di interesse.

3.4. TPR e sensibilità media

Come si può vedere nella tabella 5, il valore predittivo positivo (PPV) per ogni strumento, ad eccezione di CUSHAW3, varia dal 91% al 99,9% per gli SNV, ma la sensibilità media è molto bassa (circa il 50%). Questa discrepanza potrebbe essere dovuta a una serie di ragioni, ma la più probabile è la profondità variabile di copertura tra gli esoni. Possiamo vedere che, oltre alla bassa sensibilità SNV, la sensibilità indel è bassa (circa il 30%); tuttavia il PPV per indel è notevolmente inferiore (35,86% a 52,95%). Questo potrebbe essere dovuto a una delle seguenti ragioni: gli indel molto corti sono difficili da rilevare da NGS convenzionale, la rappresentazione degli indel da parte di diversi chiamanti di varianti può causare strumenti per sostenere erroneamente che due indel sono diversi, o strumenti di allineamento producono rappresentazioni diverse dello stesso indel.

Forse la spiegazione più probabile per entrambi i tipi di mutazioni è la questione della profondità. Come nel caso di qualsiasi studio di analisi delle varianti, un aumento della profondità di copertura porta a un aumento della sensibilità, ma è impossibile garantire una buona profondità di copertura a causa dell’incapacità dei kit di cattura dell’esoma di tirare giù uniformemente gli esoni. Inoltre, nessun singolo kit di cattura dell’esoma copre ogni esone. Infatti, è stato dimostrato che l’analisi delle varianti del sequenziamento dell’intero genoma a una profondità media inferiore a un esoma esegue meglio a causa dell’uniformità di detta profondità. Quindi, è probabile che manchi un gran numero di varianti a causa del fatto che la lista NIST-GiaB è stata creata da una compilazione di dati di sequenziamento esomico e genomico. In definitiva, per ottenere una sensibilità adeguata, sarà necessario eseguire il sequenziamento dell’intero genoma, ma questo è attualmente proibitivo in termini di costi per la maggior parte dei laboratori. Fortunatamente, questo costo sta continuando a scendere, e presto vedremo un graduale passaggio dall’analisi degli esomi alla più completa analisi del genoma intero.

3.5. Sensibilità come funzione della profondità

Poiché la sensibilità riflette una delle più importanti metriche di performance di uno strumento e la maggior parte degli strumenti lotta per raggiungere una sensibilità superiore al 50%, vorremmo esplorare ulteriormente come la profondità ha influenzato la sensibilità della chiamata delle varianti. Abbiamo esaminato una serie di diverse combinazioni di strumenti per determinare quali fossero le migliori pipeline, i variant caller e gli allineatori. Per la figura 2, abbiamo preso le cinque migliori combinazioni di variant caller e allineatori come determinato dalla loro sensibilità e dal tasso di falsi positivi (FPR). Cioè, abbiamo selezionato quelli che avevano il maggior numero di TP SNVs chiamati oltre al minor numero di FP SNVs. All’ispezione, la cosa che salta subito all’occhio è che la sensibilità è più bassa del previsto. Tutte le pipeline si comportano più o meno allo stesso livello: identificano la maggior parte delle loro varianti quando viene raggiunta una profondità di circa 150x, il che indica che questa profondità è probabilmente sufficiente e che il numero di varianti mancanti è probabilmente dovuto a certi esoni che hanno una copertura inferiore alla media. Si noti che quattro dei cinque migliori pipeline hanno GATK UnifiedGenotyper come loro variante chiamante, dimostrando le sue prestazioni superiori indipendentemente dall’allineatore utilizzato come mostrato nella Figura 3 (b).

Figura 2

Sensibilità come funzione di profondità per i primi cinque pipeline. Le prime cinque pipeline sono mostrate qui con la profondità di ogni SNV tracciata rispetto alla sensibilità.

(a)

(b)

(a)
(b)

Figura 3

Sensibilità in funzione della profondità per il top aligner e il top variant caller. (a) I risultati per la profondità di ogni SNV tracciati contro la sensibilità per l’allineatore superiore, BWA mem, accoppiato con ogni variant caller. (b) Risultati per la profondità di ogni SNV tracciati contro la sensibilità per il miglior chiamante di varianti, GATK UnifiedGenotyper, accoppiato con ogni allineatore.

Oltre a guardare le prime cinque pipeline, abbiamo determinato che sarebbe utile eseguire la stessa analisi sul miglior allineatore accoppiato con ogni chiamante di varianti (Figura 3(a)), così come il miglior chiamante di varianti accoppiato con ogni allineatore (Figura 3(b)). Come per le pipeline, il miglior allineatore è stato identificato come quello che ha prodotto il maggior numero di TP SNVs e il minor numero di FP SNVs – in questo caso BWA mem. Pur avendo il miglior allineamento con cui lavorare, vediamo ancora una differenza abbastanza grande tra i chiamanti di varianti, che è probabilmente attribuibile ai diversi algoritmi che impiegano (Figura 3(a)). Tuttavia, nel caso del variant caller più performante, GATK UnifiedGenotyper, sembra esserci meno variazione tra i primi quattro allineatori, indicando che si comporta abbastanza bene nella maggior parte delle situazioni, con le eccezioni di CUSAHW3 e MOSAIK.

3.6. Varianti condivise tra le migliori pipeline

Infine, abbiamo voluto sapere quanto fossero unici i set di varianti chiamate tra le diverse pipeline. Per fare questo, ci siamo di nuovo concentrati sulle cinque principali pipeline di chiamata delle varianti: Bowtie2 più UnifiedGenotyper, BWA mem più UnifiedGenotyper, BWA sampe più HaplotypeCaller, BWA sampe più UnifiedGenotyper, e Novoalign più UnifiedGenotyper. Come si può vedere nella Figura 4, c’è una grande quantità di sovrapposizione tra le cinque diverse pipeline in questione, con 15.489 SNV (70%) condivisi su un totale di 22.324 varianti distinte. Tuttavia, si potrebbe anche sostenere che questo è in gran parte dovuto a quattro delle cinque pipeline che utilizzano UnifiedGenotyper come loro variante caller. Questa nozione è corroborata dal fatto che il maggior numero di varianti uniche per una pipeline, 367, appartiene alla combinazione BWA sampe più HaplotypeCaller. Vale anche la pena notare che il secondo maggior numero di SNV unici appartiene anch’esso all’allineatore BWA sampe, quindi è possibile che l’alto numero di SNV unici sia meglio attribuito all’allineatore che al variant caller.

Figura 4

L’intersezione degli SNV identificati dalle prime cinque pipeline.

4. Conclusioni

Abbiamo scoperto che tra le trenta diverse pipeline testate Novoalign più GATK UnifiedGenotyper ha esibito la più alta sensibilità mantenendo un basso numero di FP per gli SNV. Tra gli allineatori, BWA mem ha sempre dato il meglio, ma i risultati variavano ancora notevolmente a seconda del variant caller utilizzato. Naturalmente, ne consegue che il miglior variant caller, GATK UnifiedGenotyper, ha prodotto per lo più risultati simili indipendentemente dall’allineatore utilizzato. Tuttavia, è subito evidente che gli indel sono ancora difficili da gestire per qualsiasi pipeline con nessuna delle pipeline che raggiunge una sensibilità media superiore al 33% o un PPV superiore al 53%. Oltre alla bassa performance complessiva che vediamo nel rilevamento degli indel, la sensibilità, indipendentemente dal tipo di mutazione, è un problema per ogni pipeline descritta in questo articolo. Il numero previsto di SNV per l’esoma di NA12878 è 34.886, ma anche quando si usa l’unione di tutte le varianti identificate dalle cinque migliori pipeline, il maggior numero identificato è molto basso (22.324). Sembra che, pur essendo ancora molto utile, l’analisi dell’esoma abbia i suoi limiti anche quando si tratta di qualcosa di apparentemente semplice come il rilevamento di SNV.

Discorso

Adam Cornish è uno studente laureato nel laboratorio di Chittibabu Guda con formazione in informatica e genomica. Chittibabu Guda (professore associato) ha un background interdisciplinare in biologia molecolare e computazionale. Ha pubblicato una serie di metodi computazionali con una varietà di applicazioni nella ricerca biomedica, dal 2001.

Conflitto di interessi

Gli autori non sono a conoscenza di interessi concorrenti.

Contributo degli autori

Adam Cornish ha progettato lo studio, eseguito tutte le analisi, fatto le figure, e scritto il documento. Chittibabu Guda ha fornito un feedback essenziale sui miglioramenti dell’articolo e l’input sulle analisi stesse e ha modificato completamente l’articolo. Tutti gli autori hanno letto e approvato l’articolo finale.

Riconoscimenti

Questo lavoro è stato sostenuto dai fondi di sviluppo a Chittibabu Guda dall’Università del Nebraska Medical Center (UNMC). Gli autori desiderano ringraziare i creatori di Novoalign per aver reso disponibile il software in versione di prova e la struttura Bioinformatica e Systems Biology Core dell’UNMC per il supporto infrastrutturale che ha facilitato questa ricerca.