Step 04. Post-processing (MarkDuplicates & BQSR)

정렬된 BAM에서 중복 리드를 마킹하고(MarkDuplicates), 알려진 변이 사이트를 이용해 BQSR(Base Quality Score Recalibration)을 수행합니다. 기본 쓰레드는 20을 사용합니다.

A) 왜 이 단계가 필요할까?

중복 마킹: PCR/광학적 중복은 같은 분자에서 유래한 리드로, 변이 콜링에 편향을 줍니다. MarkDuplicates는 이를 식별해 FLAG로 표시(제거 X)하고, 이후 분석에서 자동으로 제외되도록 합니다.
BQSR: 염기 품질(Q)은 기기/배치/염기 맥락에 따라 체계적인 오차가 납니다. BaseRecalibrator → ApplyBQSR로 품질 점수를 보정하여 변이 콜 품질을 향상시킵니다.

A-1) BQSR(품질 보정)의 작동 원리

알려진 변이 사이트(known sites) + 경험적 오류확률 → 새 품질 점수

목적: 기기/배치/염기맥락 때문에 생기는 체계적 오류를 보정.
오류 집계: BaseRecalibrator가 리드–레퍼런스 불일치를 세는데, 알려진 변이 위치는 “진짜 변이”로 보고 제외합니다. (진짜 변이를 오류로 세면 모델이 틀어짐)
경험적 오류확률 추정:

p_emp = (불일치 수) / (검사한 염기 수)

Q_emp = -10 · log10(p_emp)
공변량(covariates)별 보정: 리드그룹/사이클(염기 위치)/염기 맥락(k-mer)/보고된 Q 등을 계층적 가산 모델로 학습하여 ΔQ를 추정하고, ApplyBQSR 단계에서 새 품질로 바꿉니다.
직관: “기계가 Q30이라 했지만 이 조건에선 경험적으로 Q27쯤이네 → 조정”.

결론: BQSR은 경험적 빈도 → 오류확률 → Phred Q로 환산하는 확률적(통계적) 보정이며, known sites는 “오류 집계에서 빼주는 마스크” 역할을 합니다.

B) 입력 & 전제

이전 단계 산출물: aln/SAMPLE.sorted.bam (+ .bai).
레퍼런스: Homo_sapiens_assembly38.fasta (+ .fai, .dict).
Known sites(예, hg38): dbsnp_146.hg38.vcf.gz, Mills_and_1000G_gold_standard.indels.hg38.vcf.gz, 1000G_phase1.snps.high_confidence.hg38.vcf.gz (각각 .tbi 인덱스 필요).

C) 공통 변수 설정


# 스레드 & 경로
THREADS=20
REF_DIR="/home/kang/raw_data/gatk_bundle/v0"
REF="${REF_DIR}/Homo_sapiens_assembly38.fasta"

IN="/home/kang/God_Nas/WGS_2nd/aln"          # 03단계 출력 위치
OUT="/home/kang/God_Nas/WGS_2nd/postproc"    # 04단계 출력 위치
mkdir -p "${OUT}"

# Known sites (gzip + .tbi 필요)
DBSNP="${REF_DIR}/dbsnp_146.hg38.vcf.gz"
MILLS="${REF_DIR}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
G1KHC="${REF_DIR}/1000G_phase1.snps.high_confidence.hg38.vcf.gz"

D) 단일 샘플 예제


SAMPLE="NP-21"

# 1) MarkDuplicates (Spark 버전: 빠름)
gatk MarkDuplicatesSpark \
  -I "${IN}/${SAMPLE}.sorted.bam" \
  -O "${OUT}/${SAMPLE}.dup.bam" \
  --metrics-file "${OUT}/${SAMPLE}.dup.metrics.txt" \
  --conf 'spark.executor.cores='"${THREADS}" \
  --conf 'spark.local.dir=/tmp' \
  --create-output-bam-index true

# (대안) Picard MarkDuplicates (단일 스레드/호환성 높음)
# gatk MarkDuplicates \
#   -I "${IN}/${SAMPLE}.sorted.bam" \
#   -O "${OUT}/${SAMPLE}.dup.bam" \
#   -M "${OUT}/${SAMPLE}.dup.metrics.txt" \
#   --CREATE_INDEX true

# 2) BaseRecalibrator (리포트 생성)
gatk BaseRecalibrator \
  -R "${REF}" \
  -I "${OUT}/${SAMPLE}.dup.bam" \
  --known-sites "${DBSNP}" \
  --known-sites "${MILLS}" \
  --known-sites "${G1KHC}" \
  -O "${OUT}/${SAMPLE}.recal.table"

