Parte 1: Descripción del método¶
Traducción asistida por IA - más información y sugerencias
El llamado de variantes es un método de análisis genómico que tiene como objetivo identificar variaciones en una secuencia genómica en relación con un genoma de referencia. Aquí vamos a usar herramientas y métodos diseñados para llamar variantes germinales cortas, es decir SNPs e indels, en datos de secuenciación de genoma completo.

Un pipeline completo de llamado de variantes típicamente involucra muchos pasos, incluyendo el mapeo a la referencia (a veces referido como alineamiento del genoma) y el filtrado y priorización de variantes. Por simplicidad, en este curso nos vamos a enfocar solo en la parte de llamado de variantes.
Métodos¶
Vamos a mostrarte dos formas de aplicar el llamado de variantes a muestras de secuenciación de genoma completo para identificar SNPs e indels germinales. Primero comenzaremos con un enfoque simple por muestra que llama variantes independientemente de cada muestra. Luego te mostraremos un enfoque de llamado conjunto más sofisticado que analiza múltiples muestras juntas, produciendo resultados más precisos e informativos.
Antes de escribir cualquier código de workflow para cualquiera de los dos enfoques, vamos a probar los comandos manualmente en algunos datos de prueba.
Conjunto de datos¶
Proporcionamos los siguientes datos y recursos relacionados:
- Un genoma de referencia que consiste en una pequeña región del cromosoma humano 20 (de hg19/b37) y sus archivos accesorios (índice y diccionario de secuencia).
- Tres muestras de secuenciación de genoma completo correspondientes a un trío familiar (madre, padre e hijo), que han sido reducidas a una pequeña porción de datos en el cromosoma 20 para mantener los tamaños de archivo pequeños. Estos son datos de secuenciación Illumina de lecturas cortas que ya han sido mapeados al genoma de referencia, proporcionados en formato BAM (Binary Alignment Map, una versión comprimida de SAM, Sequence Alignment Map).
- Una lista de intervalos genómicos, es decir, coordenadas en el genoma donde nuestras muestras tienen datos adecuados para llamar variantes, proporcionada en formato BED.
Software¶
Las dos herramientas principales involucradas son Samtools, un conjunto de herramientas ampliamente utilizado para manipular archivos de alineamiento de secuencias, y GATK (Genome Analysis Toolkit), un conjunto de herramientas para el descubrimiento de variantes desarrollado en el Broad Institute.
Estas herramientas no están instaladas en el entorno de GitHub Codespaces, así que las usaremos a través de contenedores obtenidos mediante el servicio Seqera Containers (ver Hello Containers).
Consejo
Asegúrate de estar en el directorio nf4-science/genomics para que la última parte de la ruta mostrada cuando escribes pwd sea genomics.
1. Llamado de variantes por muestra¶
El llamado de variantes por muestra procesa cada muestra independientemente: el llamador de variantes examina los datos de secuenciación de una muestra a la vez e identifica posiciones donde la muestra difiere de la referencia.
En esta sección probamos los dos comandos que conforman el enfoque de llamado de variantes por muestra: indexar un archivo BAM con Samtools y llamar variantes con GATK HaplotypeCaller. Estos son los comandos que envolveremos en un workflow de Nextflow en la Parte 2 de este curso.
- Generar un archivo índice para un archivo de entrada BAM usando Samtools
- Ejecutar GATK HaplotypeCaller en el archivo BAM indexado para generar llamados de variantes por muestra en VCF (Variant Call Format)
Comenzamos probando los dos comandos en solo una muestra.
1.1. Indexar un archivo de entrada BAM con Samtools¶
Los archivos índice son una característica común de los formatos de archivo bioinformáticos; contienen información sobre la estructura del archivo principal que permite a herramientas como GATK acceder a un subconjunto de los datos sin tener que leer todo el archivo. Esto es importante debido a lo grandes que pueden llegar a ser estos archivos.
Los archivos BAM a menudo se proporcionan sin un índice, por lo que el primer paso en muchos workflows de análisis es generar uno usando samtools index.
Vamos a descargar un contenedor de Samtools, iniciarlo de forma interactiva y ejecutar el comando samtools index en uno de los archivos BAM.
1.1.1. Descargar el contenedor de Samtools¶
Ejecuta el comando docker pull para descargar la imagen del contenedor de Samtools:
Salida del comando
1.20--b5dfbd93de237464: Pulling from library/samtools
6360b3717211: Pull complete
2ec3f7ad9b3c: Pull complete
7716ca300600: Pull complete
4f4fb700ef54: Pull complete
8c61d418774c: Pull complete
03dae77ff45c: Pull complete
aab7f787139d: Pull complete
4f4fb700ef54: Pull complete
837d55536720: Pull complete
897362c12ca7: Pull complete
3893cbe24e91: Pull complete
d1b61e94977b: Pull complete
c72ff66fb90f: Pull complete
0e0388f29b6d: Pull complete
Digest: sha256:bbfc45b4f228975bde86cba95e303dd94ecf2fdacea5bfb2e2f34b0d7b141e41
Status: Downloaded newer image for community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
Si no has descargado esta imagen antes, puede tardar un minuto en completarse. Una vez que termine, tienes una copia local de la imagen del contenedor.
1.1.2. Iniciar el contenedor de Samtools de forma interactiva¶
Para ejecutar el contenedor de forma interactiva, usa docker run con las opciones -it.
La opción -v ./data:/data monta el directorio local data dentro del contenedor para que las herramientas puedan acceder a los archivos de entrada.
Notarás que tu prompt cambia a algo como (base) root@a1b2c3d4e5f6:/tmp#, indicando que ahora estás dentro del contenedor.
Verifica que puedes ver los archivos de datos de secuencia bajo /data/bam:
Con eso, estás listo para probar tu primer comando.
1.1.3. Ejecutar el comando de indexación¶
La documentación de Samtools nos da la línea de comando para ejecutar para indexar un archivo BAM.
Solo necesitamos proporcionar el archivo de entrada; la herramienta generará automáticamente un nombre para la salida agregando .bai al nombre del archivo de entrada.
Ejecuta el comando samtools index en un archivo de datos:
El comando no produce ninguna salida en la terminal, pero ahora deberías ver un archivo llamado reads_mother.bam.bai en el mismo directorio que el archivo de entrada BAM original.
Contenido del directorio
Eso completa la prueba del primer paso.
1.1.4. Salir del contenedor de Samtools¶
Para salir del contenedor, escribe exit.
Tu prompt ahora debería volver a lo que era antes de iniciar el contenedor.
1.2. Llamar variantes con GATK HaplotypeCaller¶
Queremos ejecutar el comando gatk HaplotypeCaller en el archivo BAM que acabamos de indexar.
1.2.1. Descargar el contenedor de GATK¶
Primero, ejecutemos el comando docker pull para descargar la imagen del contenedor de GATK:
Salida del comando
Algunas capas muestran Already exists porque son compartidas con la imagen del contenedor de Samtools que descargamos anteriormente.
4.5.0.0--730ee8817e436867: Pulling from library/gatk4
6360b3717211: Already exists
2ec3f7ad9b3c: Already exists
7716ca300600: Already exists
4f4fb700ef54: Already exists
8c61d418774c: Already exists
03dae77ff45c: Already exists
aab7f787139d: Already exists
4f4fb700ef54: Already exists
837d55536720: Already exists
897362c12ca7: Already exists
3893cbe24e91: Already exists
d1b61e94977b: Already exists
e5c558f54708: Pull complete
087cce32d294: Pull complete
Digest: sha256:e33413b9100f834fcc62fd5bc9edc1e881e820aafa606e09301eac2303d8724b
Status: Downloaded newer image for community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
Esto debería ser más rápido que la primera descarga porque las dos imágenes de contenedor comparten la mayoría de sus capas.
1.2.2. Iniciar el contenedor de GATK de forma interactiva¶
Inicia el contenedor de GATK de forma interactiva con el directorio de datos montado, tal como lo hicimos para Samtools.
Tu prompt cambia para indicar que ahora estás dentro del contenedor de GATK.
1.2.3. Ejecutar el comando de llamado de variantes¶
La documentación de GATK nos da la línea de comando para ejecutar para realizar el llamado de variantes en un archivo BAM.
Necesitamos proporcionar el archivo de entrada BAM (-I) así como el genoma de referencia (-R), un nombre para el archivo de salida (-O) y una lista de intervalos genómicos a analizar (-L).
Sin embargo, no necesitamos especificar la ruta al archivo índice; la herramienta lo buscará automáticamente en el mismo directorio, basándose en la convención establecida de nomenclatura y ubicación.
Lo mismo aplica para los archivos accesorios del genoma de referencia (archivos de índice y diccionario de secuencia, *.fai y *.dict).
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_mother.bam \
-O /data/vcf/reads_mother.vcf \
-L /data/ref/intervals.bed
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.vcf -L /data/ref/intervals.bed
00:27:50.687 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
00:27:50.854 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.858 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
00:27:50.858 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
00:27:50.858 INFO HaplotypeCaller - Executing as root@a1fe8ff42d07 on Linux v6.10.14-linuxkit amd64
00:27:50.858 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
00:27:50.859 INFO HaplotypeCaller - Start Date/Time: February 8, 2026 at 12:27:50 AM GMT
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.861 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
00:27:50.861 INFO HaplotypeCaller - Picard Version: 3.1.1
00:27:50.861 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
00:27:50.863 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
00:27:50.864 INFO HaplotypeCaller - Deflater: IntelDeflater
00:27:50.864 INFO HaplotypeCaller - Inflater: IntelInflater
00:27:50.864 INFO HaplotypeCaller - GCS max retries/reopens: 20
00:27:50.864 INFO HaplotypeCaller - Requester pays: disabled
00:27:50.865 INFO HaplotypeCaller - Initializing engine
00:27:50.991 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
00:27:51.016 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
00:27:51.029 INFO HaplotypeCaller - Done initializing engine
00:27:51.040 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
00:27:51.042 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
00:27:51.042 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
00:27:51.046 INFO HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
00:27:51.063 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
00:27:51.085 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
00:27:51.086 INFO IntelPairHmm - Available threads: 10
00:27:51.086 INFO IntelPairHmm - Requested threads: 4
00:27:51.086 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
00:27:51.128 INFO ProgressMeter - Starting traversal
00:27:51.136 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
00:27:51.882 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
00:27:52.969 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
00:27:52.971 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1145.7
00:27:52.971 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
00:27:52.976 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.003346916
00:27:52.976 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.045731709
00:27:52.977 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
00:27:52.981 INFO HaplotypeCaller - Shutting down engine
[February 8, 2026 at 12:27:52 AM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=203423744
La salida de registro es muy detallada, por lo que hemos resaltado las líneas más relevantes en el ejemplo anterior.
Los archivos de salida, reads_mother.vcf y su archivo índice, reads_mother.vcf.idx, se crean dentro de tu directorio de trabajo en el contenedor.
El archivo VCF contiene los llamados de variantes, como veremos en un minuto, y el archivo índice tiene la misma función que el archivo índice BAM, permitir que las herramientas busquen y recuperen subconjuntos de datos sin cargar todo el archivo.
Como VCF es un formato de texto y este es un archivo de prueba pequeño, puedes ejecutar cat reads_mother.vcf para abrirlo y ver su contenido.
Si te desplazas hacia el inicio del archivo, encontrarás un encabezado compuesto de muchas líneas de metadatos, seguido de una lista de llamados de variantes, uno por línea.
Contenidos del archivo (abreviado)
En la salida de ejemplo anterior, hemos resaltado la última línea del encabezado, que marca el inicio de los datos de llamados de variantes. Cada línea de datos describe una posible variante identificada en los datos de secuenciación de la muestra. Para orientación sobre cómo interpretar el formato VCF, consulta este artículo útil.
1.2.4. Mover los archivos de salida¶
Cualquier cosa que permanezca dentro del contenedor será inaccesible para trabajo futuro.
El archivo índice BAM se creó directamente en el directorio /data/bam en el sistema de archivos montado, pero no el archivo VCF y su índice, por lo que necesitamos mover esos dos manualmente.
Contenido del directorio
Una vez hecho esto, todos los archivos ahora son accesibles en tu sistema de archivos normal.
1.2.5. Salir del contenedor de GATK¶
Para salir del contenedor, escribe exit.
Tu prompt debería volver a la normalidad. Eso concluye la prueba de llamado de variantes por muestra.
¡Escríbelo como un workflow!
Siéntete libre de pasar directamente a la Parte 2 si deseas comenzar a implementar este análisis como un workflow de Nextflow. Solo necesitarás volver para completar la segunda ronda de pruebas antes de pasar a la Parte 3.
2. Llamado conjunto en una cohorte¶
El enfoque de llamado de variantes que acabamos de usar genera llamados de variantes por muestra. Eso está bien para observar variantes de cada muestra de forma aislada, pero produce información limitada. A menudo es más interesante observar cómo difieren los llamados de variantes entre múltiples muestras. GATK ofrece un método alternativo llamado llamado conjunto de variantes para este propósito.
El llamado conjunto de variantes implica generar un tipo especial de salida de variantes llamada GVCF (por Genomic VCF) para cada muestra, luego combinar los datos GVCF de todas las muestras y ejecutar un análisis estadístico de 'genotipado conjunto'.

Lo que es especial sobre el GVCF de una muestra es que contiene registros que resumen estadísticas de datos de secuencia sobre todas las posiciones en el área objetivo del genoma, no solo las posiciones donde el programa encontró evidencia de variación. Esto es crítico para el cálculo de genotipado conjunto (lectura adicional).
El GVCF es producido por GATK HaplotypeCaller, la misma herramienta que acabamos de probar, con un parámetro adicional (-ERC GVCF).
La combinación de los GVCFs se realiza con GATK GenomicsDBImport, que combina los llamados por muestra en un almacén de datos (análogo a una base de datos).
El análisis de 'genotipado conjunto' propiamente dicho se realiza entonces con GATK GenotypeGVCFs.
Aquí probamos los comandos necesarios para generar GVCFs y ejecutar el genotipado conjunto. Estos son los comandos que envolveremos en un workflow de Nextflow en la Parte 3 de este curso.
- Generar un archivo índice para cada archivo de entrada BAM usando Samtools
- Ejecutar GATK HaplotypeCaller en cada archivo de entrada BAM para generar un GVCF de llamados de variantes genómicas por muestra
- Recolectar todos los GVCFs y combinarlos en un almacén de datos GenomicsDB
- Ejecutar el genotipado conjunto en el almacén de datos GVCF combinado para producir un VCF a nivel de cohorte
Ahora necesitamos probar todos estos comandos, comenzando con la indexación de los tres archivos BAM.
2.1. Indexar archivos BAM para las tres muestras¶
En la primera sección anterior, solo indexamos un archivo BAM. Ahora necesitamos indexar las tres muestras para que GATK HaplotypeCaller pueda procesarlas.
2.1.1. Iniciar el contenedor de Samtools de forma interactiva¶
Ya descargamos la imagen del contenedor de Samtools, así que podemos iniciarlo directamente:
Tu prompt cambia para indicar que estás dentro del contenedor, con el directorio de datos montado como antes.
2.1.2. Ejecutar el comando de indexación en las tres muestras¶
Ejecuta el comando de indexación en cada uno de los tres archivos BAM:
samtools index /data/bam/reads_mother.bam
samtools index /data/bam/reads_father.bam
samtools index /data/bam/reads_son.bam
Contenido del directorio
Esto debería producir los archivos índice en el mismo directorio que los archivos BAM correspondientes.
2.1.3. Salir del contenedor de Samtools¶
Para salir del contenedor, escribe exit.
Tu prompt debería volver a la normalidad.
2.2. Generar GVCFs para las tres muestras¶
Para ejecutar el paso de genotipado conjunto, necesitamos GVCFs para las tres muestras.
2.2.1. Iniciar el contenedor de GATK de forma interactiva¶
Ya descargamos la imagen del contenedor de GATK anteriormente, así que podemos iniciarlo directamente:
Tu prompt cambia para indicar que estás dentro del contenedor de GATK.
2.2.2. Ejecutar el comando de llamado de variantes con la opción GVCF¶
Para producir un VCF genómico (GVCF), agregamos la opción -ERC GVCF al comando base, que activa el modo GVCF de HaplotypeCaller.
También cambiamos la extensión del archivo de salida de .vcf a .g.vcf.
Esto técnicamente no es un requisito, pero es una convención fuertemente recomendada.
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_mother.bam \
-O reads_mother.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.g.vcf -L /data/ref/intervals.bed -ERC GVCF
16:51:00.620 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
16:51:00.749 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.751 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
16:51:00.751 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
16:51:00.751 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
16:51:00.751 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
16:51:00.752 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 4:51:00 PM GMT
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
16:51:00.753 INFO HaplotypeCaller - Picard Version: 3.1.1
16:51:00.753 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
16:51:00.754 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
16:51:00.754 INFO HaplotypeCaller - Deflater: IntelDeflater
16:51:00.754 INFO HaplotypeCaller - Inflater: IntelInflater
16:51:00.754 INFO HaplotypeCaller - GCS max retries/reopens: 20
16:51:00.754 INFO HaplotypeCaller - Requester pays: disabled
16:51:00.755 INFO HaplotypeCaller - Initializing engine
16:51:00.893 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
16:51:00.905 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
16:51:00.910 INFO HaplotypeCaller - Done initializing engine
16:51:00.912 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
16:51:00.917 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
16:51:00.919 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
16:51:00.919 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
16:51:00.923 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
16:51:00.923 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
16:51:00.933 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
16:51:00.945 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
16:51:00.945 INFO IntelPairHmm - Available threads: 4
16:51:00.945 INFO IntelPairHmm - Requested threads: 4
16:51:00.945 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
16:51:00.984 INFO ProgressMeter - Starting traversal
16:51:00.985 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
16:51:01.452 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
16:51:02.358 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
16:51:02.359 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1529.5
16:51:02.359 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
16:51:02.361 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0022800000000000003
16:51:02.361 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.061637120000000004
16:51:02.361 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
16:51:02.362 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 4:51:02 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=257949696
Esto crea el archivo de salida GVCF reads_mother.g.vcf en el directorio de trabajo actual en el contenedor, así como su archivo índice, reads_mother.g.vcf.idx.
Si ejecutas head -200 reads_mother.g.vcf para ver las primeras 200 líneas del contenido del archivo, verás que es mucho más largo que el VCF equivalente que generamos en la primera sección, y la mayoría de las líneas se ven bastante diferentes de lo que vimos en el VCF.
Contenidos del archivo (abreviado)
| reads_mother.g.vcf | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 | |
Hemos resaltado una vez más la última línea del encabezado, así como los primeros tres llamados de variantes 'propiamente dichos' en el archivo.
Notarás que las líneas de llamados de variantes están intercaladas con muchas líneas no variantes, que representan regiones no variantes donde el llamador de variantes no encontró evidencia de variación. Como se mencionó brevemente arriba, esto es lo que es especial sobre el modo GVCF de llamado de variantes: el llamador de variantes captura algunas estadísticas que describen su nivel de confianza en la ausencia de variación. Esto hace posible distinguir entre dos cifras de casos muy diferentes: (1) hay datos de buena calidad que muestran que la muestra es homocigota-referencia, y (2) no hay suficientes datos buenos disponibles para hacer una determinación de cualquier manera.
En un GVCF como este, típicamente hay muchas de estas líneas no variantes, con un número menor de registros de variantes dispersos entre ellas.
2.2.3. Repetir el proceso en las otras dos muestras¶
Ahora generemos GVCFs para las dos muestras restantes ejecutando los comandos a continuación, uno tras otro.
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_father.bam \
-O reads_father.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_father.bam -O reads_father.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:28:30.677 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:28:30.801 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.803 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:28:30.804 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:28:30.804 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:28:30.804 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:28:30.804 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:28:30 PM GMT
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.805 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:28:30.805 INFO HaplotypeCaller - Picard Version: 3.1.1
17:28:30.805 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:28:30.806 INFO HaplotypeCaller - Deflater: IntelDeflater
17:28:30.807 INFO HaplotypeCaller - Inflater: IntelInflater
17:28:30.807 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:28:30.807 INFO HaplotypeCaller - Requester pays: disabled
17:28:30.807 INFO HaplotypeCaller - Initializing engine
17:28:30.933 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:28:30.946 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:28:30.951 INFO HaplotypeCaller - Done initializing engine
17:28:30.953 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:28:30.957 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:28:30.959 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:28:30.960 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:28:30.963 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:28:30.963 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:28:30.972 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:28:30.987 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:28:30.987 INFO IntelPairHmm - Available threads: 4
17:28:30.987 INFO IntelPairHmm - Requested threads: 4
17:28:30.987 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:28:31.034 INFO ProgressMeter - Starting traversal
17:28:31.034 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:28:31.570 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:28:32.865 INFO HaplotypeCaller - 9 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
9 total reads filtered out of 2064 reads processed
17:28:32.866 INFO ProgressMeter - 20_10037292_10066351:13338 0.0 38 1245.2
17:28:32.866 INFO ProgressMeter - Traversal complete. Processed 38 total regions in 0.0 minutes.
17:28:32.868 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0035923200000000004
17:28:32.868 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.10765202500000001
17:28:32.868 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.03 sec
17:28:32.869 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:28:32 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=299892736
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_son.bam \
-O reads_son.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_son.bam -O reads_son.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:30:10.017 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:30:10.156 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.159 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:30:10.159 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:30:10.159 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:30:10.159 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:30:10.159 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:30:09 PM GMT
17:30:10.159 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:30:10.160 INFO HaplotypeCaller - Picard Version: 3.1.1
17:30:10.161 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:30:10.161 INFO HaplotypeCaller - Deflater: IntelDeflater
17:30:10.162 INFO HaplotypeCaller - Inflater: IntelInflater
17:30:10.162 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:30:10.162 INFO HaplotypeCaller - Requester pays: disabled
17:30:10.162 INFO HaplotypeCaller - Initializing engine
17:30:10.277 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:30:10.290 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:30:10.296 INFO HaplotypeCaller - Done initializing engine
17:30:10.298 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:30:10.302 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:30:10.303 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:30:10.304 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:30:10.307 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:30:10.307 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:30:10.315 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:30:10.328 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:30:10.329 INFO IntelPairHmm - Available threads: 4
17:30:10.329 INFO IntelPairHmm - Requested threads: 4
17:30:10.329 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:30:10.368 INFO ProgressMeter - Starting traversal
17:30:10.369 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:30:10.875 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:30:11.980 INFO HaplotypeCaller - 14 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
14 total reads filtered out of 1981 reads processed
17:30:11.981 INFO ProgressMeter - 20_10037292_10066351:13223 0.0 35 1302.7
17:30:11.981 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
17:30:11.983 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0034843710000000004
17:30:11.983 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.048108363
17:30:11.983 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
17:30:11.984 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:30:11 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=226492416
Una vez que esto se complete, deberías tener tres archivos que terminan en .g.vcf en tu directorio actual (uno por muestra) y sus respectivos archivos índice que terminan en .g.vcf.idx.
Contenido del directorio
En este punto, hemos llamado variantes en modo GVCF para cada una de nuestras muestras de entrada. Es hora de pasar al llamado conjunto.
¡Pero no salgas del contenedor! Vamos a usar el mismo en el siguiente paso.
2.3. Ejecutar el genotipado conjunto¶
Ahora que tenemos todos los GVCFs, podemos probar el enfoque de genotipado conjunto para generar llamados de variantes para una cohorte de muestras. Es un método de dos pasos que consiste en combinar los datos de todos los GVCFs en un almacén de datos, luego ejecutar el análisis de genotipado conjunto propiamente dicho para generar el VCF final de variantes llamadas conjuntamente.
2.3.1. Combinar todos los GVCFs por muestra¶
Este primer paso usa otra herramienta de GATK, llamada GenomicsDBImport, para combinar los datos de todos los GVCFs en un almacén de datos GenomicsDB. El almacén de datos GenomicsDB es una especie de formato de base de datos que sirve como almacenamiento intermedio para la información de variantes.
gatk GenomicsDBImport \
-V reads_mother.g.vcf \
-V reads_father.g.vcf \
-V reads_son.g.vcf \
-L /data/ref/intervals.bed \
--genomicsdb-workspace-path family_trio_gdb
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar GenomicsDBImport -V reads_mother.g.vcf -V reads_father.g.vcf -V reads_son.g.vcf -L /data/ref/intervals.bed --genomicsdb-workspace-path family_trio_gdb
17:37:07.569 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:37:07.699 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.702 INFO GenomicsDBImport - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:37:07.702 INFO GenomicsDBImport - For support and documentation go to https://software.broadinstitute.org/gatk/
17:37:07.703 INFO GenomicsDBImport - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:37:07.703 INFO GenomicsDBImport - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:37:07.704 INFO GenomicsDBImport - Start Date/Time: February 11, 2026 at 5:37:07 PM GMT
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.706 INFO GenomicsDBImport - HTSJDK Version: 4.1.0
17:37:07.706 INFO GenomicsDBImport - Picard Version: 3.1.1
17:37:07.707 INFO GenomicsDBImport - Built for Spark Version: 3.5.0
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:37:07.710 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:37:07.710 INFO GenomicsDBImport - Deflater: IntelDeflater
17:37:07.711 INFO GenomicsDBImport - Inflater: IntelInflater
17:37:07.711 INFO GenomicsDBImport - GCS max retries/reopens: 20
17:37:07.711 INFO GenomicsDBImport - Requester pays: disabled
17:37:07.712 INFO GenomicsDBImport - Initializing engine
17:37:07.883 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:37:07.886 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:37:07.889 INFO GenomicsDBImport - Done initializing engine
17:37:08.560 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:37:08.561 INFO GenomicsDBImport - Vid Map JSON file will be written to /tmp/family_trio_gdb/vidmap.json
17:37:08.561 INFO GenomicsDBImport - Callset Map JSON file will be written to /tmp/family_trio_gdb/callset.json
17:37:08.561 INFO GenomicsDBImport - Complete VCF Header will be written to /tmp/family_trio_gdb/vcfheader.vcf
17:37:08.561 INFO GenomicsDBImport - Importing to workspace - /tmp/family_trio_gdb
17:37:08.878 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.359 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.487 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.591 INFO GenomicsDBImport - Done importing batch 1/1
17:37:09.592 INFO GenomicsDBImport - Import completed!
17:37:09.592 INFO GenomicsDBImport - Shutting down engine
[February 11, 2026 at 5:37:09 PM GMT] org.broadinstitute.hellbender.tools.genomicsdb.GenomicsDBImport done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=113246208
Tool returned:
true
La salida de este paso es efectivamente un directorio que contiene un conjunto de directorios anidados adicionales que contienen los datos de variantes combinados en forma de múltiples archivos diferentes. Puedes explorarlo pero rápidamente verás que este formato de almacén de datos no está diseñado para ser leído directamente por humanos.
Consejo
GATK incluye herramientas que hacen posible inspeccionar y extraer datos de llamados de variantes del almacén de datos según sea necesario.
2.3.2. Ejecutar el análisis de genotipado conjunto propiamente dicho¶
Este segundo paso usa otra herramienta de GATK, llamada GenotypeGVCFs, para recalcular las estadísticas de variantes y los genotipos individuales a la luz de los datos disponibles en todas las muestras de la cohorte.
Salida del comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar GenotypeGVCFs -R /data/ref/ref.fasta -V gendb://family_trio_gdb -O family_trio.vcf
17:38:45.084 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:38:45.217 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.220 INFO GenotypeGVCFs - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:38:45.220 INFO GenotypeGVCFs - For support and documentation go to https://software.broadinstitute.org/gatk/
17:38:45.220 INFO GenotypeGVCFs - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:38:45.220 INFO GenotypeGVCFs - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:38:45.221 INFO GenotypeGVCFs - Start Date/Time: February 11, 2026 at 5:38:45 PM GMT
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - HTSJDK Version: 4.1.0
17:38:45.222 INFO GenotypeGVCFs - Picard Version: 3.1.1
17:38:45.222 INFO GenotypeGVCFs - Built for Spark Version: 3.5.0
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:38:45.223 INFO GenotypeGVCFs - Deflater: IntelDeflater
17:38:45.223 INFO GenotypeGVCFs - Inflater: IntelInflater
17:38:45.223 INFO GenotypeGVCFs - GCS max retries/reopens: 20
17:38:45.223 INFO GenotypeGVCFs - Requester pays: disabled
17:38:45.223 INFO GenotypeGVCFs - Initializing engine
17:38:45.544 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field InbreedingCoeff - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAC - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAF - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.577 INFO GenotypeGVCFs - Done initializing engine
17:38:45.615 INFO ProgressMeter - Starting traversal
17:38:45.615 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:38:45.903 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
GENOMICSDB_TIMER,GenomicsDB iterator next() timer,Wall-clock time(s),0.07757032800000006,Cpu time(s),0.07253379200000037
17:38:46.421 INFO ProgressMeter - 20_10037292_10066351:13953 0.0 3390 252357.3
17:38:46.422 INFO ProgressMeter - Traversal complete. Processed 3390 total variants in 0.0 minutes.
17:38:46.423 INFO GenotypeGVCFs - Shutting down engine
[February 11, 2026 at 5:38:46 PM GMT] org.broadinstitute.hellbender.tools.walkers.GenotypeGVCFs done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=203423744
Esto crea el archivo de salida VCF family_trio.vcf en el directorio de trabajo actual en el contenedor, así como su índice, family_trio.vcf.idx.
Es otro archivo razonablemente pequeño, así que puedes ejecutar cat family_trio.vcf para ver el contenido del archivo, y desplazarte hacia abajo para encontrar las primeras líneas de variantes.
Contenidos del archivo (abreviado)
Hemos resaltado una vez más la última línea del encabezado, que marca el inicio de los datos de llamados de variantes.
Esto se ve similar al VCF que generamos anteriormente, excepto que esta vez tenemos información a nivel de genotipo para las tres muestras. Las últimas tres columnas en el archivo son los bloques de genotipo para las muestras, listados en orden alfabético de su campo ID, como se muestra en la línea de encabezado resaltada.
Si observamos los genotipos llamados para nuestro trío familiar de prueba para la primera variante, vemos que el padre es heterocigoto-variante (0/1), y la madre y el hijo son ambos homocigotos-variante (1/1).
¡Esa es en última instancia la información que estamos buscando extraer del conjunto de datos!
2.3.3. Mover los archivos de salida¶
Como se mencionó anteriormente, cualquier cosa que permanezca dentro del contenedor será inaccesible para trabajo futuro. Antes de salir del contenedor, vamos a mover los archivos GVCF, el VCF final multi-muestra y todos sus archivos índice manualmente al sistema de archivos fuera del contenedor. De esa manera, tendremos algo con qué comparar cuando construyamos nuestro workflow para automatizar todo este trabajo.
Contenidos del directorio" hl_lines="14-19 22-23
data
├── bam
│ ├── reads_father.bam
│ ├── reads_father.bam.bai
│ ├── reads_mother.bam
│ ├── reads_mother.bam.bai
│ ├── reads_son.bam
│ └── reads_son.bam.bai
├── ref
│ ├── intervals.bed
│ ├── ref.dict
│ ├── ref.fasta
│ └── ref.fasta.fai
├── samplesheet.csv
└── vcf
├── family_trio.vcf
├── family_trio.vcf.idx
├── reads_father.g.vcf
├── reads_father.g.vcf.idx
├── reads_mother.g.vcf
├── reads_mother.g.vcf.idx
├── reads_mother.vcf
├── reads_mother.vcf.idx
├── reads_son.g.vcf
└── reads_son.g.vcf.idx
Una vez hecho esto, todos los archivos ahora son accesibles en tu sistema de archivos normal.
2.3.4. Salir del contenedor de GATK¶
Para salir del contenedor, escribe exit.
Tu prompt debería volver a la normalidad. Eso concluye las pruebas manuales de los comandos de llamado conjunto de variantes.
Conclusión¶
Sabes cómo probar los comandos de indexación de Samtools y llamado de variantes de GATK en sus respectivos contenedores, incluyendo cómo generar GVCFs y ejecutar el genotipado conjunto en múltiples muestras.
¿Qué sigue?¶
Toma un descanso, luego dirígete a la Parte 2 para aprender cómo envolver esos mismos comandos en workflows que usan contenedores para ejecutar el trabajo.