Correct Pysam Workflow Patterns

These patterns target pysam 0.24.0 and make filtering and coordinate semantics explicit. Adapt thresholds to the assay rather than treating them as universal defaults.

Preflight an Analysis

Before combining files, verify:

Reference assembly and contig naming agree (chr1 versus 1).
Numeric coordinates use 0-based, half-open intervals.
Alignment and variant inputs are sorted as expected.
Random-access inputs have valid indexes.
CRAM has the exact reference FASTA available.
Read-group/sample metadata identifies the intended samples.
Duplicate, secondary, supplementary, QC-fail, and quality policies are stated.

Use the bundled inspectors:

bash

python scripts/inspect_hts.py sample.bam
python scripts/inspect_hts.py cohort.vcf.gz
python scripts/inspect_hts.py reference.fa

Streaming Alignment QC

The bundled script scans in file order and does not require an index:

bash

python scripts/alignment_qc.py sample.bam --output sample.qc.json
python scripts/alignment_qc.py sample.cram \
  --reference reference.fa \
  --max-records 100000

The report counts alignment records, not unique templates or fragments. Secondary and supplementary records are reported separately. Use --max-records for a sampling pass; omit it for a complete scan.

For index-level counts without reading every record:

python

import pysam

with pysam.AlignmentFile("sample.bam", "rb", require_index=True) as bam:
    for stats in bam.get_index_statistics():
        print(stats.contig, stats.mapped, stats.unmapped, stats.total)

Index statistics are fast but cannot replace custom record-level QC.

Zero-Aware Coverage

pileup() omits positions with no columns. Use count_coverage() when zeros must be represented:

python

import pysam


def base_depth(
    bam: pysam.AlignmentFile,
    contig: str,
    start: int,
    stop: int,
    *,
    min_base_quality: int = 20,
) -> list[int]:
    a, c, g, t = bam.count_coverage(
        contig,
        start,
        stop,
        quality_threshold=min_base_quality,
        read_callback="all",
    )
    return [sum(counts) for counts in zip(a, c, g, t)]


with pysam.AlignmentFile("sample.bam", "rb") as bam:
    depth = base_depth(bam, "chr1", 1_000, 2_000)

read_callback="all" excludes unmapped, secondary, QC-fail, and duplicate records, but not supplementary records. Use a custom callback when supplementary or low-MAPQ reads must also be excluded:

python

def usable_read(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
    )


a, c, g, t = bam.count_coverage(
    "chr1",
    1_000,
    2_000,
    quality_threshold=20,
    read_callback=usable_read,
)

Only A/C/G/T query bases contribute.

Convert low-depth bases into intervals

python

def below_threshold_intervals(
    depths: list[int],
    start: int,
    threshold: int,
):
    interval_start = None

    for offset, value in enumerate(depths):
        position = start + offset
        if value < threshold and interval_start is None:
            interval_start = position
        elif value >= threshold and interval_start is not None:
            yield interval_start, position
            interval_start = None

    if interval_start is not None:
        yield interval_start, start + len(depths)

The yielded intervals remain 0-based, half-open and correctly include regions with no aligned columns.

Exact Pileup at a Position

Use explicit pileup options and retain the iterator:

python

