S42 Barreras de la GWAS ultrarrápida — Explosión computacional de big data con 500.000 personas

En S42 (PLINK2 y modelos lineales mixtos), abordamos el modelo lineal mixto (LMM) que controla los falsos positivos causados por la estructura poblacional mediante la matriz de parentesco genético (GRM) y la descomposición en valores propios. En cohortes del orden de miles de individuos, las herramientas LMM tradicionales como GEMMA o GCTA funcionan bien.

Sin embargo, con la llegada de la era del big data de cohortes de cientos de miles de personas, como UK Biobank y FinnGen, los LMM tradicionales se han topado frontalmente con un muro computacional. Cuando el número de muestras aumenta a $N = 500.000$, solo la capacidad de memoria para almacenar la matriz GRM de tamaño $N \times N$ alcanza aproximadamente $500.000 \times 500.000 \times 8 \text{ bytes} \approx 2 \text{ TB}$. Sin 2 TB de memoria no es posible cargar la matriz, y la complejidad computacional de la descomposición en valores propios $O(N^3)$ multiplica el tiempo de cálculo hasta varios meses.

Además, en el análisis de enfermedades raras con una proporción desequilibrada entre casos y controles inferior a 1:100, el supuesto de normalidad asintótica del test de Wald tradicional comienza a debilitarse, lo que puede inflar significativamente los valores P falsos positivos. En esta sección examinamos la metodología de aprendizaje en bloque Ridge de dos pasos de REGENIE, surgida para superar estas limitaciones de big data, y el principio de aproximación del punto de silla (SPA, Saddlepoint Approximation) adoptado por SAIGE. Ambas herramientas apuntan al mismo problema (datos grandes e desequilibrados), pero emplean algoritmos distintos para la corrección.

Derivación de los algoritmos clave de REGENIE y SAIGE

1. REGENIE: Regresión Ridge de dos pasos y partición por bloques (Block Binning)

REGENIE modela los efectos poligénicos del genoma completo mediante regresión Ridge por bloques, sin necesidad de cargar toda la matriz GRM $N \times N$.

Paso 1 (Step 1): Cálculo del predictor genómico Ridge

Dividimos los $M_{QC}$ SNPs que han pasado el control de calidad en bloques consecutivos de 1.000 unidades: $B_1, \dots, B_K$. Dentro de cada bloque $B_k$, ajustamos modelos de regresión Ridge para múltiples parámetros de regularización $\lambda_l$ y obtenemos los predictores locales $\hat{\mathbf{w}}_{k, l}$.

$\hat{\mathbf{y}}_{k, l} = \mathbf{X}_{B_k} (\mathbf{X}_{B_k}^T \mathbf{X}_{B_k} + \lambda_l \mathbf{I})^{-1} \mathbf{X}_{B_k}^T y$

Aquí, $\hat{\mathbf{y}}_{k, l}$ es el predictor ajustado de dimensión $N$ construido a partir de la combinación del bloque $B_k$ y el parámetro de regularización $\lambda_l$, y la parte entre paréntesis $(\mathbf{X}_{B_k}^T \mathbf{X}_{B_k} + \lambda_l \mathbf{I})^{-1} \mathbf{X}_{B_k}^T y$ corresponde a los coeficientes de la regresión Ridge. Se aplica una regresión Ridge de Nivel 2 a los predictores de bloque acumulados para componer el valor final de predicción poligénica $\hat{y}_{\text{poly}}$. La cantidad que debe cargarse en la memoria a la vez se reduce a $O(N \cdot \text{bsize})$ por unidad de tamaño de bloque (por defecto, 1.000) en lugar del total de SNPs $M_{QC}$.

Paso 2: Prueba de un solo SNP con aplicación LOCO

Para cada SNP $j$, se realiza una regresión simple introduciendo como covariable el predictor $\hat{y}_{\text{poly}, -c}$ que excluye el cromosoma correspondiente (LOCO: Leave-One-Chromosome-Out).

$y = \mathbf{x}_j \beta_j + \hat{y}_{\text{poly}, -c} \gamma + \mathbf{W} \boldsymbol{\alpha} + \epsilon$

