Version Compatibility

Reference examples tested with: pysam 0.22+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

SAM/BAM/CRAM Basics

"Read a BAM file" → Open a binary alignment file and iterate over aligned reads with their mapping coordinates, flags, and quality scores.

• Python: pysam.AlignmentFile() (pysam)
• CLI: samtools view (samtools)
• R: scanBam() (Rsamtools)

View and convert alignment files using samtools and pysam.

Format Overview

SAM Format Structure

@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:248956422
@RG ID:sample1 SM:sample1
@PG ID:bwa PN:bwa VN:0.7.17
read1  0   chr1  100  60  50M  *  0  0  ACGT...  FFFF...  NM:i:0

Header lines start with @:

• @HD - Header metadata (version, sort order)
• @SQ - Reference sequence dictionary
• @RG - Read group information
• @PG - Program used to create file

Alignment fields (tab-separated):

1. QNAME - Read name
2. FLAG - Bitwise flag
3. RNAME - Reference name
4. POS - 1-based position
5. MAPQ - Mapping quality
6. CIGAR - Alignment description
7. RNEXT - Mate reference
8. PNEXT - Mate position
9. TLEN - Template length
10. SEQ - Read sequence
11. QUAL - Base qualities
12. Optional tags (NM:i:0, MD:Z:50, etc.)

samtools view

View BAM as SAM

samtools view input.bam | head

View with Header

samtools view -h input.bam | head -100

View Header Only

samtools view -H input.bam

View Specific Region

samtools view input.bam chr1:1000-2000

Count Alignments

samtools view -c input.bam

Format Conversion

Goal: Convert between SAM (text), BAM (binary), and CRAM (reference-compressed) alignment formats.

Approach: Use samtools view with format flags (-b for BAM, -C for CRAM, -h for SAM with header). CRAM requires a reference FASTA with -T.

BAM to SAM

samtools view -h -o output.sam input.bam

SAM to BAM

samtools view -b -o output.bam input.sam

BAM to CRAM

samtools view -C -T reference.fa -o output.cram input.bam

CRAM to BAM

samtools view -b -T reference.fa -o output.bam input.cram

Pipe Conversion

samtools view -b input.sam > output.bam

Common Flags

Decode Flags (Bidirectional)

# Number to mnemonics
samtools flags 147
# 0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2

# Mnemonics to number
samtools flags PAIRED,PROPER_PAIR,REVERSE,READ2   # 147

Secondary vs Supplementary (Different Semantics)

Two different concepts that are routinely conflated:

-F 2304 removes both. Strip supplementary only when downstream is small-variant calling; keep supplementary for SV calling, fusion detection, or any analysis that follows split-reads.

MAPQ Is Not Portable Across Aligners

samtools view -q 30 does different things depending on what produced the BAM. MAPQ is an aligner-specific scale, not a universal probability:

Verify the actual scale of any unfamiliar BAM:

samtools view input.bam | awk '{print $5}' | sort -un | head
samtools view -H input.bam | grep '^@PG' | head -1   # which aligner produced this BAM

0-Based vs 1-Based Coordinates (Footgun)

samtools view bam chr1:100-200 and bam.fetch('chr1', 100, 200) return different read sets at boundaries.

CIGAR Operations

Example: 50M2I30M = 50 bases match, 2 base insertion, 30 bases match

CIGAR M is overloaded -- it is the union of = and X. Some aligners (minimap2, BWA with -Y) emit =/X directly; bcftools / Picard often need M and rebuild MD/NM with samtools calmd. N operations break naive coverage calculations: a 1000 bp RNA-seq read with one 50 kb intron does not cover 50 kb. Distinguish soft-clip (S, bases retained) from hard-clip (H, bases discarded -- irreversible).

Context-Specific Tags

Beyond the standard fields, downstream tools depend on optional tags whose presence depends on aligner and assay. Inspect with samtools view input.bam | head -1 | tr '\t' '\n' or pysam read.get_tag('XX').

Missing tags fail in two modes: silently wrong (featureCounts ignoring multimappers without NH; markdup marking nothing without MC/MS) or loudly (consensus tools rejecting input without MD).

Provenance: @PG Chain

The @PG lines record every tool that touched the BAM, linked through PP (previous program) tags. This is the audit trail.

