GATK BQSR 완전 이해 — Phred 점수를 왜 재보정하는지, 알려진 사이트가 왜 필요한지

Phred 점수는 왜 부정확한가요

지난 편에서 PCR 중복을 잡았으니 이제 각 리드가 독립 관찰이라는 가정이 어느 정도 성립합니다. 다음 단계는 각 관찰이 얼마나 정확한가에 대한 재평가 — BQSR (Base Quality Score Recalibration)입니다.

시퀀서가 리드 하나하나에 붙여서 뱉는 Phred 점수는 광학·화학 신호에서 유도한 사전 확률입니다. 완벽하지 않습니다. 특히 다음 세 요인에 대해 시스템적으로 편향됩니다.

리드 사이클 위치: 앞부분보다 뒷부분이 오류율이 급증하는데 Phred는 이걸 충분히 반영 못 함.
염기 문맥 (dinucleotide context): 특정 두 염기 조합(예: GG, CG) 뒤에서 오류율이 다름.
기존 Phred bin: Phred가 보고한 Q값과 실제 오류율의 매핑이 어긋남.

BQSR은 이 편향을 dbSNP를 참 값으로 삼아 데이터 주도로 보정합니다. 어떻게 하는지 통계로 유도해봅시다.

핵심 아이디어 — 참·거짓의 근사

각 리드 위치에서 "이 염기가 참조와 다르다"는 이벤트가 관찰됩니다. 이 이벤트가 두 유형입니다.

진짜 변이 (variant): 이 개인의 실제 SNP. dbSNP나 알려진 변이 카탈로그에 있음.
시퀀싱 오류 (error): 랜덤 노이즈.

dbSNP에 알려진 위치에서 참조와 다른 것은 변이로 간주하고, 알려지지 않은 위치에서 참조와 다른 것을 오류로 간주하는 근사가 BQSR의 출발입니다. 근사이지만 dbSNP가 매우 촘촘해서 이 근사가 대체로 잘 맞습니다.

\text{empirical error rate} = \frac{\text{alt at unknown sites}}{\text{total bases}}

이 경험적 오류율에서 Phred 점수 Q_empirical을 역산합니다.

Q_{\text{empirical}} = -10 \cdot \log_{10}(\text{error rate})

다변량 층화

경험적 Q를 그냥 하나의 값으로만 계산하면 정보를 잃습니다. BQSR은 네 공변량으로 층화합니다.

원 Phred bin (reported Q): Q10, Q20, Q30 각 bin별로 따로.
리드 사이클 위치: 1～150bp 사이 위치별로.
염기 문맥: 앞 염기 X 현재 염기 조합 (16개).
머신·플로우셀 그룹 (선택).

즉 (reported Q, cycle, context) 조합의 각 셀에서 경험적 오류율을 따로 계산합니다. 그 결과 Q_empirical이 셀마다 다른 값이 나옵니다.

재보정 함수

각 리드 위치의 새 Phred는 이렇게 결정됩니다.

Q_{\text{recal}} = Q_{\text{reported}} + \Delta Q_{\text{global}} + \Delta Q_{\text{cycle}} + \Delta Q_{\text{context}}

각 Δ는 해당 층에서의 경험적 오차와 원 Phred의 차이.
로그 스케일이라 덧셈으로 표현.

이 함수를 리드 파일 전체에 적용해 새 BAM을 만듭니다.

GATK4 명령 — 두 단계

BQSR은 두 단계로 실행합니다.

1단계: 재보정 테이블 학습

bash

gatk BaseRecalibrator \
    -I S1.dedup.bam \
    -R hg38.fa \
    --known-sites resources/dbsnp_146.hg38.vcf.gz \
    --known-sites resources/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    -O S1.recal.table

--known-sites: dbSNP + 알려진 인델 (Mills, 1000G).
이 단계가 무거워도 인터벌(-L intervals.bed)로 스코프 좁히면 빨라집니다.

2단계: 학습된 테이블로 BAM 재작성

bash

gatk ApplyBQSR \
    -I S1.dedup.bam \
    -R hg38.fa \
    --bqsr-recal-file S1.recal.table \
    -O S1.recal.bam
samtools index S1.recal.bam

