S44 인과 변이는 찾았는데, 그 다음은?

S44(Fine-mapping)에서 우리는 GWAS 신호 뭉치 안에서 통계적으로 가장 그럴듯한 인과 변이를 credible set으로 좁혔습니다. 하지만 "이 SNP가 인과 변이다"라는 결론만으로는 부족합니다. 그 변이가 구체적으로 무엇을 바꿔서 질병 위험에 영향을 미치는지 알아야 합니다.

인간 유전 변이의 상당수는 단백질 서열을 바꾸는 것이 아니라, 프로모터·인핸서 같은 조절 영역에 위치해 유전자가 얼마나 많이 발현되는지를 바꿉니다. 이 연결고리를 정량화하는 것이 eQTL(expression Quantitative Trait Loci, 발현 수량형질좌위) 매핑입니다. 그런데 유전자 조절은 조직마다 다르게 작동합니다 — 같은 변이라도 간에서는 발현을 올리고 뇌에서는 아무 영향이 없을 수 있습니다. GTEx(Genotype-Tissue Expression) 프로젝트는 49개 조직에서 이 조직 특이성을 대규모로 지도화했습니다.

이 편에서는 cis-eQTL 회귀 모델의 원리, 조직 간 신호 공유를 측정하는 법, 그리고 GWAS 신호와 eQTL 신호가 "같은 인과 변이를 공유하는지" 검정하는 **colocalization(공동정위)**의 베이즈 통계를 유도합니다.

cis-eQTL 회귀 모델과 조직 간 공동추론

기본 모델: 유전자형이 발현량을 설명하는 선형 회귀

유전자 $g$ 근처(보통 TSS ± 1Mb, cis-window)의 SNP $j$에 대해, 개체 $i$의 정규화된 발현량 $y_{ig}$를 유전자형 $x_{ij} \in \{0,1,2\}$로 설명하는 회귀식은 다음과 같습니다.

$y_{ig} = \mu + x_{ij}\beta_{jg} + \mathbf{C}_i \boldsymbol{\alpha} + \epsilon_i$

여기서 $\mathbf{C}_i$는 성별, 유전체 PCA, 그리고 PEER factor(발현량에 숨어있는 배치 효과·기술적 변동뿐 아니라 세포 조성·생물학적 상태 같은 요인까지 데이터 기반으로 함께 추출되는 잠재 공변량)를 포함하는 공변량 벡터입니다. 발현량 $y_{ig}$는 원본 read count가 아니라, 한 유전자 안에서 개체 간 순위를 표준정규분포에 맞추는 **순위 기반 역정규변환(rank-based inverse normal transform, INT)**을 거친 값을 씁니다 — 발현량 분포의 왜도(skewness)·이상치가 회귀의 정규성 가정을 깨뜨리지 않도록 개체 간 분포를 정규화하는 것이 목적이며, 이 값 자체로 서로 다른 유전자의 절대 발현 수준을 직접 비교할 수 있게 되는 것은 아닙니다.

계층적 다중검정: gene-SNP 쌍의 폭발을 다스리기

유전자 하나당 cis-window 안에 SNP가 수천 개씩 있고, 유전자는 2만 개가 넘습니다. 모든 gene-SNP 쌍에 단순 Bonferroni를 적용하면 지나치게 보수적이라 진짜 신호도 놓칩니다. GTEx가 쓰는 방식은 2단계입니다.

1. 유전자 내부 보정: 각 유전자마다 cis-window 안 최고 SNP의 p-value에 대해, SNP 간 LD를 반영한 순열검정(permutation) 또는 그 근사인 beta-분포 근사로 유전자 수준 p-value $p_{g}$를 구합니다.
2. 유전자 간 보정: $p_g$들에 Storey의 q-value(FDR) 방법을 적용해, 원하는 FDR 수준(보통 5%)에서 유의미한 eGene(eQTL을 가진 유전자)을 정합니다.

조직 간 신호 공유: $\pi_1$ 통계량

