Alignment Files: SAM, BAM, and CRAM

This reference targets pysam 0.24.0. All numeric coordinates shown here are 0-based, half-open.

Open Modes and Handles

import pysam

with pysam.AlignmentFile("input.bam", "rb", threads=4) as alignments:
    print(alignments.references)

For CRAM, supply the exact reference where possible:

with pysam.AlignmentFile(
    "input.cram",
    "rc",
    reference_filename="GRCh38.fa",
    threads=4,
) as alignments:
    ...

AlignmentFile can use a path or a real file object that exposes fileno(). In-memory objects such as io.BytesIO are not supported by HTSlib. Use "-" for stdin/stdout. When an existing file object is accepted, duplicate_filehandle=True (the default) prevents pysam from closing the caller's descriptor.

Useful constructor options:

• index_filename=: nonstandard, remote, or separately named index
• require_index=True: fail early if random access is required
• reference_filename=: CRAM reference FASTA
• threads=: compression/decompression threads
• format_options=["key=value"]: HTSlib format options
• ignore_truncation=True: downgrade a missing BGZF EOF marker to a warning; do not combine with threads > 1

Headers and Index State

AlignmentFile.header is an AlignmentHeader, not a plain dictionary.

with pysam.AlignmentFile("input.bam", "rb") as bam:
    header_dict = bam.header.to_dict()
    header_text = str(bam.header)
    contigs = dict(zip(bam.references, bam.lengths))
    has_index = bam.has_index()

Use check_index() when a missing index should be an error. It raises for SAM, closed files, or unusable indexes. get_index_statistics() exposes per-contig mapped/unmapped counts recorded in an available index; these are index statistics, not a fresh scan of every record.

Iteration Choices

Indexed region query

with pysam.AlignmentFile("input.bam", "rb") as bam:
    for read in bam.fetch("chr1", 1_000, 2_000):
        ...

• Requires BAI/CSI for BAM or CRAI for CRAM.
• Returns reads overlapping the interval, including reads that start before it or end after it.
• Records are returned in coordinate/index order.
• fetch() with no region still requires an index and returns mapped records.

Sequential scan

with pysam.AlignmentFile("input.bam", "rb") as bam:
    for read in bam.fetch(until_eof=True):
        ...

• Does not require an index.
• Starts at the current file position.
• Preserves file order.
• Includes unplaced unmapped records.

fetch("*") requests only unplaced unmapped records at the end of a coordinate-sorted alignment file.

Multiple active iterators

multiple_iterators belongs to fetch(), not the AlignmentFile constructor:

with pysam.AlignmentFile("input.bam", "rb") as bam:
    chr1_reads = bam.fetch("chr1", multiple_iterators=True)
    chr2_reads = bam.fetch("chr2", multiple_iterators=True)

Each such iterator reopens the file and has overhead. Prefer one ordered pass when possible.

AlignedSegment Essentials

Identity and sequence

• query_name
• query_sequence
• query_qualities: numeric Phred scores, or None
• query_length
• query_alignment_sequence: query bases participating in the alignment
• query_alignment_qualities
• get_forward_sequence() / get_forward_qualities(): original sequencer orientation

Assigning query_sequence invalidates query_qualities. Save and reassign qualities after changing sequence.

Reference placement

• reference_name
• reference_id
• reference_start: 0-based inclusive
• reference_end: 0-based exclusive; derived from CIGAR
• mapping_quality
• next_reference_name, next_reference_start
• template_length

Many placement-derived properties are None or sentinel values for unmapped or CIGAR-less records. Check is_unmapped before using them.

Flags

Common boolean properties:

• pairing: is_paired, is_proper_pair, is_read1, is_read2
• orientation: is_reverse, mate_is_reverse
• mapping: is_unmapped, mate_is_unmapped
• status: is_secondary, is_supplementary, is_qcfail, is_duplicate

For a typical primary mapped-read filter:

def keep_primary(read: pysam.AlignedSegment) -> bool:
    return (
        not read.is_unmapped
        and not read.is_secondary
        and not read.is_supplementary
        and not read.is_qcfail
        and not read.is_duplicate
        and read.mapping_quality >= 20
    )

State whether supplementary alignments and duplicates are intentionally excluded; there is no universal filter for every analysis.

CIGAR Operations and Alignment Geometry

cigartuples stores (operation, length), the reverse of textual SAM CIGAR notation.

Prefer enum members in new code. Top-level aliases such as pysam.CMATCH remain in 0.24 for compatibility but are expected to be removed in a future release.

Useful geometry methods:

blocks = read.get_blocks()  # aligned reference blocks; gaps at D/N
pairs = read.get_aligned_pairs(matches_only=False, with_cigar=True)
reference_positions = read.get_reference_positions(full_length=True)

get_aligned_pairs(with_seq=True) needs an MD tag and returns reference bases derived from that tag. It does not consult a separately opened FASTA.

Optional Tags and Modified Bases

