BioPlayground

🧬
목록으로

GATK MarkDuplicates: PCR 중복은 어떻게 판정되고 왜 표시만 하나요

PCR 증폭 편향을 통계적으로 왜 통제해야 하는지, MarkDuplicates가 어떤 규칙으로 duplicate를 판정하는지, 왜 실제로 삭제하지 않고 플래그로 표시만 하는지 정리합니다.

중급
|
18분
|
검증 완료 (2026-07-22)
PCR duplicatesMarkDuplicateslibrary complexity
진행률0/52 (0%)

PCR 중복이 왜 문제인가요

지난 편 지도에서 MarkDuplicates가 파이프라인 3번째에 위치한 이유를 짧게 짚었습니다. 이번 편에서 그 이유를 통계적으로 풀어봅니다.

라이브러리 준비 시 원본 DNA 조각(insert)을 시퀀싱할 만큼 증폭시키기 위해 PCR을 20～30 사이클 돌립니다. 이 증폭은 원리적으로 편향적입니다.

  • GC 편향: GC가 극단인 구간이 덜 증폭됩니다.
  • 길이 편향: 짧은 인서트일수록 더 잘 증폭됩니다.
  • 초기 카피 수 편향: 초반 몇 사이클에 어떤 원본이 우연히 더 증폭되면 그 편향이 지수적으로 증폭됩니다.

결과적으로 시퀀싱 리드 안에는 원본 인서트 하나가 수백 번 복제된 사본들이 섞여 있습니다. 이 사본들을 각각 독립 관찰로 세면 커버리지 통계가 왜곡되고, 변이 콜링이 잘못된 확신을 얻습니다. MarkDuplicates가 이 편향을 통제하는 첫 관문입니다.

어떤 리드가 duplicate로 판정되나요

MarkDuplicates 규칙은 놀랍도록 단순합니다.

같은 참조 좌표(같은 위치·같은 스트랜드·같은 CIGAR 시작점)를 가진 리드들은 하나의 원본에서 유래한 duplicate로 간주.

Paired-end이면 R1의 시작점 + R2의 시작점이 같아야 합니다. R1 좌표만 같고 R2 좌표가 다르면 다른 원본입니다.

대표 선택

같은 좌표 그룹에서 한 리드만 원본으로 남기고 나머지는 duplicate로 표시합니다. 대표를 뽑는 기준은 총 base quality 합계가 가장 높은 리드입니다. 즉 품질이 가장 좋은 사본을 살리고 나머지에 플래그 1024를 붙입니다.

  • 코드 안 규칙: flag |= 0x400 (=1024).
  • 리드 자체는 파일에서 삭제되지 않고 그대로 남습니다.
  • HaplotypeCaller와 같은 뒷단이 이 플래그를 보고 통계에서 제외합니다.

왜 삭제 안 하고 표시만 하나요

세 가지 이유가 있습니다.

  1. 되돌릴 수 있음: 나중에 다른 알고리즘으로 재평가하고 싶으면 플래그만 지우면 됩니다.
  2. QC 지표 계산: duplication rate를 계산하려면 duplicate 리드도 파일에 남아야 합니다.
  3. 일부 분석은 duplicate 활용: RNA-seq이나 아주 저커버리지 분석에서는 duplicate를 오히려 신호로 씁니다.

이 "표시만 하기" 관례가 실무 파이프라인 유연성의 기반입니다.

Optical duplicates — 광학 중복의 정체

같은 좌표에 붙는 리드들 중, 플로우셀 상 물리적으로 근접한 리드들이 있습니다. 이건 사실 하나의 클러스터가 시퀀서가 두 클러스터로 잘못 판독한 것입니다. PCR 중복이 아니라 광학 중복(optical duplicate)입니다.

MarkDuplicates는 read name(예: HWI-ST177:290:C0TECACXX:1:1101:1225:2130)에서 x, y 좌표를 추출해 두 리드가 100픽셀 이내면 optical duplicate로 별도 분류합니다.

text
riD_TILE : 1101, x: 1225, y: 2130
riD_TILE : 1101, x: 1230, y: 2135   # 근접 → optical dup 후보

OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 옵션이 그 임계입니다. NovaSeq X 같은 최신 플랫폼은 2500 정도로 늘려야 합니다. 이 값을 잘못 두면 optical dup 카운트가 이상해집니다.

실행 명령

bash
gatk MarkDuplicates \
    -I S1.sorted.bam \
    -O S1.dedup.bam \
    -M S1.metrics.txt \
    --VALIDATION_STRINGENCY LENIENT \
    --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
    --TMP_DIR /data/tmp
samtools index S1.dedup.bam

옵션 세 개 정도만 잡으면 됩니다.

  • -M metrics.txt: 리포트 파일 (라이브러리별 duplication rate 등).
  • --VALIDATION_STRINGENCY LENIENT: SAM 규격 위반 리드가 있어도 stop 안 함.
  • --OPTICAL_DUPLICATE_PIXEL_DISTANCE: 플랫폼별 값.
  • --TMP_DIR: 큰 파티션 지정.

metrics.txt 읽는 법

