¿Por qué GATK4 Best Practices?

La respuesta a la pregunta "¿cuál es el pipeline estándar?" en el análisis del genoma humano ha sido siempre la misma durante los últimos 12 años: GATK Best Practices del Broad Institute. Los laboratorios de diagnóstico clínico, los consorcios a gran escala (gnomAD, All of Us, UK Biobank) y la mayoría de las reproducciones de artículos científicos utilizan este pipeline tal cual. Cada vez que aparece una nueva herramienta, se realiza un benchmarking "en comparación con GATK".

En esta sección, trazamos el panorama general de GATK4 Best Practices en forma de un mapa único. En las siguientes cinco secciones (S07～S11), profundizaremos paso a paso en las opciones y la derivación estadística de cada etapa. Por ahora, veamos primero el mapa.

Germline vs Somatic — El pipeline se divide en dos caminos

El primer concepto que debemos captar es que el pipeline se divide en dos propósitos.

• Germline (germinales): Genotipos de las células normales. Variantes heredadas de los padres. Un individuo = un genotipo.
• Somatic (somáticas): Variantes que aparecen, por ejemplo en el cáncer, en células alteradas respecto de las células normales. El tejido tumoral es una mezcla de múltiples subclones.

Los pipelines para ambos propósitos son completamente idénticos en la parte inicial (alineación y control de calidad), pero se dividen en la etapa de llamada de variantes. Germline utiliza HaplotypeCaller, mientras que Somatic utiliza Mutect2. Esto se debe a que sus supuestos estadísticos son diferentes.

Mapa de 6 etapas

Las 6 etapas estándar del pipeline Germline GATK4 Best Practices son las siguientes:

[1] Raw FASTQ
     │
[2] Alineamiento BWA-MEM  ────────  etiquetado de grupos de reads
     │
[3] MarkDuplicates (S07)  ─  marcado de duplicados de PCR
     │
[4] BQSR (S08)  ────────────  recalibración de la calidad de bases
     │
[5] HaplotypeCaller (S09) ─  detección de variantes candidatas (GVCF por muestra)
     │
[6] Joint Genotyping + VQSR (S11) ─ integración de múltiples muestras + filtrado
     │
     ↓
   Final VCF

Veamos por qué es necesario cada paso, una oración a la vez.

[3] MarkDuplicates — Control del sesgo de PCR

Durante la preparación de la librería, la etapa de amplificación por PCR replica algunas lecturas 100 veces y otras solo una vez. Ignorar este sesgo provocaría una representación excesiva en la cobertura de las posiciones de las lecturas más replicadas. MarkDuplicates marca como duplicados los pares de lecturas que comparten la misma coordenada inicial y el mismo CIGAR (estableciendo la bandera 1024). Al no eliminarlos físicamente, sino solo marcarlos, pueden ser ignorados estadísticamente en etapas posteriores.

[4] BQSR — Recalibración posterior de los puntajes Phred

Los puntajes Phred reportados por el secuenciador no son perfectos; presentan sesgos sistemáticos dependiendo del contexto de dinucleótidos o la posición del ciclo. BQSR (Base Quality Score Recalibration) utiliza sitios de SNP conocidos (dbSNP) como valores verdaderos para aprender la tasa de error real en cada posición de la lectura y recalcular los puntajes Phred. Este recálculo satisface las suposiciones de HaplotypeCaller, reduciendo significativamente los falsos positivos.

[5] HaplotypeCaller — Ensamblaje local + llamada bayesiana

Los llamadores de variantes anteriores (como UnifiedGenotyper en GATK3) analizaban cada posición de forma independiente. HaplotypeCaller extrae una ventana pequeña (por ejemplo, 300 pb), realiza un ensamblaje local con las lecturas de ese intervalo y determina el genotipo basándose en qué tan bien se ajustan las lecturas a los haplotipos candidatos. Este enfoque es dramáticamente más preciso cerca de los indels. Genera GVCFs por muestra individual, lo que permite su integración posterior en múltiples muestras.

[6] Genotipado conjunto + VQSR — Integración y filtrado de múltiples muestras

• Genotipado conjunto: Integra los GVCFs de varias muestras → reevalúa el genotipo de cada muestra utilizando las frecuencias alélicas a nivel poblacional.
• VQSR (Variant Quality Score Recalibration): Entrena un GaussianMixtureModel tomando conjuntos conocidos de SNP e indels como valores verdaderos, asignando una confianza probabilística a los nuevos candidatos de variantes.

VQSR requiere al menos 30 muestras para ser estable. Para una sola muestra o tamaños pequeños (<30), se sustituye por Hard Filtering (filtrado con umbrales simples). Analizamos esta decisión en detalle en S11.

Rama somática — Mutect2 + PoN

Si desea extraer únicamente las variantes somáticas del tejido tumoral, reemplace la rama germinativa de la siguiente manera.

[3] MarkDuplicates
     │
[4] BQSR
     │
[5-somatic] Mutect2 (S10)  ─  par tumor-normal o solo tumor (requiere PoN)
     │
[6-somatic] FilterMutectCalls  ─  filtrar solo variantes somáticas

Dos diferencias clave.

• Detección de frecuencia alélica (VAF) baja: Se deben identificar variantes presentes en menos del 5% de las células dentro del tejido tumoral.
• Panel of Normals (PoN): Cuando no hay un control normal disponible, se utiliza la combinación de ruido de secuenciación de múltiples muestras normales como base para filtrar los falsos positivos.

Si captas intuitivamente dónde divergen las vías Germline y Somatic, te resultará mucho más fácil abordar los capítulos S09/S10 posteriormente.

¿Por qué este orden? — La cadena de supuestos estadísticos

Capturad intuitivamente por qué cada paso de las Buenas Prácticas sigue ese orden específico.

