De la especificidad tisular a la estructura poblacional: volviendo a los fundamentos y estimando directamente

En S42, abordamos la "estratificación poblacional" como un problema que genera falsos positivos en los estudios de asociación del genoma completo (GWAS) y aprendimos cómo corregir su influencia mediante PCA y LMM. Sin embargo, no tratamos realmente "cómo estimar la propia estructura poblacional". En este capítulo, profundizamos un nivel más para verificar directamente, utilizando únicamente datos de genotipado, qué tan diferentes son los orígenes ancestrales entre los individuos y si cada SNP se hereda normalmente dentro de esa población.

Abordaremos dos aspectos. Primero, la prueba de equilibrio de Hardy-Weinberg (HWE) para determinar si un SNP individual cumple con el supuesto de apareamiento aleatorio. Segundo, el modelo de máxima verosimilitud de ADMIXTURE para estimar cómo el genoma completo de un individuo se compone de una mezcla de varias poblaciones ancestrales.

Deducción del modelo de verosimilitud de HWE y ADMIXTURE

Equilibrio de Hardy-Weinberg: la proporción esperada de genotipos resultante del apareamiento aleatorio

Si una población se aparea al azar y no hay selección natural, migración ni mutación, la proporción de genotipos queda determinada exactamente por las frecuencias alélicas $p$ (alelo de referencia) y $q = 1-p$ (alelo alternativo) de la siguiente manera:

$P(\text{AA}) = p^2, \quad P(\text{Aa}) = 2pq, \quad P(\text{aa}) = q^2$

Durante la etapa de control de calidad (QC), los SNPs que se desvían significativamente de esta proporción pueden deberse a errores en la determinación del genotipo, pérdida de alelos, contaminación de muestras, inclusión de parientes, el efecto Wahlund (déficit de heterocigotos artificial causado por combinar subpoblaciones heterogéneas como una sola), endogamia y, con menor frecuencia, señales de selección natural. Entre las causas relacionadas con el genotipado, las fuentes prácticas son las que suelen requerir una verificación prioritaria. Se cuantifica la discrepancia entre los valores observados y los esperados mediante una prueba de bondad de ajuste chi-cuadrada (aunque para fines educativos es simple, para variantes raras con recuentos esperados bajos, la prueba exacta es más adecuada).

$\chi^2 = \sum \frac{(O_i - E_i)^2}{E_i}, \quad df = 1$

(Dado que estimamos una frecuencia alélica a partir de los datos entre las tres categorías, los grados de libertad son $3-1-1 = 1$).

Ejemplo de cálculo manual: encontrar SNPs que violan HWE en una cohorte de 100 individuos

Se observaron 50 individuos con genotipo AA, 30 con Aa y 20 con aa para un SNP específico. La frecuencia alélica es $p = \frac{2(50)+30}{200} = 0.65$ y $q = 0.35$.

Número esperado de genotipos: $E_{AA} = 100 \times 0.65^2 = 42.25$, $E_{Aa} = 100 \times 2(0.65)(0.35) = 45.5$, $E_{aa} = 100 \times 0.35^2 = 12.25$.

$\chi^2 = \frac{(50-42.25)^2}{42.25} + \frac{(30-45.5)^2}{45.5} + \frac{(20-12.25)^2}{12.25} \approx 1.42 + 5.28 + 4.91 = 11.61$

El valor crítico para un nivel de significancia de 0.05 en la distribución chi-cuadrado con 1 grado de libertad es 3.84. Dado que $11.61 > 3.84$, este SNP se desvía significativamente del equilibrio de Hardy-Weinberg (HWE) ($p < 0.001$); el patrón de tener muchos menos heterocigotos (Aa) y más homocigotos de lo esperado es, comúnmente, una señal de errores en la determinación del genotipo o de mezcla poblacional oculta.

ADMIXTURE: Estimación simultánea de las proporciones ancestrales y las frecuencias alélicas por grupo ancestral

Supongamos que el individuo $i$ tiene una mezcla de proporciones $q_{i1}, \dots, q_{iK}$ ($\sum_k q_{ik}=1$) de $K$ grupos ancestrales, y que la frecuencia del alelo en el SNP $j$ dentro del grupo ancestral $k$ es $p_{kj}$. La frecuencia alélica esperada para cada individuo para el SNP $j$ en el individuo $i$ es el promedio ponderado de cada componente ancestral.

$\pi_{ij} = \sum_{k=1}^{K} q_{ik} \, p_{kj}$

