Una comparación de las tuberías de llamada de variantes utilizando Genome in a Bottle como referencia
Abstract
La secuenciación de alto rendimiento, especialmente de exomas, es una herramienta de diagnóstico popular, pero es difícil determinar qué herramientas son las mejores para analizar estos datos. En este estudio, utilizamos los resultados de NIST Genome in a Bottle como un recurso novedoso para la validación de nuestra tubería de análisis de exomas. Utilizamos seis alineadores diferentes y cinco llamadores de variantes diferentes para determinar qué línea de trabajo, de un total de 30, funciona mejor en un exoma humano que se utilizó para ayudar a generar la lista de variantes detectadas por el Consorcio Genome in a Bottle. De estos 30 pipelines, descubrimos que Novoalign, junto con GATK UnifiedGenotyper, mostró la mayor sensibilidad, manteniendo un bajo número de falsos positivos para SNVs. Sin embargo, es evidente que los indels siguen siendo difíciles de manejar para cualquier canalización, ya que ninguna de las herramientas logró una sensibilidad media superior al 33% o un valor predictivo positivo (VPP) superior al 53%. Por último, como se esperaba, se descubrió que los alineadores pueden desempeñar un papel tan vital en la detección de variantes como los propios llamadores de variantes.
1. Antecedentes
En los últimos años se han producido muchos avances en las tecnologías de secuenciación de alto rendimiento. Debido a estos avances, ahora es posible detectar un gran número de variantes potencialmente causantes de enfermedades , y, en unos pocos casos, los datos de secuenciación de próxima generación (NGS) se han utilizado incluso con fines de diagnóstico . Esto se debe en parte a la evolución de las tecnologías de secuenciación en los últimos años, pero también a la cantidad de mejoras introducidas en las diversas herramientas bioinformáticas utilizadas para analizar las montañas de datos producidos por los instrumentos de NGS .
Cuando se buscan mutaciones en un paciente, un flujo de trabajo típico consiste en secuenciar su exoma con un secuenciador Illumina, alinear los datos brutos con el genoma humano de referencia y, a continuación, identificar las variantes de un solo nucleótido (SNV) o las inserciones y deleciones cortas (indels) que podrían causar o influir en el fenotipo de interés . Aunque esto es bastante sencillo, no lo es decidir cuáles son las mejores herramientas para cada etapa del proceso de análisis. Hay un gran número de herramientas que se utilizan en varios pasos intermedios, pero los dos pasos más importantes de todo el proceso son la alineación de las lecturas brutas con el genoma y la búsqueda de variantes (es decir, SNVs e indels). En este estudio, pretendemos ayudar al bioinformático de hoy en día dilucidando la combinación correcta de la herramienta de alineación de lecturas cortas y la herramienta de búsqueda de variantes para procesar los datos de secuenciación del exoma producidos por los instrumentos de NGS.
En el pasado se han realizado varios estudios de este tipo, pero todos tenían inconvenientes de una forma u otra. Lo ideal sería tener una lista de todas las variantes conocidas contenidas en una muestra, de manera que cuando se ejecute un pipeline de herramientas de análisis, se pueda probar para saber con certeza que está funcionando correctamente. Sin embargo, en el pasado no existía tal lista, por lo que la validación tenía que realizarse mediante métodos menos completos. En algunos casos, la validación se realizaba generando datos simulados para crear un conjunto de verdaderos positivos (TP) y verdaderos negativos (TN) conocidos. Aunque esto proporciona convenientemente una lista de todos los TP y TN en el conjunto de datos, no hace un buen trabajo para representar con precisión la biología. Otros métodos para validar las líneas de llamadas de variantes incluyen el uso de matrices de genotipado o la secuenciación Sanger para obtener una lista de TP y falsos positivos (FP). Estos métodos tienen la ventaja de proporcionar resultados biológicamente validados, pero también tienen la desventaja de no ser exhaustivos debido al número limitado de puntos en las matrices de genotipado y el coste prohibitivo de la validación de Sanger cuando se realiza miles de veces. Por último, ninguno de estos estudios pretendía analizar el efecto del alineador de lecturas cortas en la determinación de variantes. En consecuencia, el efecto del rendimiento del alineador no pudo ser evaluado de forma independiente.
En este estudio, tenemos la ventaja de una lista de variantes para una mujer anónima de Utah (ID del sujeto: NA12878, originalmente secuenciado para el proyecto 1000 Genomes ) que fue validado experimentalmente por el Consorcio Genome in a Bottle (GiaB) dirigido por el NIST. Esta lista de variantes se creó integrando 14 conjuntos de datos diferentes procedentes de cinco secuenciadores distintos, y nos permite validar cualquier lista de variantes generada por nuestros pipelines de análisis del exoma . La novedad de este trabajo es validar la combinación correcta de alineadores y llamadores de variantes frente a un conjunto de datos de variantes completo y determinado experimentalmente: NIST-GiaB.
Para realizar nuestro análisis utilizaremos uno de los conjuntos de datos de exomas utilizados originalmente para crear la lista NIST-GiaB. Elegimos sólo uno de los exomas originales generados por Illumina TruSeq porque queríamos proporcionar un escenario de uso estándar para alguien que desee realizar un análisis NGS, y mientras la secuenciación del genoma completo sigue bajando de precio, la secuenciación del exoma sigue siendo una alternativa popular y viable . También es importante señalar que, según Bamshad et al., actualmente el número esperado de SNVs por exoma europeo-americano es de 20.283 ± 523 . A pesar de esto, el número total de SNVs encontrados en la lista NIST-GiaB con el potencial de existir en el conjunto de datos del exoma TruSeq fue de 34.886, lo que es significativamente mayor de lo esperado. Esto se debe probablemente al hecho de que, aunque el kit de exoma se utilizó para generar datos NIST-GiaB, también se complementó con la secuenciación del genoma completo.
Por último, consideramos un gran número de alineadores y llamadores de variantes, pero finalmente elegimos las 11 herramientas basándonos en la prevalencia, la popularidad y la relevancia para nuestro conjunto de datos (por ejemplo, SNVMix, VarScan2 y MuTect no se utilizaron ya que están destinados a ser utilizados en muestras derivadas de tumores). Nuestro análisis en sí implica la comparación de seis alineadores (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK , y Novoalign) y cinco llamadores de variantes (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup , y SNPSVM ). En este estudio también intentamos determinar el efecto, si es que hay alguno, del alineador en la llamada de variantes y qué alineadores funcionan mejor cuando se utiliza una muestra normal del exoma de Illumina. Hasta donde sabemos, este es el primer informe que valida todas las combinaciones posibles (un total de 30 pipelines) de una amplia gama de alineadores y llamadores de variantes.
2. Métodos
2.1. Datasets
El genoma humano de referencia hg19 se descargó del navegador de la UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) y se utilizó para realizar los alineamientos. El exoma humano, SRR098401, se descargó del Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Para fines de anotación y calibración, se utilizaron dbSNP137 sin sitios después de la versión 129, HapMap 3.3, Human Omni 2.5 BeadChip, y las listas de conjuntos de indeles estándar de oro de Mills y 1000 G (todas de ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).
2.2. La figura 1 muestra el flujo de trabajo utilizado en este estudio, que es similar al descrito en la guía de mejores prácticas elaborada por el Instituto Broad. Esto implica una serie de pasos para asegurar que los archivos de alineación producidos son de la más alta calidad, así como varios más para garantizar que las variantes se llaman correctamente. En primer lugar, las lecturas en bruto se alinearon con hg19 y, a continuación, se eliminaron del alineamiento los duplicados de PCR. A continuación, para ayudar a la identificación de indels más adelante en el proceso, se realizó una realineación de lecturas alrededor de los indels. El último paso del proceso de alineación fue realizar un paso de recalibración de la puntuación de calidad de las bases, que ayuda a mejorar el sesgo inherente y las inexactitudes de las puntuaciones emitidas por los secuenciadores. Lamentablemente, a pesar de estos pasos, la tasa de alineación de cada alineador fue significativamente inferior a la esperada, por lo que para compensar esto, se utilizó el kit de herramientas fastx para filtrar las lecturas de baja calidad (Tabla 1). Las lecturas de baja calidad se definieron como aquellas que tenían al menos la mitad de sus puntuaciones de calidad por debajo de 30. Tras el procesamiento de la alineación, se realizó la llamada de variantes y el filtrado de variantes.
Aligner
% lecturas alineadas (sin filtrar)
% de lecturas alineadas (filtradas)
Profundidad media de cobertura
Bowtie2
89.73
98.73
47.97
BWA mem
92.91
99.85
46.89
BWA sampe
85,95
97,49
46,67
CUSHAW3
85.00
99,81
47,69
MOSAIK
85,68
96,22
45.14
Novoalign
82,21
94,20
45,62
Tabla 1
Porcentajes de alineación para lecturas filtradas y sin filtrar. La profundidad media de cobertura corresponde a los archivos de alineación creados con las lecturas filtradas.
Figura 1
Esquema de la línea de análisis de datos utilizada. Para garantizar la creación de alineaciones de máxima calidad, las lecturas se filtran primero y se alinean con el genoma humano de referencia, hg19. A continuación, se eliminan los duplicados de la PCR, se alinean las lecturas alrededor de los indels putativos y se recalibran las puntuaciones de calidad de las bases. Por último, las variantes se llaman y se validan con la lista de variantes NIST-GiaB.
Las seis herramientas utilizadas para generar alineaciones fueron Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK y Novoalign, y las cinco herramientas utilizadas para generar variantes fueron FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup y SNPSVM, como puede verse en la Tabla 2.
Herramienta
Tipo
Versión
Referencia
Bowtie2
Alineador
2.1.0
BWA sampe
Aligner
0.7.5a
BWA mem
Aligner
0.7.5a
CUSHAW3
Aligner
3.0.3
MOSAIK
Aligner
2.2.3
Novoalign
Aligner
3.02.07
N/A
FreeBayes
Genotyper
v9.9.2-19-g011561f
GATK HaplotypeCaller
Genotyper
2.7-2
N/A
GATK UnifiedGenotyper
Genotyper
2.7-2
SAMtools mpileup
Genotyper
0.1.19
SNPSVM
Genotyper
0.01
Tabla 2
Estas son las 11 herramientas diferentes utilizadas que componen los 30 (seis alineadores cinco llamadores de variantes) pipelines diferentes. También se incluyen las versiones del software para garantizar la reproducibilidad.
2.3. Filtrado
Los datos crudos se adquirieron del SRA (SRR098401), se dividieron con fastq-dump y se filtraron utilizando el kit de herramientas fastx. En concreto, fastq-dump utilizó los indicadores -split_files y -split_spot, y fastq_quality_filter se ejecutó con los siguientes argumentos: -Q 33 -q 30 -p 50. A continuación, las lecturas se emparejaron correctamente con un script personalizado.
2.4. Alineación
Los alineadores utilizaron los argumentos por defecto excepto cuando se utilizó un argumento de hilos cuando estaba disponible. Los comandos utilizados son los siguientes.
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. Cálculo de la profundidad de cobertura de los alineamientos
Para garantizar un cálculo adecuado de la profundidad de cobertura, se utilizó el módulo CalculateHsMetrics de Picard Tools con los siguientes argumentos:(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.bedEs importante tener en cuenta que el archivo TruSeq exome bed debe tener la cabecera del archivo de alineación SAM como anexo para que este módulo funcione. Además, la columna 6 debe trasladarse a la columna 4, y la columna 5 debe eliminarse del archivo de lecho de TruSeq.
2.6. Procesamiento de archivos de alineación
El procesamiento de los archivos de alineación (archivos SAM/BAM) requiere los siguientes pasos para todos los alineadores:(1)Conversión de SAM a BAM con SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Clasificación de archivos BAM con el módulo de 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=coordinada(3)Eliminación de duplicados de PCR mediante el módulo de herramientas Picard, MarkDuplicates:(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=alineaciones/NA12878.ALN.sorted.bam O=alineaciones/NA12878.ALN.dups_removed.bam REMOVE_DUPLICATES=true M=alignments/metrics(4)Grupo de lectura añadido a los archivos de alineación mediante el módulo de herramientas Picard, AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=coordinación RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Realineación en torno a los indeles utilizando los módulos GATK RealignerTargetCreator e 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)Recalibración de bases utilizando los módulos GATK BaseRecalibrator y PrintReads:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T BaseRecalibrator -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -sitiosconocidosgenomedbsnp_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. Llamada de variantes
Se utilizaron los argumentos por defecto para cada llamador de variantes a menos que contuviera una bandera «threads» o «parallel» en cuyo caso también se utilizó. Además, los indels se llamaron por separado de los SNV cuando fue posible. En concreto, los comandos utilizados son los siguientes:
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.vcfDebido a la inexistencia de las banderas CIGAR requeridas en el archivo de alineación, SNPSVM no pudo llamar a las variantes para CUSHAW3, y SAMtools mpileup no pudo llamar a las variantes en las alineaciones MOSAIK por la misma razón. Además, debido al hecho de que SNPSVM sólo detecta SNVs, no se reportaron indels para este programa.
2.8. Filtración de variantes
La filtración varió dependiendo del llamador de variantes que se utilizó. En los casos de GATK HaplotypeCaller y GATK UnifiedGenotyper, se utilizaron los módulos de GATK, VariantRecalibrator y ApplyRecalibration, para filtrar SNVs utilizando HapMap 3.3, el Omni 2.5 SNP BeadChip, y dbSNP 137 sin datos de 1000 Genome como conjuntos de entrenamiento. Para SNPSVM, se utilizaron puntuaciones QUAL ≥ 4 y valores DP ≥ 6. Para FreeBayes y SAMtools, se utilizaron puntuaciones QUAL ≥ 20 y valores DP ≥ 6.
2.9. Comparación de variantes
Para la comparación de variantes, se utilizó USeq 8.8.1 para comparar SNVs compartidos entre todos los conjuntos de datos. Para comparar indels, se utilizó la herramienta vcflib vcfintersect. El archivo TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, obtenido el 11 de diciembre de 2013, se utilizó para restringir las comparaciones a las ubicaciones que podían ser capturadas por el kit pull down del exoma utilizado en la secuenciación de SRR098401. Este archivo se puede obtener de Illumina aquí: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Para garantizar que las variantes estuvieran representadas de forma idéntica entre diferentes conjuntos de llamadas, se utilizó la herramienta vcflib vcfallelicprimitives para preprocesar los archivos vcf.
2.10. Cálculos estadísticos
Verdaderos positivos (TP). Es una mutación que fue detectada por el pipeline que se está probando y es una que existe en la lista NIST-GiaB.
Falso Positivo (FP). Es una mutación que fue detectada por la tubería que se está probando pero es una que no existe en la lista NIST-GiaB.
Verdadero negativo (TN). Es una mutación que no fue detectada por la tubería que se está probando y es una que no existe en la lista NIST-GiaB.
Falso negativo (FN). Es una mutación que no fue detectada por la tubería que se está probando pero que sí existe en la lista NIST-GiaB:
3. Resultados y Discusión
3.1. Prefiltrado de variantes
Cuando se realiza un análisis de variantes, uno de los muchos escollos que hay que tener en cuenta es el espacio de la secuencia del exoma (tal y como lo define el kit de captura del exoma) y cómo puede afectar a los resultados del análisis. En este caso, teníamos un único exoma (SRR098401) que fue extraído utilizando el kit de exoma TruSeq de Illumina y secuenciado en un HiSeq 2000. Con esto en mente, queríamos asegurarnos de que estábamos midiendo la capacidad de las herramientas bioinformáticas para hacer su trabajo y no lo bien que funcionaba el kit de captura del exoma TruSeq de Illumina. Es decir, sólo queremos tratar de llamar a las variantes que se supone que están presentes en los exones como se define por el kit de extracción. Por esta razón, utilizamos el archivo de cama proporcionado por Illumina, no un archivo de cama de anotación genérico, por ejemplo, RefSeq para hg19. Encontramos que para este individuo en particular, según la lista NIST-GiaB, debería haber un total de 34.886 SNVs y 1.473 indels dentro de las regiones definidas por el archivo de cama TruSeq.
Una vez que filtramos las variantes que no estaban localizadas en las regiones definidas por el archivo de cama del exoma TruSeq de Illumina, pasamos de cientos de miles de variantes putativas (datos no mostrados) a, en promedio, unas 23.000 variantes (SNVs e indels) por pipeline (Tabla 3). Este es un paso importante para los investigadores para empezar, ya que reduce significativamente el espacio de búsqueda de variantes potencialmente interesantes.
Aligner
Genotyper
Raw TP SNVs
Raw FP SNVs
Raw TP indels
Raw FP indels
Bowtie2
FreeBayes
23,985
73,473
806
2,482
Bowtie2
GATK HC
21,631
273
771
1,103
Bowtie2
GATK UG
25,136
2,276
418
420
Bowtie2
mpileup
21,930
1,030
734
1,414
Bowtie2
SNPSVM
17,613
47
–
–
BWA mem
FreeBayes
23,857
18.256
785
2.088
BWA mem
GATK HC
21,707
367
779
1,348
BWA mem
GATK UG
21,925
213
402
408
BWA mem
mpileup
25,081
2,129
761
1,772
BWA mem
SNPSVM
17,920
65
–
–
BWA sampe
FreeBayes
23,789
27,143
737
1,872
BWA sampe
GATK HC
21,878
263
758
1,161
BWA sampe
GATK UG
22,153
321
394
385
BWA sampe
mpileup
25,206
2,205
684
1,401
BWA sampe
SNPSVM
18,017
78
–
–
CUSHAW3
FreeBayes
23,191
53,525
624
3,310
CUSHAW3
GATK HC
19,673
14,814
751
4,727
CUSHAW3
GATK UG
19,113
13,184
360
1,005
CUSHAW3
mpileup
22,171
9,694
681
1,983
CUSHAW3
SNPSVM
–
–
MOSAIK
FreeBayes
23,373
39,203
783
3,359
MOSAIK
GATK HC
13,528
111
500
458
MOSAIK
GATK UG
17,147
76
392
284
MOSAIK
mpileup
–
–
–
–
MOSAIK
SNPSVM
14,586
8
–
–
Novoalign
FreeBayes
22,794
2,970
678
1,554
Novoalign
GATK HC
21,407
473
779
1,370
Novoalign
GATK UG
21,113
144
387
365
Novoalign
mpileup
24,512
1,861
773
1,781
Novoalign
SNPSVM
17,109
164
–
–
Tabla 3
Estadísticas de variantes en bruto para las 30 tuberías, incluyendo SNVs e indels.
3.2. Resultados de variantes crudas comparados con GiaB
Un aspecto que queríamos entender al hacer esta comparación era la importancia de filtrar las variantes detectadas por estas herramientas. La razón es que lo ideal es tener un nivel de sensibilidad lo más alto posible para que las mutaciones de interés no se pierdan en el proceso de filtrado. Por lo tanto, nos corresponde determinar si este paso es necesario o no y hasta qué punto lo es, ya que está claro, a partir de los resultados del NIST-GiaB y de la revisión de Bamshad et al., que la sensibilidad podría ser un problema.
Como podemos ver en la Tabla 3, el filtrado es más necesario para algunos llamadores de variantes que para otros cuando se trata de SNVs, y es absolutamente necesario para los indels. En la mayoría de los casos, el número de variantes TP es cercano o superior a nuestro número esperado de alrededor de 20.000 , pero, por otro lado, en algunos casos el número de FPs es muy alto.
Claramente hay mucha variación en los números generados por cada pipeline. Sin embargo, uno puede encontrar algunos puntos en común en los números que probablemente provienen de los orígenes algorítmicos de cada herramienta. FreeBayes produce tanto el mayor número de variantes sin filtrar como el mayor número de PFs. Es probable que sólo veamos este tipo de rendimiento en esta herramienta debido al hecho de que, aunque no es el único llamador de variantes basado en la inferencia bayesiana, es único en su interpretación de los alineamientos. Es decir, es un llamador basado en el haplotipo que identifica las variantes basándose en la secuencia de las propias lecturas en lugar de la alineación, que es como opera el UnifiedGenotyper de GATK.
Además, vemos que los alineadores basados en Burrows-Wheeler tienen un rendimiento muy similar entre sí: Bowtie2, BWA mem, y BWA sampe logran resultados similares en todos los ámbitos. Se podría suponer que esto se debe al hecho de que todas estas herramientas utilizan algoritmos similares al realizar su tarea designada. Esta observación se ve apoyada por el hecho de que MOSAIK (alineaciones con huecos utilizando el algoritmo Smith-Waterman) y CUSHAW3 (un enfoque híbrido de siembra) tienen algoritmos subyacentes muy diferentes y, por lo tanto, producen resultados muy diferentes.
Esta diferencia en los resultados que se correlacionan con diferentes algoritmos se ve mejor en los resultados de SNPSVM. De los llamadores de variantes, es el único que utiliza máquinas de vectores de apoyo y construcción de modelos para generar llamadas SNV. Parece que, aunque tiene la desventaja de no ser tan sensible como otros métodos, se beneficia de ser extremadamente preciso independientemente del alineador que se utilice. Esto sugiere que se puede omitir el paso de filtrado cuando se utiliza este llamador de variantes.
Con respecto a los indels, ningún alineador parece destacar entre los demás como uno que maneja bien este tipo de mutación. De hecho, cuando se observa el número de PFs, está claro que es el llamador de variantes el que juega el mayor papel en la precisión de la identificación de indels. Además, no hay datos para los pipelines CUSHAW3 más SNPSVM ni MOSAIK más SAMtools mpileup debido a que los archivos de alineación no contienen las cadenas CIGAR necesarias para que los llamadores de variantes funcionen aguas abajo. Por último, la razón por la que no hay datos de indel para SNPSVM es porque esta herramienta se utiliza únicamente para la identificación de SNVs.
3.3. Resultados filtrados en comparación con GiaVM Resultados filtrados comparados con GiaB
Como en la Tabla 4, las prácticas de filtrado estándar consiguen eliminar un gran número de SNVs FP para cada tubería; sin embargo, parece que estos filtros son un poco demasiado agresivos en la mayoría de los casos para la detección de SNVs, pero no lo suficientemente estrictos para los indels. Esto se hace evidente cuando se observan las diferencias en el número de FPs reportados en cada conjunto de datos. Por ejemplo, Bowtie2 con Freebayes ve una eliminación de 72.570 SNVs FP (una reducción del 98%) pero sólo una eliminación de 1.736 indels FP (una reducción del 70%). También hay que tener en cuenta que los filtros utilizados dependían de la tubería y, en su mayor parte, dentro de cada tubería produjeron reducciones similares en SNV y FP indel. La única excepción es el número de variantes identificadas a partir de los alineamientos CUSHAW3 en comparación con otros alineamientos: en general el número de SNVs TP es menor, el número de SNVs FP es mayor, y es el único alineador que produce más de 1.000 indels FP después del filtrado.
Alineador
Genotipador
SNVs TP filtrados
Filtrados SNVs
Indels TP filtrados
Indels FP filtrados
Bowtie2
FreeBayes
17,504
903
481
746
Bowtie2
GATK HC
17,330
29
648
687
Bowtie2
GATK UG
19,937
49
395
338
Bowtie2
mpileup
17,049
153
402
541
Bowtie2
SNPSVM
13,983
8
–
–
BWA mem
FreeBayes
17,376
347
461
739
BWA mem
GATK HC
19,388
302
689
860
BWA mem
GATK UG
20,000
48
397
355
BWA mem
mpileup
17,070
57
403
606
BWA mem
SNPSVM
15,060
10
–
–
BWA sampe
FreeBayes
17,435
450
443
647
BWA sampe
GATK HC
19,438
214
630
725
BWA sampe
GATK UG
19,557
27
384
336
BWA sampe
mpileup
17,049
111
387
518
BWA sampe
SNPSVM
15,218
10
–
CUSHAW3
FreeBayes
16,620
7,627
362
1,294
CUSHAW3
GATK HC
16,590
2,195
665
1,551
CUSHAW3
GATK UG
17,939
2.202
357
545
CUSHAW3
multiplicación
15,942
4,029
368
796
CUSHAW3
SNPSVM
–
–
–
–
MOSAIK
FreeBayes
17,177
679
458
645
MOSAIK
GATK HC
11,616
33
426
255
MOSAIK
GATK UG
16,423
42
381
224
MOSAIK
mpileup
–
–
–
–
MOSAIK
SNPSVM
4,727
3
–
–
Novoalign
FreeBayes
16,658
219
384
559
Novoalign
GATK HC
19,406
385
702
872
Novoalign
GATK UG
20,521
46
386
315
Novoalign
mpileup
16,493
62
396
579
Novoalign
SNPSVM
14,451
18
–
–
Tabla 4
Estadísticas de variantes filtradas para los 30 pipelines, incluyendo SNVs e indels.
Dado el hecho de que el filtrado reduce significativamente el número de variantes TP, podría ser prudente, con la excepción de los pipelines que utilizan CUSHAW3 y FreeBayes, omitir este paso cuando se buscan variantes raras y de alto impacto. En su lugar, se podría dedicar más tiempo a un proceso de filtrado que se base en la biología y no en la estadística. Por ejemplo, puede tener más sentido invertir tiempo en la identificación de una pequeña lista de variantes que probablemente sean de alto impacto: mutaciones en el sitio de empalme, indels que causan cambios de marco, mutaciones de truncamiento, mutaciones de pérdida de parada o mutaciones en genes que se sabe que son biológicamente relevantes para el fenotipo de interés. TPR y sensibilidad media
Como puede verse en la Tabla 5, el valor predictivo positivo (PPV) de cada herramienta, con la excepción de CUSHAW3, oscila entre el 91% y el 99,9% para las SNV, pero la sensibilidad media es muy baja (alrededor del 50%). Esta discrepancia podría deberse a varias razones, pero la más probable es la profundidad variable de la cobertura en los exones. Podemos ver que, además de la baja sensibilidad de los SNV, la sensibilidad de los indels es baja (alrededor del 30%); sin embargo, el VPP de los indels es considerablemente menor (35,86% a 52,95%). Esto podría deberse a cualquiera de las siguientes razones: los indels muy cortos son difíciles de detectar por NGS convencional , la representación de los indels por diferentes llamadores de variantes puede hacer que las herramientas afirmen incorrectamente que dos indels son diferentes, o las herramientas de alineación producen diferentes representaciones del mismo indel .
Herramienta
Promedio SNV PPV
Promedio SNV sensibilidad
Promedio de PPV de indel
Sensibilidad media de indel
Bowtie2
98.69%
49,19%
45,45%
32,69%
BWA mem
99.15%
50,96%
43,24%
33,10%
BWA sampe
99.09%
50,85%
45,31%
31,30%
CUSHAW3
80.69%
48,08%
29,50%
29,74%
MOSAIK
98.51%
35,79%
52,95%
28,63%
Novoalign
99,17%
50.18%
44,55%
31,70%
FreeBayes
90,95%
51,00%
35.86%
32,79%
GATK HaplotypeCaller
97,05%
51,03%
43,17%
31.79%
GATK UnifiedGenotyper
97,93%
50,77%
52,12%
31.57%
SAMtools mpileup
94,99%
50,76%
39,15%
31,30%
SNPSVM
99,92%
50.85%
N/A
N/A
Tabla 5
Valor predictivo positivo (VPP) medio y sensibilidad para cada herramienta.
Tal vez la explicación más probable para ambos tipos de mutaciones sea la cuestión de la profundidad. Como es el caso de cualquier estudio de análisis de variantes, un aumento de la profundidad de la cobertura conduce a un aumento de la sensibilidad, pero es imposible garantizar una buena profundidad de la cobertura debido a la incapacidad de los kits de captura del exoma para tirar uniformemente de los exones . Además, ningún kit de captura del exoma cubre todos los exones. De hecho, se ha demostrado que el análisis de variantes de la secuenciación del genoma completo con una profundidad media inferior a la de un exoma rinde más debido a la uniformidad de dicha profundidad. Por lo tanto, es probable que falte un gran número de variantes debido a que la lista NIST-GiaB se creó a partir de una compilación de datos de secuenciación exómica y genómica. En última instancia, para lograr una sensibilidad adecuada habrá que realizar una secuenciación del genoma completo, pero actualmente su coste es prohibitivo para la mayoría de los laboratorios. Afortunadamente, este coste sigue bajando, y pronto veremos un cambio gradual del análisis del exoma al análisis más completo del genoma completo.
3.5. La sensibilidad como función de la profundidad
Debido a que la sensibilidad refleja una de las métricas de rendimiento más importantes de una herramienta y a que la mayoría de las herramientas luchan por alcanzar una sensibilidad superior al 50%, nos gustaría seguir explorando cómo la profundidad afectaba a la sensibilidad de la llamada de variantes. Hemos analizado diferentes combinaciones de herramientas para determinar cuáles eran los mejores pipelines, llamadores de variantes y alineadores. Para la Figura 2, tomamos las cinco mejores combinaciones de llamadores de variantes y alineadores según su sensibilidad y tasa de falsos positivos (FPR). Es decir, seleccionamos los que tenían el mayor número de SNVs TP llamados además del menor número de SNVs FP. Tras la inspección, lo que salta a la vista es que la sensibilidad es menor de lo esperado. Todos los pipelines rinden más o menos al mismo nivel: identifican la mayoría de sus variantes cuando se alcanza una profundidad de aproximadamente 150x, lo que indica que esta profundidad es probablemente suficiente y que el número de variantes perdidas se debe probablemente a que ciertos exones tienen una cobertura inferior a la media. Obsérvese que cuatro de las cinco tuberías con mejor rendimiento tienen GATK UnifiedGenotyper como llamador de variantes, lo que demuestra su rendimiento superior independientemente del alineador utilizado, como se muestra en la Figura 3(b).
Figura 2
Sensibilidad en función de la profundidad para las cinco tuberías principales. Los cinco principales pipelines se muestran aquí con la profundidad de cada SNV trazada contra la sensibilidad.
(a)
(b)
(a)
(b)
Figura 3
Sensibilidad en función de la profundidad para el alineador superior y el llamador de variantes superior. (a) Resultados para la profundidad de cada SNV trazados contra la sensibilidad para el alineador superior, BWA mem, emparejado con cada llamador de variantes. (b) Resultados para la profundidad de cada SNV trazados contra la sensibilidad para el mejor llamador de variantes, GATK UnifiedGenotyper, emparejado con cada alineador.
Además de mirar los cinco mejores pipelines, determinamos que sería útil realizar el mismo análisis en el mejor alineador acoplado con cada llamador de variantes (Figura 3(a)), así como el mejor llamador de variantes acoplado con cada alineador (Figura 3(b)). Al igual que con los pipelines, el mejor alineador se identificó como aquel que produjo el mayor número de SNVs TP y el menor número de SNVs FP, en este caso BWA mem. A pesar de tener el mejor alineador para trabajar, todavía vemos una diferencia bastante grande entre los llamadores de variantes, lo que es probablemente atribuible a los diferentes algoritmos que emplean (Figura 3(a)). Sin embargo, en el caso del mejor llamador de variantes, GATK UnifiedGenotyper, parece haber menos variación entre los cuatro mejores alineadores, lo que indica que funciona bastante bien en la mayoría de las situaciones, con las excepciones de CUSAHW3 y MOSAIK.
3.6. Variantes compartidas entre los mejores alineadores
Por último, queríamos saber hasta qué punto los conjuntos de llamadas de variantes eran únicos entre los diferentes alineadores. Para ello, nos centramos de nuevo en las cinco principales tuberías de llamadas de variantes: Bowtie2 más UnifiedGenotyper, BWA mem más UnifiedGenotyper, BWA sampe más HaplotypeCaller, BWA sampe más UnifiedGenotyper y Novoalign más UnifiedGenotyper. Como se puede ver en la Figura 4, hay una gran cantidad de solapamiento entre los cinco pipelines diferentes en cuestión, con 15.489 SNVs (70%) compartidos de un total de 22.324 variantes distintas. Sin embargo, también se podría argumentar que esto se debe en gran medida a que cuatro de los cinco pipelines utilizan el UnifiedGenotyper como su llamador de variantes. Esta idea se ve corroborada por el hecho de que el mayor número de variantes exclusivas de un pipeline, 367, pertenece a la combinación BWA sampe más HaplotypeCaller. También vale la pena señalar que el segundo número más alto de SNVs únicos también pertenece al alineador BWA sampe, por lo que es posible que el alto número de SNVs únicos se atribuya mejor al alineador que al llamador de variantes.
Figura 4
La intersección de los SNVs identificados por los cinco mejores pipelines.
4. Conclusiones
Encontramos que entre los treinta pipelines diferentes probados Novoalign más GATK UnifiedGenotyper mostraron la mayor sensibilidad manteniendo un bajo número de FP para SNVs. De los alineadores, BWA mem fue el que mejor funcionó, pero los resultados variaron mucho dependiendo del llamador de variantes utilizado. Naturalmente, se deduce que el mejor llamador de variantes, GATK UnifiedGenotyper, produjo en su mayoría resultados similares independientemente del alineador utilizado. Sin embargo, es evidente que los indels siguen siendo difíciles de manejar para cualquier línea de producción, ya que ninguna de las líneas de producción logró una sensibilidad media superior al 33% o un VPP superior al 53%. Además del bajo rendimiento general que vemos en la detección de indels, la sensibilidad, independientemente del tipo de mutación, es un problema para todas las canalizaciones descritas en este documento. El número esperado de SNVs para el exoma de NA12878 es de 34.886, pero incluso cuando se utiliza la unión de todas las variantes identificadas por los cinco mejores pipelines, el mayor número identificado fue muy bajo (22.324). Parece que, aunque sigue siendo muy útil, el análisis del exoma tiene sus limitaciones incluso cuando se trata de algo tan aparentemente sencillo como la detección de SNV.
Divulgación
Adam Cornish es un estudiante graduado en el laboratorio de Chittibabu Guda con formación en informática y genómica. Chittibabu Guda (profesor asociado) tiene una formación interdisciplinar en biología molecular y computacional. Ha publicado una serie de métodos computacionales con una variedad de aplicaciones en la investigación biomédica, desde 2001.
Conflicto de intereses
Los autores no tienen conocimiento de ningún interés en competencia.
Contribución de los autores
Adam Cornish diseñó el estudio, realizó todos los análisis, hizo las figuras y escribió el artículo. Chittibabu Guda aportó comentarios esenciales para mejorar el artículo, así como información sobre los propios análisis, y editó minuciosamente el artículo. Todos los autores leyeron y aprobaron el documento final.
Agradecimientos
Este trabajo fue apoyado por los fondos de desarrollo a Chittibabu Guda del Centro Médico de la Universidad de Nebraska (UNMC). Los autores desean dar las gracias a los creadores de Novoalign por poner a disposición el software como versión de prueba y a las instalaciones del Núcleo de Bioinformática y Biología de Sistemas del UNMC por el apoyo a la infraestructura que facilitó esta investigación.
|
||||||||||||||||||||||||||||||||||||||||
Las seis herramientas utilizadas para generar alineaciones fueron Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK y Novoalign, y las cinco herramientas utilizadas para generar variantes fueron FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup y SNPSVM, como puede verse en la Tabla 2.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2.3. Filtrado
Los datos crudos se adquirieron del SRA (SRR098401), se dividieron con fastq-dump y se filtraron utilizando el kit de herramientas fastx. En concreto, fastq-dump utilizó los indicadores -split_files y -split_spot, y fastq_quality_filter se ejecutó con los siguientes argumentos: -Q 33 -q 30 -p 50. A continuación, las lecturas se emparejaron correctamente con un script personalizado.
2.4. Alineación
Los alineadores utilizaron los argumentos por defecto excepto cuando se utilizó un argumento de hilos cuando estaba disponible. Los comandos utilizados son los siguientes.
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. Cálculo de la profundidad de cobertura de los alineamientos
Para garantizar un cálculo adecuado de la profundidad de cobertura, se utilizó el módulo CalculateHsMetrics de Picard Tools con los siguientes argumentos:(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.bedEs importante tener en cuenta que el archivo TruSeq exome bed debe tener la cabecera del archivo de alineación SAM como anexo para que este módulo funcione. Además, la columna 6 debe trasladarse a la columna 4, y la columna 5 debe eliminarse del archivo de lecho de TruSeq.
2.6. Procesamiento de archivos de alineación
El procesamiento de los archivos de alineación (archivos SAM/BAM) requiere los siguientes pasos para todos los alineadores:(1)Conversión de SAM a BAM con SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Clasificación de archivos BAM con el módulo de 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=coordinada(3)Eliminación de duplicados de PCR mediante el módulo de herramientas Picard, MarkDuplicates:(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=alineaciones/NA12878.ALN.sorted.bam O=alineaciones/NA12878.ALN.dups_removed.bam REMOVE_DUPLICATES=true M=alignments/metrics(4)Grupo de lectura añadido a los archivos de alineación mediante el módulo de herramientas Picard, AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=coordinación RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Realineación en torno a los indeles utilizando los módulos GATK RealignerTargetCreator e 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)Recalibración de bases utilizando los módulos GATK BaseRecalibrator y PrintReads:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T BaseRecalibrator -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -sitiosconocidosgenomedbsnp_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. Llamada de variantes
Se utilizaron los argumentos por defecto para cada llamador de variantes a menos que contuviera una bandera «threads» o «parallel» en cuyo caso también se utilizó. Además, los indels se llamaron por separado de los SNV cuando fue posible. En concreto, los comandos utilizados son los siguientes:
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.vcfDebido a la inexistencia de las banderas CIGAR requeridas en el archivo de alineación, SNPSVM no pudo llamar a las variantes para CUSHAW3, y SAMtools mpileup no pudo llamar a las variantes en las alineaciones MOSAIK por la misma razón. Además, debido al hecho de que SNPSVM sólo detecta SNVs, no se reportaron indels para este programa.
2.8. Filtración de variantes
La filtración varió dependiendo del llamador de variantes que se utilizó. En los casos de GATK HaplotypeCaller y GATK UnifiedGenotyper, se utilizaron los módulos de GATK, VariantRecalibrator y ApplyRecalibration, para filtrar SNVs utilizando HapMap 3.3, el Omni 2.5 SNP BeadChip, y dbSNP 137 sin datos de 1000 Genome como conjuntos de entrenamiento. Para SNPSVM, se utilizaron puntuaciones QUAL ≥ 4 y valores DP ≥ 6. Para FreeBayes y SAMtools, se utilizaron puntuaciones QUAL ≥ 20 y valores DP ≥ 6.
2.9. Comparación de variantes
Para la comparación de variantes, se utilizó USeq 8.8.1 para comparar SNVs compartidos entre todos los conjuntos de datos. Para comparar indels, se utilizó la herramienta vcflib vcfintersect. El archivo TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, obtenido el 11 de diciembre de 2013, se utilizó para restringir las comparaciones a las ubicaciones que podían ser capturadas por el kit pull down del exoma utilizado en la secuenciación de SRR098401. Este archivo se puede obtener de Illumina aquí: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Para garantizar que las variantes estuvieran representadas de forma idéntica entre diferentes conjuntos de llamadas, se utilizó la herramienta vcflib vcfallelicprimitives para preprocesar los archivos vcf.
2.10. Cálculos estadísticos
Verdaderos positivos (TP). Es una mutación que fue detectada por el pipeline que se está probando y es una que existe en la lista NIST-GiaB.
Falso Positivo (FP). Es una mutación que fue detectada por la tubería que se está probando pero es una que no existe en la lista NIST-GiaB.
Verdadero negativo (TN). Es una mutación que no fue detectada por la tubería que se está probando y es una que no existe en la lista NIST-GiaB.
Falso negativo (FN). Es una mutación que no fue detectada por la tubería que se está probando pero que sí existe en la lista NIST-GiaB:
3. Resultados y Discusión
3.1. Prefiltrado de variantes
Cuando se realiza un análisis de variantes, uno de los muchos escollos que hay que tener en cuenta es el espacio de la secuencia del exoma (tal y como lo define el kit de captura del exoma) y cómo puede afectar a los resultados del análisis. En este caso, teníamos un único exoma (SRR098401) que fue extraído utilizando el kit de exoma TruSeq de Illumina y secuenciado en un HiSeq 2000. Con esto en mente, queríamos asegurarnos de que estábamos midiendo la capacidad de las herramientas bioinformáticas para hacer su trabajo y no lo bien que funcionaba el kit de captura del exoma TruSeq de Illumina. Es decir, sólo queremos tratar de llamar a las variantes que se supone que están presentes en los exones como se define por el kit de extracción. Por esta razón, utilizamos el archivo de cama proporcionado por Illumina, no un archivo de cama de anotación genérico, por ejemplo, RefSeq para hg19. Encontramos que para este individuo en particular, según la lista NIST-GiaB, debería haber un total de 34.886 SNVs y 1.473 indels dentro de las regiones definidas por el archivo de cama TruSeq.
Una vez que filtramos las variantes que no estaban localizadas en las regiones definidas por el archivo de cama del exoma TruSeq de Illumina, pasamos de cientos de miles de variantes putativas (datos no mostrados) a, en promedio, unas 23.000 variantes (SNVs e indels) por pipeline (Tabla 3). Este es un paso importante para los investigadores para empezar, ya que reduce significativamente el espacio de búsqueda de variantes potencialmente interesantes.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3.2. Resultados de variantes crudas comparados con GiaB
Un aspecto que queríamos entender al hacer esta comparación era la importancia de filtrar las variantes detectadas por estas herramientas. La razón es que lo ideal es tener un nivel de sensibilidad lo más alto posible para que las mutaciones de interés no se pierdan en el proceso de filtrado. Por lo tanto, nos corresponde determinar si este paso es necesario o no y hasta qué punto lo es, ya que está claro, a partir de los resultados del NIST-GiaB y de la revisión de Bamshad et al., que la sensibilidad podría ser un problema.
Como podemos ver en la Tabla 3, el filtrado es más necesario para algunos llamadores de variantes que para otros cuando se trata de SNVs, y es absolutamente necesario para los indels. En la mayoría de los casos, el número de variantes TP es cercano o superior a nuestro número esperado de alrededor de 20.000 , pero, por otro lado, en algunos casos el número de FPs es muy alto.
Claramente hay mucha variación en los números generados por cada pipeline. Sin embargo, uno puede encontrar algunos puntos en común en los números que probablemente provienen de los orígenes algorítmicos de cada herramienta. FreeBayes produce tanto el mayor número de variantes sin filtrar como el mayor número de PFs. Es probable que sólo veamos este tipo de rendimiento en esta herramienta debido al hecho de que, aunque no es el único llamador de variantes basado en la inferencia bayesiana, es único en su interpretación de los alineamientos. Es decir, es un llamador basado en el haplotipo que identifica las variantes basándose en la secuencia de las propias lecturas en lugar de la alineación, que es como opera el UnifiedGenotyper de GATK.
Además, vemos que los alineadores basados en Burrows-Wheeler tienen un rendimiento muy similar entre sí: Bowtie2, BWA mem, y BWA sampe logran resultados similares en todos los ámbitos. Se podría suponer que esto se debe al hecho de que todas estas herramientas utilizan algoritmos similares al realizar su tarea designada. Esta observación se ve apoyada por el hecho de que MOSAIK (alineaciones con huecos utilizando el algoritmo Smith-Waterman) y CUSHAW3 (un enfoque híbrido de siembra) tienen algoritmos subyacentes muy diferentes y, por lo tanto, producen resultados muy diferentes.
Esta diferencia en los resultados que se correlacionan con diferentes algoritmos se ve mejor en los resultados de SNPSVM. De los llamadores de variantes, es el único que utiliza máquinas de vectores de apoyo y construcción de modelos para generar llamadas SNV. Parece que, aunque tiene la desventaja de no ser tan sensible como otros métodos, se beneficia de ser extremadamente preciso independientemente del alineador que se utilice. Esto sugiere que se puede omitir el paso de filtrado cuando se utiliza este llamador de variantes.
Con respecto a los indels, ningún alineador parece destacar entre los demás como uno que maneja bien este tipo de mutación. De hecho, cuando se observa el número de PFs, está claro que es el llamador de variantes el que juega el mayor papel en la precisión de la identificación de indels. Además, no hay datos para los pipelines CUSHAW3 más SNPSVM ni MOSAIK más SAMtools mpileup debido a que los archivos de alineación no contienen las cadenas CIGAR necesarias para que los llamadores de variantes funcionen aguas abajo. Por último, la razón por la que no hay datos de indel para SNPSVM es porque esta herramienta se utiliza únicamente para la identificación de SNVs.
3.3. Resultados filtrados en comparación con GiaVM Resultados filtrados comparados con GiaB
Como en la Tabla 4, las prácticas de filtrado estándar consiguen eliminar un gran número de SNVs FP para cada tubería; sin embargo, parece que estos filtros son un poco demasiado agresivos en la mayoría de los casos para la detección de SNVs, pero no lo suficientemente estrictos para los indels. Esto se hace evidente cuando se observan las diferencias en el número de FPs reportados en cada conjunto de datos. Por ejemplo, Bowtie2 con Freebayes ve una eliminación de 72.570 SNVs FP (una reducción del 98%) pero sólo una eliminación de 1.736 indels FP (una reducción del 70%). También hay que tener en cuenta que los filtros utilizados dependían de la tubería y, en su mayor parte, dentro de cada tubería produjeron reducciones similares en SNV y FP indel. La única excepción es el número de variantes identificadas a partir de los alineamientos CUSHAW3 en comparación con otros alineamientos: en general el número de SNVs TP es menor, el número de SNVs FP es mayor, y es el único alineador que produce más de 1.000 indels FP después del filtrado.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Dado el hecho de que el filtrado reduce significativamente el número de variantes TP, podría ser prudente, con la excepción de los pipelines que utilizan CUSHAW3 y FreeBayes, omitir este paso cuando se buscan variantes raras y de alto impacto. En su lugar, se podría dedicar más tiempo a un proceso de filtrado que se base en la biología y no en la estadística. Por ejemplo, puede tener más sentido invertir tiempo en la identificación de una pequeña lista de variantes que probablemente sean de alto impacto: mutaciones en el sitio de empalme, indels que causan cambios de marco, mutaciones de truncamiento, mutaciones de pérdida de parada o mutaciones en genes que se sabe que son biológicamente relevantes para el fenotipo de interés. TPR y sensibilidad media
Como puede verse en la Tabla 5, el valor predictivo positivo (PPV) de cada herramienta, con la excepción de CUSHAW3, oscila entre el 91% y el 99,9% para las SNV, pero la sensibilidad media es muy baja (alrededor del 50%). Esta discrepancia podría deberse a varias razones, pero la más probable es la profundidad variable de la cobertura en los exones. Podemos ver que, además de la baja sensibilidad de los SNV, la sensibilidad de los indels es baja (alrededor del 30%); sin embargo, el VPP de los indels es considerablemente menor (35,86% a 52,95%). Esto podría deberse a cualquiera de las siguientes razones: los indels muy cortos son difíciles de detectar por NGS convencional , la representación de los indels por diferentes llamadores de variantes puede hacer que las herramientas afirmen incorrectamente que dos indels son diferentes, o las herramientas de alineación producen diferentes representaciones del mismo indel .
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Tal vez la explicación más probable para ambos tipos de mutaciones sea la cuestión de la profundidad. Como es el caso de cualquier estudio de análisis de variantes, un aumento de la profundidad de la cobertura conduce a un aumento de la sensibilidad, pero es imposible garantizar una buena profundidad de la cobertura debido a la incapacidad de los kits de captura del exoma para tirar uniformemente de los exones . Además, ningún kit de captura del exoma cubre todos los exones. De hecho, se ha demostrado que el análisis de variantes de la secuenciación del genoma completo con una profundidad media inferior a la de un exoma rinde más debido a la uniformidad de dicha profundidad. Por lo tanto, es probable que falte un gran número de variantes debido a que la lista NIST-GiaB se creó a partir de una compilación de datos de secuenciación exómica y genómica. En última instancia, para lograr una sensibilidad adecuada habrá que realizar una secuenciación del genoma completo, pero actualmente su coste es prohibitivo para la mayoría de los laboratorios. Afortunadamente, este coste sigue bajando, y pronto veremos un cambio gradual del análisis del exoma al análisis más completo del genoma completo.
3.5. La sensibilidad como función de la profundidad
Debido a que la sensibilidad refleja una de las métricas de rendimiento más importantes de una herramienta y a que la mayoría de las herramientas luchan por alcanzar una sensibilidad superior al 50%, nos gustaría seguir explorando cómo la profundidad afectaba a la sensibilidad de la llamada de variantes. Hemos analizado diferentes combinaciones de herramientas para determinar cuáles eran los mejores pipelines, llamadores de variantes y alineadores. Para la Figura 2, tomamos las cinco mejores combinaciones de llamadores de variantes y alineadores según su sensibilidad y tasa de falsos positivos (FPR). Es decir, seleccionamos los que tenían el mayor número de SNVs TP llamados además del menor número de SNVs FP. Tras la inspección, lo que salta a la vista es que la sensibilidad es menor de lo esperado. Todos los pipelines rinden más o menos al mismo nivel: identifican la mayoría de sus variantes cuando se alcanza una profundidad de aproximadamente 150x, lo que indica que esta profundidad es probablemente suficiente y que el número de variantes perdidas se debe probablemente a que ciertos exones tienen una cobertura inferior a la media. Obsérvese que cuatro de las cinco tuberías con mejor rendimiento tienen GATK UnifiedGenotyper como llamador de variantes, lo que demuestra su rendimiento superior independientemente del alineador utilizado, como se muestra en la Figura 3(b).
(a)
(b)
(a)
(b)
Además de mirar los cinco mejores pipelines, determinamos que sería útil realizar el mismo análisis en el mejor alineador acoplado con cada llamador de variantes (Figura 3(a)), así como el mejor llamador de variantes acoplado con cada alineador (Figura 3(b)). Al igual que con los pipelines, el mejor alineador se identificó como aquel que produjo el mayor número de SNVs TP y el menor número de SNVs FP, en este caso BWA mem. A pesar de tener el mejor alineador para trabajar, todavía vemos una diferencia bastante grande entre los llamadores de variantes, lo que es probablemente atribuible a los diferentes algoritmos que emplean (Figura 3(a)). Sin embargo, en el caso del mejor llamador de variantes, GATK UnifiedGenotyper, parece haber menos variación entre los cuatro mejores alineadores, lo que indica que funciona bastante bien en la mayoría de las situaciones, con las excepciones de CUSAHW3 y MOSAIK.
3.6. Variantes compartidas entre los mejores alineadores
Por último, queríamos saber hasta qué punto los conjuntos de llamadas de variantes eran únicos entre los diferentes alineadores. Para ello, nos centramos de nuevo en las cinco principales tuberías de llamadas de variantes: Bowtie2 más UnifiedGenotyper, BWA mem más UnifiedGenotyper, BWA sampe más HaplotypeCaller, BWA sampe más UnifiedGenotyper y Novoalign más UnifiedGenotyper. Como se puede ver en la Figura 4, hay una gran cantidad de solapamiento entre los cinco pipelines diferentes en cuestión, con 15.489 SNVs (70%) compartidos de un total de 22.324 variantes distintas. Sin embargo, también se podría argumentar que esto se debe en gran medida a que cuatro de los cinco pipelines utilizan el UnifiedGenotyper como su llamador de variantes. Esta idea se ve corroborada por el hecho de que el mayor número de variantes exclusivas de un pipeline, 367, pertenece a la combinación BWA sampe más HaplotypeCaller. También vale la pena señalar que el segundo número más alto de SNVs únicos también pertenece al alineador BWA sampe, por lo que es posible que el alto número de SNVs únicos se atribuya mejor al alineador que al llamador de variantes.
4. Conclusiones
Encontramos que entre los treinta pipelines diferentes probados Novoalign más GATK UnifiedGenotyper mostraron la mayor sensibilidad manteniendo un bajo número de FP para SNVs. De los alineadores, BWA mem fue el que mejor funcionó, pero los resultados variaron mucho dependiendo del llamador de variantes utilizado. Naturalmente, se deduce que el mejor llamador de variantes, GATK UnifiedGenotyper, produjo en su mayoría resultados similares independientemente del alineador utilizado. Sin embargo, es evidente que los indels siguen siendo difíciles de manejar para cualquier línea de producción, ya que ninguna de las líneas de producción logró una sensibilidad media superior al 33% o un VPP superior al 53%. Además del bajo rendimiento general que vemos en la detección de indels, la sensibilidad, independientemente del tipo de mutación, es un problema para todas las canalizaciones descritas en este documento. El número esperado de SNVs para el exoma de NA12878 es de 34.886, pero incluso cuando se utiliza la unión de todas las variantes identificadas por los cinco mejores pipelines, el mayor número identificado fue muy bajo (22.324). Parece que, aunque sigue siendo muy útil, el análisis del exoma tiene sus limitaciones incluso cuando se trata de algo tan aparentemente sencillo como la detección de SNV.
Divulgación
Adam Cornish es un estudiante graduado en el laboratorio de Chittibabu Guda con formación en informática y genómica. Chittibabu Guda (profesor asociado) tiene una formación interdisciplinar en biología molecular y computacional. Ha publicado una serie de métodos computacionales con una variedad de aplicaciones en la investigación biomédica, desde 2001.
Conflicto de intereses
Los autores no tienen conocimiento de ningún interés en competencia.
Contribución de los autores
Adam Cornish diseñó el estudio, realizó todos los análisis, hizo las figuras y escribió el artículo. Chittibabu Guda aportó comentarios esenciales para mejorar el artículo, así como información sobre los propios análisis, y editó minuciosamente el artículo. Todos los autores leyeron y aprobaron el documento final.
Agradecimientos
Este trabajo fue apoyado por los fondos de desarrollo a Chittibabu Guda del Centro Médico de la Universidad de Nebraska (UNMC). Los autores desean dar las gracias a los creadores de Novoalign por poner a disposición el software como versión de prueba y a las instalaciones del Núcleo de Bioinformática y Biología de Sistemas del UNMC por el apoyo a la infraestructura que facilitó esta investigación.