[WGS 분석] 변이 호출 (GATK HaplotypeCaller: Variant Calling)

변이 호출을 시작하기 전에 GATK 가 있는지부터 확인.

gatk --version

 

GATK가 깔려 있지 않다면,

conda install -c bioconda gatk4

---

변이 호출 전 준비(GATK가 요구하는 사전 조건들을 맞추는 과정)

변이 호출 전체 흐름

1. GATK용 reference 인덱스 준비 (.dict 생성)
2. BAM에 Read Group 있는지 확인 
3. (필요시) Read Group 추가
4. HaplotypeCaller 실행 ← 진짜 변이 호출
5. VCF 결과 확인

---

GATK용 reference 인덱스 준비

GATK를 성공적으로 설치하고 나서 GATK가 reference에 필요로 하는 추가 인덱스를 생성해줄 필요가 있음.

gatk CreateSequenceDictionary -R ~/bioinfo/bi_study/wgs_practice/ref/chr20.fa

 

GATK는 .dict 파일이 따로 필요하고, 파일이 잘 생성된 것을 확인.

---

BAM에 Read Group 있는지 확인

samtools view -H ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634.bam | grep "@RG"

 

위 명령어를 통해 BAM에 Read Group 있는지 확인했지만, 출력이 따로 없어 Read Group이 없는 것을 확인.

 

---

Read Group 추가

BAM에 Read Group이 없어 samtools로 Read Group 추가.

gatk AddOrReplaceReadGroups \
  -I ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634.bam \
  -O ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634_RG.bam \
  -RGID SRR062634 \
  -RGLB lib1 \
  -RGPL ILLUMINA \
  -RGPU unit1 \
  -RGSM NA12878

 

옵션 설명

-RGID — Read Group ID. 이 시퀀싱 run을 구분하는 고유 식별자. 보통 flowcell+lane 조합을 쓰는데 지금은 간단히 SRR062634로 설정.

-RGLB — Library 이름. DNA 라이브러리를 구분하는 식별자. 같은 샘플이라도 라이브러리를 다르게 준비하면 다른 값을 줌.

-RGPL — Platform. 어떤 시퀀싱 장비로 만들었는지. ILLUMINA, PACBIO, ONT 등.

-RGPU — Platform Unit. flowcell-barcode-lane 조합. 더 세밀한 구분용인데 지금은 간단히 unit1.

-RGSM — Sample 이름. 가장 중요한 값. 사용하는 샘플이 NA12878이니까 이걸로 설정.

 

GATK가 여러 샘플을 합쳐서 분석할 때 이 RGSM 기준으로 "이 read가 누구 거다"를 구분.

Read Group를 추가한 BAM 파일에 인덱스 생성

새롭게 생성한 BAM 파일(SRR062634_RG.bam)에 인덱스 다시 만듬.

samtools index ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634_RG.bam

 

위 명령어를 실행 후, 다시 인덱스 생성되었는지 확인.

samtools view -H ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634_RG.bam | grep "@RG"

 

잘 생성된 것을 확인.

---

변이 호출 : HaplotypeCaller 실행

gatk HaplotypeCaller \
  -R ~/bioinfo/bi_study/wgs_practice/ref/chr20.fa \
  -I ~/bioinfo/bi_study/wgs_practice/aligned/SRR062634_RG.bam \
  -O ~/bioinfo/bi_study/wgs_practice/variants/SRR062634.vcf.gz

 

변이 호출 과정은 파일이 작은 편이어도 시간이 다소 걸림.

옵션 설명

 

-R : Reference genome. (chr20.fa)

-I : Input BAM. Read Group 추가된 정렬 결과.

-O : Output VCF. 변이 호출 결과가 저장될 파일.

 

HaplotypeCaller가 하는 일은, BAM에 있는 read들을 보면서 reference랑 다른 부분(SNP, Indel)을 찾아내는 것.

단순히 한 위치씩 보는 게 아니라, 주변 영역을 다시 조립(local assembly)해서 더 정확하게 변이를 잡아내는 방식.

그래서 다른 단순 caller보다 정확도가 높음.

 