1. Sin alineación, no hay posición de referencia y, por tanto, la variante no puede definirse.
2. Sin MarkDuplicates, el sesgo de PCR distorsiona las estadísticas de cobertura.
3. Sin BQSR, HaplotypeCaller viola los supuestos sobre la probabilidad de error.
4. Sin HaplotypeCaller, no hay variantes candidatas.
5. Sin VQSR/Filtros duros, los falsos positivos permanecen sin filtrar.

Estas cinco condiciones forman una cadena lógica. Cambiar el orden rompe los supuestos estadísticos de las etapas posteriores. En los próximos cinco capítulos profundizaremos en esta cadena a nivel estadístico.

Unificación del entorno GATK con Docker/Singularity

La configuración del entorno para GATK4 suele fallar frecuentemente debido a la interdependencia de Java, Python, R y otras herramientas. Broad Institute distribuye imágenes oficiales de Docker.

docker pull broadinstitute/gatk:latest
docker run --rm -v $(pwd):/data broadinstitute/gatk:latest \
    gatk HaplotypeCaller -R /data/hg38.fa -I /data/S1.bam -O /data/S1.vcf.gz

O bien, en un HPC con Singularity/Apptainer:

apptainer pull docker://broadinstitute/gatk:latest
apptainer exec gatk_latest.sif gatk HaplotypeCaller -R hg38.fa -I S1.bam -O S1.vcf.gz

Estandarizar el trabajo en un solo contenedor mejora considerablemente la reproducibilidad del equipo. En los siguientes capítulos, los ejemplos de CLI mostrarán únicamente los comandos ejecutados dentro del contenedor, no el comando del contenedor.

Archivos de recursos: referencias y sitios conocidos

GATK requiere varios archivos auxiliares. Conviene descargarlos previamente del GATK Resource Bundle al configurar la canalización por primera vez.

• FASTA de referencia + .fai + .dict
• VCF de dbSNP (para BQSR)
• VCF de indels conocidos (Mills, indels de 1000G)
• HapMap · Omni · 1000G Phase1 SNP (VQSR training)

Todos están disponibles públicamente en la página de GATK Resource Bundle. Para hg38, el total ronda los 100 GB.

Tres consejos prácticos

1. Asigne un partición grande a BP083P0002

Las herramientas de GATK crean numerosos archivos temporales. Si /tmp es pequeño, la ejecución puede fallar. Especifique un volumen con espacio suficiente, por ejemplo --tmp-dir /data/gatk_tmp.

2. Reduzca el alcance con un archivo de intervalo

Para WES, limite el análisis a los intervalos del exoma; para genes concretos, solo a esas regiones. Reducir el ámbito con -L intervals.bed disminuye drásticamente el tiempo de ejecución.

3. Adquiera el hábito de registrar

gatk HaplotypeCaller ... 2>&1 | tee S1_hc.log

El 90% de la depuración posterior a la ejecución del pipeline consiste en revisar los registros (logs). Ejecutar sin logs convierte el rastreo de causas de fallo en un infierno.

Presets de GATK en Galaxy EU

Si la línea de comandos (CLI) resulta abrumadora, puede ejecutar directamente el flujo de trabajo de GATK4 en Galaxy EU.

1. Busque Whole Genome Sequencing (WGS) - Germline en Shared Data → Workflows.
2. Impórtelo y ejecútelo desde el History.
3. Entrada: FASTQ de extremo emparejado (paired-end), referencia hg38.

Galaxy ejecuta internamente exactamente estos 6 pasos en orden secuencial. Al tener las opciones expuestas en la interfaz de usuario, es ideal para familiarizarse antes de aprender la CLI.

Mapeo CS

• Tubería multicapa: La salida de cada paso se convierte en la entrada del siguiente. Es una extensión práctica de la filosofía de las tuberías de Unix.
• Precondicionamiento estadístico: MarkDuplicates y BQSR son preprocesamientos necesarios para satisfacer los supuestos estadísticos de las etapas posteriores.
• Política de reproducibilidad: Bloqueo de versiones de la imagen Docker y del Resource Bundle garantiza la reproducibilidad de los resultados.

Conclusión

Ya tenemos el mapa. En las próximas 5 entregas, profundizaremos en cada uno de estos 6 pasos.

• S07: El algoritmo de detección de duplicados de PCR en MarkDuplicates.
• S08: La deducción estadística de la recalibración Phred en BQSR.
• S09: El ensamblaje local y la llamada bayesiana en HaplotypeCaller.
• S10: Las estadísticas somáticas de Mutect2 + PoN.
• S11: Las reglas de decisión entre VQSR y el filtrado estricto (Hard Filtering).

Manipular las opciones del pipeline sin un mapa es perderse. Memorice la figura de estos 6 pasos de esta entrega y pasemos a la siguiente.

Para profundizar más

• Van der Auwera G. & O'Connor B. (2020), Genomics in the Cloud: Using Docker, GATK, and WDL in Terra. O'Reilly. El manual práctico oficial de los autores de GATK4.
• DePristo M. et al. (2011), A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43:491. El artículo original de GATK.
• Broad Institute BroadE YouTube — GATK Best Practices Overview Conferencia de 30 minutos. Con subtítulos completos. Excelente para comprender el contexto detrás de las decisiones de cada paso.
• GATK Official Walkthrough (https://gatk.broadinstitute.org/hc/en-us/articles/360035531192): Diagrama de flujo paso a paso.
• Sanger Genome Structure Analysis YouTube Channel: Demo práctica de GATK4.

Los próximos 5 episodios son el corazón del análisis genómico práctico. Si crea una cuenta en Galaxy EU y ejecuta los 6 pasos con un conjunto de datos pequeño antes de pasar al siguiente episodio, cada etapa será mucho más clara.