BioPlayground

🧬
목록으로

HaplotypeCaller: 로컬 어셈블리와 Bayesian 콜링으로 Germline 변이 부르기

왜 HaplotypeCaller가 위치별 콜러보다 정확할까요. 로컬 De Bruijn 어셈블리로 후보 하플로타입을 만들고, PairHMM으로 리드를 매핑하고, Bayes로 유전자형을 확정하는 세 단계를 유도합니다.

중급
|
22분
|
검증 완료 (2026-07-22)
germline variant callingHaplotypeCallerGVCF
진행률0/52 (0%)

지도 위 위치 다시 확인

앞선 세 편에서 우리는 라이브러리 편향(MarkDuplicates)과 시퀀서 편향(BQSR)을 통제했습니다. 이제 남은 리드가 그 자체로 "이 개인의 유전체가 참조와 어떻게 다른지"에 대한 정직한 관찰이 됩니다. 이 관찰에서 실제 변이를 뽑는 게 이번 편의 HaplotypeCaller입니다.

Germline 파이프라인의 심장부입니다. 서두르지 맙시다.

왜 로컬 어셈블리인가

기존 방법(GATK3 UnifiedGenotyper 이전)은 각 위치를 독립적으로 봤습니다. "이 좌표에 리드가 몇 개 참조와 다른 염기를 보고했나?" 이 접근이 SNP는 잘 잡았지만 인델 근처에서 크게 약했습니다. 이유는 이렇습니다.

인델이 있으면 그 주위 정렬이 여러 방식으로 가능해서 리드마다 다르게 붙습니다. 위치별 관찰이 흩어져서 "이 좌표에 SNP가 있음"이 아니라 "이 구간에 정체불명의 노이즈가 많음"으로 나옵니다.

HaplotypeCaller의 해법은 이렇습니다.

  1. 관심 구간(수백 bp)의 리드를 다 모읍니다.
  2. 로컬 어셈블리로 그 구간에서 가능한 하플로타입 후보를 만듭니다.
  3. 각 후보 하플로타입에 리드를 다시 붙여봅니다. 어느 후보가 리드 세트를 가장 잘 설명하나요?
  4. Bayes 규칙으로 유전자형을 확정합니다.

이 접근이 위치별 접근에 비해 인델 근처에서 극적으로 정확합니다.

네 단계로 나눠서 유도

1단계: Active Region 검출

전 게놈을 다 어셈블할 수 없으니 "뭔가 다른 게 있는 구간"만 골라야 합니다.

  • 리드가 참조와 다르게 붙은 지점, 소프트 클립이 많은 지점, MAPQ가 낮은 지점을 스캔.
  • 이 지점을 중심으로 ±수백 bp 창을 active region으로 지정.
  • 나머지 구간은 참조와 동일하다고 간주 (연산 절약).

전 게놈 3×10⁹bp 중 실제로 조사해야 할 active region은 대략 1% 미만입니다.

2단계: 로컬 어셈블리 — De Bruijn 그래프

active region의 리드들로 소형 De Bruijn 그래프를 만듭니다 (M23 참조). 그래프에서 참조 경로와 대안 경로들을 뽑아냅니다.

  • 참조 경로: 참조 그대로.
  • 대안 경로 1: 참조와 한 위치에 SNP.
  • 대안 경로 2: 참조에 인델.
  • 대안 경로 3: 참조와 완전히 다른 구조 (구조변이 후보).

이 후보들이 후보 하플로타입입니다. HaplotypeCaller가 어셈블리를 통해 만들기 때문에 인델·복잡변이도 자연스럽게 후보에 포함됩니다.

3단계: PairHMM으로 각 리드가 각 하플로타입에 얼마나 잘 붙나 계산

각 (리드 r, 하플로타입 h) 쌍에 대해 P(r|h)를 계산합니다. 이 계산이 PairHMM입니다.

  • 상태: 매치 · 인서션 · 딜리션.
  • 방출 확률: 리드의 Phred (BQSR 후) 기반.
  • 전이 확률: 갭 오픈·연장 확률.
  • Forward 알고리즘 (M19 참조)으로 전체 정렬 확률을 로그 합.

