Post-Alignment Processing

intermediate

Post-Alignment Processing Overview

GATK Best Practices BAM Processing

📋Sorted BAM

🔖MarkDuplicates

📐BaseRecalibrator

✅ApplyBQSR

🎯Recalibrated BAM

MarkDuplicates: Flags PCR duplicates (FLAG 1024)

BaseRecalibrator: Build recalibration model vs dbSNP/Mills

ApplyBQSR: Corrected quality scores

Recalibrated BAM: Ready for HaplotypeCaller

Mark Duplicates

PCR duplicates arise from amplification of the same DNA fragment. They inflate apparent coverage and create false variant evidence. Must be marked before variant calling.

code

# GATK MarkDuplicates (recommended - faster than Picard)
gatk MarkDuplicates \
  -I sample.sorted.bam \
  -O dedup.bam \
  -M dedup_metrics.txt \
  --REMOVE_DUPLICATES false  # mark only, do not remove

samtools index dedup.bam

# Check duplication rate
grep "ESTIMATED_LIBRARY_SIZE" -A1 dedup_metrics.txt

Optical duplicates: from same cluster on flow cell; flagged separately
PCR duplicates: from amplification; same start/end coordinates on both reads
UMIs (Unique Molecular Identifiers): barcodes added pre-PCR; allow true duplicate collapse
Duplication rate: <15% ideal for WES; >30% indicates library quality issues
REMOVE_DUPLICATES false: mark with FLAG 1024 but keep; allows GATK to ignore them

Base Quality Score Recalibration (BQSR)

BQSR detects systematic patterns in base quality scores using machine learning. Sequencers mis-estimate quality for specific positions, cycles, and base contexts. BQSR corrects this, improving variant calling accuracy by ~5–10%.

code

# Step 1: Build recalibration model
gatk BaseRecalibrator \
  -I dedup.bam \
  -R hg38.fa \
  --known-sites dbsnp_138.hg38.vcf.gz \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
  -O recal_data.table

# Step 2: Apply recalibration
gatk ApplyBQSR \
  -I dedup.bam \
  -R hg38.fa \
  --bqsr-recal-file recal_data.table \
  -O recal.bam

# Optional: Evaluate improvement
gatk AnalyzeCovariates \
  -before recal_data.table \
  -after post_recal.table \
  -plots recal_plots.pdf

Interval Padding for WES

code

# Most GATK tools accept interval lists for WES
# Add 100bp padding to capture exon flanks (splice sites)
gatk PreprocessIntervals \
  -R hg38.fa \
  -L capture.bed \
  --padding 100 \
  --bin-length 0 \
  -O intervals.interval_list

# Use in BaseRecalibrator
gatk BaseRecalibrator \
  -I dedup.bam \
  -R hg38.fa \
  -L intervals.interval_list \
  --known-sites dbsnp.vcf.gz \
  -O recal_data.table

Sequence Alignment Variant Calling