2. SAIGE: Derivación de la aproximación del punto de silla (Saddlepoint Approximation, SPA)

SAIGE resuelve el problema en el que se rompe el supuesto de normalidad de la estadística $S = \mathbf{x}^T (y - \hat{\mu})$ en datos de enfermedades raras o con proporciones extremas de controles (por ejemplo, 500 casos frente a 499.500 controles) mediante la aproximación del punto de silla (SPA).

La siguiente función generadora de cumulantes (CGF) es una forma simplificada que omite las covariables y la corrección del modelo mixto para ilustrar el principio. En la implementación real, SAIGE construye la CGF incluyendo los residuos genotípicos corregidos por el modelo mixto y el procedimiento de razón de varianzas. La CGF $K_S(t) = \ln \mathbb{E}[e^{t S}]$ y sus derivadas primera y segunda son las siguientes:

$K_S(t) = \sum_{i=1}^N \ln \left( 1 - \hat{\mu}_i + \hat{\mu}_i e^{t x_i} \right) - t \sum_{i=1}^N x_i \hat{\mu}_i$

$K_S'(t) = \sum_{i=1}^N \frac{x_i \hat{\mu}_i e^{t x_i}}{1 - \hat{\mu}_i + \hat{\mu}_i e^{t x_i}} - \sum_{i=1}^N x_i \hat{\mu}_i, \quad K_S''(t) = \sum_{i=1}^N \frac{x_i^2 \hat{\mu}_i (1 - \hat{\mu}_i) e^{t x_i}}{\left( 1 - \hat{\mu}_i + \hat{\mu}_i e^{t x_i} \right)^2}$

Dado un valor observado de la estadística, $s$, si se encuentra la solución $\hat{\zeta}$ de la ecuación de punto de silla $K_S'(\hat{\zeta}) = s$, entonces, a través de la fórmula de Lugannani-Rice, se puede aproximar la probabilidad de cola exacta $P(S \ge s)$.

$P(S \ge s) \approx 1 - \Phi(\hat{w}) + \phi(\hat{w}) \left( \frac{1}{\hat{w}} - \frac{1}{\hat{u}} \right)$

$\hat{w} = \text{sgn}(\hat{\zeta}) \sqrt{2 \left( \hat{\zeta} s - K_S(\hat{\zeta}) \right)}, \quad \hat{u} = \hat{\zeta} \sqrt{K_S''(\hat{\zeta})}$

Ejemplo de cálculo a pequeña escala: por qué la aproximación normal falla en datos extremadamente desequilibrados

Veamos un ejemplo muy pequeño en el que podemos calcular exactamente para comprobar cuánto puede fallar la aproximación normal. Supongamos que seleccionamos solo 5 individuos que poseen la variante rara ($x_i = 1$). Bajo la hipótesis nula, la probabilidad de que cada uno desarrolle la enfermedad es la misma, $\hat{\mu}_i = 0.05$, y son independientes. La estadística de puntuación es $S = \sum_i (y_i - \hat{\mu}_i)$, y el valor observado es $s = 3 - 5 \times 0.05 = 2.75$ si 3 de estos 5 individuos son realmente pacientes.

El número de pacientes $K$ sigue una distribución binomial $K \sim \text{Binomial}(5, 0.05)$, por lo que la probabilidad de obtener un valor observado o superior puede calcularse exactamente.

$P(K \ge 3) = \sum_{k=3}^{5} \binom{5}{k} 0.05^k \, 0.95^{5-k} \approx 0.001128 + 0.0000297 + 0.0000003 \approx 1.16 \times 10^{-3}$

Por otro lado, si aplicamos la aproximación normal a la misma situación, tenemos $\text{Var}(S) = 5 \times 0.05 \times 0.95 = 0.2375$, por lo que

$Z = \frac{2.75}{\sqrt{0.2375}} \approx 5.64, \quad P(Z \ge 5.64) \approx 8.6 \times 10^{-9}$

Entre la probabilidad exacta (aprox. $1.16 \times 10^{-3}$) y la aproximación normal (aprox. $8.6 \times 10^{-9}$) existe una diferencia de unas 130,000 veces. Si se adopta la aproximación normal tal cual como el p-valor de GWAS, se informará erróneamente una significancia absurdamente fuerte para una señal que, en realidad, podría aparecer por azar.