Si asumimos que el HWE se cumple dentro de cada componente ancestral, la probabilidad de observar el genotipo del SNP $j$ en el individuo $i$ es equivalente a la fórmula de HWE, utilizando $\pi_{ij}$ como frecuencia alélica.

$P(x_{ij}=0) = (1-\pi_{ij})^2, \quad P(x_{ij}=1) = 2\pi_{ij}(1-\pi_{ij}), \quad P(x_{ij}=2) = \pi_{ij}^2$

El objetivo de ADMIXTURE es maximizar la verosimilitud total $\sum_i \sum_j \log P(x_{ij} \mid q_i, p)$ con respecto a $\{q_{ik}\}$ y $\{p_{kj}\}$. Mientras que STRUCTURE (Pritchard et al., 2000) resolvió esto mediante muestreo MCMC bayesiano, ADMIXTURE (Alexander et al., 2009) optimiza directamente la misma función de verosimilitud utilizando relajación por bloques: se actualizan secuencialmente los bloques de $p$ mediante programación cuadrática secuencial, manteniendo $q$ fijo, y luego se actualiza $q$ manteniendo $p$ fijo, repitiendo este proceso junto con una aceleración cuasi-Newton, lo que permite una convergencia decenas o cientos de veces más rápida. Aunque comparte una intuición similar al algoritmo EM al tratar iterativamente la asignación de ancestros latentes, es importante tener en cuenta que el propio algoritmo de optimización de ADMIXTURE no es EM.

Ejemplo práctico: comparación de la verosimilitud de genotipos para un SNP en dos poblaciones ancestrales

Supongamos que las frecuencias alélicas de un SNP en $K=2$ poblaciones ancestrales son $p_1 = 0.8$ (población 1) y $p_2 = 0.2$ (población 2). Para dos individuos, A y B, observados como heterocigotos ($x=1$), con proporciones de ancestro $q_A = (0.5, 0.5)$ y $q_B = (0.9, 0.1)$ respectivamente:

$\pi_A = 0.5(0.8) + 0.5(0.2) = 0.5 \implies P(x=1 \mid q_A) = 2(0.5)(0.5) = 0.5$

$\pi_B = 0.9(0.8) + 0.1(0.2) = 0.74 \implies P(x=1 \mid q_B) = 2(0.74)(0.26) = 0.3848$

Aunque la observación de heterocigosidad es la misma, la verosimilitud varía según la proporción de ancestro del individuo. ADMIXTURE encuentra los $q_i$ que maximizan el producto (suma logarítmica) de estas verosimilitudes a lo largo de cientos de miles de SNPs, utilizando el método de relajación por bloques con aceleración cuasi-Newton descrito anteriormente.

Práctica con PLINK2 + R

#!/usr/bin/env bash
# Pipeline de filtrado HWE + ejecución de ADMIXTURE
set -e

# 1. Eliminar SNPs con valor p de HWE < 1e-6 (valor umbral de ejemplo; ajustar según tamaño de muestra y diseño del estudio)
# --hwe de PLINK2 utiliza por defecto la prueba exacta (diferente de la aproximación de chi-cuadrado usada con fines educativos)
# Si el diseño es caso-control, se debe extraer solo el grupo control con --keep antes de aplicar este filtrado
plink2 --bfile raw_genotypes --hwe 1e-6 --make-bed --out hwe_filtered

# 2. Poda de LD y conversión a binario PLINK 1 (requisito de entrada para ADMIXTURE)
plink2 --bfile hwe_filtered --indep-pairwise 50 5 0.2 --out ld_pruned
plink2 --bfile hwe_filtered --extract ld_pruned.prune.in --make-bed --out admixture_input

# 3. Ejecutar ADMIXTURE variando el número de grupos ancestrales desde K=2 hasta K=5 (seleccionar K óptimo mediante validación cruzada)
for K in 2 3 4 5; do
  admixture --cv admixture_input.bed $K | tee log_K${K}.out
done

# Reproducir en código el ejemplo de cálculo de las dos primeras componentes principales
p1 <- 0.8; p2 <- 0.2

genotype_likelihood <- function(q, p1, p2, x) {
  pi_val <- q[1] * p1 + q[2] * p2
  switch(x + 1, (1 - pi_val)^2, 2 * pi_val * (1 - pi_val), pi_val^2)
}

cat("Verosimilitud de heterocigosis del individuo A (q=0.5,0.5):", genotype_likelihood(c(0.5, 0.5), p1, p2, 1), "\n")
cat("Verosimilitud del genotipo heterocigoto del individuo B (q=0.9,0.1):", genotype_likelihood(c(0.9, 0.1), p1, p2, 1), "\n")

