BioPlayground

🧬
목록으로

REGENIE와 SAIGE: 50만 명의 UK Biobank를 삼키는 대규모 초고속 GWAS

50만 명 수준의 바이오뱅크 빅데이터에서 O(N^3) LMM 병목을 극복합니다. REGENIE의 2단계 Ridge 블록 분할과 SAIGE의 안장점 근사(SPA) 원리를 각자의 방식으로 살펴보고 실습합니다.

심화
|
25분
|
검증 완료 (2026-07-29)
GWASREGENIESAIGEsaddlepoint approximationcase-control imbalance
진행률0/86 (0%)

S42 초고속 GWAS의 장벽 — 50만 명 빅데이터의 계산 폭발

S42(PLINK2와 선형 혼합 모델)에서 우리는 유전적 친족 행렬(GRM)과 고유값 분해를 이용해 집단 구조화가 일으키는 거짓 양성을 제어하는 선형 혼합 모델(LMM)을 다루었습니다. 수천 명 규모의 코호트에서는 GEMMA나 GCTA 같은 전통적 LMM 도구가 잘 작동합니다.

하지만 UK Biobank, FinnGen 등 수십만 명 코호트 빅데이터 시대가 열리면서 기존 LMM은 정면으로 계산의 벽에 부딪혔습니다. 표본 수가 N=500,000N = 500,000명으로 늘어나면, N×NN \times N 크기의 GRM 행렬 메모리 용량만 500,000×500,000×8 bytes≈2 TB500,000 \times 500,000 \times 8 \text{ bytes} \approx 2 \text{ TB}에 달합니다. 2TB 메모리가 아니면 행렬을 올릴 수 없으며, 고유값 분해 연산 복잡도 O(N3)O(N^3)은 계산 시간을 수 개월 단위로 늘려버립니다.

또한 환자군 대 대조군 비율이 1:100 이하로 불균형한 희귀 질환 분석에서는 기존 Wald 검정의 점근 정규성 가정이 흔들리기 시작해, 거짓 양성 P-값이 크게 부풀려질 수 있습니다. 이 편에서는 이러한 빅데이터 한계를 극복하기 위해 등장한 REGENIE의 2단계 Ridge 블록 학습과, SAIGE가 독자적으로 채택한 안장점 근사(SPA, Saddlepoint Approximation) 원리를 각각 살펴봅시다. 두 도구는 같은 문제(대규모·불균형 데이터)를 겨냥하지만 보정 방식은 서로 다른 별개의 알고리즘입니다.

REGENIE와 SAIGE의 핵심 알고리즘 유도

1. REGENIE: 2단계 Ridge 회귀와 블록 분할(Block Binning)

REGENIE는 N×NN \times N GRM 행렬 전체를 올리지 않고, 전게놈 다유전자 효과를 블록 단위 분할 Ridge 회귀로 모형화합니다.

1단계 (Step 1): 전게놈 Ridge 예측자 산출

전게놈 품질 관리를 통과한 MQCM_{QC}개의 SNP를 연속된 1,000개 단위의 블록 B1,…,BKB_1, \dots, B_K로 분할합니다. 블록 BkB_k 내에서 여러 정규화 파라미터 λl\lambda_l에 대해 Ridge 회귀 모델을 학습시켜 국소 예측자 w^k,l\hat{\mathbf{w}}_{k, l}를 구합니다.

y^k,l=XBk(XBkTXBk+λlI)−1XBkTy\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

여기서 y^k,l\hat{\mathbf{y}}_{k, l}은 블록 BkB_k·정규화 파라미터 λl\lambda_l 조합으로 만든 N차원 예측값(fitted predictor)이며, 괄호 안 (XBkTXBk+λlI)−1XBkTy(\mathbf{X}_{B_k}^T \mathbf{X}_{B_k} + \lambda_l \mathbf{I})^{-1} \mathbf{X}_{B_k}^T y 부분이 Ridge 회귀 계수입니다. 축적된 블록 예측값들에 Level 2 Ridge 회귀를 적용하여 최종 다유전자 예측값 y^poly\hat{y}_{\text{poly}}를 구성합니다. 한 번에 메모리에 올려야 하는 양은 전체 SNP MQCM_{QC}가 아니라 블록 크기(기본 1,000개) 단위 O(N⋅bsize)O(N \cdot \text{bsize})로 줄어듭니다.