def aligned_depth_at(
    bam: pysam.AlignmentFile,
    fasta: pysam.FastaFile,
    contig: str,
    position: int,
) -> int:
    iterator = bam.pileup(
        contig,
        position,
        position + 1,
        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:
        if column.reference_pos == position:
            return column.get_num_aligned()
    return 0

This definition is not identical to VCF INFO/DP from a caller. Name custom annotations so their provenance and filters remain clear.

SNP Base Support

This helper is deliberately limited to single-nucleotide REF/ALT alleles:

python

from collections import Counter


def snp_base_counts(
    bam: pysam.AlignmentFile,
    fasta: pysam.FastaFile,
    record: pysam.VariantRecord,
) -> Counter[str]:
    if (
        len(record.ref) != 1
        or not record.alts
        or any(len(alt) != 1 for alt in record.alts)
    ):
        raise ValueError("snp_base_counts only supports simple SNP records")

    counts: Counter[str] = Counter()
    iterator = bam.pileup(
        record.contig,
        record.start,
        record.start + 1,
        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
            if query_position is None:
                continue
            base = pileup_read.alignment.query_sequence[query_position]
            counts[base.upper()] += 1
    return counts

For indels, inspect PileupRead.indel, CIGAR, and normalized alleles or use a dedicated variant caller. This SNP method must not be generalized to symbolic or breakend alleles.

Annotate a VCF with BAM-Derived Depth

Copy the header, declare a new field, translate copied records, and use record.start directly:

python

import pysam


def annotate_depth(
    input_vcf: str,
    input_bam: str,
    reference_fasta: str,
    output_vcf: str,
) -> None:
    with pysam.FastaFile(reference_fasta) as fasta, pysam.AlignmentFile(
        input_bam,
        "rb",
        reference_filename=reference_fasta,
    ) as bam, pysam.VariantFile(input_vcf) as source:
        header = source.header.copy()
        if "BAM_BASE_DP" in header.info:
            raise ValueError("BAM_BASE_DP already exists in input header")
        header.info.add(
            "BAM_BASE_DP",
            number=1,
            type="Integer",
            description=(
                "Aligned base depth at POS; MAPQ>=20, baseQ>=20, "
                "samtools pileup filters, paired overlaps collapsed"
            ),
        )

        with pysam.VariantFile(output_vcf, "w", header=header) as output:
            for source_record in source:
                record = source_record.copy()
                record.translate(header)
                record.info["BAM_BASE_DP"] = aligned_depth_at(
                    bam,
                    fasta,
                    record.contig,
                    record.start,
                )
                output.write(record)

For CRAM input, reference_filename is essential. For a large unindexed VCF, this pattern can issue many BAM seeks; process by contig or sorted windows to improve locality.

Validate VCF REF Alleles Against FASTA

python

def reference_matches(
    record: pysam.VariantRecord,
    fasta: pysam.FastaFile,
) -> bool:
    observed = fasta.fetch(
        record.contig,
        record.start,
        record.start + len(record.ref),
    )
    return observed.upper() == record.ref.upper()


with pysam.FastaFile("reference.fa") as fasta, pysam.VariantFile(
    "variants.vcf.gz"
) as variants:
    mismatches = [
        (record.contig, record.pos, record.ref)
        for record in variants
        if not reference_matches(record, fasta)
    ]

Do not "fix" mismatches automatically. Investigate assembly version, contig aliases, left normalization, and strand/representation errors.

Filter an Alignment File

Use the bundled script for common primary-read filtering:

bash

python scripts/alignment_qc.py input.bam --max-records 10000
python scripts/filter_alignments.py input.bam filtered.bam \
  --min-mapq 20 --exclude-duplicates --exclude-supplementary --index

If implementing custom filtering, preserve the source header and stream records:

python

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

Filtering preserves input order; it does not sort. Index only coordinate-sorted output. If the header's sort-order declaration is wrong, fix the workflow rather than trusting it.

Extract Strand-Aware BED Sequences

python

IUPAC_COMPLEMENT = str.maketrans(
    "ACGTRYMKBDHVNacgtrymkbdhvn",
    "TGCAYRKMVHDBNtgcayrkmvhdbn",
)


with pysam.TabixFile(
    "genes.bed.gz", parser=pysam.asBed()
) as genes, pysam.FastaFile("reference.fa") as fasta, open(
    "genes.fa", "x", encoding="utf-8"
) as output:
    for gene in genes.fetch():
        sequence = fasta.fetch(gene.contig, gene.start, gene.end)
        if gene.strand == "-":
            sequence = sequence.translate(IUPAC_COMPLEMENT)[::-1]
        output.write(f">{gene.name}\n{sequence}\n")

BED parser coordinates are already 0-based. Sanitize or encode record names if they will be consumed by strict downstream FASTA parsers.

Count RNA Splice Junctions

find_introns() counts N CIGAR operations:

python

with pysam.AlignmentFile("rna.bam", "rb") as bam:
    primary_reads = (
        read
        for read in bam.fetch("chr1")
        if not read.is_secondary
        and not read.is_supplementary
        and not read.is_duplicate
        and read.mapping_quality >= 20
    )
    junction_counts = bam.find_introns(primary_reads)

Keys are (start, stop) 0-based splice intervals. Filter strand and library orientation according to the assay.

Bulk Operations via Wrapped Tools

For sort/index and normalization, mature command implementations are usually preferable to Python record loops:

python

import pysam.samtools
import pysam.bcftools

pysam.samtools.sort(
    "-@", "4",
    "-o", "sorted.bam",
    "input.bam",
    catch_stdout=False,
)
pysam.samtools.index(
    "-@", "4",
    "sorted.bam",
    catch_stdout=False,
)

pysam.bcftools.norm(
    "-f", "reference.fa",
    "-m", "-any",
    "-Oz",
    "-o", "normalized.vcf.gz",
    "input.vcf.gz",
    catch_stdout=False,
)
pysam.bcftools.index(
    "--csi",
    "normalized.vcf.gz",
    catch_stdout=False,
)

Pass each argument as its own string. Do not split or evaluate an untrusted shell command. Use catch_stdout=False when -o writes large or binary data.

Do Not Hand-Roll Complex VCF Merges

Combining records by (contig, pos, ref, alts) is insufficient because inputs can differ in:

allele normalization and multiallelic decomposition
contig order and metadata
INFO/FORMAT Number and Type definitions
FILTER definitions
sample names and ploidy
duplicate positions and phasing

Normalize and validate inputs, then use bcftools merge for samples or bcftools concat for disjoint genomic partitions as appropriate.

Output Validation

Alignment output:

python

pysam.samtools.quickcheck("-v", "filtered.bam")

Variant output:

python

with pysam.VariantFile("annotated.vcf.gz") as variants:
    assert "BAM_BASE_DP" in variants.header.info

Then index final sorted output and fetch known intervals at contig starts, interval boundaries, and high-coordinate regions. Format validity does not prove biological validity; compare counts and selected records to an independent tool when results matter.