omics.co.in logo omics.co.in

5: SAM/BAM/CRAM Processing & Duplication Metrics

Tested Stacks: GATK 4.6.1 Samtools 1.21 Nextflow 24.10 Seurat v5.1 Python 3.12
Prerequisites: Basic Linux CLI & Sequencing Formats
📚 Chapter 5 of 13 0% Course Completed

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 vs. Optical Duplicates

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 scores

CIGAR 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.bam

3. 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.txt

4. 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.cram
⚠️

Common 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 chr1 or 1; divergent naming causes silent alignment drops.
Was this guide helpful to your research?
Help us maintain 100% accurate, peer-reviewed bioinformatics standards.
📝 Propose Edit on GitHub ↗