2단계 (Step 2): LOCO 적용 단일 SNP 검정

각 SNP jj에 대해 해당 염색체를 뺀 LOCO(Leave-One-Chromosome-Out) 예측자 y^poly,−c\hat{y}_{\text{poly}, -c}를 공변량으로 투입하여 단일 회귀를 수행합니다.

y=xjβj+y^poly,−cγ+Wα+ϵy = \mathbf{x}_j \beta_j + \hat{y}_{\text{poly}, -c} \gamma + \mathbf{W} \boldsymbol{\alpha} + \epsilon

2. SAIGE: 안장점 근사(Saddlepoint Approximation, SPA) 유도

희귀 질환이나 극단적 대조군 비율(예: Case 500명 vs Control 499,500명) 데이터에서 통계량 S=xT(y−μ^)S = \mathbf{x}^T (y - \hat{\mu})의 정규분포 가정이 깨지는 문제를 SAIGE는 **안장점 근사(SPA)**로 해결합니다.

아래 CGF는 원리를 보여주기 위해 공변량과 혼합모형 보정을 뺀 단순화한 형태입니다. 실제 SAIGE는 혼합모형으로 보정된 genotype residual과 variance-ratio 절차까지 포함해 CGF를 구성합니다. CGF KS(t)=ln⁡E[etS]K_S(t) = \ln \mathbb{E}[e^{t S}]와 그 1차, 2차 도함수는 다음과 같습니다.

KS(t)=∑i=1Nln⁡(1−μ^i+μ^ietxi)−t∑i=1Nxiμ^iK_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

KS′(t)=∑i=1Nxiμ^ietxi1−μ^i+μ^ietxi−∑i=1Nxiμ^i,KS′′(t)=∑i=1Nxi2μ^i(1−μ^i)etxi(1−μ^i+μ^ietxi)2K_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}

관측된 통계량 값 s에 대해 안장점 방정식 KS′(ζ^)=sK_S'(\hat{\zeta}) = s의 해 ζ^\hat{\zeta}를 구하면, Lugannani-Rice 공식을 통해 정확한 꼬리 확률 P(S≥s)P(S \ge s)를 근사합니다.

P(S≥s)≈1−Φ(w^)+ϕ(w^)(1w^−1u^)P(S \ge s) \approx 1 - \Phi(\hat{w}) + \phi(\hat{w}) \left( \frac{1}{\hat{w}} - \frac{1}{\hat{u}} \right)

w^=sgn(ζ^)2(ζ^s−KS(ζ^)),u^=ζ^KS′′(ζ^)\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})}

작은 손 계산 예제: 극단적 불균형 데이터에서 정규 근사가 틀리는 이유

정규 근사가 얼마나 크게 틀릴 수 있는지, 직접 전수 계산이 가능한 아주 작은 예제로 확인해봅시다. 희귀 변이 보유자 5명(xi=1x_i=1)만 골라낸 상황을 가정합니다. 귀무가설 하에서 각자의 발병 확률은 μ^i=0.05\hat{\mu}_i = 0.05로 동일하고 서로 독립입니다. 스코어 통계량은 S=∑i(yi−μ^i)S = \sum_i (y_i - \hat{\mu}_i)이며, 관측값은 이 5명 중 3명이 실제 환자인 경우 s=3−5×0.05=2.75s = 3 - 5 \times 0.05 = 2.75입니다.

발병자 수 KK는 이항분포 K∼Binomial(5,0.05)K \sim \text{Binomial}(5, 0.05)를 따르므로, 관측값 이상이 나올 확률은 전수 계산으로 정확히 구할 수 있습니다.

