GATK으로 변이 호출하기 — BAM에서 VCF까지 (NGS Pipelines 3)

TL;DR — 정렬이 끝난 BAM을 손에 쥐었다면, 다음 질문은 "이 사람의 유전체에서 참조와 다른 자리(변이)는 어디인가"입니다. 업계 표준 도구가 GATK (Genome Analysis Toolkit)이고, Broad Institute가 정리한 GATK Best Practices 워크플로를 따르면 됩니다.

흐름은 MarkDuplicates (PCR 중복 표시) → BaseRecalibrator·ApplyBQSR (염기품질 재보정) → HaplotypeCaller (변이 호출, GVCF 모드) → GenotypeGVCFs (유전형 결정) → VariantFiltration (하드필터)입니다. 결과물은 SNP과 INDEL이 담긴 VCF (변이 포맷) 파일입니다.

핵심 함정 세 가지만 미리 짚습니다. 참조 유전체엔 .fai·.dict가 같이 있어야 하고, BAM의 모든 read엔 리드그룹(read group)이 붙어 있어야 하며, BQSR엔 알려진 변이 자리(known-sites) 파일이 필요합니다. 이 셋 중 하나라도 빠지면 GATK는 첫 단계에서 멈춥니다.

🔗 관련 글  ·  변이 호출이 무엇인지 개념부터 보려면 변이 호출(variant calling)이란?  ·  결과물 VCF 포맷의 구조는 VCF 파일이란?  ·  그 앞 단계인 정렬은 서열 정렬(read alignment)이란?

이 글에서 만들 것

서열 정렬 글에서 read를 참조 유전체에 맞춰 BAM을 만들었습니다. BAM은 "각 read가 유전체 어디에 붙었나"를 담고 있지만, 그 자체로는 변이 목록이 아닙니다. 같은 자리에 붙은 수십 개의 read를 모아 "이 자리는 참조가 A인데 이 샘플은 G가 많네 → 변이일 가능성이 높다"를 통계적으로 판정하는 단계가 변이 호출(variant calling)입니다.

이 글은 conda로 GATK4와 samtools를 깔고, 참조 유전체를 준비한 다음, GATK Best Practices의 생식세포(germline) 단변이 호출 워크플로를 단계별로 돌립니다. 중복 표시부터 염기품질 재보정, 변이 호출, 유전형 결정, 필터링까지 거쳐 최종 VCF를 얻습니다. 마지막에 자주 만나는 에러를 정리하고, GVCF 기반 joint genotyping과 딥러닝 대안인 DeepVariant까지 한 번에 살펴봅니다.

# 핵심 흐름 한눈에 (GATK Best Practices · 생식세포 단변이)
gatk MarkDuplicates      -I aln.bam   -O dedup.bam   -M metrics.txt   # ① PCR 중복 표시
gatk BaseRecalibrator    -I dedup.bam -R ref.fa --known-sites known.vcf.gz -O recal.table  # ② 재보정 표 학습
gatk ApplyBQSR           -I dedup.bam -R ref.fa --bqsr-recal-file recal.table -O bqsr.bam   #    재보정 적용
gatk HaplotypeCaller     -I bqsr.bam  -R ref.fa -O sample.g.vcf.gz -ERC GVCF               # ③ 변이 호출(GVCF)
gatk GenotypeGVCFs       -V sample.g.vcf.gz -R ref.fa -O raw.vcf.gz                          # ④ 유전형 결정
gatk VariantFiltration   -V raw.vcf.gz -R ref.fa --filter-name ... -O filtered.vcf.gz        # ⑤ 하드필터

그림 1. GATK Best Practices 생식세포 단변이 워크플로 — BAM을 중복 표시·재보정한 뒤 변이를 호출하고, 유전형 결정과 필터링을 거쳐 최종 VCF를 만든다.

1. 사전 준비 — 설치·참조 준비·입력 BAM

설치 (conda)

GATK4는 bioconda로 설치합니다. 2025년 4월 기준 최신은 4.6.2.0입니다. GATK는 자바(Java) 위에서 도는데, conda 패키지가 의존성으로 OpenJDK까지 같이 깔아 주므로 따로 자바를 맞출 필요는 없습니다. 보조 도구인 samtools는 같은 환경에 함께 넣습니다.

