pysam-genomics

pysam — BAM/CRAM Python Interface

Hard rules

1. No fabricated citations. Every cited work must resolve to a verifiable
2. No claim without provenance. Every quantitative or factual claim
3. No silent failure. Every script invocation, API call, or tool use must declare its exit status and what to do on non-zero. A skill that silently swallows errors is a violation.

When to use

• Programmatically iterate BAM alignments with Python logic.
• Extract variant or coverage data from a BAM.
• Build a custom pipeline that needs alignment-level access.
• Query regions from indexed BAM/CRAM without loading a whole file.

When NOT to use

• You need a BAM stats summary → use samtools stats or mosdepth.
• You need variant calling → use GATK/HaplotypeCaller or freebayes.
• Visual → IGV, IGV.js, or JBrowse.

Prerequisites

• pysam ≥ 0.22 (PyPI or conda)
• htslib ≥ 1.19 (linked automatically)

Core workflow

1. Open an alignment file with pysam.AlignmentFile.
2. Iterate and filter with Python (no shell calls).
3. Query by region with pileup / fetch.
4. Write optional BAM/CRAM output.
5. Close explicitly or use a context manager.

Code patterns

Open a BAM (read)

import pysam
"
bam = pysam.AlignmentFile("input.bam", "rb")

Open a SAM (auto-detect gzipped)

sam = pysam.AlignmentFile("input.sam", "r")

Open a indexed region

Access a small region without loading the whole BAM:

bam = pysam.AlignmentFile("input.bam", "rb")
for read in bam.fetch("chr1", 100_000, 200_000):
    print(read.query_name, read.reference_name, read.pos, read.cigartuples)

Iterate through all alignments

bam = pysam.AlignmentFile("input.bam", "rb")
for read in bam:
    # read is a pysam AlignedSegment
    if read.is_proper_pair:
        print(read.query_name, read.tlen)

Filter by FLAG in Python

def proper_paired(read):
    return read.is_paired and read.is_proper_pair

bam = pysam.AlignmentFile("input.bam", "rb")
proper = [r for r in bam if proper_paired(r)]

Extract MAPQ, CIGAR, cigarstring

for read in bam:
    print(read.query_name,
          read.mapq,
          read.cigartuples,   # [(0, 100)] = 100M
          read.cigarstring)   # "100M"

Extract the aligned sequence

for read in bam:
    if not read.is_unmapped:
        print(read.query_sequence[:50])  # First 50 bp

Access mate information

for read in bam:
    if read.is_paired and read.mate_is_unmapped:
        print(f"{read.query_name} mate is unmapped")

Pileup at a position

bam = pysam.AlignmentFile("input.bam", "rb")
for column in bam.pileup("chr1", 100_000, 100_001):
    for read in column.pileups:
        print(read.alignment.query_name)

Write a filtered BAM

with pysam.AlignmentFile("filtered.bam", "wb", header=bam.header) as out:
    for read in bam:
        if read.mapq >= 30 and read.is_proper_pair:
            out.write(read)

Iterating over a tabix-indexed BED/VCF

vcf = pysam.TabixFile("variants.vcf.gz")
for record in vcf.contigs:
    pass  # Use region query for actual records

Read from CRAM (requires reference)

cram = pysam.AlignmentFile("input.cram", "rc", reference_fasta="genome.fa")

Access the header as Python dict

header = bam.header.to_dict()
print(header["RG"])  # List of read groups

Convert SAM to Python objects

with pysam.AlignmentFile("input.sam", "r") as sam:
    for read in sam:
        # read is an AlignedSegment
        ...

Get index statistics without loading

index = pysam.IndexedReads(pysam.AlignmentFile("input.bam"))

The Index object is for random-access.

Create custom tags

bam = pysam.AlignmentFile("input.bam", "rb", fields=["NM", "MD", "RG"])
# Fields are automatically loaded if they exist

Check if a position is covered

bam = pysam.AlignmentFile("input.bam", "rb")
for col in bam.pileup("chr1", 100_000, 100_001):
    coverage = col.nsegments
    if coverage >= 10:
        print(f"Position covered: {coverage}")

Extract all reads overlapping a gene

# Use fetch for gene coordinates
bam = pysam.AlignmentFile("input.bam", "rb")
gene_reads = [r for r in bam.fetch("chr1", start, end)]

Read group filter

rg_bam = pysam.AlignmentFile("input.bam", "rb", header=header)
for read in rg_bam:
    # read.get_tag("RG") - read group from optional tags
    pass

Common pitfalls

• Forgetting to close the file. Use a context manager (with ...) or close explicitly.
• CIGAR interpretation. read.cigartuples returns (operation, length) tuples. 0 = M, 1 = I, 2 = D.
• MAPQ zero means unmapped in this alignment, not zero quality.
• Reference fasta for CRAM. Must be the exact FASTA used to create the CRAM.
• Duplicate reads being filtered. samtools markdup -r is the modern way.
• Pileup object is lazy. If you need it, materialize it (e.g., with list()) otherwise it's consumed by iteration.

Validation

• After filtering, samtools flagstat filtered.bam should show the expected read count.
• read_mapq should be >= your threshold.
• read.is_proper_pair checks both ends.

Open alternatives

References

• pysam documentation: https://pysam.readthedocs.io/
• htslib: http://www.htslib.org/
• SAM spec: https://samtools.github.io/hts-specs/SAMv1.pdf
• Companion: ors-bioinformatics-sequence-sam-bam-basics, ors-bioinformatics-sequence-bcftools-variant-manipulation.

Changelog

• 1.0.0 (2026-06-10): Initial adaptation by Pradyumna Jayaram from pysam-genomic-files (SciAgent-Skills).

Cross-references

Other skills in this category:

• batch-processing
• bowtie2-alignment
• bwa-alignment
• bwa-mem2-alignment
• codon-usage
• compressed-sequence-files
• fastq-quality-scores
• filter-sequences
• format-conversion
• hisat2-alignment
• motif-search
• paired-end-fastq
• read-write-sequences
• reverse-complement
• sam-bam-basics
• samtools-bam-processing
• seq-objects
• sequence-properties
• sequence-slicing
• sequence-statistics
• star-alignment
• transcription-translation