feat: customizable homozygous reference sample selection

[dlopez-bioinfo]

Motivación

Actualmente, la contabilización de homocigoto referencia (N_HOM_REF) en afquery se realiza de forma implícita: para una posición dada, se consideran hom-ref todas las muestras que (a) están cubiertas por el BED de captura de su tecnología en esa posición y (b) no aparecen con una variante en el VCF, o aparecen pero no son ni heterocigotas ni homocigotas alternativas ni FAIL.

Esto puede generar falsos hom-ref: el VCF puede no contener la posición simplemente porque ninguna muestra tenía variante allí, sin que ello garantice que la posición haya sido efectivamente interrogada en la cohorte con calidad suficiente.

Ver src/afquery/query.py:276 (y análogos en :454, :532):

N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL

donde eligible proviene de _compute_eligible (query.py:178), que solo aplica el BED de captura y el filtro de muestras.

Propuesta

Hacer configurable el criterio con el que una posición del BED de captura se considera "apta" para contabilizar hom-ref. Las dos alternativas operan al mismo nivel conceptual: refinan el BED efectivo de una tecnología a partir de evidencia observada en las variantes de la cohorte con ese mismo BED. Una posición que no pasa el criterio no cuenta en AN para las muestras de esa tecnología — es equivalente a que esté fuera del BED de captura.

