De la epigenética a la genómica de poblaciones en S41: La lucha contra las señales falsas

Hasta S41 (secuenciación bisulfito), hemos tratado análisis de biología molecular que leen los cambios en la metilación del ADN y el transcriptoma a nivel de célula única. Ahora ampliamos nuestra perspectiva hacia poblaciones de miles o cientos de miles de individuos. El estudio de asociación del genoma completo (GWAS, Genome-Wide Association Study) es una técnica central de la bioinformática que busca variantes genéticas estadísticamente asociadas con una enfermedad o rasgo específico entre millones de SNPs.

Sin embargo, ejecutar un análisis de regresión lineal simple sin más conduce a resultados terribles. Si los linajes (ancestry) del grupo de pacientes y del grupo de control difieren sutilmente, o si hay parientes lejanos mezclados en la muestra, se desbordarán valores P falsamente positivos para decenas de miles de SNPs irrelevantes para la enfermedad. A esto se le denomina problema de estructura poblacional (Population Stratification) y parentesco oculto (Cryptic Relatedness).

En esta sección, confirmaremos por qué falla el modelo de regresión simple, derivaremos los principios matemáticos del modelo mixto lineal (LMM, Linear Mixed Model), que corrige en gran medida este problema, así como el método de cálculo de la matriz de parentesco genético (GRM). Luego, practicaremos el control de calidad (QC) y el análisis de componentes principales (PCA) con PLINK2, y verificaremos directamente en R el efecto de corrección de la inflación Lambda_GC (la ejecución real de GEMMA/GCTA se tratará en la siguiente sección práctica).

Los límites de la regresión simple y derivación de la fórmula del LMM

Colapso del modelo de regresión lineal simple

Dado un conjunto de N muestras y M SNPs, el efecto del SNP $j$-ésimo $x_j \in \{0, 1, 2\}$ sobre el rasgo $y$ puede modelarse mediante una regresión simple como sigue:

$y_i = \mu + x_{ij} \beta_j + e_i, \quad e_i \sim \mathcal{N}(0, \sigma_e^2)$

El supuesto crítico de este modelo es que los residuos $e_i$ son independientes entre todos los individuos ($i.i.d.$). Sin embargo, en las poblaciones humanas reales, debido a compartir ancestros comunes, los genotipos de todo el genoma se parecen sistemáticamente. Cuando se rompe el supuesto de independencia de los residuos, los errores se acumulan e inflan la puntuación z de las estadísticas. El grado de inflación de las estadísticas se mide mediante el coeficiente de inflación genómica $\lambda_{GC}$.

$\lambda_{GC} = \frac{\text{median}(\chi^2_{\text{observed}})}{0.455}$

En una población aleatoria, $\lambda_{GC}$ debería ser 1.0; sin embargo, al ejecutar una regresión simple con mezcla poblacional, $\lambda_{GC}$ se eleva por encima de 1.5.

Introducción del modelo mixto lineal (LMM)

El modelo de efectos mixtos lineales trata el efecto $\beta_j$ del SNP $j$-ésimo como un efecto fijo, y separa el efecto acumulativo del resto de los SNPs genómicos en un efecto aleatorio poligénico $u$.

$y = \mathbf{X} \boldsymbol{\beta} + \mathbf{x}_j \beta_j + \mathbf{u} + \mathbf{e}$

$\mathbf{u} \sim \mathcal{N}(0, \sigma_g^2 \mathbf{K}), \quad \mathbf{e} \sim \mathcal{N}(0, \sigma_e^2 \mathbf{I}_N)$

La matriz de varianza-covarianza completa del rasgo $y$, $\mathbf{V}$, es la siguiente:

$\mathbf{V} = \text{Var}(y) = \sigma_g^2 \mathbf{K} + \sigma_e^2 \mathbf{I}_N$

donde $\sigma_g^2$ es la varianza genética, $\sigma_e^2$ es la varianza ambiental, y $\mathbf{K}$ es la Matriz de Relación Genómica (GRM, por sus siglas en inglés) de tamaño N×N.

Derivación de la fórmula de la Matriz de Relación Genómica (GRM)

Para los datos de SNPs estandarizados de $M$ loci, el genotipo del individuo $i$, $x_{im}$, se estandariza de la siguiente manera:

$z_{im} = \frac{x_{im} - 2p_m}{\sqrt{2p_m(1 - p_m)}}$

El elemento $K_{ik}$ de la matriz GRM $\mathbf{K}$ se define como:

$K_{ik} = \frac{1}{M} \sum_{m=1}^M z_{im} z_{km} = \frac{1}{M} \sum_{m=1}^M \frac{(x_{im} - 2p_m)(x_{km} - 2p_m)}{2p_m(1 - p_m)}$

El elemento diagonal $K_{ii}$ contiene información sobre el coeficiente de endogamia del individuo $i$, mientras que los elementos fuera de la diagonal $K_{ik}$ representan la similitud genética entre el individuo $i$ y el individuo $k$.

Ejemplo de cálculo manual: Obtención directa de la GRM para 3 individuos y 2 SNPs