각 리드에 대해 P(r|h1), P(r|h2), ... 벡터가 나옵니다.

4단계: Bayes 유전자형 결정

두 하플로타입 조합 (h_i, h_j)이 개인의 유전자형이라면, 관찰된 리드 세트의 우도는

P(reads∣hi,hj)=∏rP(r∣hi)+P(r∣hj)2P(\text{reads}|h_i, h_j) = \prod_{r} \frac{P(r|h_i) + P(r|h_j)}{2}

각 리드는 두 하플로타입 중 하나에서 왔을 확률이 반반이라는 가정입니다. 이 우도에 사전 확률(인구 알렐 빈도)을 곱해 사후 확률을 얻고, argmax 유전자형이 최종 콜입니다.

GT=arg⁡max⁡(hi,hj)P(hi,hj∣reads)\text{GT} = \arg\max_{(h_i, h_j)} P(h_i, h_j | \text{reads})

이게 HaplotypeCaller의 핵심 통계입니다.

GVCF — 다중 샘플 통합의 관문

한 샘플만 콜하면 그냥 VCF를 뱉으면 됩니다. 그런데 30명, 300명, 30,000명을 통합 분석하고 싶으면?

각 샘플을 독립적으로 VCF로 뱉고 나중에 합치면 문제가 생깁니다. A 샘플에서만 콜된 변이가 다른 샘플에는 아예 기록이 없어 (missing data) 결측인지 참조인지 알 수 없어집니다.

**GVCF (Genomic VCF)**가 이 문제를 해결합니다.

  • 각 샘플에 대해 모든 위치를 기록합니다.
  • 변이가 있는 위치는 상세 정보와 함께.
  • 변이가 없는 위치는 "이 구간은 참조와 같다"고 블록으로 기록.

이후 여러 GVCF를 통합할 때 각 위치의 참조/변이 여부가 명시적입니다. 통합 후 GenotypeGVCFs로 다중 샘플 콜을 진행합니다 (S11에서 다룹니다).

실행 명령

Germline WGS 표준입니다.

bash
gatk HaplotypeCaller \
    -I S1.recal.bam \
    -R hg38.fa \
    -O S1.g.vcf.gz \
    -ERC GVCF \
    --tmp-dir /data/tmp
  • -ERC GVCF: GVCF 모드 (다중 샘플 통합 준비).
  • -ERC BP_RESOLUTION: 위치별 상세 (정밀 진단, 파일 크기 큼).
  • 단일 샘플 최종 VCF만 원하면 -ERC NONE (기본).

옵션 두어 개만 더 알아둡시다.

  • --emit-ref-confidence: 참조 신뢰도 표기.
  • --stand-call-conf 30: 최소 콜 신뢰도 (기본).
  • -L intervals.bed: 스코프 좁히기.

실무 팁 다섯 개

1. WES는 intervals bed 필수

WES는 엑솜만 관심. -L exome.bed로 스코프 좁히면 실행 시간이 20배 줄어듭니다.

2. WGS 병렬화

WGS는 크니까 염색체별로 나눠 병렬 실행 후 통합.

bash
for chr in chr{1..22} chrX chrY chrM; do
    gatk HaplotypeCaller ... -L $chr -O S1_$chr.g.vcf.gz &
done
wait
gatk MergeVcfs -I S1_chr1.g.vcf.gz -I S1_chr2.g.vcf.gz ... -O S1.g.vcf.gz

Snakemake/Nextflow에서 자동 처리.

3. --pair-hmm-implementation

기본 로직스라 GATK4는 AVX-512 SIMD 구현 사용. 컨테이너가 CPU 지원 안 하면 스칼라 fallback (~5배 느림). 서버 CPU 확인.

4. GVCF의 크기

