5: SAM/BAM/CRAM Processing & Duplication Metrics
Chapter 5: SAM/BAM/CRAM Processing & Duplication Metrics
Sequence Alignment/Map (SAM) format and its binary counterpart (BAM) represent the universal container format for aligned genomic reads. Mastering the SAM specification, decoding bitwise flags, and properly handling duplicate reads are indispensable skills for computational genomics.
PCR Duplicates arise when a single genomic fragment is amplified disproportionately during library prep PCR cycles, producing identical molecules that map to the exact same 5' coordinates. Optical Duplicates occur when a single sequencing cluster on a patterned flow cell is falsely identified as two adjacent clusters by the optical camera. If duplicates are not marked, variant callers treat them as independent biological observations, generating false-positive heterozygous/homozygous calls.
1. The SAM Format & CIGAR Strings
Every alignment record contains 11 mandatory tab-delimited columns:
# Col 1: QNAME - Query template name (read ID)
# Col 2: FLAG - Bitwise condition flag (e.g. 99, 147)
# Col 3: RNAME - Reference sequence name (chr1, chr2)
# Col 4: POS - 1-based leftmost mapping position
# Col 5: MAPQ - Mapping quality (-10 * log10 P(wrong mapping))
# Col 6: CIGAR - Alignment operations (e.g. 150M, 75M2I73M, 30S120M)
# Col 7: RNEXT - Reference name of mate (= if same chromosome)
# Col 8: PNEXT - 1-based leftmost mapping position of mate
# Col 9: TLEN - Observed template length (insert size)
# Col 10: SEQ - Segment sequence
# Col 11: QUAL - ASCII Phred quality scoresCIGAR Operations: M (alignment match or mismatch), I (insertion into reference), D (deletion from reference), N (skipped region/intron in RNA-seq), S (soft clipping: bases present in read but not aligned to reference), and H (hard clipping: bases stripped from record).
2. Bitwise FLAG Arithmetic
# Inspect reads matching specific bitwise conditions:
# Flag 1 (0x1): Template having multiple segments in sequencing (paired)
# Flag 2 (0x2): Each segment properly aligned according to aligner
# Flag 4 (0x4): Segment unmapped
# Flag 16 (0x10): SEQ being reverse complemented
# Flag 1024 (0x400): Read is PCR or optical duplicate
# Filter for properly paired reads only (exclude unmapped, retain flag 2)
samtools view -f 2 -F 3852 -b sample1.bam > sample1.properly_paired.bam3. Duplicate Marking with Picard / Sambamba
# Mark duplicates using Sambamba (high-speed multithreaded alternative to Picard)
sambamba markdup \
-t 8 \
--tmpdir=/tmp/sambamba \
sample1_aligned.sorted.bam \
sample1_aligned.marked_dup.bam
# Generate comprehensive alignment metrics with samtools
samtools flagstat sample1_aligned.marked_dup.bam
samtools stats sample1_aligned.marked_dup.bam > sample1.stats.txt4. CRAM Reference-Based Archival Compression
# Convert heavy BAM (100 GB) to space-saving CRAM (~35 GB)
samtools view -C -T GRCh38_noalt_decoy.fasta -@ 8 -o sample1.cram sample1_aligned.marked_dup.bam
samtools index sample1.cramCommon Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.