# gatk4 전용 conda 환경 생성 (samtools 동봉)
conda create -n gatk4 -c conda-forge -c bioconda gatk4 samtools
conda activate gatk4

# 설치 확인
gatk --version          # The Genome Analysis Toolkit (GATK) v4.6.2.0
samtools --version | head -1

gatk은 내부적으로 Picard·GATK 엔진을 한 명령으로 묶은 래퍼(wrapper)입니다. 예전에 Picard를 따로 받아 java -jar picard.jar로 돌리던 MarkDuplicates 같은 도구도 GATK4에서는 gatk MarkDuplicates로 똑같이 부릅니다.

참조 유전체 준비 — .fai와 .dict는 필수

GATK는 참조 유전체 FASTA 하나만으로는 동작하지 않습니다. 같은 폴더에 인덱스(.fai)와 사전(.dict) 두 파일이 더 있어야 합니다. .fai는 samtools가, .dict는 GATK(CreateSequenceDictionary)가 만듭니다. 이 둘을 빠뜨리는 게 GATK 입문자가 가장 자주 막히는 지점입니다.

# 참조 FASTA (예: ref.fa) 가 준비됐다고 가정

# ① FASTA 인덱스 → ref.fa.fai 생성
samtools faidx ref.fa

# ② 서열 사전 → ref.dict 생성 (ref.fa.dict 가 아니라 ref.dict)
gatk CreateSequenceDictionary -R ref.fa

ls ref.*
# ref.fa  ref.fa.fai  ref.dict

입력 BAM — 정렬·정렬순서·리드그룹·인덱스

입력 BAM은 네 가지 조건을 갖춰야 합니다. ⓐ 참조에 정렬돼 있고, ⓑ 좌표 순으로 정렬(coordinate-sorted)됐으며, ⓒ 모든 read에 리드그룹(read group, @RG)이 붙어 있고, ⓓ BAM 인덱스(.bai)가 있어야 합니다. 이 중 리드그룹 누락이 특히 흔한 사고인데, HaplotypeCaller가 샘플을 식별하는 데 @RG의 SM (sample) 태그를 쓰기 때문입니다.

# 좌표 정렬 + 인덱스 (이미 정렬돼 있으면 인덱스만)
samtools sort aln.bam -o aln.sorted.bam
samtools index aln.sorted.bam

# 리드그룹이 없다면 추가 (SM= 샘플명, 이 값이 VCF 컬럼명이 됨)
gatk AddOrReplaceReadGroups \
  -I aln.sorted.bam \
  -O aln.rg.bam \
  -RGID sample1 -RGLB lib1 -RGPL ILLUMINA -RGPU unit1 -RGSM sample1
samtools index aln.rg.bam

💡 재현용 소형 데이터 — 직접 돌려 보고 싶다면 GATK가 배포하는 테스트 번들이나 samtools의 작은 예제로 시작하면 됩니다. 대표적으로 효모(S. cerevisiae)나 대장균(E. coli) 같은 작은 유전체 + 공개 SRA read 일부를 BWA로 정렬하면 노트북에서도 몇 분 안에 전 과정을 돌릴 수 있습니다. 사람 데이터를 쓴다면 BQSR용 known-sites는 GATK Resource Bundle (Broad가 공개)의 dbSNP·Mills indel VCF를 받습니다.

2. 1단계 — MarkDuplicates (PCR 중복 표시)

라이브러리를 만들 때 PCR로 DNA를 증폭하면, 원래 한 조각이던 분자가 여러 벌로 복제됩니다. 시퀀싱하면 이 복제본들이 같은 자리에 붙은 똑같은 read로 나타나는데, 이걸 진짜 독립 증거로 세면 변이 신호가 부풀려집니다. MarkDuplicates는 이런 중복 read를 찾아 플래그(flag)만 달아 둡니다 — 지우는 게 아니라 표시만 해서, 뒤 단계가 "이건 중복이니 한 번만 세라"를 알게 합니다.

# ① PCR/광학 중복 표시 (제거가 아니라 마킹)
gatk MarkDuplicates \
  -I aln.rg.bam \
  -O dedup.bam \
  -M dedup_metrics.txt          # 중복률 등 통계 리포트

