Chapter 5: Step 3 — Post-Alignment Processing: Deduplication & GATK BQSR

Chapter 5: Step 3 — Post-Alignment Processing: Deduplication & GATK BQSR

Raw BAM files generated directly by aligners cannot be immediately subjected to variant calling. Library amplification introduces duplicate copies of identical DNA molecules, while sequencers exhibit systematic cycle-dependent quality estimation errors. In this chapter, we master coordinate sorting, PCR and optical duplicate marking, and GATK Base Quality Score Recalibration (BQSR).

1. Duplicate Marking: Optical vs PCR Duplicates

During library preparation, PCR amplification cycles can produce thousands of identical copies of the same original genomic fragment. If a mutation arose during PCR Cycle 1, it could be present across 50 duplicate reads, creating a devastating false-positive call.

• PCR Duplicates: Reads that have the exact same 5' alignment coordinates and orientation for both forward and reverse reads. Picard retains the single read pair with the highest sum of base quality scores and sets the 0x400 bitwise SAM flag on all duplicate copies.
• Optical Duplicates: Clusters on the flow cell that physically originate from a single cluster that spread into adjacent nanowells. Picard uses pixel distance thresholds to distinguish optical from library duplicates.

2. Base Quality Score Recalibration (BQSR)

Illumina instruments assign Phred quality scores based on physical optical properties. However, real-world accuracy drifts systematically based on:

• Cycle Bias: Phred scores typically decline toward the 3' end of the read due to phasing/pre-phasing decay.
• Dinucleotide Context Bias: Errors occur disproportionately after specific motifs (e.g. after GG or CC homopolymers).

GATK BQSR analyzes all aligned bases that do not overlap known polymorphic sites (masking known SNPs in dbSNP and 1000 Genomes), builds a multidimensional empirical error covariate table, and recalculates accurate Phred values.

# Step 3a: Mark Duplicates with GATK MarkDuplicatesSpark (multithreaded)
gatk MarkDuplicatesSpark     -I patient_sorted.bam     -O patient_dedup.bam     -M dedup_metrics.txt     --conf 'spark.executor.cores=8'

# Step 3b: Generate BQSR Recalibration Covariate Table
gatk BaseRecalibrator     -I patient_dedup.bam     -R GRCh38.primary_assembly.fa     --known-sites dbsnp_146.hg38.vcf.gz     --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz     -O recal_data.table

# Step 3c: Apply Recalibration to Generate Final Analysis-Ready BAM
gatk ApplyBQSR     -I patient_dedup.bam     -R GRCh38.primary_assembly.fa     --bqsr-recal-file recal_data.table     -O patient_recalibrated.bam

# Index the final recalibrated BAM
samtools index patient_recalibrated.bam