집단 유전학 완벽 이해 — 하디-바인베르크 평형과 ADMIXTURE 조상 추정

S45 조직 특이성에서 집단 구조로 — 처음으로 돌아가서 직접 추정해보기

S42에서 우리는 "집단 구조화(Population Stratification)"를 GWAS의 거짓 양성을 만드는 골칫거리로 다뤘고, PCA와 LMM으로 그 영향을 교정하는 법을 배웠습니다. 그런데 정작 "집단 구조 자체를 어떻게 추정하는가"는 다루지 않았습니다. 이 편에서는 한 단계 더 근본으로 내려가, 유전자형 데이터 하나로 개체가 얼마나 조상 집단이 다른지, 그리고 각 SNP가 그 집단 안에서 정상적으로 유전되고 있는지를 직접 검정합니다.

두 가지를 다룹니다. 첫째, 개별 SNP가 무작위 교배 가정을 만족하는지 보는 하디-바인베르크 평형(HWE, Hardy-Weinberg Equilibrium) 검정. 둘째, 개체 전체 유전체가 여러 조상 집단의 혼합으로 어떻게 구성되어 있는지 추정하는 ADMIXTURE의 최대우도 모델입니다.

HWE 검정과 ADMIXTURE 우도 모델 유도

하디-바인베르크 평형: 무작위 교배가 만드는 기대 유전자형 비율

한 집단이 무작위로 교배하고, 자연선택·이주·돌연변이가 없다면, 대립유전자 빈도 $p$ (참조 대립유전자)와 $q=1-p$ (대립 유전자)로부터 유전자형 비율이 다음과 같이 정확히 정해집니다.

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

QC 단계에서 이 비율이 크게 깨진 SNP는 유전형 결정(genotyping) 오류, allele dropout, 표본 오염, 친족 포함, Wahlund 효과(이질적 하위집단을 하나로 합쳐 생기는 인위적 이형접합 결핍), 근친교배, 그리고 드물게는 자연선택의 신호일 수 있습니다 — genotyping 관련 원인이 실무에서는 더 흔한 우선 점검 대상입니다. 카이제곱 적합도 검정(교육용으로는 단순하지만, 기대도수가 작은 희귀변이에는 정확검정(exact test)이 더 적절합니다)으로 관측치와 기대치의 괴리를 정량화합니다.

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

(범주 3개에서 대립유전자 빈도 1개를 데이터로 추정했으므로 자유도는 $3-1-1=1$ 입니다.)

손 계산 예제: 100명 코호트에서 HWE 위반 SNP 찾기

한 SNP에서 AA 50명, Aa 30명, aa 20명이 관측됐습니다. 대립유전자 빈도는 $p = \frac{2(50)+30}{200} = 0.65$ , $q = 0.35$ 입니다.

기대 유전자형 수: $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$

자유도 1의 카이제곱 분포에서 유의수준 0.05 임계값은 3.84입니다. $11.61 > 3.84$ 이므로 이 SNP는 HWE를 유의미하게 벗어났습니다( $p < 0.001$ ) — 이형접합자(Aa)가 기대보다 훨씬 적고 동형접합자가 많은 패턴은 흔히 유전형 결정 오류나 숨은 집단 혼합의 신호입니다.

ADMIXTURE: 조상 비율과 조상 집단별 대립유전자 빈도의 동시 추정

개체 $i$ 가 $K$ 개의 조상 집단으로부터 비율 $q_{i1}, \dots, q_{iK}$ ( $\sum_k q_{ik}=1$ )만큼 혼합되어 있고, 조상 집단 $k$ 에서 SNP $j$ 의 대립유전자 빈도가 $p_{kj}$ 라고 합시다. 개체 $i$ 가 SNP $j$ 에서 갖는 개체별 기대 대립유전자 빈도는 각 조상 성분의 가중 평균입니다.

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

각 조상 성분 내부에서는 HWE가 성립한다고 가정하면, 개체 $i$ 의 SNP $j$ 유전자형이 관측될 우도는 $\pi_{ij}$ 를 대립유전자 빈도로 쓴 HWE 식과 같습니다.

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

전체 로그우도 $\sum_i \sum_j \log P(x_{ij} \mid q_i, p)$ 를 $\{q_{ik}\}$ 와 $\{p_{kj}\}$ 에 대해 최대화하는 것이 ADMIXTURE의 목표입니다. STRUCTURE(Pritchard et al., 2000)가 이를 베이즈 MCMC 샘플링으로 풀었다면, ADMIXTURE(Alexander et al., 2009)는 같은 우도 함수를 블록 좌표 상승법(block relaxation) — $q$ 를 고정하고 $p$ 블록을 순차 이차계획법(sequential quadratic programming)으로 갱신, 다시 $p$ 를 고정하고 $q$ 를 갱신하는 과정을 준-뉴턴(quasi-Newton) 가속과 함께 반복 — 으로 직접 최적화해 수십～수백 배 빠르게 수렴시킵니다. 잠재 조상 배정을 반복적으로 다루는 생성 모델이라는 점에서 EM과 직관을 공유하지만, ADMIXTURE 자체의 최적화 알고리즘은 EM이 아니라는 점에 주의해야 합니다.

손 계산 예제: 두 조상 집단, 한 SNP에서의 유전자형 우도 비교

$K=2$ 조상 집단의 한 SNP 대립유전자 빈도가 $p_1 = 0.8$ (집단 1), $p_2 = 0.2$ (집단 2)라고 합시다. 이형접합자( $x=1$ )로 관측된 두 개체 A, B의 조상 비율이 각각 $q_A = (0.5, 0.5)$ , $q_B = (0.9, 0.1)$ 일 때,

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