# 3) ApplyBQSR (보정 적용 → 최종 BAM)
gatk ApplyBQSR \
  -R "${REF}" \
  -I "${OUT}/${SAMPLE}.dup.bam" \
  --bqsr-recal-file "${OUT}/${SAMPLE}.recal.table" \
  -O "${OUT}/${SAMPLE}.final.bam"

# 인덱스 & QC
samtools index -@ ${THREADS} "${OUT}/${SAMPLE}.final.bam"
samtools flagstat -@ ${THREADS} "${OUT}/${SAMPLE}.final.bam" > "${OUT}/${SAMPLE}.final.flagstat.txt"

E) 주요 옵션 설명

MarkDuplicatesSpark: Spark 로컬 모드로 병렬 수행. --metrics-file에 중복 통계가 기록되며, --create-output-bam-index로 자동 색인 생성.
BaseRecalibrator: 알려진 변이 사이트에서 실제 변이는 오류로 학습하지 않도록 제외하고 시스템적 오차 모델을 학습해 recal.table 생성.
ApplyBQSR: 학습된 보정 테이블을 적용하여 염기 품질(Q) 재계산.
Known-sites: dbSNP(다양한 SNP), Mills(인델), 1000G high-conf SNP 등(레퍼런스 빌드 일치·.tbi 필수).

F) 배치 실행 스크립트 (`run_postproc.sh`)

03단계에서 생성한 aln/*.sorted.bam을 모두 찾아 일괄 처리합니다.


#!/usr/bin/env bash
set -euo pipefail

THREADS=20
REF_DIR="/home/kang/raw_data/gatk_bundle/v0"
REF="${REF_DIR}/Homo_sapiens_assembly38.fasta"

IN="/home/kang/God_Nas/WGS_2nd/aln"
OUT="/home/kang/God_Nas/WGS_2nd/postproc"
mkdir -p "${OUT}"

DBSNP="${REF_DIR}/dbsnp_146.hg38.vcf.gz"
MILLS="${REF_DIR}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
G1KHC="${REF_DIR}/1000G_phase1.snps.high_confidence.hg38.vcf.gz"

shopt -s nullglob
for BAM in "${IN}"/*.sorted.bam; do
  SAMPLE="$(basename "${BAM}" .sorted.bam)"
  echo "[POSTPROC] ${SAMPLE}"

  gatk MarkDuplicatesSpark \
    -I "${BAM}" \
    -O "${OUT}/${SAMPLE}.dup.bam" \
    --metrics-file "${OUT}/${SAMPLE}.dup.metrics.txt" \
    --conf 'spark.executor.cores='"${THREADS}" \
    --conf 'spark.local.dir=/tmp' \
    --create-output-bam-index true

  gatk BaseRecalibrator \
    -R "${REF}" \
    -I "${OUT}/${SAMPLE}.dup.bam" \
    --known-sites "${DBSNP}" \
    --known-sites "${MILLS}" \
    --known-sites "${G1KHC}" \
    -O "${OUT}/${SAMPLE}.recal.table"

  gatk ApplyBQSR \
    -R "${REF}" \
    -I "${OUT}/${SAMPLE}.dup.bam" \
    --bqsr-recal-file "${OUT}/${SAMPLE}.recal.table" \
    -O "${OUT}/${SAMPLE}.final.bam"

  samtools index -@ ${THREADS} "${OUT}/${SAMPLE}.final.bam"
  samtools flagstat -@ ${THREADS} "${OUT}/${SAMPLE}.final.bam" > "${OUT}/${SAMPLE}.final.flagstat.txt"
done

echo "All done. Final BAMs in: ${OUT}"

요약

MarkDuplicates로 PCR/광학적 중복을 마킹, BQSR로 체계적 품질 오차를 보정.
입력: sorted.bam → 출력: dup.bam → final.bam(+.bai).
Known-sites는 레퍼런스 빌드 일치와 .tbi 인덱스 필요.
BQSR은 p_emp, Q_emp 기반의 확률적 보정이며 known sites는 오류 집계에서 제외하는 마스크.
QC: *.dup.metrics.txt·flagstat로 중복률/매핑 요약 확인.

← Back: Alignment Next: Variant Calling → GATK Resource Bundle ↗ GATK Tools Docs ↗