[서열분석개론] 서열 분석 - 변이 호출 파이프라인

정렬이 끝난 BAM이 있다고 바로 변이를 부를 수 있는 건 아닙니다.

전처리를 거쳐야 신뢰도 높은 변이 목록(VCF)을 얻을 수 있고, 그렇게 얻은 VCF도 표기를 통일해야 다른 DB와 비교가 가능해집니다.
이 단계를 순서대로 정리합니다.

---

GATK Best Practices 흐름

sorted BAM → MarkDuplicates → BQSR → Variant Caller → VCF

 

GATK(Genome Analysis Toolkit)는 Broad Institute에서 만든 변이 검출 툴킷으로, 여기서 다루는 MarkDuplicates, BaseRecalibrator, HaplotypeCaller/Mutect2가 모두 이 안에 포함되어 있습니다.
변이 호출 분야의 golden standard로 통합니다.

MarkDuplicates — PCR 중복 표시

라이브러리를 만드는 과정에서 PCR 증폭을 거치는데, 이때 같은 원본 DNA 조각이 여러 개의 read로 복제됩니다.
문제는 이 중복 read들이 독립적인 관측이 아니라는 점입니다.
그대로 세면 특정 변이의 depth나 대립빈도(AF)가 실제보다 부풀려져서, 가짜 변이가 진짜처럼 보일 수 있습니다.

 

MarkDuplicates는 같은 시작 좌표를 가진 read들을 찾아서 duplicate 플래그(FLAG 0x400/1024)로 표시합니다.

여기서 중요한 건 "제거"가 아니라 "표시"라는 점입니다.

이후 변이 검출기가 이 플래그를 보고 중복 read를 한 번만 세도록 만드는 게 목적입니다.

BQSR — 염기 품질 재보정

시퀀서가 자체적으로 매기는 Phred 품질 점수(Q)는 완벽하지 않습니다.

BQSR(Base Quality Score Recalibration)은 실제 오류율을 다시 계산해서 이 점수를 보정하는 과정입니다.

1. read와 레퍼런스가 다른 위치를 셉니다.
   단, 알려진 SNP(mills, known_indels 등)는 오류가 아니라 진짜 변이일 수 있으니 카운트에서 제외합니다.
2. 예를 들어 "Q30(1/1000 확률로 틀림)"이라고 표시된 염기가 실제로는 1/500 확률로 틀렸다면, 진짜 실력은 Q27인 셈이라 점수를 낮춥니다.
3. read 내 위치(cycle), 앞뒤 염기 배열(context) 같은 조건별로 그룹을 나눠서 각각 따로 보정합니다.

정리하면 BQSR은 "시퀀서가 자기 실력을 과대평가한 부분을 실측 데이터로 바로잡는" 과정입니다.

---

Germline vs Somatic

어떤 변이 검출기를 쓸지는 찾으려는 변이의 성격에 따라 달라집니다.

구분	Germline (생식세포)	Somatic (체세포)
변이 특성	유전됨, 모든 세포에 공통	후천적으로 생김, 일부 세포에만 존재 (예: 암)
대립빈도(AF)	~0.5(이형접합) / ~1.0(동형접합)	다양하고 낮을 수 있음
대표 도구	HaplotypeCaller	Mutect2

 

희귀질환처럼 유전되는 변이를 찾을 때는 germline 표준 도구인 HaplotypeCaller를 씁니다.

정식 germline 파이프라인은 GVCF를 만들고 joint genotyping을 거치는 과정까지 포함하는 게 일반적입니다.

 

Mutect2는 원래 종양 조직에서 정상 조직 대비 후천적으로 생긴 변이를 찾는 somatic caller입니다.

다만 tumor-only 모드로 쓰면 스텝이 단순해져서(Mutect2 → FilterMutectCalls 2단계), 교육이나 실습 목적으로 germline 변이를 빠르게 훑는 용도로 쓰기도 합니다.

실제 진단용 germline 분석에서는 HaplotypeCaller를 쓰는 게 표준이라는 점은 구분해서 알아둘 필요가 있습니다.

---

VCF 구조

VCF(Variant Call Format)는 변이를 담는 표준 텍스트 포맷입니다.

##fileformat=VCFv4.2
##INFO=<ID=DP,...>
#CHROM  POS       ID  REF  ALT  QUAL  FILTER  INFO       FORMAT       HG002
chr21   44287133  .   G    A    .     PASS    DP=73;...  GT:AF:DP     0/1:0.468:72

• ##로 시작하는 줄은 헤더로, 포맷 버전과 INFO/FORMAT 필드 정의, 레퍼런스 정보 등이 담깁니다
• POS는 변이 위치(1-based), REF→ALT가 실제 서열 변화(G→A)
• FILTER가 PASS면 필터를 통과했다는 뜻. INFO는 변이 수준의 주석 정보
• FORMAT은 샘플별로 어떤 값이 들어있는지 정의. GT(유전형, 0/1은 이형접합), AF(대립빈도), DP(depth) 순서로 값이 붙습니다.

---

왜 정규화가 필요한가

같은 변이도 표기 방식이 여러 가지일 수 있습니다.

표기를 통일해야 ClinVar, gnomAD 같은 DB와 정확히 매칭됩니다.

이 작업을 vt라는 도구가 담당하고, 크게 두 단계로 나뉩니다.

decompose

다중대립(multiallelic) 변이를 개별 줄로 분해합니다.

# 분해 전 (한 줄에 ALT가 2개)
chr1  100  A  G,T  GT:AD  1/2:0,15,20

# 분해 후 (각각 별도 줄)
chr1  100  A  G  GT:AD  ...
chr1  100  A  T  GT:AD  ...

normalize

indel을 left-align하고 가장 간결한 형태로 통일합니다.

# 정규화 전 (레퍼런스: ...G A A A T... A가 반복되는 구간)
chr1  102  AA  A

# 정규화 후 (같은 1bp 결실이지만 가장 왼쪽·최소 표현으로 통일)
chr1  100  GA  G

 

이렇게 표기를 통일해야 같은 변이인데도 표기 방식이 달라서 DB 매칭에서 놓치는 일을 막을 수 있습니다.