La aproximación de punto de silla (SPA) es un método para reducir esta distorsión. Sin embargo, la ecuación CGF anterior $K_S'(\hat{\zeta}) = s$ generalmente no tiene una solución en forma cerrada, por lo que en la práctica se debe encontrar la raíz numéricamente (SAIGE también utiliza el mismo método internamente). En este ejemplo, la raíz es $\hat{\zeta} = \ln(28.5) \approx 3.3499$, y al sustituirla en la fórmula de Lugannani-Rice, el valor aproximado de SPA es aproximadamente $3.9 \times 10^{-4}$. Aunque sigue siendo unas 3 veces menor que el valor del cálculo exhaustivo ($1.16 \times 10^{-3}$), es incomparablemente más cercano a la aproximación normal (aprox. $8.6 \times 10^{-9}$). Debido a que es una situación extrema con solo 5 muestras, persiste un error propio de la distribución de red (discreta); a medida que el tamaño de la muestra aumenta, la SPA se aproxima más al cálculo exhaustivo. Ejecute el siguiente código para comparar los tres valores (cálculo exhaustivo, aproximación normal y SPA).

Práctica de REGENIE & R SPA

Pipeline de ejecución de la etapa 2 de REGENIE y código de implementación de la aproximación de punto de silla (SPA) en R.

Etapa 1: Pipeline de REGENIE CLI (Step 1 & Step 2)

#!/usr/bin/env bash
# Canalización ultrarrápida de GWAS en dos pasos con REGENIE
set -e

# [Paso 1] Obtener el predictor poligénico mediante regresión ridge de todo el genoma (--bsize 1000)
regenie \
  --step 1 --bed ukb_qc_array \
  --phenoFile phenotypes.txt --covarFile covariates.txt \
  --bsize 1000 --bt --lowmem --out regenie_step1_out

# [Paso 2] Operación por SNP sobre variantes imputadas + corrección de Firth aproximada propia de REGENIE
# (--firth --approx: opción de REGENIE para corregir el desequilibrio, distinta de la SPA de SAIGE)
regenie \
  --step 2 --pgen ukb_imputed_chr1 \
  --phenoFile phenotypes.txt --covarFile covariates.txt \
  --pred regenie_step1_out_pred.list \
  --bsize 400 --bt --firth --approx --out regenie_step2_chr1_out

Paso 2: Implementación de la corrección mediante aproximación de puntos de silla (SPA) en R

# Implementación de la aproximación de punto de silla (SPA) en R — reproduce el ejemplo manual de cinco personas del texto
cgf_k <- function(t, mu, x) sum(log(1 - mu + mu * exp(t * x))) - t * sum(x * mu)
cgf_k1 <- function(t, mu, x) sum((x * mu * exp(t * x)) / (1 - mu + mu * exp(t * x))) - sum(x * mu)
cgf_k2 <- function(t, mu, x) sum((x^2 * mu * (1 - mu) * exp(t * x)) / (1 - mu + mu * exp(t * x))^2)

# Como K1(t)-s_obs crece de forma monótona con t, Newton-Raphson sin restricciones puede
# divergir fácilmente según el valor inicial (en este ejemplo también diverge a NaN al comenzar en t=0).
# En la práctica se usa una bisección segura que va estrechando el intervalo.
solve_zeta <- function(s_obs, mu, x, lower = -50, upper = 50) {
  f <- function(t) cgf_k1(t, mu, x) - s_obs
  uniroot(f, c(lower, upper), tol = 1e-9)$root
}

# Ejemplo del texto: las cinco personas (x=1) tienen mu=0.05; s observado = 3 - 5*0.05 = 2.75
mu <- rep(0.05, 5); x <- rep(1, 5); s_obs <- 2.75

# 1) Cálculo exhaustivo (valor verdadero)
exact_p <- sum(dbinom(3:5, size = 5, prob = 0.05))

# 2) Aproximación normal (Wald)
var_s <- sum(x^2 * mu * (1 - mu))
normal_p <- 1 - pnorm(s_obs / sqrt(var_s))

