Volver a verificar la ubicación en el mapa

En las tres secciones anteriores, controlamos el sesgo de la biblioteca (MarkDuplicates) y el sesgo del secuenciador (BQSR). Ahora, las lecturas restantes se convierten en una observación honesta de "cómo difiere el genoma de este individuo de la referencia". Extraer las variantes reales a partir de esta observación es el objetivo de HaplotypeCaller en esta sección.

Es el núcleo del pipeline de germline. No nos apresuremos.

¿Por qué ensamblaje local?

Los métodos anteriores (antes de UnifiedGenotyper de GATK3) analizaban cada posición de forma independiente: "¿Cuántas lecturas informan de una base diferente a la referencia en esta coordenada?". Este enfoque capturaba bien los SNP, pero era muy débil cerca de los indels. La razón es la siguiente.

Cuando hay un indel, el alineamiento alrededor de este puede ocurrir de múltiples maneras, por lo que las lecturas se alinean de forma diferente en cada caso. Las observaciones por posición se dispersan, resultando en "hay mucho ruido desconocido en esta región" en lugar de "hay un SNP en esta coordenada".

La solución de HaplotypeCaller es la siguiente:

1. Recopilar todas las lecturas de una región de interés (de cientos de pb).
2. Crear candidatos de haplotipos posibles dentro de esa región mediante ensamblaje local.
3. Realinear cada lectura candidata contra cada haplotipo candidato. ¿Cuál candidato explica mejor el conjunto de lecturas?
4. Determinar el genotipo utilizando la regla de Bayes.

Este enfoque es dramáticamente más preciso cerca de los indels en comparación con el enfoque por posición.

Desglose en cuatro pasos para guiar el proceso

Paso 1: Detección de regiones activas

No podemos ensamblar todo el genoma, por lo que debemos seleccionar solo las regiones "donde hay algo diferente".

• Escanear puntos donde las lecturas se alinean de forma diferente a la referencia, con muchas supresiones suaves (soft clips) o con MAPQ bajo.
• Definir una ventana de ±cientos de pb centrada en estos puntos como región activa.
• Asumir que el resto de las regiones son idénticas a la referencia (ahorro computacional).

De los 3×10⁹ pb del genoma completo, menos del 1% corresponde realmente a regiones activas que deben investigarse.

Paso 2: Ensamblaje local — Grafo de De Bruijn

Construimos un grafo de De Bruijn pequeño con las lecturas de la región activa (ver M23). Del grafo, extraemos la ruta de referencia y las rutas alternativas.

• Ruta de referencia: La referencia tal cual.
• Ruta alternativa 1: Un SNP en una posición respecto a la referencia.
• Ruta alternativa 2: Un indel en la referencia.
• Ruta alternativa 3: Una estructura completamente diferente (candidato a variante estructural).

Estas rutas son los candidatos de haplotipos. Dado que HaplotypeCaller los genera mediante ensamblaje, las variantes complejas e indels se incluyen naturalmente como candidatos.

Paso 3: Cálculo de la alineación de cada lectura con cada haplotipo usando PairHMM

Para cada par (lectura r, haplotipo h), calculamos P(r|h). Este cálculo corresponde a PairHMM.

• Estado: coincidencia, inserción y deleción.
• Probabilidad de emisión: basada en Phred (después de BQSR) de las lecturas.
• Probabilidad de transición: probabilidades de apertura y extensión de huecos.
• El algoritmo Forward (ver M19) calcula la suma logarítmica de la probabilidad de alineamiento global.

Para cada lectura se obtiene un vector con P(r|h1), P(r|h2), ...

Paso 4: Determinación del genotipo mediante Bayes

Si la combinación de dos haplotipos (h_i, h_j) constituye el genotipo de un individuo, la verosimilitud del conjunto de lecturas observadas es

Se asume que cada lectura proviene con igual probabilidad de uno u otro haplotipo. Multiplicando esta verosimilitud por la probabilidad a priori (frecuencia alélica poblacional) se obtiene la probabilidad a posteriori, y el genotipo con argmax es la llamada final.

Esta es la estadística central de HaplotypeCaller.

GVCF — La puerta de entrada a la integración multimuestra

Si solo se realiza una llamada para una muestra, basta con generar un VCF. Pero ¿qué ocurre cuando se desea integrar el análisis de 30, 300 o 30.000 muestras?

Generar un VCF independiente para cada muestra y combinarlos después genera problemas. Las variantes llamadas solo en la muestra A no están registradas en las otras muestras, por lo que resulta imposible distinguir si se trata de datos faltantes (missing data) o de referencia.

GVCF (Genomic VCF) resuelve este problema.

• Registra todas las posiciones para cada muestra.
• En las posiciones con variantes, se incluye la información detallada.
• En las posiciones sin variantes, se registra en bloque que "este segmento es idéntico a la referencia".

Al integrar múltiples GVCFs posteriormente, la condición de referencia o variante en cada posición queda explícita. Tras la integración, se procede a la llamada multimuestra con GenotypeGVCFs (se aborda en S11).

Comando de ejecución

Estándar para WGS germinales.

gatk HaplotypeCaller \
    -I S1.recal.bam \
    -R hg38.fa \
    -O S1.g.vcf.gz \
    -ERC GVCF \
    --tmp-dir /data/tmp

• -ERC GVCF: Modo GVCF (preparación para la integración de múltiples muestras).
• -ERC BP_RESOLUTION: Detalles por posición (diagnóstico preciso, tamaño de archivo grande).
• Si solo se desea el VCF final de una sola muestra, utilizar -ERC NONE (predeterminado).