GVCF가 VCF보다 크지만, 압축률이 좋아서 인간 WGS 하나가 대략 500MB～1GB. 30샘플이면 15～30GB.

5. HaplotypeCaller가 놓치는 경우

  • 리드 커버리지가 극히 낮은 구간 (예: 5× 이하).
  • 반복 서열 (short tandem repeat) 근처 (S14 SV에서 별도 도구).
  • 구조변이 (100bp 이상 대형 변이) — HaplotypeCaller는 인델 ~40bp까지가 실용 한계.

이런 경우 DeepVariant(S12)나 SV 도구(S14)로 보완.

손 계산 — Bayes 유전자형 판단

active region에서 리드 20개가 관찰됐고, 후보 하플로타입이 두 개 있다고 합시다. h1은 참조, h2는 SNP 하나.

각 리드 r에 대해 PairHMM이 이렇게 나왔다고 합시다 (평균값).

  • 12개 리드: P(r|h2) ≫ P(r|h1) → h2에 유리.
  • 8개 리드: P(r|h1) ≈ P(r|h2) → 판단 애매.

세 유전자형 후보(h1/h1, h1/h2, h2/h2) 우도를 계산:

  • h1/h1 (참조 동형): 12개 리드 관찰이 불리 (log P 크게 감소).
  • h1/h2 (헤테로): 반반 가정으로 모든 리드 설명 가능.
  • h2/h2 (변이 동형): 애매한 8개 리드 관찰이 불리.

Bayes로는 h1/h2가 argmax. 헤테로 SNP 콜.

이 로직이 진짜 HaplotypeCaller의 판단입니다.

VCF 열어보기

bash
bcftools view S1.g.vcf.gz | head -50

주요 필드.

  • CHROM · POS · REF · ALT: 위치와 참조/대안 염기.
  • QUAL: 변이 신뢰도 (Phred).
  • INFO: DP(depth), MQ(mean MAPQ), FS(strand bias) 등 통계.
  • FORMAT · Sample: GT(유전자형), AD(각 대립유전자 depth), DP, GQ(유전자형 신뢰도).

GT의 대표 값.

  • 0/0: 참조 동형.
  • 0/1: 헤테로.
  • 1/1: 대안 동형.
  • ./.: 결측.

CS 매핑

  • De Bruijn 로컬 어셈블리: M23의 그 알고리즘을 active region 단위로 재활용.
  • PairHMM: HMM의 M18～M19에서 유도한 Forward 알고리즘을 리드-하플로타입 쌍에 적용.
  • Bayes 사후 확률 최대화: 사전(알렐 빈도) + 우도(PairHMM) → 사후.

DryBench "HMM Forward 실무" 편(hmm-forward-in-practice) 참고.

마무리

HaplotypeCaller는 로컬 어셈블리·PairHMM·Bayes 정리 세 개념의 조립체입니다. 다음 편(S10)에서는 이 프레임워크가 어떻게 종양-정상 관계로 확장되어 Mutect2가 되는지 봅니다. Germline과 Somatic이 어디서 갈리는지 통계 수준에서 확인합시다.

더 깊게 파고 싶다면

  • Poplin R. et al. (2018), Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. HaplotypeCaller의 로컬 어셈블리와 GVCF 확장성.
  • GATK 공식 HaplotypeCaller 가이드 (https://gatk.broadinstitute.org/hc/en-us/articles/360036721051): 옵션 원전.
  • Broad Institute BroadE — HaplotypeCaller deep dive 유튜브 강의. active region · PairHMM 시각화가 좋습니다.
  • Rausch T. et al. (2019), Alfred: interactive multi-sample BAM alignment statistics. Bioinformatics 35:2489. 결과 검증 도구.
  • 1000 Genomes Project 공식 데이터 (https://www.internationalgenome.org/): 벤치마크 참 값.

로컬 어셈블리가 왜 인델을 잘 잡는지 감으로 잡히면 다음 편의 Somatic 콜링도 훨씬 자연스럽게 이해됩니다. 다음 편에서 Mutect2로 갑시다.