조직 A에서 유의미한 eQTL이 조직 B에서도 진짜 신호인지(단지 검정력 부족으로 안 보이는 것인지) 측정하려면, 조직 B의 해당 p-value 분포에서 귀무가설이 아닌 비율 $\pi_1 = 1 - \pi_0$을 추정합니다($\pi_0$은 Storey 방법으로 추정하는 "진짜 귀무가설 비율"). $\pi_1$이 1에 가까울수록 조직 A의 신호 대부분이 조직 B에서도 비귀무(non-null)로 나타난다는 뜻이고, 낮으면 조직 특이적 조절일 가능성이 큽니다. 다만 $\pi_1$이 높다고 해서 두 조직에서 정확히 같은 인과 SNP·같은 효과 방향까지 보장되는 것은 아니며, 이는 어디까지나 "이 신호가 우연이 아니다"의 공유 정도를 요약한 값입니다. GTEx 실측에서 근육·피부처럼 유사한 조직 쌍은 $\pi_1$이 높고, 뇌와 혈액처럼 이질적인 조직 쌍은 낮게 나타나는 경향이 있습니다.

Colocalization: GWAS 신호와 eQTL 신호가 "같은 원인"을 공유하는가

S44의 베이즈 인수(ABF) 프레임워크를 그대로 확장하면, 두 개의 서로 다른 형질(GWAS 형질과 유전자 발현량)이 같은 인과 변이를 공유하는지 검정할 수 있습니다. coloc(Giambartolomei et al., 2014)은 locus 안 각 SNP $j$에 대해 형질 1의 ABF($a_j$)와 형질 2의 ABF($b_j$)를 각각 구한 뒤, 5가지 배타적 가설에 사후확률을 매깁니다.

$\text{Sum}_1 = \sum_j a_j, \quad \text{Sum}_2 = \sum_j b_j, \quad \text{Sum}_4 = \sum_j a_j b_j, \quad \text{Sum}_3 = \text{Sum}_1 \text{Sum}_2 - \text{Sum}_4$

사전확률 $p_1, p_2, p_{12}$(SNP 하나가 각각 형질1/형질2/두 형질 모두의 인과 변이일 사전확률, coloc 기본값은 $p_1=p_2=10^{-4}$, $p_{12}=10^{-5}$)를 곱해 정규화하면 $\text{PP.H4}$(공유 인과 변이 사후확률)를 얻습니다.

손 계산 예제: 두 SNP locus에서 PP.H4 구하기

아래 숫자는 실제 데이터가 아니라 계산 과정을 보여주기 위한 단순 산술 예시입니다. 2개 SNP locus에서 형질 1(GWAS)의 ABF가 $a = [366.3, 366.3]$(두 SNP의 GWAS 신호가 통계적으로 구분되지 않는 상황), 형질 2(eQTL)의 ABF가 $b = [500, 2]$(eQTL 쪽에서는 SNP1을 뚜렷하게 지목하는 상황)라고 합시다. 실제 coloc 분석에서는 두 형질의 SNP 집합·LD·ancestry가 서로 맞춰져 있는지(harmonization)를 먼저 확인해야 하며, 같은 locus인데 한쪽만 완전 LD로 보이고 다른 쪽은 아닌 경우는 변이 조합이나 QC 문제일 수 있어 그 자체로 별도 점검 대상입니다.

$\text{Sum}_1 = 732.6, \quad \text{Sum}_2 = 502, \quad \text{Sum}_4 = 366.3 \times 500 + 366.3 \times 2 = 183{,}882.6$

$\text{Sum}_3 = 732.6 \times 502 - 183{,}882.6 = 183{,}882.6$

기본 사전확률로 각 가설의 비정규화 값을 구하면:

$H0{:}~1,\quad H1{:}~10^{-4}\times732.6=0.0733,\quad H2{:}~10^{-4}\times502=0.0502$ $H3{:}~10^{-8}\times183{,}882.6=0.00184,\quad H4{:}~10^{-5}\times183{,}882.6=1.8388$

총합 $2.9641$로 나누면 $\text{PP.H4} \approx 0.62$, $\text{PP.H0} \approx 0.34$입니다. Sum3와 Sum4가 우연히 같은 값인데도 $\text{PP.H4}$가 $\text{PP.H3}$보다 압도적으로 큰 이유는, 사전확률 자체가 $p_{12}(10^{-5}) \gg p_1 p_2(10^{-8})$이기 때문입니다 — "같은 SNP가 두 형질 모두의 원인일 사전 확률"이 "서로 다른 두 SNP가 각각 원인일 사전 확률"보다 애초에 훨씬 높게 설정되어 있는 것입니다. 이 예제에서는 PP.H4가 0.62로 공유 인과 변이 쪽에 무게가 실리지만, 0.34의 PP.H0도 남아있어 확정적 결론은 아닙니다.