출력되는 VCF는 압축 형식(.vcf.gz)으로 나오고, SNP랑 Indel이 한 파일에 같이 들어감.

---

 

variant calling 잘 이루어졌는지 확인

ls -lh ~/bioinfo/bi_study/wgs_practice/variants/

 

 

SRR062634.vcf.gz랑 .vcf.gz.tbi 두 파일 있는지 확인.

gzcat ~/bioinfo/bi_study/wgs_practice/variants/SRR062634.vcf.gz | head -30

 

SRR062634.vcf.gz 파일 내용 확인.

##fileformat=VCFv4.2
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller --output /Users/kofe/bioinfo/bi_study/wgs_practice/variants/SRR062634.vcf.gz --input /Users/kofe/bioinfo/bi_study/wgs_practice/aligned/SRR062634_RG.bam --reference /Users/kofe/bioinfo/bi_study/wgs_practice/ref/chr20.fa --use-posteriors-to-calculate-qual false --dont-use-dragstr-priors false --use-new-qual-calculator true --annotate-with-num-discovered-alleles false --heterozygosity 0.001 --indel-heterozygosity 1.25E-4 --heterozygosity-stdev 0.01 --standard-min-confidence-threshold-for-calling 30.0 --max-alternate-alleles 6 --max-genotype-count 1024 --sample-ploidy 2 --num-reference-samples-if-no-call 0 --genotype-assignment-method USE_PLS_TO_ASSIGN --contamination-fraction-to-filter 0.0 --output-mode EMIT_VARIANTS_ONLY --all-site-pls false --flow-likelihood-parallel-threads 0 --flow-likelihood-optimized-comp false --trim-to-haplotype true --exact-matching false --flow-use-t0-tag false --flow-remove-non-single-base-pair-indels false --flow-remove-one-zero-probs false --flow-quantization-bins 121 --flow-fill-empty-bins-value 0.001 --flow-symmetric-indel-probs false --flow-report-insertion-or-deletion false --flow-disallow-probs-larger-than-call false --flow-lump-probs false --flow-retain-max-n-probs-base-format false --flow-probability-scaling-factor 10 --flow-order-cycle-length 4 --keep-boundary-flows false --gvcf-gq-bands 1 --gvcf-gq-bands 2 --gvcf-gq-bands 3 --gvcf-gq-bands 4 --gvcf-gq-bands 5 --gvcf-gq-bands 6 --gvcf-gq-bands 7 --gvcf-gq-bands 8 --gvcf-gq-bands 9 --gvcf-gq-bands 10 --gvcf-gq-bands 11 --gvcf-gq-bands 12 --gvcf-gq-bands 13 --gvcf-gq-bands 14 --gvcf-gq-bands 15 --gvcf-gq-bands 16 --gvcf-gq-bands 17 --gvcf-gq-bands 18 --gvcf-gq-bands 19 --gvcf-gq-bands 20 --gvcf-gq-bands 21 --gvcf-gq-bands 22 --gvcf-gq-bands 23 --gvcf-gq-bands 24 --gvcf-gq-bands 25 --gvcf-gq-bands 26 --gvcf-gq-bands 27 --gvcf-gq-bands 28 --gvcf-gq-bands 29 --gvcf-gq-bands 30 --gvcf-gq-bands 31 --gvcf-gq-bands 32 --gvcf-gq-bands 33 --gvcf-gq-bands 34 --gvcf-gq-bands 35 --gvcf-gq-bands 36 --gvcf-gq-bands 37 --gvcf-gq-bands 38 --gvcf-gq-bands 39 --gvcf-gq-bands 40 --gvcf-gq-bands 41 --gvcf-gq-bands 42 --gvcf-gq-bands 43 --gvcf-gq-bands 44 --gvcf-gq-bands 45 --gvcf-gq-bands 46 --gvcf-gq-bands 47 --gvcf-gq-bands 48 --gvcf-gq-bands 49 --gvcf-gq-bands 50 --gvcf-gq-bands 51 --gvcf-gq-bands 52 --gvcf-gq-bands 53 --gvcf-gq-bands 54 --gvcf-gq-bands 55 --gvcf-gq-bands 56 --gvcf-gq-bands 57 --gvcf-gq-bands 58 --gvcf-gq-bands 59 --gvcf-gq-bands 60 --gvcf-gq-bands 70 --gvcf-gq-bands 80 --gvcf-gq-bands 90 --gvcf-gq-bands 99 --floor-blocks false --indel-size-to-eliminate-in-ref-model 10 --disable-optimizations false --dragen-mode false --dragen-378-concordance-mode false --flow-mode NONE --apply-bqd false --apply-frd false --disable-spanning-event-genotyping false --transform-dragen-mapping-quality false --mapping-quality-threshold-for-genotyping 20 --max-effective-depth-adjustment-for-frd 0 --just-determine-active-regions false --dont-genotype false --do-not-run-physical-phasing false --do-not-correct-overlapping-quality false --use-filtered-reads-for-annotations false --use-flow-aligner-for-stepwise-hc-filtering false --adaptive-pruning false --do-not-recover-dangling-branches false --recover-dangling-heads false --kmer-size 10 --kmer-size 25 --dont-increase-kmer-sizes-for-cycles false --allow-non-unique-kmers-in-ref false --num-pruning-samples 1 --min-dangling-branch-length 4 --recover-all-dangling-branches false --max-num-haplotypes-in-population 128 --min-pruning 2 --adaptive-pruning-initial-error-rate 0.001 --pruning-lod-threshold 2.302585092994046 --pruning-seeding-lod-threshold 9.210340371976184 --max-unpruned-variants 100 --linked-de-bruijn-graph false --disable-artificial-haplotype-recovery false --enable-legacy-graph-cycle-detection false --debug-assembly false --debug-graph-transformations false --capture-assembly-failure-bam false --num-matching-bases-in-dangling-end-to-recover -1 --error-correction-log-odds -Infinity --error-correct-reads false --kmer-length-for-read-error-correction 25 --min-observations-for-kmer-to-be-solid 20 --likelihood-calculation-engine PairHMM --base-quality-score-threshold 18 --dragstr-het-hom-ratio 2 --dont-use-dragstr-pair-hmm-scores false --pair-hmm-gap-continuation-penalty 10 --expected-mismatch-rate-for-read-disqualification 0.02 --pair-hmm-implementation FASTEST_AVAILABLE --pcr-indel-model CONSERVATIVE --phred-scaled-global-read-mismapping-rate 45 --disable-symmetric-hmm-normalizing false --disable-cap-base-qualities-to-map-quality false --enable-dynamic-read-disqualification-for-genotyping false --dynamic-read-disqualification-threshold 1.0 --native-pair-hmm-threads 4 --native-pair-hmm-use-double-precision false --flow-hmm-engine-min-indel-adjust 6 --flow-hmm-engine-flat-insertion-penatly 45 --flow-hmm-engine-flat-deletion-penatly 45 --pileup-detection false --use-pdhmm false --use-pdhmm-overlap-optimization false --make-determined-haps-from-pd-code false --print-pileupcalling-status false --fallback-gga-if-pdhmm-fails true --pileup-detection-enable-indel-pileup-calling false --pileup-detection-active-region-phred-threshold 0.0 --num-artificial-haplotypes-to-add-per-allele 5 --artifical-haplotype-filtering-kmer-size 10 --pileup-detection-snp-alt-threshold 0.1 --pileup-detection-indel-alt-threshold 0.1 --pileup-detection-absolute-alt-depth 0.0 --pileup-detection-snp-adjacent-to-assembled-indel-range 5 --pileup-detection-snp-basequality-filter 12 --pileup-detection-bad-read-tolerance 0.0 --pileup-detection-proper-pair-read-badness true --pileup-detection-edit-distance-read-badness-threshold 0.08 --pileup-detection-chimeric-read-badness true --pileup-detection-template-mean-badness-threshold 0.0 --pileup-detection-template-std-badness-threshold 0.0 --pileup-detection-filter-assembly-alt-bad-read-tolerance 0.0 --pileup-detection-edit-distance-read-badness-for-assembly-filtering-threshold 0.12 --bam-writer-type CALLED_HAPLOTYPES --dont-use-soft-clipped-bases false --override-fragment-softclip-check false --min-base-quality-score 10 --smith-waterman FASTEST_AVAILABLE --emit-ref-confidence NONE --max-mnp-distance 0 --force-call-filtered-alleles false --reference-model-deletion-quality 30 --soft-clip-low-quality-ends false --allele-informative-reads-overlap-margin 2 --smith-waterman-dangling-end-match-value 25 --smith-waterman-dangling-end-mismatch-penalty -50 --smith-waterman-dangling-end-gap-open-penalty -110 --smith-waterman-dangling-end-gap-extend-penalty -6 --smith-waterman-haplotype-to-reference-match-value 200 --smith-waterman-haplotype-to-reference-mismatch-penalty -150 --smith-waterman-haplotype-to-reference-gap-open-penalty -260 --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 --smith-waterman-read-to-haplotype-match-value 10 --smith-waterman-read-to-haplotype-mismatch-penalty -15 --smith-waterman-read-to-haplotype-gap-open-penalty -30 --smith-waterman-read-to-haplotype-gap-extend-penalty -5 --flow-assembly-collapse-hmer-size 0 --flow-assembly-collapse-partial-mode false --flow-filter-alleles false --flow-filter-alleles-qual-threshold 30.0 --flow-filter-alleles-sor-threshold 3.0 --flow-filter-lone-alleles false --flow-filter-alleles-debug-graphs false --min-assembly-region-size 50 --max-assembly-region-size 300 --active-probability-threshold 0.002 --max-prob-propagation-distance 50 --force-active false --assembly-region-padding 100 --padding-around-indels 75 --padding-around-snps 20 --padding-around-strs 75 --max-extension-into-assembly-region-padding-legacy 25 --max-reads-per-alignment-start 50 --enable-legacy-assembly-region-trimming false --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false --minimum-mapping-quality 20 --disable-tool-default-annotations false --enable-all-annotations false --allow-old-rms-mapping-quality-annotation-data false",Version="4.6.2.0",Date="2026? 7? 1? ?? 2? 30? 29? KST">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
##contig=<ID=chr20,length=64444167>
##source=HaplotypeCaller
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA12878
chr20	70882	.	T	C	47.32	.	AC=2;AF=1.00;AN=2;DP=3;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=41.73;QD=23.66;SOR=2.303	GT:AD:DP:GQ:PL	1/1:0,2:2:6:59,6,0
chr20	70968	.	C	CTAA	35.44	.	AC=2;AF=1.00;AN=2;DP=1;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=40.00;QD=25.36;SOR=1.609	GT:AD:DP:GQ:PL	1/1:0,1:1:3:45,3,0
chr20	70975	.	T	A	35.48	.	AC=2;AF=1.00;AN=2;DP=1;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=40.00;QD=28.73;SOR=1.609	GT:AD:DP:GQ:PL	1/1:0,1:1:3:45,3,0
chr20	70976	.	C	A	35.47	.	AC=2;AF=1.00;AN=2;DP=1;ExcessHet=0.0000;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=40.00;QD=30.97;SOR=1.609	GT:AD:DP:GQ:PL	1/1:0,1:1:3:45,3,0

 

성공적으로 변이 호출이 됨. 

VCF 파일이 정상적으로 생성됐고 실제 변이 데이터가 들어있음.

 

VCF 구조를 간단히 보면,

헤더 (##로 시작) : 각 필드가 뭘 의미하는지 정의해둔 부분. INFO, FORMAT 필드 설명이 다 여기 있음.

 

변이 레코드 (마지막 줄들)

chr20  70882  .  T  C  47.32  .  AC=2;AF=1.00;...  GT:AD:DP:GQ:PL  1/1:0,2:2:6:59,6,0

 

필드	값	의미
CHROM	chr20	어느 염색체
POS	70882	위치
REF	T	reference 염기
ALT	C	변이 염기
QUAL	47.32	변이 신뢰도 점수
GT	1/1	유전자형 (homozygous alt)

 

1/1은 양쪽 대립유전자 모두 변이형이라는 뜻, 즉 동형접합 변이.