이제 S1.recal.bam이 다음 단계 HaplotypeCaller의 입력이 됩니다.

3단계 (선택): 재보정 효과 시각화

bash

gatk BaseRecalibrator \
    -I S1.recal.bam \
    -R hg38.fa \
    --known-sites resources/dbsnp_146.hg38.vcf.gz \
    -O S1.recal_after.table

gatk AnalyzeCovariates \
    -before S1.recal.table \
    -after S1.recal_after.table \
    -plots S1.bqsr_plots.pdf

PDF에 재보정 전후의 Q값 정렬도가 시각화됩니다. 이 그래프가 우측 대각선에 가까이 붙어야 정상. 붙지 않으면 known sites 파일 문제 또는 라이브러리 결함 신호입니다.

dbSNP · 알려진 인델 — 왜 이 두 세트인가

BQSR의 known-sites 인자에 넣는 표준 파일이 셋입니다.

파일	내용	출처
dbSNP 146+	알려진 SNP 위치	NCBI
Mills · 1000G Gold Standard Indels	검증된 인델	Broad + 1000G
1000G Phase 1 High Confidence SNPs	인구 SNP	1000G

Mills가 왜 별도인지: dbSNP는 SNP 위주, 인델은 검증 어려워 소규모. BQSR은 인델 근처가 특히 편향적이라 별도 파일이 필요합니다.

중요: 참조 게놈 버전(hg19 vs hg38)과 known-sites 파일 버전이 반드시 일치해야 합니다. hg19 dbSNP를 hg38 BAM에 쓰면 좌표가 밀려서 전체 재보정이 무너집니다. Resource Bundle에서 참조 버전에 맞는 세트를 받읍시다.

손 계산 — 왜 BQSR이 필요한가

가상의 시나리오를 봅시다. Illumina가 보고한 Q20 bin의 리드 위치 30번(리드 뒷부분)에서, dbSNP 미등록 위치에 참조 불일치가 10만 위치에 발생. 그 위치의 총 관찰 수 = 10,000,000.

관찰 오류율: 100,000 / 10,000,000 = 0.01
실제 Q_empirical: −10 log(0.01) = Q20

이 경우 Phred가 이미 잘 맞으므로 Δ가 거의 0. 재보정이 필요 없습니다.

이제 다른 시나리오. 같은 Q20 bin이지만 리드 위치 149번(리드 끝)에서 참조 불일치가 100만 위치. 총 관찰 = 10,000,000.

관찰 오류율: 1,000,000 / 10,000,000 = 0.1
실제 Q_empirical: −10 log(0.1) = Q10

Phred는 Q20이라 했는데 실제로는 Q10입니다. BQSR은 이 리드 위치의 Phred를 Q10으로 낮춥니다. 그 결과 HaplotypeCaller가 이 위치의 관찰을 덜 신뢰하고, 위양성 변이를 줄입니다.

이게 BQSR의 실무 효과입니다.

결과 검증 — 세 지표

BQSR 실행 후 아래 세 개를 확인합니다.

1. AnalyzeCovariates 그래프

Q_reported vs Q_empirical 산점도가 대각선에 붙어야 정상. 재보정 후 대각선에 훨씬 가까워졌다면 성공.

2. 재보정 전후 파일 크기

BQSR은 Phred만 갈아주므로 파일 크기가 크게 안 변합니다. 갑자기 절반으로 줄면 리드 유실 사고 의심.

3. 재보정 후 HaplotypeCaller 결과의 Ti/Tv 비율

Ti/Tv (transition/transversion ratio): 인간 WGS에서 2.0～2.1이 정상.
BQSR 없이 콜하면 오류가 랜덤에 가까워 Ti/Tv가 1.5 근처로 떨어짐.
BQSR 후 콜은 Ti/Tv가 2.0에 붙어야 함.

Ti/Tv가 왜 지표인가

DNA 돌연변이는 두 유형이 있습니다.

Transition (Ti): purine ↔ purine (A↔G) 또는 pyrimidine ↔ pyrimidine (C↔T). 화학적으로 잘 일어남.
Transversion (Tv): purine ↔ pyrimidine (A↔C, A↔T, G↔C, G↔T). 상대적으로 덜 일어남.