Conozcamos también algunas opciones adicionales.

• --emit-ref-confidence: Notación de confianza de referencia.
• --stand-call-conf 30: Confianza mínima de llamada (predeterminada).
• -L intervals.bed: Reducción del alcance.

Cinco consejos prácticos

1. WES requiere un archivo intervals.bed

En WES solo interesan los exones. Utilizar -L exome.bed para reducir el alcance disminuye el tiempo de ejecución en un factor de 20.

2. Paralelización de WGS

Dado que WGS es grande, dividir por cromosomas para ejecutar en paralelo y luego integrar los resultados.

for chr in chr{1..22} chrX chrY chrM; do
    gatk HaplotypeCaller ... -L $chr -O S1_$chr.g.vcf.gz &
done
wait
gatk MergeVcfs -I S1_chr1.g.vcf.gz -I S1_chr2.g.vcf.gz ... -O S1.g.vcf.gz

Procesamiento automático en Snakemake/Nextflow.

3. --pair-hmm-implementation

La lógica básica indica que GATK4 utiliza la implementación SIMD AVX-512. Si el contenedor no admite la CPU, se produce una degradación a cálculo escalar (~5 veces más lento). Verificar la CPU del servidor.

4. Tamaño de GVCF

Aunque GVCF es más grande que VCF, su tasa de compresión es alta, por lo que un WGS humano completo tiene aproximadamente entre 500 MB y 1 GB. Para 30 muestras, esto equivale a 15～30 GB.

5. Casos en los que HaplotypeCaller puede pasar por alto variantes

• Regiones con cobertura de lecturas extremadamente baja (por ejemplo, ≤5×).
• Cerca de secuencias repetitivas (repetición en tándem corta) (herramienta separada en S14 SV).
• Variantes estructurales (variantes grandes de ≥100 pb): el límite práctico de HaplotypeCaller para indels es de hasta ~40 pb.

En estos casos, se puede complementar con DeepVariant (S12) o herramientas de SV (S14).

Cálculo manual — Determinación del genotipo por Bayes

Supongamos que en la región activa se observaron 20 lecturas y hay dos haplotipos candidatos: h1 (referencia) y h2 (un SNP).

Para cada lectura r, supongamos que PairHMM produce los siguientes resultados (valores promedio):

• 12 lecturas: P(r|h2) ≫ P(r|h1) → favorece a h2.
• 8 lecturas: P(r|h1) ≈ P(r|h2) → decisión ambigua.

Se calculan las verosimilitudes para los tres genotipos candidatos (h1/h1, h1/h2, h2/h2):

• h1/h1 (homocigoto de referencia): La observación de 12 lecturas es desfavorable (log P disminuye significativamente).
• h1/h2 (heterocigoto): Bajo la suposición de proporción 50/50, todas las lecturas pueden explicarse adecuadamente.
• h2/h2 (homocigoto variante): La observación ambigua de 8 lecturas es desfavorable.

Según Bayes, h1/h2 es el argmax. Se llama a la variante como SNP heterocigoto.

Esta lógica corresponde exactamente al criterio de decisión de HaplotypeCaller.

Abrir VCF

bcftools view S1.g.vcf.gz | head -50

Campos principales:

• CHROM · POS · REF · ALT: posición y bases de referencia/alternativa.
• QUAL: confianza de la variante (Phred).
• INFO: estadísticas como DP (profundidad), MQ (MAPQ medio) y FS (sesgo de hebra).
• FORMAT · Sample: GT (genotipo), AD (profundidad de cada alelo), DP y GQ (confianza del genotipo).

Valores habituales de GT:

• 0/0: homocigoto de referencia.
• 0/1: heterocigoto.
• 1/1: homocigoto alternativo.
• ./.: dato ausente.

Mapeo CS

• Ensamblaje local de De Bruijn: reutiliza el algoritmo de M23 en cada región activa.
• PairHMM: aplica el algoritmo Forward derivado en M18～M19 de HMM a cada par lectura-haplotipo.
• Maximización de la probabilidad posterior bayesiana: previa (frecuencia alélica) + verosimilitud (PairHMM) → posterior.

Consulte el artículo de DryBench "HMM Forward en la práctica" (hmm-forward-in-practice).

Conclusión

HaplotypeCaller integra tres conceptos: ensamblaje local, PairHMM y el teorema de Bayes. En el próximo artículo (S10) veremos cómo este marco se amplía a la relación tumor-normal para convertirse en Mutect2. Examinaremos a nivel estadístico dónde se separan los análisis germinal y somático.

Para profundizar

• Poplin R. et al. (2018), Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. Ensamblaje local de HaplotypeCaller y escalabilidad de GVCF.
• Guía oficial de GATK sobre HaplotypeCaller (https://gatk.broadinstitute.org/hc/en-us/articles/360036721051): fuente original de las opciones.
• Broad Institute BroadE — HaplotypeCaller deep dive, conferencia en YouTube. Ofrece buenas visualizaciones de las regiones activas y PairHMM.
• Rausch T. et al. (2019), Alfred: interactive multi-sample BAM alignment statistics. Bioinformatics 35:2489. Herramienta para validar los resultados.
• Datos oficiales del Proyecto 1000 Genomes (https://www.internationalgenome.org/): Valor de referencia (ground truth).

Si logra comprender intuitivamente por qué el ensamblaje local detecta bien los indels, la detección de variantes somáticas (Somatic calling) del próximo episodio se entenderá de forma mucho más natural. En el próximo episodio, pasaremos a Mutect2.