Calculemos directamente la GRM $\mathbf{K}$ basándonos en los datos de 3 individuos (A, B, C) y 2 SNPs. Para simplificar el cálculo, asumimos que la frecuencia alélica de referencia estimada previamente en una población más grande para ambos SNPs es $p_1 = 0.5$ y $p_2 = 0.5$ (esta no es la frecuencia de la muestra de estos tres individuos; en la práctica, se utilizan frecuencias estimadas a partir de toda la cohorte que ha pasado el control de calidad). En este caso, el denominador de estandarización es $\sqrt{2(0.5)(0.5)} \approx 0.7071$.

Calcule los valores estandarizados $z_{im} = \frac{x_{im} - 1}{0.7071}$ para cada SNP.

• SNP 1: $z_{A1} = -1.4142$, $z_{B1} = 0.0$, $z_{C1} = 1.4142$
• SNP 2: $z_{A2} = -1.4142$, $z_{B2} = 1.4142$, $z_{C2} = 1.4142$

Aplique la fórmula $K_{ik} = \frac{1}{2} (z_{i1}z_{k1} + z_{i2}z_{k2})$.

• $K_{AA} = \frac{1}{2} (2 + 2) = 2.0$, $K_{BB} = \frac{1}{2} (0 + 2) = 1.0$, $K_{CC} = \frac{1}{2} (2 + 2) = 2.0$
• $K_{AB} = \frac{1}{2} (0 - 2) = -1.0$, $K_{AC} = \frac{1}{2} (-2 - 2) = -2.0$, $K_{BC} = \frac{1}{2} (0 + 2) = 1.0$

$\mathbf{K} = \begin{bmatrix} 2.0 & -1.0 & -2.0 \\ -1.0 & 1.0 & 1.0 \\ -2.0 & 1.0 & 2.0 \end{bmatrix}$

De esta manera, la matriz de parentesco definida por los productos internos normalizados entre individuos produce la covarianza entre individuos que presentan correlación.

Cálculo rápido mediante transformación de rotación (pre-blancamiento)

Aplicemos una descomposición en valores propios a la matriz de covarianza $\mathbf{V}$.

$\mathbf{K} = \mathbf{U} \mathbf{S} \mathbf{U}^T$

Al multiplicar ambos lados del modelo por $\mathbf{U}^T$, la covarianza entre individuos entrelazada se convierte en una matriz diagonal, transformando cada componente en observaciones independientes (para lograr un blanqueamiento completo y una normalización de varianza unitaria, sería necesario multiplicar adicionalmente por $(\sigma_g^2\mathbf{S}+\sigma_e^2\mathbf{I})^{-1/2}$, pero con solo la rotación ya se logra suficientemente el objetivo de evitar operaciones iterativas de inversión matricial).

$\mathbf{U}^T y = \mathbf{U}^T \mathbf{X} \boldsymbol{\beta} + \mathbf{U}^T \mathbf{x}_j \beta_j + \mathbf{U}^T \mathbf{u} + \mathbf{U}^T \mathbf{e}$

Dado que la varianza de cada componente en el espacio rotado se simplifica a $\sigma_g^2 s_i + \sigma_e^2$, es posible realizar pruebas mucho más rápidas para millones de SNPs sin necesidad de operaciones iterativas de inversión matricial.

Práctica de GWAS basada en PLINK2 y R

Realicemos la gestión de calidad (QC) y la extracción de covariables de PCA mediante PLINK2, así como la validación de la corrección lambda con LMM en R.

Paso 1: Gestión de calidad de genotipos y extracción de PCA utilizando PLINK2

#!/usr/bin/env bash
# Práctica de control de calidad (QC) de datos PLINK2 y análisis de componentes principales (PCA)
set -e

# 1. Control de calidad de SNP y muestras (MAF > 0.01, p HWE > 1e-6)
plink2 --bfile raw_genotypes \
  --maf 0.01 --hwe 1e-6 --geno 0.02 --mind 0.02 \
  --make-bed --out qc_filtered

# 2. Poda LD y cálculo de los 10 primeros componentes principales (PCA)
plink2 --bfile qc_filtered --indep-pairwise 50 5 0.2 --out ld_pruned
plink2 --bfile qc_filtered --extract ld_pruned.prune.in --pca 10 --out pca_results

Paso 2: Corrección de lambda LMM y visualización del gráfico Q-Q utilizando R

# Análisis de resultados de cálculos LMM en entorno R y cálculo de Lambda_GC
library(ggplot2)

set.seed(42)
n_snps <- 100000
true_lambda <- 1.4  # Factor de inflación observado cuando queda estructura poblacional (valor hipotético)

# Reproducir el fenómeno de "aumento de SNPs aparentemente significativos bajo la hipótesis nula" mediante regresión simple (inflación):
# Reproducir el fenómeno de "aumento de SNPs aparentemente significativos bajo la hipótesis nula" mediante regresión simple (inflación):
chisq_raw <- true_lambda * rchisq(n_snps, df = 1)
p_raw <- pchisq(chisq_raw, df = 1, lower.tail = FALSE)