samtools index dedup.bam

-M으로 받는 메트릭 파일에는 라이브러리별 중복 비율(PERCENT_DUPLICATION)이 적힙니다. 이 값이 비정상적으로 높으면(예: 50%↑) 라이브러리 복잡도가 낮거나 input DNA가 적었다는 신호라, 호출 결과를 해석할 때 참고합니다.

3. 2단계 — BQSR (염기품질 재보정)

시퀀서가 매기는 염기품질점수(base quality score)에는 기계·화학에서 오는 계통오차(systematic error)가 섞여 있습니다. 예를 들어 특정 염기서열 맥락이나 read 위치에서 품질을 일관되게 과대·과소평가하는 식입니다. 염기품질점수 재보정(Base Quality Score Recalibration, BQSR)은 이미 알려진 변이 자리(known-sites)를 빼고 남은 불일치를 "진짜 오류"로 간주해 오차 패턴을 학습하고, 그 패턴으로 품질점수를 다시 매깁니다.

BQSR은 두 단계입니다. BaseRecalibrator가 재보정 모델을 학습해 표(table)로 내보내고, ApplyBQSR가 그 표로 BAM의 품질점수를 갈아 끼웁니다.

# ② 재보정 모델 학습 → recal.table
#    --known-sites 는 "여기 불일치는 진짜 변이지 오류 아님"을 알려주는 정답지
gatk BaseRecalibrator \
  -I dedup.bam \
  -R ref.fa \
  --known-sites dbsnp.vcf.gz \
  --known-sites Mills_and_1000G_gold_standard.indels.vcf.gz \
  -O recal.table

# ② 학습한 표를 적용 → 품질점수 재보정된 BAM
gatk ApplyBQSR \
  -I dedup.bam \
  -R ref.fa \
  --bqsr-recal-file recal.table \
  -O bqsr.bam

--known-sites가 없으면 BQSR는 모든 불일치를 오류로 착각합니다 — 실제 변이까지 오류로 학습해 품질을 망가뜨리므로, 이 파일은 사실상 필수입니다. 사람이라면 dbSNP과 Mills/1000G gold standard indel VCF를 GATK Resource Bundle에서 받아 씁니다. 모델 생물이나 known-sites가 없는 비표준 종이라면, 한 번 변이를 호출해 고신뢰 세트를 만든 뒤 그걸 known-sites로 되먹이는 부트스트랩(bootstrap) 방식을 씁니다.

4. 3단계 — HaplotypeCaller (변이 호출, GVCF 모드)

이제 변이를 부릅니다. HaplotypeCaller (변이 호출기)는 단순히 자리마다 염기를 세지 않고, 변이가 많아 보이는 활성 구간(active region)을 잡아 그 부위의 read를 국소 재조립(local de novo assembly)합니다. 후보 하플로타입(haplotype, 한쪽 염색체의 서열 조합)을 여러 개 세우고, 각 read가 어느 하플로타입에서 나왔을 확률을 계산해 유전형을 정합니다. 단순 카운팅보다 INDEL 주변에서 특히 정확합니다.

# ③ 변이 호출 — GVCF 모드 (-ERC GVCF)
gatk HaplotypeCaller \
  -I bqsr.bam \
  -R ref.fa \
  -O sample1.g.vcf.gz \
  -ERC GVCF                      # 모든 자리의 신뢰도를 기록하는 GVCF로 출력

-ERC GVCF가 핵심입니다. 이걸 빼면 변이가 있는 자리만 적힌 평범한 VCF가 나오지만, 붙이면 모든 자리(변이가 없는 자리 포함)의 신뢰도를 블록으로 기록한 GVCF (genomic VCF)가 나옵니다. GVCF는 그 자체로 끝이 아니라, 여러 샘플을 함께 유전형 결정할 때 쓰는 중간 산출물입니다.

GVCF와 joint genotyping이 왜 중요한가

