S41 후성유전학에서 집단 유전체학으로 — 거짓 신호와의 싸움

S41(바이설파이트 시퀀싱)까지 우리는 단일 세포 수준에서 DNA 메틸화와 전사체 변화를 읽어내는 분자생물학적 분석을 다루었습니다. 이제 시선을 수천 내지 수십만 명의 **집단(population)**으로 확장합니다. 전게놈 연관성 분석(GWAS, Genome-Wide Association Study)은 수백만 개의 SNP 중 특정 질환이나 형질과 통계적으로 연관된 유전 변이를 탐색하는 바이오인포매틱스의 핵심 기법입니다.

하지만 무턱대고 단순 선형 회귀 분석을 돌리면 끔찍한 결과를 맞이합니다. 환자군과 대조군의 혈통(ancestry)이 미세하게 다르거나 표본 내에 사촌 관계가 섞여 있으면, 질병과 무관한 수만 개의 SNP에서 거짓 양성 P-값이 쏟아져 나옵니다. 이를 집단 구조화(Population Stratification) 및 숨은 친족 관계(Cryptic Relatedness) 문제라 부릅니다.

이 편에서는 단순 회귀 모델이 무너지는 이유를 확인하고, 이를 상당 부분 교정하는 **선형 혼합 모델(LMM, Linear Mixed Model)**의 수학적 원리와 유전적 친족 행렬(GRM) 계산법을 유도한 뒤, PLINK2로 QC·PCA를 실습하고 R로 Lambda_GC 인플레이션 교정 효과를 직접 확인해봅시다(실제 GEMMA/GCTA 실행은 다음 실습 편에서 다룹니다).

단순 회귀의 한계와 LMM 수식 유도

단순 선형 회귀 모델의 붕괴

N명의 샘플과 M개의 SNP가 주어졌을 때, j번째 SNP $x_j \in \{0, 1, 2\}$가 형질 $y$에 미치는 영향은 다음과 같이 단순 회귀로 모형화할 수 있습니다.

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

이 모형의 치명적 가정은 잔차 $e_i$가 모든 개체 사이에서 독립이라는 점($i.i.d.$)입니다. 그러나 실제 인간 집단에서는 공통 조상을 공유함에 따라 유전자형 전체가 체계적으로 닮아 있습니다. 잔차가 서로 독립이라는 가정이 깨지면 오차가 축적되어 통계량 z-score가 부풀려집니다. 통계량의 인플레이션 정도는 **유전체 인플레이션 계수 $\lambda_{GC}$**로 측정합니다.

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

무작위 집단에서는 $\lambda_{GC} = 1.0$이어야 하지만, 집단 혼재 시 단순 회귀를 실행하면 $\lambda_{GC}$가 1.5 이상으로 솟구치게 됩니다.

선형 혼합 모델(LMM)의 도입

선형 혼합 모델은 j번째 SNP의 효과 $\beta_j$는 고정 효과(fixed effect)로 두고, 나머지 전게놈 SNP들의 누적 효과를 다유전자 무작위 효과(polygenic random effect) $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)$

형질 $y$의 전체 분산-공분산 행렬 $\mathbf{V}$는 다음과 같습니다.

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

여기서 $\sigma_g^2$는 유전적 분산, $\sigma_e^2$는 환경적 분산이며, $\mathbf{K}$는 N×N 크기의 **유전적 친족 행렬(GRM, Genomic Relationship Matrix)**입니다.

유전적 친족 행렬(GRM) 수식 유도

M개의 표준화된 SNP 데이터에 대해, 개체 i의 유전자형 $x_{im}$을 다음과 같이 표준화합니다.

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

GRM 행렬 $\mathbf{K}$의 원소 $K_{ik}$는 다음과 같이 정의됩니다.

$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)}$

대각선 원소 $K_{ii}$는 개체 i의 근친교배 계수 정보를 담고 있으며, 비대각선 원소 $K_{ik}$는 개체 i와 k 사이의 유전적 유사도를 나타냅니다.