자연 돌연변이는 Ti/Tv ≈ 2.0. 랜덤 오류는 Ti/Tv ≈ 0.5. 이 비율이 콜된 변이 세트의 자연성을 진단합니다.

실무 팁

BQSR을 스킵하는 경우

매우 작은 표적 패널 (dbSNP 커버리지 부족).
인간 이외 종에서 known-sites가 부재.
이 경우 GATK가 --skip-annotations로 스킵 가능하나 Best Practices 이탈.

CRAM에서 BQSR

CRAM은 참조 기반 압축이라 BQSR 결과를 재저장할 때 참조 버전 유지 필수. 실무는 BAM으로 처리 후 마지막에 CRAM 변환.

인터벌 파일 활용

전체 게놈이 아닌 표적 영역만 재보정하려면 -L intervals.bed. 학습이 훨씬 빠릅니다.

Colab에서 소형 실습

Colab에서 chr22만 대상으로 GATK Docker를 부르면 소형 재보정이 가능합니다.

bash

!docker pull broadinstitute/gatk:latest
# chr22 BAM + chr22 dbSNP 슬라이스로 축소해서 실행
!docker run --rm -v $(pwd):/data broadinstitute/gatk:latest \
    gatk BaseRecalibrator \
    -I /data/chr22.bam \
    -R /data/hg38_chr22.fa \
    --known-sites /data/dbsnp_chr22.vcf.gz \
    -O /data/chr22.recal.table

용량이 크면 콜랩에서는 안 될 수 있어 GCP Batch나 Slurm으로 옮겨가는 게 실무.

CS 매핑

경험적 오류율 학습: 데이터로부터 실제 확률을 재계산하는 empirical Bayes 정신.
다변량 층화: 여러 공변량 조합으로 세분화하여 편향 캡처.
로그 확률 표기: Phred = −10 log P의 그 로그 스케일 그대로.

DryBench "경험적 오류 모델링" 편(empirical-error-modeling) 참고.

마무리

BQSR은 눈에 보이는 결과보다는 뒷단 콜링의 통계 가정을 만족시키는 목적의 도구입니다. AnalyzeCovariates 그래프가 대각선에 붙는 순간, 우리 파이프라인이 HaplotypeCaller의 자연 지형에 서게 됩니다. 다음 편(S09)에서 그 HaplotypeCaller가 어떻게 로컬 어셈블리와 Bayesian 콜링을 조합하는지 파고들어 봅시다.

더 깊게 파고 싶다면

DePristo M. et al. (2011), A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43:491. BQSR 원 방법.
GATK 공식 BQSR 가이드 (https://gatk.broadinstitute.org/hc/en-us/articles/360035890531): 옵션 원전.
Broad Institute BroadE — Base Quality Score Recalibration deep dive 유튜브 강의. AnalyzeCovariates 그래프의 임상 해석 예시.
Wall J.D. et al. (2014), Estimating genotype error rates from high-coverage next-generation sequence data. Genome Research 24:1734. 실제 오류율 추정 연구.
NCBI dbSNP 공식 (https://www.ncbi.nlm.nih.gov/snp/): known-sites 원전.

Ti/Tv 지표가 뜨는 순간 파이프라인 진단이 훨씬 편해집니다. 다음 편에서 HaplotypeCaller의 로컬 어셈블리로 갑시다.

BioPlayground

BQSR: Phred 점수 재보정의 통계 유도와 실무

Phred 점수는 왜 부정확한가요

핵심 아이디어 — 참·거짓의 근사

다변량 층화

재보정 함수

GATK4 명령 — 두 단계

1단계: 재보정 테이블 학습

2단계: 학습된 테이블로 BAM 재작성

3단계 (선택): 재보정 효과 시각화

dbSNP · 알려진 인델 — 왜 이 두 세트인가

손 계산 — 왜 BQSR이 필요한가

결과 검증 — 세 지표

1. AnalyzeCovariates 그래프

2. 재보정 전후 파일 크기

3. 재보정 후 HaplotypeCaller 결과의 Ti/Tv 비율

Ti/Tv가 왜 지표인가

실무 팁

BQSR을 스킵하는 경우

CRAM에서 BQSR

인터벌 파일 활용

Colab에서 소형 실습

CS 매핑

마무리

더 깊게 파고 싶다면