P(K≥3)=∑k=35(5k)0.05k 0.955−k≈0.001128+0.0000297+0.0000003≈1.16×10−3P(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}

반면 같은 상황에 정규 근사를 적용하면, Var(S)=5×0.05×0.95=0.2375\text{Var}(S) = 5 \times 0.05 \times 0.95 = 0.2375이므로

Z=2.750.2375≈5.64,P(Z≥5.64)≈8.6×10−9Z = \frac{2.75}{\sqrt{0.2375}} \approx 5.64, \quad P(Z \ge 5.64) \approx 8.6 \times 10^{-9}

정확한 확률(약 1.16×10−31.16 \times 10^{-3})과 정규 근사(약 8.6×10−98.6 \times 10^{-9}) 사이에는 약 13만 배의 차이가 있습니다. 정규 근사를 그대로 GWAS p-value로 채택하면, 실제로는 우연히 나올 수 있는 수준의 신호를 터무니없이 강한 유의성으로 잘못 보고하게 됩니다.

안장점 근사(SPA)는 이 왜곡을 줄이기 위한 방법입니다. 다만 앞의 CGF 방정식 KS′(ζ^)=sK_S'(\hat{\zeta}) = s는 일반적으로 닫힌 형태 해가 없어, 실제로는 수치적으로 근을 찾아야 합니다(SAIGE도 내부적으로 같은 방식을 씁니다). 이 예제에서 근은 ζ^=ln⁡(28.5)≈3.3499\hat{\zeta} = \ln(28.5) \approx 3.3499이고, 이를 Lugannani-Rice 식에 대입하면 SPA 근사값은 약 3.9×10−43.9 \times 10^{-4}입니다. 전수 계산값(1.16×10−31.16 \times 10^{-3})보다는 여전히 3배가량 작지만, 정규 근사(약 8.6×10−98.6 \times 10^{-9})와는 비교가 안 될 만큼 훨씬 가깝습니다. 표본이 5명뿐인 극단적 상황이라 격자(이산) 분포 특유의 오차가 남아있는 것이며, 표본이 커질수록 SPA는 전수 계산에 더 가까워집니다. 아래 코드를 직접 실행해 세 값(전수 계산·정규 근사·SPA)을 비교해봅시다.

REGENIE & R SPA 실습

REGENIE 2단계 실행 파이프라인과 R에서의 안장점 근사(SPA) 구현 코드입니다.

1단계: REGENIE CLI 파이프라인 (Step 1 & Step 2)

bash
#!/usr/bin/env bash
# REGENIE 2단계 초고속 GWAS 파이프라인
set -e


# [Step 1] 전게놈 Ridge 회귀로 다유전자 예측자 산출 (--bsize 1000)
regenie \
  --step 1 --bed ukb_qc_array \
  --phenoFile phenotypes.txt --covarFile covariates.txt \
  --bsize 1000 --bt --lowmem --out regenie_step1_out


# [Step 2] Imputed 변이 대상 단일 SNP 연산 + REGENIE 자체의 근사 Firth 보정
# (--firth --approx: SAIGE의 SPA와는 별개인 REGENIE 고유의 불균형 보정 옵션)
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

2단계: R 기반 안장점 근사(SPA) 보정 구현

r
# R 환경에서 안장점 근사(SPA) 함수 구현 — 본문의 5명 손 계산 예제 재현
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)

# K1(t)-s_obs는 t에 대해 단조증가이므로, 무제약 뉴턴-랩슨은 초기값에 따라
# 쉽게 발산합니다(이 예제에서도 t=0 시작 시 NaN으로 발산 확인). 실무에서는
# 구간을 좁혀가는 안전한 이분법(bisection)을 씁니다.
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
}

# 본문 예제: 5명(x=1)이 모두 mu=0.05, 관측 s = 3 - 5*0.05 = 2.75
mu <- rep(0.05, 5); x <- rep(1, 5); s_obs <- 2.75