손 계산 예제: 3명 개체와 2개 SNP의 GRM 직접 구하기

3명의 개체(A, B, C)와 2개의 SNP 데이터를 바탕으로 GRM $\mathbf{K}$를 직접 계산해봅시다. 계산을 단순하게 유지하기 위해, 두 SNP 모두 더 큰 참조 집단에서 미리 추정된 대립유전자 빈도가 $p_1 = 0.5$, $p_2 = 0.5$라고 가정합니다(이 3명만의 표본 빈도가 아닙니다 — 실무에서는 QC를 통과한 전체 코호트에서 추정한 빈도를 씁니다). 이때 표준화 분모는 $\sqrt{2(0.5)(0.5)} \approx 0.7071$입니다.

각 SNP별 표준화 값 $z_{im} = \frac{x_{im} - 1}{0.7071}$을 산출합니다.

• 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$

$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}$

이처럼 개체 간 정규화 내적으로 정의된 친족 행렬은 상관관계가 존재하는 개체 간 공분산을 산출해냅니다.

회전 변환(Pre-whitening)을 통한 고속 연산

공분산 행렬 $\mathbf{V}$에 고유값 분해(Eigen-decomposition)를 적용해봅시다.

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

모델 양변에 $\mathbf{U}^T$를 곱하면, 얽혀 있던 개체 간 공분산이 대각행렬로 바뀌면서 각 성분이 서로 독립인 관측값으로 변환됩니다(완전한 백색화·단위분산 정규화까지 하려면 여기에 $(\sigma_g^2\mathbf{S}+\sigma_e^2\mathbf{I})^{-1/2}$ 스케일링을 한 번 더 곱해야 하지만, 회전만으로도 반복 역행렬 연산을 피하는 목적은 충분히 달성됩니다).

$\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}$

회전된 공간에서 각 성분 분산은 $\sigma_g^2 s_i + \sigma_e^2$로 단순 대각화되므로, 수백만 개의 SNP를 반복 역행렬 연산 없이 훨씬 빠르게 검정할 수 있습니다.

PLINK2 및 R 기반 GWAS 실습

PLINK2를 이용한 품질 관리 및 PCA 공변량 추출, 그리고 R 기반 LMM 람다 교정 검증을 수행해봅시다.

1단계: PLINK2를 활용한 유전자형 QC 및 PCA 추출

#!/usr/bin/env bash
# PLINK2 데이터 QC 및 주성분 분석(PCA) 실습
set -e

# 1. SNP 및 샘플 품질 관리 (MAF > 0.01, HWE p > 1e-6)
plink2 --bfile raw_genotypes \
  --maf 0.01 --hwe 1e-6 --geno 0.02 --mind 0.02 \
  --make-bed --out qc_filtered

# 2. LD 프루닝 및 상위 10개 주성분(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

2단계: R을 활용한 LMM 람다 교정 및 Q-Q Plot 시각화

# R 환경에서 LMM 연산 결과 분석 및 Lambda_GC 계산
library(ggplot2)

set.seed(42)
n_snps <- 100000
true_lambda <- 1.4  # 집단 구조가 남아있을 때 관측되는 인플레이션 배율(가정값)

# 단순 회귀(인플레이션): 카이제곱 통계량 자체를 true_lambda배 부풀려
# "귀무가설 하인데도 유의미해 보이는 SNP가 늘어나는" 현상을 재현합니다.
chisq_raw <- true_lambda * rchisq(n_snps, df = 1)
p_raw <- pchisq(chisq_raw, df = 1, lower.tail = FALSE)

# LMM 교정 후: 귀무가설이 성립하므로 p-value는 균등분포를 따라야 합니다.
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("단순 선형 회귀 Lambda_GC: %.3f\n", calc_lambda(p_raw)))
cat(sprintf("선형 혼합 모델(LMM) Lambda_GC: %.3f\n", calc_lambda(p_lmm)))

# Q-Q plot 시각화
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()

