¿Por qué el filtrado es una etapa separada?

En las dos entregas anteriores, llamamos variantes candidatas con HaplotypeCaller y Mutect2. Estas llamadas en bruto contienen una cantidad considerable de falsos positivos. El llamador captura tantas variantes como puede, pero pospone la decisión sobre cuáles son reales a la siguiente etapa.

Esta entrega corresponde a esa siguiente etapa. Las Mejores Prácticas de Broad proponen dos enfoques:

• VQSR (Recalibración de Puntuación de Calidad de Variantes): Un método estadísticamente sofisticado. Requiere un número mínimo de muestras.
• Filtrado duro: Un filtrado simple basado en umbrales. Se utiliza en cualquier contexto.

Aquí se resumen las reglas de selección entre ambos métodos y las opciones prácticas.

VQSR — La estética del GaussianMixtureModel

La idea central de VQSR se puede resumir en una frase:

Utilizar un conjunto conocido de variantes verdaderas (HapMap, 1000G) como datos de entrenamiento para aprender mediante estadística multivariante la probabilidad de que cada variante sea real.

Vectores de características

Cada candidato a variante se representa mediante las siguientes estadísticas:

• QD (Calidad por Profundidad): Confianza normalizada por la cobertura.
• FS (Fisher Strand): Sesgo de hebra.
• SOR (Relación de Probabilidades de Hebra): Sesgo hacia la hebra alternativa.
• MQ (MAPQ medio): Promedio de la calidad de mapeo de las lecturas.
• MQRankSum: Diferencia en el MAPQ entre las lecturas de referencia y las alternativas.
• ReadPosRankSum: Sesgo de la posición dentro de la lectura.

En este vector de 6 a 7 dimensiones, las variantes verdaderas tienden a agruparse en una región específica, mientras que los falsos positivos se agrupan en otra región distinta.

Entrenamiento del GMM

VQSR entrena dos Modelos de Mezcla Gaussiana (Gaussian Mixture Model, GMM) en este espacio de vectores:

• GMM positivo: Entrenado con conjuntos verdaderos de HapMap y 1000G.
• GMM negativo: Entrenado con artefactos repetitivos como PoN.

Para cada llamada, la razón de verosimilitud logarítmica entre los dos GMM (VQSLOD) se utiliza como medida de confianza. Un VQSLOD más alto indica una mayor probabilidad de que la variante sea real.

Tranche de sensibilidad

La única decisión que debe tomar el usuario de VQSR es: "¿En qué punto cortamos la sensibilidad respecto al conjunto de referencia?"

• Tranche del 99,7 %: Mantiene el 99,7 % del conjunto verdadero. Tiene relativamente más falsos positivos.
• Tranche del 99,0 %: Mantiene el 99,0 % del conjunto verdadero. Reduce los falsos positivos.
• Tranche del 90,0 %: Muy estricto.

Por lo general, se utiliza el 99,5 % para SNP y el 99,0 % para INDEL. Consulte la curva de tranche (SensitivityVSFalsePositive) y decida en función de las políticas del equipo.

Comando de ejecución

Se ejecuta en dos pasos.

# 1) Entrenar
gatk VariantRecalibrator \
    -R hg38.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=true,prior=12 omni.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10 1000G.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2 dbsnp.vcf.gz \
    -an QD -an FS -an SOR -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.snp.recal \
    --tranches-file cohort.snp.tranches

# 2) Aplicar
gatk ApplyVQSR \
    -R hg38.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.snp.recal \
    --tranches-file cohort.snp.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.snp.recalibrated.vcf.gz

Las INDEL se ejecutan por separado (modo INDEL).

Cuando no se usa VQSR — Filtrado duro

Por convención, VQSR requiere una cohorte de muestras de 30 o más para que el entrenamiento del GMM sea estable. Por debajo de este umbral, el GMM puede sobreajustarse al ruido, lo que empeora la calidad del filtrado.

• Diagnóstico clínico de una sola muestra.
• Cohortes pequeñas con 20 muestras o menos.
• Especies no humanas (falta un conjunto de referencia).

En estos casos, el filtro duro es el estándar.

Umbrales recomendados por GATK para el filtro duro

SNP:

gatk VariantFiltration \
    -V raw.snp.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    --filter-expression "SOR > 3.0" --filter-name "SOR3" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
    --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
    -O snp.filtered.vcf.gz

INDEL:

gatk VariantFiltration \
    -V raw.indel.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 200.0" --filter-name "FS200" \
    --filter-expression "SOR > 10.0" --filter-name "SOR10" \
    --filter-expression "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20" \
    -O indel.filtered.vcf.gz

Este es el umbral derivado por Broad a partir de una cohorte de más de 30 individuos. Estos valores se basan en la referencia humana hg38.

VQSR vs Filtro duro — Reglas de decisión

Cálculo manual — ¿Por qué GMM?

¿Por qué no se puede tratar VQSR simplemente como un filtro basado en umbrales en lugar de usar GMM? La razón es que las características están correlacionadas.

• Las regiones con alta cobertura tienden a tener QD alto y FS bajo.
• Las regiones repetitivas tienden a tener MQ bajo y SOR alto.