text
LIBRARY   UNPAIRED_READS  READ_PAIRS  UNMAPPED  UNPAIRED_DUP  READ_PAIR_DUP  READ_PAIR_OPT_DUP  PERCENT_DUPLICATION  ESTIMATED_LIBRARY_SIZE
lib1      100             49,999,900  200,000   50            2,500,000      15,000             0.05                 950,000,000

두 지표를 매번 봅니다.

  • PERCENT_DUPLICATION: 0.05 = 5%. WGS면 3～10% 정상 범위. 20% 이상이면 라이브러리 재제작 결재감.
  • ESTIMATED_LIBRARY_SIZE: 유일 인서트 추정치. 라이브러리 복잡도. 클수록 좋음.

이 두 수치가 라이브러리 품질의 실무 성적표입니다.

손 계산 — 두 시나리오

같은 좌표에 리드 10개가 붙었다고 합시다.

시나리오 A: 커버리지 30배 WGS

  • 인간 게놈 3×10⁹ bp × 30× / 리드 150bp = 6억 리드 pair.
  • 같은 좌표에 리드 10개가 붙을 확률은 통계적으로 낮음 (평균 커버리지 30, 표준편차 5).
  • 실제로는 광학·PCR 중복 신호.

시나리오 B: 표적 시퀀싱 (앰플리콘)

  • PCR로 특정 유전자 부분만 증폭한 라이브러리.
  • 같은 좌표에 리드가 수천 개 붙는 게 정상 — 원본 리드가 진짜로 그 좌표에서 나온 거지 PCR 사고가 아님.
  • 이 경우 MarkDuplicates를 스킵하거나 UMI 기반 중복 판정을 씁니다.

즉 실험 설계에 따라 규칙이 완전히 달라집니다.

UMI — 라이브러리 복잡도 문제의 근본 해법

PCR 중복을 좌표만으로 판정하면 오판이 남습니다. 특히 저용량 샘플이나 표적 시퀀싱에서 그렇습니다. 근본 해법은 UMI (Unique Molecular Identifier) 입니다.

각 원본 인서트에 라이브러리 준비 단계에서 8～12bp 랜덤 바코드를 붙입니다. 리드 헤더에 UMI가 함께 들어오면 좌표 + UMI 조합으로 원본을 판정합니다. 같은 좌표에 리드가 10개 있어도 UMI가 다 다르면 진짜 서로 다른 원본입니다.

  • fgbio: Fulcrum Genomics의 UMI 처리 도구.
  • umi_tools: cgat에서 만든 파이썬 도구.
  • GATK4에도 MarkDuplicatesWithMateCigar나 UmiAwareMarkDuplicatesWithMateCigar가 있습니다.

RNA-seq 저발현 유전자 정량이나 액체생검 저VAF 검출에서는 UMI가 이제 표준입니다. GATK Best Practices 기본 파이프라인은 UMI 없는 관례라는 걸 알아두시면 됩니다.

Colab에서 소형 실습

bash
!apt-get install -y samtools default-jre >/dev/null
!wget -q https://github.com/broadinstitute/picard/releases/download/3.1.0/picard.jar
!wget -q https://github.com/samtools/samtools/raw/master/examples/toy.bam
!java -jar picard.jar MarkDuplicates \
    I=toy.bam \
    O=toy.dedup.bam \
    M=toy.metrics.txt \
    VALIDATION_STRINGENCY=LENIENT
!cat toy.metrics.txt | head -15

toy.bam은 아주 작은 예제라 통계 지표가 인공적입니다. 그래도 metrics.txt 열어보는 감을 잡을 수 있습니다.

CS 매핑

  • 해시 기반 그룹핑: 리드들을 (참조좌표, 스트랜드) 키로 해시맵에 그룹핑 → 그룹별로 대표 선택.
  • 최댓값 선택: 그룹 안에서 base quality 합계 argmax.
  • 플래그 마킹 vs 삭제: 되돌릴 수 있는 소프트 삭제 관례.

DryBench "해시 기반 중복 제거" 편(hash-based-deduplication) 참고.

마무리

MarkDuplicates는 겉보기엔 단순한 도구지만, PCR 편향 통제라는 통계적 목적을 가진 필수 관문입니다. 다음 편(S08)에서 BQSR로 넘어갑니다. Phred 점수가 왜 재보정이 필요한지, 어떻게 학습되는지 통계 수준으로 파고들어 봅시다.

더 깊게 파고 싶다면

  • Picard MarkDuplicates 공식 문서 (https://broadinstitute.github.io/picard/): 옵션 원전.
  • Ebbert M. et al. (2016), Evaluating the necessity of PCR duplicates removal from next-generation sequencing data and a comparison of approaches. BMC Bioinformatics 17:239. Duplicate 제거 필요성의 실증 연구.
  • Smith T. et al. (2017), UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Research 27:491. UMI 도구의 원 논문.
  • Broad Institute BroadE — MarkDuplicates deep dive 유튜브 강의. optical dup의 광학 원리가 시각적으로 나옵니다.
  • fgbio 공식 문서 (https://github.com/fulcrumgenomics/fgbio): UMI 파이프라인 실무 표준.

플래그 1024가 왜 삭제가 아닌지 감으로 잡히면 뒷단 파이프라인 옵션이 훨씬 명료해집니다. 다음 편에서 BQSR로 갑시다.