위 코드에서 단순 회귀 쪽 Lambda_GC는 가정한 인플레이션 배율(1.4) 부근으로 나오고, LMM 교정 후 Lambda_GC는 1.0 근처로 낮아지는 것을 확인할 수 있습니다. 실제 코호트에서 "Lambda_GC가 몇 이하면 통과"라는 고정된 절대 기준은 없습니다 — 표본 크기와 형질의 다유전자성(polygenicity)에 따라 정상 범위 자체가 달라지므로, 1.0 대비 얼마나 벗어났는지를 다른 QC 지표(친족·PCA 잔차)와 함께 판단해야 합니다.

CS 매핑

• 공간 백색화(Pre-whitening): 다변량 데이터의 공분산 행렬을 대각화·정규화하여 무색 신호(white noise) 상태로 직교 변환하는 기법입니다. LMM의 친족 행렬 고유값 분해는 비독립 데이터를 독립 성분으로 회전시키는 첫 단계로, 여기에 분산 스케일링을 더하면 완전한 프리화이트닝이 됩니다.
• 대각화(Eigen-decomposition): 대칭 행렬 $\mathbf{K}$를 고유벡터와 고유값의 곱 $\mathbf{U} \mathbf{S} \mathbf{U}^T$로 분해하여 $O(N^3)$ 기하학적 연산을 대각성분 $O(N)$ 연산으로 압축하는 선형대수학 개념입니다.
• 제한 최대우도법(REML): 고정 효과 계수 추정에 의한 자유도 손실을 보정해, 일반 최대우도법보다 편향을 줄여 변량 성분을 추정하는 최적화 알고리즘입니다.

자주 만나는 결함

• LOCO(Leave-One-Chromosome-Out) 미적용: GRM 산출 시 검정 대상 SNP가 속한 염색체를 빼지 않으면, 검정 SNP 신호를 무작위 효과가 미리 흡수해 진짜 신호가 깎여 나가는 현상(proximal contamination)이 발생합니다.
• 근친 개체 미선별: KING coefficient 0.0884는 통상 2촌(second-degree)과 3촌(third-degree) 관계를 가르는 경계값입니다. 이 값 이상으로 가까운 친족 쌍을 사전 QC에서 제거하지 않고 그대로 GLM에 넣으면 통계량이 왜곡됩니다.
• 공변량 누락: 성별, 연령, PCA 상위 성분(PC1-10)을 고정 효과 공변량으로 투입하지 않으면 잔차 분산이 줄어들지 않습니다.

더 깊게 파고 싶다면

본문은 BPD 연구진이 직접 재구성한 서술입니다. 더욱 깊은 이론적 탐구를 위해 다음 자료를 참고합시다.

• StatQuest: Joshua Starmer 교수의 통계/유전학 시각화 강의 목록
  • StatQuest — Linear Models and Linear Mixed Models
• GCTA 원 논문: Yang et al. (2011), GCTA: A Tool for Genome-wide Complex Trait Analysis, AJHG, 88(1):76-82.
• GEMMA 원 논문: Zhou & Stephens (2012), Genome-wide efficient mixed-model analysis for association studies, Nature Genetics, 44(7):821-824.
• PLINK2 공식 문서: Purcell & Chang, PLINK 2.0 Reference Manual (cog-genomics.org/plink/2.0).

PLINK2와 LMM을 통해 수천 명 규모의 GWAS에서 집단 구조 문제를 제어하는 법을 배웠습니다. 하지만 표본 수가 50만 명에 달하는 UK Biobank 빅데이터 시대에는 $N=500,000$의 $N \times N$ GRM 행렬을 고유값 분해하는 $O(N^3)$ LMM은 메모리 폭발로 멈춰 서게 됩니다.

다음 편 S43에서는 50만 명 거대 바이오뱅크 데이터를 분할 정복(Ridge Block) 및 안장점 근사(SPA)로 수시간 만에 완파하는 REGENIE와 SAIGE의 초고속 확장 알고리즘으로 넘어가겠습니다.