Es decir, si se establecen umbrales individuales por separado, se pueden pasar por alto "variantes normales de baja cobertura" o no capturar correctamente los "artefactos en regiones repetitivas". El GMM aprende automáticamente estas correlaciones.

Ejemplo hipotético: si las variantes verdaderas se agrupan cerca de un vector de 6 dimensiones (QD=15, FS=1, SOR=1, MQ=60, MQRankSum=0, ReadPosRankSum=0) y los falsos positivos se agrupan cerca de otro cluster (QD=3, FS=50, SOR=8, MQ=35, MQRankSum=-10, ReadPosRankSum=-5), el GMM encuentra el límite de decisión entre estos dos clusters. Las nuevas llamadas se clasifican según este límite.

El filtro duro aproxima este límite mediante cajas alineadas a los ejes. Esto reduce la sensibilidad, pero hace que la interpretación sea más clara.

Genotipado conjunto — Paso previo obligatorio antes de VQSR

VQSR requiere un VCF integrado de múltiples muestras, no un VCF de una sola muestra. Esta integración se denomina genotipado conjunto (Joint Genotyping).

# 1) Integrar los GVCF en GenomicsDB
gatk GenomicsDBImport \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.bed \
    -V S1.g.vcf.gz -V S2.g.vcf.gz ... -V S30.g.vcf.gz

# 2) Realizar la llamada conjunta de varias muestras
gatk GenotypeGVCFs \
    -R hg38.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz

En este proceso, la información de todas las muestras se combina por cada posición y se calcula la frecuencia alélica poblacional. El VCF resultante sirve como entrada para VQSR.

Diagnóstico cuando falla VQSR

Razones comunes por las que falla VQSR.

1. Falta de variantes de entrenamiento: Si hay menos de 500 llamadas correspondientes al conjunto de referencia (HapMap), el GMM no se entrena.
2. Anomalía en la curva de tranche: Si la curva de sensibilidad-fp es escalonada, indica un problema con los datos de entrenamiento.
3. Los INDEL son mucho peores que los SNP: Los INDEL siempre son más difíciles de procesar que los SNP. Reduzca el tranche al 99.0%.

El estándar cuando falla VQSR es recurrir a filtros duros (Hard Filter).

Pipeline de procesamiento final del VCF

Después de aplicar VQSR o filtros duros, se ajusta el VCF final de la siguiente manera:

# 1) Conservar solo las variantes que superan el filtro
bcftools view -f PASS cohort.vcf.gz -Oz -o cohort.pass.vcf.gz

# 2) Normalizar (alineamiento a la izquierda y representación mínima)
bcftools norm -f hg38.fa cohort.pass.vcf.gz -Oz -o cohort.norm.vcf.gz

# 3) Indexar
bcftools index cohort.norm.vcf.gz

# 4) Obtener estadísticas
bcftools stats cohort.norm.vcf.gz > cohort.stats.txt

bcftools norm es una etapa que se omite con frecuencia. Como cada herramienta puede representar los indels de manera distinta, especialmente en el alineamiento a la izquierda, la normalización es imprescindible. Sin ella, al comparar dos VCF una misma variante puede parecer diferente.

Mapeo CS

• GMM (modelo de mezcla gaussiana): consulte el episodio de DryBench «Fundamentos de GMM» (gmm-basics), una técnica central de estimación de densidad multivariante.
• Tratamiento de correlaciones multivariantes: captura automáticamente la correlación entre características mediante una matriz de covarianza.
• Regla de decisión por umbral: una tranche de sensibilidad corresponde a seleccionar un punto sobre la curva ROC.

Conclusión

VQSR y los filtros duros son herramientas que deben elegirse según la situación. Recordar el umbral orientativo de unas 30 muestras facilita las decisiones iniciales de un proyecto. En el siguiente episodio (S12) saldremos del marco de GATK para estudiar el enfoque de DeepVariant, que llama variantes mediante aprendizaje profundo.

Este episodio completa el mapa del pipeline GATK4 (S06–S11). S12–S16 continuarán con llamadores alternativos, variantes estructurales, CNV e imputación; S17–S21, con cuantificación de RNA-seq y DGEA. Aquí concluye la primera parte de la serie «FASTQ to Paper», dedicada a GATK4.

Si desea profundizar

• DePristo M. et al. (2011), A framework for variation discovery and genotyping. Nature Genetics 43:491. Método original de VQSR.
• Guía oficial de VQSR de GATK (https://gatk.broadinstitute.org/hc/en-us/articles/360035531612): fuente primaria de opciones y recursos.
• Guía de filtros duros de GATK (https://gatk.broadinstitute.org/hc/en-us/articles/360035531112): origen empírico de los umbrales.
• Broad Institute BroadE — Variant Filtering, clase en YouTube con ejemplos de interpretación de curvas por tranche.
• Documentación oficial de bcftools (https://samtools.github.io/bcftools/): fuente primaria de las opciones de normalización y estadísticas.

Memorizar las reglas de selección entre VQSR y los filtros estrictos (Hard Filters) ya cubre la mitad del proceso de decisión en un flujo clínico. En el próximo capítulo, pasemos a DeepVariant.