# Reproducción de la prueba chi-cuadrado de HWE
obs <- c(AA = 50, Aa = 30, aa = 20); n <- sum(obs)
p_allele <- (2 * obs["AA"] + obs["Aa"]) / (2 * n)
exp_counts <- n * c(p_allele^2, 2 * p_allele * (1 - p_allele), (1 - p_allele)^2)
chisq <- sum((obs - exp_counts)^2 / exp_counts)
cat(sprintf("Estadístico chi-cuadrado de HWE: %.2f (p-valor: %.2e)\n", chisq, pchisq(chisq, df = 1, lower.tail = FALSE)))

Mapeo de CS

• Ascenso de coordenadas de bloques (Block Coordinate Ascent): Es el método de optimización que utiliza realmente ADMIXTURE. En lugar de optimizar todas las variables a la vez, se fija $q$ y se optimiza el bloque $p$ mediante programación cuadrática secuencial, luego se fija $p$ y se optimiza $q$, alternando entre ambos. Este enfoque de descenso/ascenso de coordenadas es común en problemas de optimización a gran escala; al añadir aceleración cuasi-newtoniana, converge mucho más rápido que el muestreo MCMC de STRUCTURE.
• Modelo de variable latente (Latent Variable Model): Asumir variables latentes no observadas (la asignación de la proporción ancestral del individuo) y modelar a partir de ellas el proceso generador de los datos observados es una estructura similar a la del agrupamiento k-medias o el modelo de mezcla gaussiana (GMM). Sin embargo, a diferencia de que el GMM se entrena habitualmente con EM, ADMIXTURE entrena este mismo tipo de modelo de variable latente mediante optimización numérica directa en lugar de EM.
• Modelo de mezcla (Mixture Model): La expresión $\pi_{ij} = \sum_k q_{ik}p_{kj}$ de ADMIXTURE corresponde a la estructura típica de un modelo de mezcla, que considera que cada observación se genera como una suma ponderada de varias componentes latentes.

Problemas comunes

• Aplicar el filtro de HWE a toda la cohorte: En un diseño caso-control, la prueba de HWE debe realizarse únicamente en el grupo control. La violación de HWE en el grupo de casos puede ser una señal de asociación con la enfermedad; aplicar el filtro aquí eliminaría señales reales.
• Omitir el ajuste de LD y ejecutar ADMIXTURE: Si se incluyen directamente SNPs fuertemente ligados, unas pocas regiones genómicas (especialmente bloques de LD a larga distancia y regiones inversas) pueden influir excesivamente en los resultados de la estimación ancestral. Dado que ADMIXTURE asume independencia entre marcadores, se recomienda encarecidamente el ajuste de LD como estándar, y es seguro verificar qué tan sensibles son los resultados a la intensidad del ajuste.
• Fijar arbitrariamente la K óptima: El "número de grupos ancestrales K" no tiene un valor verdadero biológicamente determinado. Se debe elegir el K con el error de validación cruzada (CV error) más bajo basado en los datos, pero es seguro presentar e interpretar los resultados para varios valores de K simultáneamente.

Para profundizar más

El texto ha sido reestructurado directamente por el equipo de investigación de BPD. Para profundizar, consulte el artículo original y la documentación oficial.

• Artículo original de ADMIXTURE: Alexander, Novembre & Lange (2009), Fast model-based estimation of ancestry in unrelated individuals, Genome Research, 19:1655–1664.
• ARTÍCULO ORIGINAL SOBRE ESTRUCTURA: Pritchard, Stephens y Donnelly (2000), Inferencia de la estructura poblacional mediante datos genotípicos multilocus, Genetics, 155(2): 945–959.
• ARTÍCULO ORIGINAL SOBRE PLINK: Purcell y colaboradores (2007), PLINK: un conjunto de herramientas para análisis de asociación a nivel del genoma y análisis de ligamiento basados en la población, American Journal of Human Genetics, 81(3): 559–575.
• DOCUMENTACIÓN OFICIAL DE PLINK2: cog-genomics.org/plink/2.0.

Hemos aprendido a estimar directamente la estructura poblacional. Ahora pasaremos a una de las regiones más polimórficas y clínicamente relevantes del genoma humano: los genes HLA, que el sistema inmunitario utiliza para distinguir entre lo propio y lo ajeno.

En el siguiente capítulo, S47, trataremos los principios de la determinación de tipos HLA y la imputación.