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

-b

-C

-h

-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

samtools flags 147
# 0x93 147 PAIRED,PROPER_PAIR,REVERSE,READ2

CIGAR Operations

Example:

50M2I30M

pysam Python Alternative

Goal: Read and manipulate alignment data programmatically in Python.

Approach: Use

pysam.AlignmentFile

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

samtools view file.bam

AlignmentFile('file.bam', 'rb')

samtools view -H file.bam

bam.header

samtools view -c file.bam

sum(1 for _ in bam)

samtools view file.bam chr1:1-1000

bam.fetch('chr1', 0, 1000)

samtools view -h -o out.sam in.bam

samtools view -b -o out.bam in.sam

samtools view -C -T ref.fa -o out.cram in.bam

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
• bam-statistics - Generate statistics from alignment files
• sequence-io/read-sequences - Parse FASTA/FASTQ input files