if read.has_tag("NM"):
    edit_distance = read.get_tag("NM")

read.set_tag("XX", 7, value_type="i")
all_tags = read.get_tags(with_value_type=True)

Use standard tags according to the SAM tags specification. Avoid changing alignment-derived tags such as NM/MD without recomputing them.

For base modifications encoded by MM/ML:

for (canonical_base, strand, modification), calls in (
    read.modified_bases or {}
).items():
    for query_position, quality in calls:
        probability = None if quality < 0 else quality / 256.0

The key is (canonical base, strand, modification), where strand is 0 forward or 1 reverse. pysam 0.24 removed the earlier five-modification-type limit and fixed crashes on degenerate empty MM calls.

Counting and Coverage

Record count

with pysam.AlignmentFile("input.bam", "rb") as bam:
    raw_overlap_count = bam.count("chr1", 1_000, 2_000)
    filtered_count = bam.count(
        "chr1",
        1_000,
        2_000,
        read_callback="all",
    )

count() defaults to read_callback="nofilter". "all" excludes reads with unmapped, secondary, QC-fail, or duplicate flags. It does not accept a quality= argument. Use a callback for custom record filtering:

count = bam.count(
    "chr1",
    1_000,
    2_000,
    read_callback=lambda read: keep_primary(read),
)

A/C/G/T base coverage

a, c, g, t = bam.count_coverage(
    "chr1",
    1_000,
    2_000,
    quality_threshold=20,
    read_callback="all",
)
depth = [sum(values) for values in zip(a, c, g, t)]

The result has exactly stop - start positions and therefore represents zero coverage. Only A/C/G/T bases are counted; ambiguous query bases do not contribute.

Pileup Semantics

with pysam.FastaFile("reference.fa") as fasta, pysam.AlignmentFile(
    "input.bam", "rb"
) as bam:
    iterator = bam.pileup(
        "chr1",
        1_000,
        2_000,
        truncate=True,
        stepper="samtools",
        fastafile=fasta,
        min_mapping_quality=20,
        min_base_quality=20,
        max_depth=100_000,
        ignore_overlaps=True,
        ignore_orphans=True,
    )
    for column in iterator:
        for pileup_read in column.pileups:
            if pileup_read.is_del or pileup_read.is_refskip:
                continue
            query_position = pileup_read.query_position
            base = pileup_read.alignment.query_sequence[query_position]

Key defaults and behaviors:

• Without truncate=True, columns outside the requested interval can appear when overlapping reads extend beyond the interval.
• stepper="all" filters unmapped, secondary, QC-fail, and duplicate reads.
• stepper="nofilter" disables read filtering.
• stepper="samtools" applies samtools-style processing; provide fastafile for full BAQ/reference behavior.
• Default min_base_quality is 13.
• Default max_depth is 8000.
• Paired overlap detection and orphan filtering are enabled by default.
• nsegments counts reads in the pileup column before base-level exclusions; get_num_aligned() is often the clearer aligned-base depth.

PileupColumn and PileupRead proxy objects are valid only while their iterator remains alive. Do not retain a column after iteration ends.

For SNP support, inspect base and quality. For insertions/deletions, use PileupRead.indel, deletion/refskip state, CIGAR, and normalized alleles. Single-base counting is not an indel caller.

Writing Alignments

Preserve an input header

with pysam.AlignmentFile("input.bam", "rb") as source, pysam.AlignmentFile(
    "filtered.bam", "wb", template=source, threads=4
) as destination:
    for read in source.fetch(until_eof=True):
        if keep_primary(read):
            destination.write(read)

The output retains input order. Only index it if that order is coordinate sorted and unmapped records remain in a valid location.

Construct a new record

header = pysam.AlignmentHeader.from_dict(
    {
        "HD": {"VN": "1.6", "SO": "coordinate"},
        "SQ": [{"SN": "chr1", "LN": 248_956_422}],
    }
)

with pysam.AlignmentFile("new.bam", "wb", header=header) as output:
    read = pysam.AlignedSegment(output.header)
    read.query_name = "read001"
    read.query_sequence = "ACGTACGTAA"
    read.flag = 0
    read.reference_id = output.get_tid("chr1")
    read.reference_start = 100
    read.mapping_quality = 60
    read.cigartuples = [(pysam.CIGAR_OPS.CMATCH, 10)]
    read.query_qualities = pysam.qualitystring_to_array("IIIIIIIIII")
    output.write(read)

Construct records against the destination header so numeric reference IDs map correctly. Sequence length, qualities, and query-consuming CIGAR operations must agree.

Validation

For BAM/CRAM:

pysam.samtools.quickcheck("-v", "filtered.bam")
pysam.samtools.index("filtered.bam", catch_stdout=False)

quickcheck checks basic headers and EOF markers, not biological correctness. Reopen the file, verify header/reference compatibility, and inspect expected regions. Read cram_and_performance.md for CRAM-specific validation.