샘플을 한 명씩 따로 끝까지 호출하면, 어떤 샘플에서 "이 자리는 참조와 같음"인지 "데이터가 없어 모름"인지 구별이 안 됩니다. GVCF는 변이가 없는 자리도 "참조와 같다고 이만큼 확신한다"를 적어 두기 때문에, 여러 샘플의 GVCF를 합쳐 한꺼번에 유전형을 정하는 joint genotyping (공동 유전형 결정)이 가능해집니다. 이렇게 하면 한 샘플에선 애매하던 변이가 다른 샘플의 증거로 보강돼 호출 정확도가 올라갑니다. 이것이 GATK가 단일 샘플조차 GVCF로 뽑기를 권하는 이유입니다.

# (여러 샘플) 각 샘플을 GVCF로 뽑은 뒤 합치기
#   - 적은 수의 GVCF 병합: CombineGVCFs
#   - 대규모 코호트: GenomicsDBImport (DB로 통합, 확장성↑)
gatk CombineGVCFs \
  -R ref.fa \
  -V sample1.g.vcf.gz -V sample2.g.vcf.gz -V sample3.g.vcf.gz \
  -O cohort.g.vcf.gz

5. 4단계 — GenotypeGVCFs (유전형 결정)

GVCF는 아직 "최종 변이표"가 아닙니다. GenotypeGVCFs가 GVCF(단일 샘플이면 그 GVCF, 다중이면 합친 GVCF)를 받아, 자리마다 실제 유전형(0/0, 0/1, 1/1)을 매기고 변이 자리만 추려 일반 VCF로 내보냅니다.

# ④ GVCF → 유전형이 매겨진 VCF
gatk GenotypeGVCFs \
  -R ref.fa \
  -V sample1.g.vcf.gz \
  -O raw.vcf.gz                  # 여기서 나온 게 '날것(raw)' 변이 — 아직 필터 전

여기서 나온 raw.vcf.gz는 SNP과 INDEL이 모두 들어 있지만 아직 필터를 거치지 않은 날것입니다. 위양성(false positive)이 섞여 있으므로, 그대로 쓰지 말고 다음 단계에서 걸러야 합니다.

\

그림 2. SNP과 INDEL — read 더미에서 참조와의 차이를 모아 한 자리의 변이 유형(치환·삽입·결실)을 판정한다.

6. 5단계 — 변이 필터링 (하드필터 vs VQSR)

날것 VCF에서 신뢰할 변이만 남기는 방법은 크게 둘입니다. 하드필터(hard-filtering)는 QD·FS·MQ 같은 변이 주석(annotation) 값에 직접 임계값을 걸어 거르고, VQSR (Variant Quality Score Recalibration)은 알려진 고신뢰 변이로 기계학습 모델(가우시안 혼합)을 학습해 변이마다 신뢰점수를 매깁니다.

구분	하드필터(VariantFiltration)	VQSR (VariantRecalibrator)
방식	주석 값에 고정 임계값 적용	알려진 변이로 학습한 확률 모델
데이터 요구	적음 — 한 샘플로도 가능	많음 — 대규모 샘플·exome/WGS 권장
학습 자원	불필요	known/training 자원 세트 필요(HapMap·omni·1000G·dbSNP)
장점	단순·투명·소규모에 적합	더 정밀, 대규모에서 위양성↓
한계	임계값이 거칠다	표본 적으면 모델이 불안정
언제	샘플 수가 적거나 비모델 종	대규모 코호트, 사람 WGS/WES

소규모 프로젝트나 표본이 적을 때는 하드필터가 현실적입니다. GATK가 권하는 방식은 SNP과 INDEL을 나눠 각각 다른 임계값을 거는 것입니다 — 둘은 오류 특성이 다르기 때문입니다. 먼저 SelectVariants로 타입을 분리한 뒤 VariantFiltration을 겁니다.

# 6-1. SNP 만 분리 → 하드필터
gatk SelectVariants -R ref.fa -V raw.vcf.gz \
  --select-type-to-include SNP -O raw.snp.vcf.gz

gatk VariantFiltration -R ref.fa -V raw.snp.vcf.gz -O filtered.snp.vcf.gz \
  --filter-name "QD2"     --filter-expression "QD < 2.0" \
  --filter-name "FS60"    --filter-expression "FS > 60.0" \
  --filter-name "MQ40"    --filter-expression "MQ < 40.0" \
  --filter-name "MQRS-12.5" --filter-expression "MQRankSum < -12.5" \
  --filter-name "RPRS-8"  --filter-expression "ReadPosRankSum < -8.0" \
  --filter-name "SOR3"    --filter-expression "SOR > 3.0"