R 실습: coloc 방식 PP.H4 계산

# 앞의 두 SNP 손 계산 예제를 코드로 재현
a <- c(366.3, 366.3)  # 형질 1(GWAS) ABF
b <- c(500, 2)        # 형질 2(eQTL) ABF
p1 <- 1e-4; p2 <- 1e-4; p12 <- 1e-5

sum1 <- sum(a); sum2 <- sum(b); sum4 <- sum(a * b); sum3 <- sum1 * sum2 - sum4

unnorm <- c(
  H0 = 1,
  H1 = p1 * sum1,
  H2 = p2 * sum2,
  H3 = p1 * p2 * sum3,
  H4 = p12 * sum4
)

pp <- unnorm / sum(unnorm)
print(round(pp, 4))

CS 매핑

• 베이즈 공동추론(Colocalization): 두 개의 독립적인 증거(GWAS·eQTL)가 하나의 잠재 원인을 공유하는지 검정하는 것은, 서로 다른 센서 로그가 같은 근본 이벤트에서 나온 것인지 판별하는 데이터 융합(sensor fusion) 문제와 동일한 구조입니다.
• 잠재 변수 추출(PEER factor): 관측 데이터에서 알려지지 않은 배치 효과·기술적 변동을 비지도 학습으로 뽑아내는 것은 PCA·인자분석(factor analysis)이 잠재 변수를 추정하는 것과 같은 원리입니다.
• 계층적 다중검정: 유전자 내부 보정 후 유전자 간 FDR을 적용하는 2단계 구조는, 트리 구조 데이터에서 하위 그룹 내 보정과 그룹 간 보정을 분리하는 계층적 다중비교 문제와 동일합니다.

자주 만나는 결함

• PEER factor 과다 적용: 공변량을 지나치게 많이 넣으면 진짜 eQTL 신호까지 함께 제거(overcorrection)될 수 있습니다. GTEx는 표본 수에 비례해 PEER factor 개수를 정하는 가이드라인을 둡니다.
• trans-eQTL과 cis-eQTL 혼동: cis-eQTL(근처 조절)과 trans-eQTL(먼 거리·다른 염색체 조절)은 통계적 검정력과 다중검정 부담이 완전히 다릅니다. trans-eQTL은 검정 쌍의 수가 훨씬 많아 훨씬 엄격한 보정이 필요합니다.
• colocalization 사전확률을 데이터에 안 맞게 그대로 사용: $p_{12}=10^{-5}$는 일반적인 기본값이지 모든 locus에 맞는 절대 기준은 아닙니다. 유전자 밀도가 높은 지역이나 강한 LD 블록에서는 민감도 분석(sensitivity analysis)으로 사전확률 변화에 따른 결론 안정성을 확인해야 합니다.

더 깊게 파고 싶다면

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

• GTEx Consortium 대표 논문: GTEx Consortium (2020), The GTEx Consortium atlas of genetic regulatory effects across human tissues, Science, 369(6509):1318–1330.
• coloc 원 논문: Giambartolomei et al. (2014), Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics, PLOS Genetics, 10(5):e1004383. (원조 coloc.abf는 각 형질이 분석 구간 안에서 인과 변이를 최대 1개만 갖는다고 가정합니다 — 신호가 여러 개인 locus에는 coloc.susie처럼 fine-mapping과 결합한 확장판이 필요합니다.)
• METASOFT 원 논문: Han & Eskin (2011), Random-Effects Model Aimed at Discovering Associations in Meta-Analysis of Genome-wide Association Studies, American Journal of Human Genetics, 88(5):586–598.
• GTEx Portal: gtexportal.org — eQTL 요약 통계량과 조직별 발현 시각화 무료 제공.

인과 변이가 어떤 유전자의 발현을 어느 조직에서 바꾸는지까지 연결했습니다. 이제 개별 SNP의 관점에서 벗어나, 유전체 전체 수준에서 집단(population) 간 조상 구조와 혼합(admixture)이 이 모든 분석에 어떤 영향을 주는지 살펴볼 차례입니다.

다음 편 S46에서는 PLINK와 ADMIXTURE로 집단 유전학의 기본기 — 조상 구성 추정과 하디-바인베르크 평형 검정을 다룹니다.