samtools view -H input.bam | grep '^@PG'

A clean germline pipeline:

@PG ID:bwa-mem PN:bwa VN:0.7.17
@PG ID:samtools.1 PN:samtools VN:1.20 PP:bwa-mem CL:samtools sort
@PG ID:samtools.2 PN:samtools VN:1.20 PP:samtools.1 CL:samtools fixmate
@PG ID:samtools.3 PN:samtools VN:1.20 PP:samtools.2 CL:samtools markdup

A broken/missing chain (no PP, unknown tools, gaps) means the BAM cannot be reliably reproduced. Production pipelines often reject inputs without a complete chain.

CRAM Reference Resolution (Critical)

CRAM stores reads relative to a reference; without it, the file is unreadable. htslib resolves the reference in this order:

1. Command-line -T ref.fa / --reference
2. REF_CACHE env var (local MD5-named cache)
3. REF_PATH env var (colon-separated; can include URLs)
4. UR: URL in SAM @SQ header
5. Last resort: EBI ENA download via MD5 in M5: tag (fails on offline HPC)

On HPC nodes without internet, populate a local cache once:

mkdir -p $HOME/cram_cache
seq_cache_populate.pl -root $HOME/cram_cache reference.fa
export REF_CACHE=$HOME/cram_cache/%2s/%2s/%s
export REF_PATH=$REF_CACHE   # disables ENA fallback

samtools quickcheck -v file.cram   # header + EOF only
samtools view -c file.cram          # forces full decode; proves reference reachable

CRAM operations can be irreversibly lossy: --output-fmt-option=archive=1 enables 8-bin Illumina quality binning (~30-50% additional size reduction; benign for >=30x germline WGS, harmful for low-coverage / somatic / forensic / archival). Convert against the exact reference the BAM was aligned to (matched by @SQ M5:); a different reference silently corrupts bases on read-back.

pysam Python Alternative

Goal: Read and manipulate alignment data programmatically in Python.

Approach: Use pysam.AlignmentFile to open BAM/CRAM files, iterate over reads, and access properties like coordinates, flags, CIGAR, and tags.

Open and Iterate

import pysam

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for read in bam:
        print(f'{read.query_name}\t{read.reference_name}:{read.reference_start}')

Access Header

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for sq in bam.header['SQ']:
        print(f'{sq["SN"]}: {sq["LN"]} bp')

Read Alignment Properties

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for read in bam:
        print(f'Name: {read.query_name}')
        print(f'Flag: {read.flag}')
        print(f'Chrom: {read.reference_name}')
        print(f'Pos: {read.reference_start}')  # 0-based
        print(f'MAPQ: {read.mapping_quality}')
        print(f'CIGAR: {read.cigarstring}')
        print(f'Seq: {read.query_sequence}')
        print(f'Qual: {read.query_qualities}')
        break

Check Flag Properties

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for read in bam:
        if read.is_paired and read.is_proper_pair:
            if read.is_reverse:
                strand = '-'
            else:
                strand = '+'
            print(f'{read.query_name} on {strand} strand')

Fetch Region

with pysam.AlignmentFile('input.bam', 'rb') as bam:
    for read in bam.fetch('chr1', 1000, 2000):
        print(read.query_name)

Convert BAM to SAM

with pysam.AlignmentFile('input.bam', 'rb') as infile:
    with pysam.AlignmentFile('output.sam', 'w', header=infile.header) as outfile:
        for read in infile:
            outfile.write(read)

Convert to CRAM

with pysam.AlignmentFile('input.bam', 'rb') as infile:
    with pysam.AlignmentFile('output.cram', 'wc', reference_filename='reference.fa', header=infile.header) as outfile:
        for read in infile:
            outfile.write(read)

Quick Reference

Related Skills

• alignment-indexing - Create indices for random access (required for fetch/region queries)
• alignment-sorting - Sort alignments by coordinate or name
• alignment-filtering - Filter alignments by flags, quality, regions
• alignment-validation - Sequence dictionary cross-validation (M5 checksums)
• bam-statistics - Generate statistics from alignment files
• reference-operations - REF_PATH/REF_CACHE setup for CRAM
• sequence-io/read-sequences - Parse FASTA/FASTQ input files