# 6-2. INDEL 만 분리 → 하드필터 (임계값이 SNP과 다름)
gatk SelectVariants -R ref.fa -V raw.vcf.gz \
  --select-type-to-include INDEL -O raw.indel.vcf.gz

gatk VariantFiltration -R ref.fa -V raw.indel.vcf.gz -O filtered.indel.vcf.gz \
  --filter-name "QD2"     --filter-expression "QD < 2.0" \
  --filter-name "FS200"   --filter-expression "FS > 200.0" \
  --filter-name "RPRS-20" --filter-expression "ReadPosRankSum < -20.0" \
  --filter-name "SOR10"   --filter-expression "SOR > 10.0"

# 6-3. 다시 합치기 (선택)
gatk MergeVcfs -I filtered.snp.vcf.gz -I filtered.indel.vcf.gz -O filtered.vcf.gz

VariantFiltration은 조건에 걸린 변이를 지우지 않고 VCF의 FILTER 칸에 이름표(예: QD2)만 답니다. 통과한 변이는 PASS로 남습니다. 그래서 나중에 PASS만 골라 쓰면 됩니다. 임계값의 의미를 간단히 보면 — QD(QualByDepth)는 깊이로 정규화한 변이 신뢰도라 낮으면 약한 증거, FS(FisherStrand)는 가닥 편향(strand bias)이라 높으면 한쪽 가닥에만 쏠린 의심스러운 신호, MQ(RMSMappingQuality)는 정렬 품질이라 낮으면 모호한 영역, SOR(StrandOddsRatio)도 가닥 편향의 다른 척도입니다.

# PASS 변이만 추출
gatk SelectVariants -R ref.fa -V filtered.vcf.gz \
  --exclude-filtered -O final.PASS.vcf.gz

그림 3. 변이 필터링 두 갈래 — 하드필터(고정 임계값)는 소규모에, VQSR(학습 모델)은 대규모 코호트에 적합하다.

7. 자주 발생하는 에러 & 해결

• A USER ERROR has occurred: Fasta dict file ... does not exist / index 없음 — 참조 FASTA 옆에 .dict나 .fai가 없습니다. gatk CreateSequenceDictionary -R ref.fa로 .dict를, samtools faidx ref.fa로 .fai를 만드세요. 셋(ref.fa·ref.fa.fai·ref.dict)이 같은 폴더에 있어야 합니다.
• Sample ... read group ... is malformed / does not have read group 오류 — BAM의 read에 리드그룹(@RG)이 없습니다. gatk AddOrReplaceReadGroups로 -RGSM(샘플명) 등을 채워 넣으세요. HaplotypeCaller는 @RG의 SM 태그로 샘플을 식별합니다.
• BQSR 결과가 이상하다 / --known-sites를 안 줬다 — known-sites 없이 BQSR를 돌리면 실제 변이까지 "오류"로 학습해 품질점수가 망가집니다. dbSNP·Mills indel 같은 known-sites VCF를 반드시 주세요(사람은 GATK Resource Bundle). 비모델 종은 부트스트랩으로 known-sites를 만듭니다.
• Java heap space / 메모리 부족(OutOfMemory) — 대용량 BAM·WGS에서 흔합니다. --java-options "-Xmx16g"로 힙을 키우거나(gatk --java-options "-Xmx16g" HaplotypeCaller ...), 인터벌(interval)로 영역을 쪼개 병렬로 돌리세요(-L chr1 등).
• BAM 인덱스(.bai) 없음 — 좌표 정렬한 뒤 samtools index dedup.bam으로 인덱스를 만드세요. 정렬되지 않은 BAM은 인덱스가 안 만들어지니, samtools sort부터 합니다.
• VCF가 너무 크고 느리다 — 단일 샘플을 끝까지 일반 VCF로 뽑기보다, GVCF(-ERC GVCF)로 뽑아 두면 나중에 샘플을 추가하거나 joint genotyping할 때 처음부터 다시 안 돌려도 됩니다.