# Tras corrección LMM: al cumplirse la hipótesis nula, los valores p deben seguir una distribución uniforme.
p_lmm <- runif(n_snps)

calc_lambda <- function(p_vals) {
  chisq <- qchisq(1 - p_vals, df = 1)
  return(median(chisq, na.rm = TRUE) / qchisq(0.5, df = 1))
}

cat(sprintf("Regresión lineal simple Lambda_GC: %.3f\n", calc_lambda(p_raw)))
cat(sprintf("Modelo mixto lineal (LMM) Lambda_GC: %.3f\n", calc_lambda(p_lmm)))

# Visualización del gráfico Q-Q
qq_df <- data.frame(
  observed = -log10(sort(p_lmm)),
  expected = -log10(ppoints(n_snps))
)

ggplot(qq_df, aes(x = expected, y = observed)) +
  geom_point(color = "#2b5c8f", alpha = 0.6, size = 1.2) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  labs(title = "GWAS Q-Q Plot (LMM Corrected)", x = "Expected -log10(P)", y = "Observed -log10(P)") +
  theme_minimal()

En el código anterior, se puede observar que la regresión simple Lambda_GC tiende a acercarse al multiplicador de inflación asumido (1.4), mientras que tras la corrección con LMM, Lambda_GC disminuye hasta cerca de 1.0. En una cohorte real, no existe un umbral absoluto fijo como "pasar si Lambda_GC es menor que X", ya que el rango normal varía según el tamaño de la muestra y la poligenicidad del rasgo; por lo tanto, debe evaluarse cuánto se desvía de 1.0 junto con otros indicadores de control de calidad (residuos de PCA y parentesco).

Mapeo CS

• Blanqueo espacial (Pre-whitening): Técnica que diagonaliza y normaliza la matriz de covarianza de datos multivariados para transformarlos ortogonalmente a un estado de ruido blanco. La descomposición espectral de la matriz de parentesco en LMM es el primer paso para rotar datos no independientes en componentes independientes; al añadir escalado de varianza, se obtiene el pre-whitening completo.
• Descomposición espectral (Eigen-decomposition): Concepto de álgebra lineal que descompone una matriz simétrica $\mathbf{K}$ en el producto de vectores y valores propios $\mathbf{U} \mathbf{S} \mathbf{U}^T$, comprimiendo las operaciones geométricas $O(N^3)$ en operaciones de componentes diagonales $O(N)$.
• Método de máxima verosimilitud restringida (REML): Algoritmo de optimización que estima componentes de varianza corrigiendo la pérdida de grados de libertad debida a la estimación de coeficientes de efectos fijos, reduciendo el sesgo en comparación con la máxima verosimilitud estándar.

Defectos comunes

• No aplicación de LOCO (Leave-One-Chromosome-Out): Si al calcular la GRM no se excluye el cromosoma al que pertenece el SNP candidato, la señal del SNP en prueba es absorbida previamente por los efectos aleatorios, lo que provoca una reducción de la señal verdadera (contaminación proximal).
• No selección de individuos emparentados: El coeficiente KING 0.0884 suele ser el valor umbral que separa relaciones de segundo y tercer grado. Si no se eliminan en el QC previo los pares de parientes más cercanos que este valor y se introducen directamente en el GLM, las estadísticas se distorsionan.
• Omisión de covariables: Si no se incluyen como covariables de efectos fijos el sexo, la edad y los componentes principales superiores (PC1-10), la varianza residual no disminuye adecuadamente.

Para profundizar más

El texto ha sido reescrito directamente por el equipo de investigación del BPD. Para una exploración teórica más profunda, consulte los siguientes recursos:

• StatQuest: Lista de clases de visualización de estadística y genética del profesor Joshua Starmer
  • StatQuest — Modelos lineales y modelos mixtos lineales

• Artículo original de GCTA: Yang et al. (2011), GCTA: A Tool for Genome-wide Complex Trait Analysis, AJHG, 88(1):76-82.
• Artículo original de GEMMA: Zhou & Stephens (2012), Genome-wide efficient mixed-model analysis for association studies, Nature Genetics, 44(7):821-824.
• Documentación oficial de PLINK2: Purcell & Chang, PLINK 2.0 Reference Manual (cog-genomics.org/plink/2.0).

Aprendimos a controlar la estructura poblacional en GWAS a gran escala (miles de individuos) mediante PLINK2 y LMM. Sin embargo, en la era de los big data como UK Biobank, con un tamaño de muestra de 500.000, el método LMM $O(N^3)$ para la descomposición de valores propios de la matriz GRM $N \times N$ con $N=500.000$ se detiene debido a una explosión de memoria.

En el siguiente capítulo S43, pasaremos a los algoritmos de escalado ultrarrápido de REGENIE y SAIGE, que resuelven completamente los datos masivos de biobancos de 500.000 individuos en pocas horas mediante la estrategia divide y vencerás (Ridge Block) y la aproximación del punto de silla (SPA).