# 3) Aproximación de punto de silla (SPA, Lugannani-Rice)
zeta <- solve_zeta(s_obs, mu, x)
w <- sign(zeta) * sqrt(2 * (zeta * s_obs - cgf_k(zeta, mu, x)))
v <- zeta * sqrt(cgf_k2(zeta, mu, x))
spa_p <- 1 - pnorm(w) + dnorm(w) * (1 / w - 1 / v)

cat(sprintf("Valor p del cálculo exhaustivo:    %.5e\n", exact_p))
cat(sprintf("Valor p de la aproximación normal: %.5e\n", normal_p))
cat(sprintf("Valor p de la aproximación de punto de silla (SPA): %.5e\n", spa_p))

Mapeo de CS

• Divide y vencerás (Divide and Conquer): Paradigma algorítmico que, en lugar de cargar cientos de miles de SNPs en la memoria, divide los datos en bloques de 1.000 para realizar una regresión Ridge de primer paso y combinar los resultados.
• Aproximación del punto de silla (Saddlepoint Approximation): Técnica numérica precisa que transforma la probabilidad de cola de la distribución de estadísticos complejos en una trayectoria de paso por el punto de silla en el plano complejo, permitiendo un cálculo rápido con un error mucho menor que la aproximación normal (si la muestra es muy pequeña, puede quedar un error residual característico de la distribución discreta).
• Optimización numérica y PCG: Técnica utilizada por SAIGE para evitar el cálculo de la inversa de la matriz de covarianza y reducir el uso de memoria a $O(N)$ mediante métodos numéricos iterativos para resolver ecuaciones lineales.

Errores frecuentes

• Inclusión de SNPs imputados en el Paso 1: El Paso 1 de REGENIE suele entrenar el predictor poligénico con decenas de miles a cientos de miles de SNPs de array que han pasado el control de calidad (el límite superior exacto depende del tamaño de la muestra y la configuración de memoria). Incluir directamente más de 10 millones de SNPs imputidos aumenta el tiempo de cálculo y el uso de memoria para los bloques Ridge hasta niveles difíciles de manejar.
• Falta de LOCO en el Paso 2: Si no se utiliza el predictor poligénico con Leave-One-Chromosome-Out (LOCO) que excluye el cromosoma correspondiente durante la prueba de variantes individuales del Paso 2, se produce un error al eliminar previamente la señal propia.
• No establecer el umbral de corrección de Firth: Si se desactiva la corrección Firth/SPA para variantes con pocas muestras de caso y se ejecuta una regresión logística simple, los valores P pueden inflarse significativamente.

Para profundizar más

El texto ha sido reescrito directamente por el equipo de investigación de BPD. Para explorar en profundidad los algoritmos GWAS a gran escala más recientes, consulte las siguientes fuentes originales:

• Artículo original de REGENIE: Mbatchou et al. (2021), Cálculo eficiente de la regresión de genoma completo para rasgos cuantitativos y binarios, Nature Genetics, 53(7):1097-1103. (PubMed 34017140)
• Artículo original de SAIGE: Zhou et al. (2018), Control eficiente del desequilibrio caso-control y la relación muestral en estudios de asociación genética a gran escala, Nature Genetics, 50(9):1335-1341.
• Documentación del pipeline Pan-UK Biobank: Broad Institute, Metodología técnica y pipeline de control de calidad de Pan-UK Biobank (pan.ukbb.broadinstitute.org).
• Artículo original de Lugannani & Rice: Lugannani & Rice (1980), Aproximaciones del punto silla para sumas de variables aleatorias, Advances in Applied Probability, 12(2):475-490.

De este modo, hemos examinado en detalle el núcleo de los pipelines de GWAS, desde la descomposición en valores propios (S42) de los modelos lineales mixtos (LMM) hasta la aceleración del análisis de biobancos a gran escala con 500.000 participantes mediante las técnicas de Ridge por bloques y las aproximaciones del punto silla (S43) implementadas en REGENIE y SAIGE.

En el siguiente capítulo, S44, pasaremos al análisis de mapeo fino (SuSiE y CAVIAR) para seleccionar entre los cientos de variantes asociadas identificadas por GWAS aquellas que son verdaderas variantes causales de la enfermedad.