# 1) 전수 계산 (참값)
exact_p <- sum(dbinom(3:5, size = 5, prob = 0.05))

# 2) 정규 근사 (Wald)
var_s <- sum(x^2 * mu * (1 - mu))
normal_p <- 1 - pnorm(s_obs / sqrt(var_s))

# 3) 안장점 근사 (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("전수 계산 p-value:    %.5e\n", exact_p))
cat(sprintf("정규 근사 p-value:    %.5e\n", normal_p))
cat(sprintf("안장점 근사(SPA) p-value: %.5e\n", spa_p))

CS 매핑

  • 분할 정복(Divide and Conquer): 수십만 개의 SNP를 메모리에 올리는 대신, 1,000개 블록 단위로 나누어 1단계 Ridge 회귀를 수행하고 결합하는 알고리즘 패러다임입니다.
  • 안장점 근사(Saddlepoint Approximation): 복잡한 통계량 분포의 꼬리 적분 확률을 복소평면 상의 안장점 통과 경로로 변환해, 정규 근사보다 훨씬 작은 오차로 고속 계산해내는 정밀 수치 기법입니다(표본이 매우 작으면 격자분포 특유의 잔여 오차가 남습니다).
  • 수치 최적화 및 PCG: SAIGE가 공분산 행렬 연산 시 역행렬 구하기를 피하고, 반복적 선형 방정식 수치 해법으로 메모리 사용량을 O(N)O(N)으로 줄이는 기법입니다.

자주 만나는 결함

  • Step 1에 Imputed SNP 투입: REGENIE Step 1은 통상 수만～십만 개 수준의 QC 통과 Array SNP로 다유전자 예측자를 학습합니다(정확한 상한은 표본 수·메모리 설정에 따라 달라집니다). 1,000만 개 넘는 Imputed SNP 전체를 그대로 투입하면 블록 Ridge 연산 시간과 메모리가 감당하기 어려운 수준으로 늘어납니다.
  • Step 2 LOCO 누락: Step 2 단일 변이 검정 시 해당 염색체를 뺀 LOCO 다유전자 예측자를 쓰지 않으면 자기 자신 신호를 미리 제거해버리는 오류가 발생합니다.
  • Firth 보정 임계값 미설정: Case 수가 적은 변이에 대해 Firth/SPA 보정을 끄고 단순 로지스틱 회귀를 돌리면 P-value가 크게 부풀려질 수 있습니다.

더 깊게 파고 싶다면

본문은 BPD 연구진이 직접 재구성한 서술입니다. 최신 대규모 GWAS 알고리즘 심화를 위해 다음 원문들을 탐구해봅시다.

  • REGENIE 원 논문: Mbatchou et al. (2021), Computationally efficient whole-genome regression for quantitative and binary traits, Nature Genetics, 53(7):1097-1103. (PubMed 34017140)
  • SAIGE 원 논문: Zhou et al. (2018), Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies, Nature Genetics, 50(9):1335-1341.
  • Pan-UK Biobank 파이프라인 문서: Broad Institute, Pan-UK Biobank Technical Methodology & QC pipeline (pan.ukbb.broadinstitute.org).
  • Lugannani & Rice 원 논문: Lugannani & Rice (1980), Saddle Point Approximations for Shares of Sums of Random Variables, Advances in Applied Probability, 12(2):475-490.

이로써 LMM의 고유값 분해(S42)부터 50만 명 대규모 바이오뱅크 분석을 가속화하는 REGENIE와 SAIGE의 블록 Ridge 및 안장점 근사(S43)까지, GWAS 파이프라인의 심장부를 모두 들여다보았습니다.

다음 편 S44에서는 GWAS가 찾아낸 수백 개의 연관 변이 중, 어떤 변이가 진짜 질병을 일으키는 인과 변이(Causal Variant)인지 골라내는 Fine-mapping (SuSiE 및 CAVIAR) 분석으로 넘어갑니다.