필터가 왜 별도 단계인가요

지난 두 편에서 HaplotypeCaller와 Mutect2로 변이 후보를 콜했습니다. 이 raw 콜은 위양성이 상당히 섞여 있습니다. 콜러가 잡을 수 있는 만큼 잡되, 어느 게 진짜인지 판단을 다음 단계로 미룬 것이죠.

이번 편이 그 다음 단계입니다. Broad Best Practices는 두 방향을 제시합니다.

• VQSR (Variant Quality Score Recalibration): 통계적으로 정교한 방법. 최소 샘플 수 조건 있음.
• Hard Filtering: 임계값 기반 단순 필터. 아무 데서나 씀.

두 방법의 선택 규칙과 실무 옵션을 정리합니다.

VQSR — GaussianMixtureModel의 미학

VQSR의 핵심 아이디어를 한 문장으로:

알려진 참 변이 세트(HapMap, 1000G)를 훈련 데이터로 삼아 각 변이가 참일 확률을 다변량 통계로 학습.

특징 벡터

각 변이 후보를 아래 통계로 표현합니다.

• QD (Quality by Depth): 커버리지로 정규화한 신뢰도.
• FS (Fisher Strand): 스트랜드 편향.
• SOR (Strand Odds Ratio): 대안 스트랜드 편향.
• MQ (Mean MAPQ): 리드 매핑 품질 평균.
• MQRankSum: 참조 vs 대안 리드의 MAPQ 차이.
• ReadPosRankSum: 리드 안 위치 편향.

이 6～7차원 벡터에서 참 변이들은 특정 지역에 몰려 있고, 위양성은 다른 지역에 몰려 있습니다.

GMM 훈련

VQSR은 이 벡터 공간에서 Gaussian Mixture Model (GMM) 두 개를 훈련합니다.

• Positive GMM: HapMap · 1000G 참 세트로 훈련.
• Negative GMM: PoN처럼 반복적 아티팩트로 훈련.

각 콜에 대해 두 GMM의 로그 우도 비율(VQSLOD)이 신뢰도가 됩니다. VQSLOD가 높을수록 참일 확률이 큰 것.

Sensitivity Tranche

VQSR 사용자가 결정할 것은 하나. "참 세트에 대한 감도를 어디로 자를까?"

• 99.7% tranche: 참 세트의 99.7%를 유지. 위양성이 상대적으로 많음.
• 99.0% tranche: 참 세트 99.0% 유지. 위양성 감소.
• 90.0% tranche: 매우 엄격.

일반적으로 SNP는 99.5%, INDEL은 99.0%를 씁니다. tranche 곡선(SensitivityVSFalsePositive)을 보고 팀 정책에 맞춰 결정.

실행 명령

두 단계로 실행.

# 1) 훈련
gatk VariantRecalibrator \
    -R hg38.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=true,prior=12 omni.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10 1000G.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2 dbsnp.vcf.gz \
    -an QD -an FS -an SOR -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.snp.recal \
    --tranches-file cohort.snp.tranches

# 2) 적용
gatk ApplyVQSR \
    -R hg38.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.snp.recal \
    --tranches-file cohort.snp.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.snp.recalibrated.vcf.gz

INDEL은 별도 실행 (mode INDEL).

VQSR을 안 쓸 때 — Hard Filtering

VQSR은 GMM 훈련이 안정적으로 이루어지려면 관례상 30명 이상 샘플 코호트가 필요합니다. 그 미만이면 GMM이 노이즈에 과적합해서 오히려 필터 품질이 떨어집니다.

• 단일 샘플 진단 임상.
• 20명 이하 소규모 코호트.
• 인간 이외 종 (참 세트 부족).

이 경우 Hard Filter가 표준입니다.

GATK 권장 Hard Filter 임계

SNP:

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

INDEL:

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

Broad가 30명 이상 코호트로 실증해서 뽑은 임계값입니다. 이 값들은 인간·hg38 기준.

VQSR vs Hard Filter — 결정 규칙

손 계산 — 왜 GMM인가

VQSR을 GMM 대신 그냥 임계 필터로 처리하면 왜 안 될까요? 이유는 특징이 상관되어 있기 때문입니다.

• 커버리지가 높은 구간은 QD도 높고 FS도 낮음.
• 반복 서열 구간은 MQ가 낮고 SOR도 높음.

즉 개별 임계를 나열하면 "커버리지가 낮은 정상 변이"가 놓치거나 "반복 서열의 아티팩트"가 잡히지 않습니다. GMM은 이 상관관계를 자동으로 학습합니다.