🎓 다음 학습
• (개념 복습) 변이 호출(variant calling)이란? — SNP·INDEL을 가려내는 원리를 그림으로
• (결과 포맷) VCF 파일이란? — 이 글이 만든 VCF의 컬럼·INFO·FORMAT 해부
• (앞 단계) 서열 정렬(read alignment)이란? — 변이 호출의 입력 BAM이 만들어지는 과정

8. 확장 — VCF 주석(annotation)과 DeepVariant

변이에 의미 붙이기 — VCF annotation

final.PASS.vcf.gz는 "어느 자리가 변이인가"까지만 말해 줍니다. 그 변이가 어느 유전자의, 어떤 결과를 낳는 변화인가(아미노산이 바뀌는지, 빈도가 얼마나 희귀한지, 질환과 연관됐는지)는 주석(annotation) 단계에서 붙입니다. 대표 도구는 SnpEff, VEP (Variant Effect Predictor), ANNOVAR입니다. 예를 들어 SnpEff는 변이가 단백질에 미치는 영향(missense·nonsense·synonymous 등)을 예측해 VCF의 INFO 칸에 적어 줍니다.

# 예: SnpEff 로 기능적 영향 주석 (genome DB 는 사전 설치)
snpEff -v GRCh38.105 final.PASS.vcf.gz > annotated.vcf
# → INFO 에 ANN= 필드가 추가되어 유전자·영향·HGVS 표기가 붙는다

주석까지 끝나야 비로소 "이 변이가 임상적으로/생물학적으로 의미 있나"를 해석할 수 있습니다. VCF 컬럼과 INFO 필드를 어떻게 읽는지는 VCF 파일 글에서 다룹니다.

딥러닝 대안 — DeepVariant

GATK가 통계 모델로 변이를 부른다면, Google의 DeepVariant는 같은 문제를 이미지 인식으로 바꿔 풉니다. 같은 자리에 쌓인 read 더미(pileup)를 이미지 텐서로 만들고, 합성곱 신경망(CNN)으로 가운데 염기가 변이인지 분류합니다. 2025년 기준 최신은 1.10 버전입니다.

가장 큰 차이는 운영 방식입니다. GATK는 위에서 본 것처럼 여러 명령을 순서대로 거치지만, DeepVariant는 한 명령으로 끝납니다(BQSR도 필요 없습니다).

# DeepVariant — 한 명령으로 BAM → VCF (도커 실행)
docker run \
  -v "$PWD":/data google/deepvariant:1.10.0 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/ref.fa \
  --reads=/data/bqsr.bam \
  --output_vcf=/data/dv.vcf.gz \
  --output_gvcf=/data/dv.g.vcf.gz \
  --num_shards=$(nproc)

벤치마크에서 DeepVariant는 특히 SNP과 어려운 영역(저복잡도·homopolymer)의 INDEL 정확도가 높은 편으로 보고됩니다. 반면 GATK는 대규모 코호트의 joint genotyping과 희귀변이 처리에서 표준 워크플로로 자리잡고 있습니다. 둘은 우열이라기보다 쓰임이 다른데, 임상·소수 샘플 정밀 분석엔 DeepVariant가, 집단유전체·다수 샘플 통합엔 GATK가 흔히 쓰입니다. DeepVariant도 GVCF를 뽑을 수 있어, 그 GVCF를 GATK의 GenotypeGVCFs로 joint genotyping에 넘기는 조합도 실무에서 씁니다.

References

1. Van der Auwera, G.A. & O'Connor, B.D. (2020). Genomics in the Cloud: Using Docker, GATK, and WDL in Terra. O'Reilly. (GATK Best Practices 공식 교재)
2. McKenna, A. et al. (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, 20, 1297–1303. doi:10.1101/gr.107524.110
3. DePristo, M.A. et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 43, 491–498. doi:10.1038/ng.806
4. Poplin, R. et al. (2018). A universal SNP and small-indel variant caller using deep neural networks. Nature Biotechnology, 36, 983–987. doi:10.1038/nbt.4235 (DeepVariant 원논문)
5. GATK Best Practices — Germline short variant discovery (SNPs + Indels). Broad Institute. gatk.broadinstitute.org
6. Hard-filtering germline short variants. GATK Documentation. gatk.broadinstitute.org

Pipette & Pipeline · A bio portfolio journal