Nota: una tercera vía — precomputar cobertura real desde BAM/CRAM en la ingesta y refinar el BED efectivo — se descarta como feature del core de afquery. Ese flujo puede realizarse externamente y subirse como BED por muestra (ver #29).

Alternativa A — Umbral de número de variantes PASS por posición

Una posición cuenta como apta para hom-ref solo si hay al menos K variantes PASS observadas en muestras de la cohorte con el mismo BED de captura (p. ej. K=1: exigir que la posición haya sido llamada con PASS al menos una vez en esa cohorte técnica). El umbral es un simple conteo.

• Coste ingesta: bajo. La información ya está en los bitmaps het/hom/fail almacenados en Parquet; solo hay que consultarlos y agregarlos por tecnología en query time (o precomputar una máscara por tecnología).
• Coste implementación: bajo-medio. Cambios aislados en _compute_eligible / query() y parámetros CLI.
• Limitación: un único PASS puede no ser suficiente garantía de calidad; no mira los valores de calidad subyacentes.

Alternativa B — Umbral sobre campos de calidad de las variantes observadas

Generalización de A: una posición cuenta como apta para hom-ref solo si existe al menos una muestra (o K muestras) de la cohorte con el mismo BED de captura que tenga una variante en esa posición con valores de campos del VCF por encima de un umbral (p. ej. DP ≥ 30, GQ ≥ 20, QUAL ≥ 50). La lógica es: si alguna muestra de la cohorte interrogó esta posición con buena calidad, entonces la posición está realmente cubierta en esa tecnología y se puede inferir hom-ref para el resto; si nadie alcanzó el umbral, la posición se excluye de AN para esa tecnología.

Ejemplo: en una posición X dentro del BED de WES_v1, si al menos una muestra WES_v1 tiene una variante con DP ≥ 30, la posición X es apta y las demás muestras WES_v1 sin variante ahí cuentan como hom-ref. Si ninguna muestra WES_v1 alcanza ese umbral en X, la posición X se descarta para WES_v1 y no contribuye a AN.

• Coste ingesta: medio. Hoy no se almacenan campos como DP, GQ, QUAL por variante durante la ingesta; habría que leerlos del VCF y guardar al menos la información mínima para poder evaluar el criterio en query time, o precomputar una máscara de posiciones aptas por tecnología.
• Coste implementación: medio. Decidir qué se guarda (p. ej. max(DP) por posición y tecnología) y dónde (columnas adicionales en Parquet o una estructura auxiliar por tecnología). Diseñar la API para configurar umbrales.
• Más preciso que A porque permite exigir calidad mínima, no solo presencia.

Nota sobre almacenamiento

Ambas alternativas pueden implementarse sin persistir datos por muestra-posición adicionales: basta con mantener, por tecnología, una máscara de posiciones aptas (calculada en ingesta o al vuelo a partir de los bitmaps existentes para A; calculada en ingesta a partir de los campos del VCF para B). Esto limita el impacto en tamaño de BD.

Estudio previo requerido

Antes de implementar, se debe producir un documento que evalúe:

1. Coste de ingesta estimado para cada alternativa (segundos/muestra, escala con N muestras).
2. Coste de almacenamiento para B (tamaño de la máscara de posiciones aptas por tecnología). A reutiliza lo ya almacenado.
3. Impacto en tiempo de query de cada alternativa (máscara precomputada vs cálculo al vuelo).
4. Impacto en las frecuencias calculadas en datos reales: cuánto cambian AF/AN/N_HOM_REF en un dataset representativo al aplicar cada criterio, para decidir si B aporta precisión suficiente sobre A para justificar su coste extra.
5. Combinabilidad: A y B se pueden combinar (p. ej. "al menos K variantes con DP ≥ 30").
6. Elección de campos del VCF para B (¿qué campo(s) son los más informativos y qué umbrales por defecto razonables?).
7. API propuesta (flags CLI, parámetros en QueryParams, compatibilidad hacia atrás con las BD actuales).

Resultado esperado del estudio

• Recomendación priorizada de qué alternativa(s) implementar.
• Plan de migración: las BD existentes no tienen la información necesaria para B (y posiblemente tampoco la máscara precomputada para A) — debe definirse el comportamiento por defecto (seguir usando el criterio actual) y cómo se habilitan los nuevos modos.

[dlopez-bioinfo]

Findings from PR #30 Review — Relevant to This Issue

During the code review of PR #30 (variant_info module), important behavior in the treatment of FILTER≠PASS samples was identified that directly impacts how this issue should be resolved.

Current Behavior

In _query_region_inner() (query.py:276):

het_elig = het_bm & eligible
hom_elig = hom_bm & eligible
AC = (len((het_elig | hom_elig) & haploid_elig) + len(het_elig & diploid_elig) + 2 * len(hom_elig & diploid_elig))
N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL

FILTER≠PASS samples contribute to AN/AC identically to PASS samples, despite their genotype being uncertain. Example:

Sample S00: 1/1, FILTER=LowQual
Sample S01: 1/1, FILTER=PASS

AN contribution: both +2 alleles
AC contribution: both +2 alleles (if diploid)

Such samples are only reported separately as N_FAIL in query output, but are not excluded from AN.

Inconsistency in variant_info

In variant_info (introduced in PR #30), a FILTER≠PASS sample is reported as:

{"genotype": "alt", "filter_pass": false}

The "alt" label creates a false impression that the sample receives special treatment, when in reality it counts toward AN identically to a PASS hom_ref.

Implications for Issue #28

This issue proposes refining which samples/positions contribute to AN based on quality criteria. FILTER≠PASS treatment is a specific case of this:

1. Alternative A (PASS variant threshold): should clarify whether FILTER≠PASS samples, which currently count, should continue to count.

2. Alternative B (quality thresholds): the genotype uncertainty of FILTER≠PASS samples is precisely the kind of uncertainty that B aims to address.

3. Current N_HOM_REF calculation: defined as len(eligible) - N_HET - N_HOM_ALT - N_FAIL. If a sample has FILTER≠PASS but is in hom_bitmap, it contributes to N_HOM_ALT (not N_FAIL), regardless of FILTER status.

Recommendations

When resolving #28, consider:

1. Explicitly document FILTER≠PASS handling in AN/AC for each proposed alternative.
2. Add a test verifying how AN/AF change when FILTER≠PASS samples are included vs. excluded.
3. Consider Alternative C: exclude FILTER≠PASS samples from AC but retain them in AN (more conservative than current behavior).
4. API design: New parameters for feat: customizable homozygous reference sample selection #28 could include a flag like include_fail_samples_in_an (default: true for backward-compatibility).

---

Note: This finding indicates that variant_info documentation should clarify that genotype="alt" means "uncertain genotype", not "excluded from AN".