가상 예: 6차원 벡터 공간에서 참 변이가 (QD=15, FS=1, SOR=1, MQ=60, MQRankSum=0, ReadPosRankSum=0) 근처에 몰려 있고, 위양성이 (QD=3, FS=50, SOR=8, MQ=35, MQRankSum=-10, ReadPosRankSum=-5) 근처에 몰려 있으면, GMM은 이 두 클러스터 사이 결정 경계를 찾습니다. 새 콜은 이 경계 기준으로 분류.

Hard Filter는 이 경계를 축 정렬 상자로 근사합니다. 감도가 낮아지지만 해석이 명료.

Joint Genotyping — VQSR 전 필수 단계

VQSR은 단일 샘플 VCF가 아니라 다중 샘플 통합 VCF를 요구합니다. 이 통합이 Joint Genotyping.

# 1) GVCF들을 GenomicsDB에 통합
gatk GenomicsDBImport \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.bed \
    -V S1.g.vcf.gz -V S2.g.vcf.gz ... -V S30.g.vcf.gz

# 2) 다중 샘플 콜
gatk GenotypeGVCFs \
    -R hg38.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz

이 과정에서 각 위치별로 모든 샘플의 정보가 합쳐지고 인구 알렐 빈도가 계산됩니다. 결과 VCF가 VQSR 입력.

VQSR 실패 시 진단

VQSR이 실패하는 흔한 이유.

1. 훈련 변이 부족: 참 세트(HapMap)에 해당하는 콜이 500 미만이면 GMM이 훈련되지 않음.
2. tranche 곡선 이상: sensitivity-fp 곡선이 계단형이면 훈련 데이터 이상.
3. INDEL이 SNP보다 훨씬 나쁨: INDEL은 항상 SNP보다 어려움. tranche를 99.0%로 낮춤.

VQSR 실패하면 Hard Filter로 fallback이 표준.

VCF 최종 처리 파이프라인

VQSR/Hard 필터 후 최종 VCF는 이렇게 다듬습니다.

# 1) 필터 통과한 변이만
bcftools view -f PASS cohort.vcf.gz -Oz -o cohort.pass.vcf.gz

# 2) 정규화 (좌측 정렬 · 최소 표현)
bcftools norm -f hg38.fa cohort.pass.vcf.gz -Oz -o cohort.norm.vcf.gz

# 3) 인덱스
bcftools index cohort.norm.vcf.gz

# 4) 통계
bcftools stats cohort.norm.vcf.gz > cohort.stats.txt

bcftools norm이 자주 놓치는 스텝인데, 인델 표기 방식(좌측 정렬)이 도구마다 다르니 표준화가 필수. 이 정규화 없이 두 VCF를 비교하면 같은 변이가 다르다고 나옵니다.

CS 매핑

• GMM (Gaussian Mixture Model): DryBench "GMM 기초" 편(gmm-basics) 참고. 다변량 밀도 추정의 정수.
• 다변량 상관관계 처리: 특징 간 상관을 covariance matrix로 자동 캡처.
• 임계 결정 규칙: sensitivity tranche = ROC curve 위의 point selection.

마무리

VQSR과 Hard Filter는 둘 다 도구이고, 상황에 맞춰 골라 씁니다. 30명 조건을 외워두면 프로젝트 초기 결정이 수월합니다. 다음 편(S12)에서는 GATK 프레임을 벗어나 DeepVariant가 딥러닝으로 변이를 콜하는 새로운 접근을 봅니다.

이 편으로 GATK4 파이프라인 지도(S06～S11)가 완성됩니다. 이후 S12～S16은 대안 콜러와 구조변이·CNV·임퓨테이션으로, S17～S21은 RNA-seq 정량과 DGEA로 이어집니다. "FASTQ to Paper" 시리즈 1부(S06～S11 GATK4)가 여기서 마무리됩니다.

더 깊게 파고 싶다면

• DePristo M. et al. (2011), A framework for variation discovery and genotyping. Nature Genetics 43:491. VQSR의 원 방법.
• GATK 공식 VQSR 가이드 (https://gatk.broadinstitute.org/hc/en-us/articles/360035531612): 옵션·resource 원전.
• GATK Hard Filter 가이드 (https://gatk.broadinstitute.org/hc/en-us/articles/360035531112): 임계값 실증 유래.
• Broad Institute BroadE — Variant Filtering 유튜브 강의. tranche 곡선 해석 예시.
• bcftools 공식 문서 (https://samtools.github.io/bcftools/): 정규화·통계 옵션 원전.

VQSR과 Hard Filter의 선택 규칙만 외워둬도 임상 파이프라인 결정이 절반은 끝납니다. 다음 편에서 DeepVariant로 갑시다.