[WGS 분석] VCF 파일 구조 이해하기

일반적인 VCF 구조

크게 세 부분으로 나뉨.

 

1. 메타정보 (##로 시작)

##fileformat=VCFv4.2
##INFO=<ID=DP,...>
##FORMAT=<ID=GT,...>
##GATKCommandLine=...

# 실제 VCF 파일의 메타 정보
##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

 

각 INFO, FORMAT 필드가 뭘 의미하는지 정의해두는 곳.

헤더 읽는 법만 알면 어떤 VCF든 해석 가능함.

 

2. 컬럼 헤더 (#CHROM으로 시작, 딱 1줄)

#CHROM  POS  ID  REF  ALT  QUAL  FILTER  INFO  FORMAT  NA12878

 

필드 이름을 정의하는 부분.
맨 끝이 샘플명.

 

3. 변이 레코드 (실제 데이터)

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

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

 

한 줄이 변이 하나를 나타냄.

필드 구성

• 고정 8개 필드 (CHROM~INFO) : 변이 자체 정보
• FORMAT : 샘플 정보 형식 정의
• 샘플 값 (NA12878) : 실제 샘플 데이터

---

variant calling 한 VCF 파일 읽기

gzcat ~/bioinfo/bi_study/wgs_practice/variants/SRR062634.vcf.gz | grep -v "^##" | head -5

 

위 명령어를 통해 밑의 결과물 출력됨. 

#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 파일 해석

일단 첫 번째 변이 chr20 70882를 기준으로 각 필드가 무엇을 뜻하는지 간단히 확인하면,

8개 고정 필드

필드	값	의미
CHROM	chr20	염색체
POS	70882	변이 위치
ID	.	dbSNP ID (없으면 .)
REF	T	reference 염기
ALT	C	변이한 염기
QUAL	47.32	변이 신뢰도 점수
FILTER	.	필터 통과 여부
INFO	AC=2;AF=1.00...	변이 요약 정보

그 뒤 두 필드

필드	값	의미
FORMAT	GT:AD:DP:GQ:PL	샘플 정보 형식 정의
NA12878	1/1:0,2:2:6:59,6,0	실제 샘플 값

 

이중에 INFO 필드랑 FORMAT/샘플 필드 두 가지가 핵심.

---

필드 상세 설명

CHROM

chr20

염색체 이름.

어느 염색체에서 변이가 발견됐는지 명시함.

POS

70882

변이 위치.

reference genome의 몇 번째 염기 위치인지.

1-based 좌표계 (1번부터 시작함).

ID

.

변이 식별자 . 이면 알려진 ID 없다는 뜻.

dbSNP에 등록된 변이면 rs123456789 같은 형태로 나옴.

나중에 어노테이션 단계에서 채워지는 필드.

REF

T

reference genome의 염기.

기준이 되는 서열.

ALT

C

실제 샘플에서 발견된 변이 염기.

 

REF랑 비교해서:

• REF랑 길이 같으면 → SNP (단일 염기 변이)
• ALT가 더 길면 → Insertion
• ALT가 더 짧으면 → Deletion

실제 데이터에서의 예시:

T → C       : SNP (염기 하나 바뀜)
C → CTAA    : Insertion (TAA 삽입)

QUAL

47.32

변이 호출 신뢰도로 Phred scale로 나타냄.

47.32면 오류일 확률이 약 0.002%.

보통 Q30 이상이면 신뢰할 만하다고 봄.

FILTER

.

품질 필터 결과.

• . → 필터 적용 안 됨
• PASS → 모든 필터 통과
• LowQual → 품질 낮음

지금 실제 데이터에서는 GATK에서 별도 필터링(VariantFiltration) 단계를 안 거쳤기 때문에 값이 . 임
나중에 필터링하면 PASS나 LowQual로 바뀜.

INFO

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

 

INFO 필드는 ;로 구분된 키=값 쌍.

 

유전자 변이(돌연변이)에 대한 추가 정보를 Key=Value(키=값) 형식으로 모아둔 8번째 데이터 열.

변이 빈도 관련

필드	값	의미
AC	2	ALT 대립유전자 개수
AF	1.00	ALT 빈도 (1.00 = 100%)
AN	2	전체 대립유전자 개수 (diploid라 2)

 

AF 1.00이면 이 샘플에서 reference 염기가 전혀 없고 변이만 있다는 뜻.

즉 동형접합 변이(homozygous).

커버리지 관련

필드	값	의미
DP	3	이 위치를 커버하는 read 수
MQ	41.73	평균 mapping quality
QD	23.66	QUAL을 DP로 나눈 값 (depth 보정 품질)

 

DP 3은 read가 3개밖에 없다는 것.

커버리지가 낮아서 신뢰도가 상대적으로 낮음.

편향(bias) 관련

필드	값	의미
FS	0.000	Strand bias. 0에 가까울수록 좋음
SOR	2.303	Strand bias 다른 측정 방식
ExcessHet	0.0000	과도한 이형접합 여부

 

FS랑 SOR은 변이가 특정 방향 read에서만 나오는지 체크함.

한쪽 방향에서만 변이가 보이면 진짜 변이가 아니라 시퀀싱 아티팩트(실제 생물학적 샘플에 존재하지 않는, 유전체 분석(시퀀싱) 과정의 오류나 기술적 한계 때문에 데이터에 가짜로 생기는 흔적)일 가능성이 높음.

FORMAT, SAMPLE

GT:AD:DP:GQ:PL	1/1:0,2:2:6:59,6,0

실제 예시에서 FORMAT 필드는 GT:AD:DP:GQ:PL이고 샘플 값은 1/1:0,2:2:6:59,6,0 인데 FORMAT 필드와 SAMPLE 필드가  1:1 대응임.

GT (Genotype) = 1/1

FORMAT, SAMPLE에서 가장 중요한 필드.

이 샘플의 유전자형을 나타냄.

값 의미

값	의미
0/0	homozygous reference (변이 없음)
0/1	heterozygous (한쪽만 변이)
1/1	homozygous alt (양쪽 다 변이)

 

1/1이면 diploid(2배체)인 사람에서 양쪽 염색체 모두 변이형이라는 뜻.

AD (Allelic Depth) = 0,2

REF read 수와 ALT read 수.

0,2면 reference 읽은 read 0개, 변이 읽은 read 2개라는 뜻.

GT 1/1과 값이 일치함.

DP (Depth) = 2

이 위치 총 read 수.
AD 합이랑 같아야 힘.

0+2=2.

GQ (Genotype Quality) = 6

유전자형 신뢰도(Phred scale)
6이면 낮은 편.
DP가 2밖에 안 되니까 당연하게 생각할 수 있음.

PL (Phred-scaled Likelihoods) = 59,6,0

각 유전자형일 가능성을 Phred scale로 표현한 것.

순서는 0/0, 0/1, 1/1.

0/0일 가능성 → 59 (매우 낮음)
0/1일 가능성 → 6  (낮음)
1/1일 가능성 → 0  (가장 높음)

 

값이 0에 가까울수록 그 유전자형일 가능성이 높음.

그래서 1/1로 판정된 것.