같은 이형접합 관측치라도 개체의 조상 비율에 따라 우도가 달라집니다. ADMIXTURE는 수십만 개의 SNP에 걸쳐 이런 우도를 모두 곱한(로그 합산한) 값을 최대화하는 $q_i$ 를, 앞서 설명한 블록 좌표 상승법 + 준-뉴턴 가속으로 찾습니다.

PLINK2 + R 실습

bash

#!/usr/bin/env bash
# HWE 필터링 + ADMIXTURE 실행 파이프라인
set -e

# 1. HWE p-value 1e-6 미만 SNP 제거 (예시 임계값 — 표본 수·연구 설계에 맞게 조정 필요)
# PLINK2의 --hwe는 기본적으로 exact test를 사용함(교육용 카이제곱 근사와 다름)
# 케이스-컨트롤 설계라면 --keep로 컨트롤군만 추출한 뒤 이 필터를 적용해야 함
plink2 --bfile raw_genotypes --hwe 1e-6 --make-bed --out hwe_filtered

# 2. LD 프루닝 후 PLINK 1 바이너리로 변환 (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. K=2부터 K=5까지 조상 집단 수를 바꿔가며 ADMIXTURE 실행 (교차검증으로 최적 K 선택)
for K in 2 3 4 5; do
  admixture --cv admixture_input.bed $K | tee log_K${K}.out
done

# 앞의 두 조상 집단 손 계산 예제를 코드로 재현
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("개체 A (q=0.5,0.5) 이형접합 우도:", genotype_likelihood(c(0.5, 0.5), p1, p2, 1), "\n")
cat("개체 B (q=0.9,0.1) 이형접합 우도:", genotype_likelihood(c(0.9, 0.1), p1, p2, 1), "\n")

# 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("HWE 카이제곱 통계량: %.2f (p-value: %.2e)\n", chisq, pchisq(chisq, df = 1, lower.tail = FALSE)))

CS 매핑

블록 좌표 상승법(Block Coordinate Ascent): ADMIXTURE가 실제로 쓰는 최적화 방식입니다. 전체 변수를 한꺼번에 최적화하는 대신 $q$ 를 고정하고 $p$ 블록을 순차 이차계획법으로 최적화, 다시 $p$ 를 고정하고 $q$ 를 최적화하는 식으로 번갈아 최적화하는 것은 대규모 최적화 문제에서 흔히 쓰는 좌표 하강/상승법이며, 여기에 준-뉴턴 가속을 더해 STRUCTURE의 MCMC 샘플링보다 훨씬 빠르게 수렴합니다.
잠재 변수 모델(Latent Variable Model): 관측되지 않은 잠재 변수(개체의 조상 성분 배정)를 가정하고 그로부터 관측 데이터의 생성 과정을 모델링하는 것은 k-평균 클러스터링·가우시안 혼합 모델(GMM)과 같은 구조입니다. 다만 GMM이 흔히 EM으로 학습되는 것과 달리, ADMIXTURE는 같은 종류의 잠재변수 모델을 EM이 아닌 직접 수치최적화로 학습한다는 차이가 있습니다.
혼합 모델(Mixture Model): ADMIXTURE의 $\pi_{ij} = \sum_k q_{ik}p_{kj}$ 는 각 관측치가 여러 잠재 성분의 가중 합으로 생성된다고 보는 혼합 모델의 전형적 구조입니다.

자주 만나는 결함

HWE 필터를 전체 코호트에 일괄 적용: 케이스-컨트롤 설계에서는 컨트롤군에서만 HWE를 검정해야 합니다. 케이스군의 HWE 위반은 오히려 질병 연관 신호일 수 있어, 여기에 필터를 걸면 진짜 신호를 지워버립니다.
LD 프루닝 생략 후 ADMIXTURE 실행: 서로 강하게 연관된 SNP를 그대로 넣으면 소수의 유전체 영역(특히 장거리 LD 블록·역위 영역)이 조상 추정 결과를 과도하게 좌우할 수 있습니다. ADMIXTURE는 marker 간 독립을 가정하므로 LD 프루닝이 표준적으로 강하게 권장되며, 프루닝 강도에 따라 결과가 얼마나 민감한지 확인하는 것이 안전합니다.
최적 K를 임의로 고정: "조상 집단 수 K"는 생물학적으로 정해진 참값이 없습니다. 교차검증 오차(CV error)가 가장 낮은 K를 데이터 기반으로 고르되, 여러 K에서의 결과를 함께 보고 해석하는 것이 안전합니다.

더 깊게 파고 싶다면

본문은 BPD 연구진이 직접 재구성한 서술입니다. 원 논문과 공식 문서로 심화해봅시다.

ADMIXTURE 원 논문: Alexander, Novembre & Lange (2009), Fast model-based estimation of ancestry in unrelated individuals, Genome Research, 19:1655–1664.
STRUCTURE 원 논문: Pritchard, Stephens & Donnelly (2000), Inference of Population Structure Using Multilocus Genotype Data, Genetics, 155(2):945–959.
PLINK 원 논문: Purcell et al. (2007), PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses, American Journal of Human Genetics, 81(3):559–575.
PLINK2 공식 문서: cog-genomics.org/plink/2.0.

집단 구조를 직접 추정하는 법을 배웠습니다. 이제 인간 유전체에서 가장 다형적(polymorphic)이면서도 임상적으로 중요한 영역 — 면역계가 자기/비자기를 구별하는 데 쓰이는 HLA 유전자로 넘어갑니다.

다음 편 S47에서는 HLA 유형 결정과 임퓨테이션의 원리를 다룹니다.

BioPlayground

집단 유전학 기본기: 하디-바인베르크 평형과 ADMIXTURE 조상 추정