Correct Pysam Workflow Patterns

These patterns target pysam 0.24.1 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:

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

Use the bundled inspectors:

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:

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:

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. CRAI has no record counts: the CRAM index-statistic zeros are unavailable data. The QC script counts mapped non-soft-clipped query bases, including insertions, only when CIGAR and SEQ exist. Missing SEQ does not contribute to mean length.

Zero-Aware Coverage

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

import pysam


def base_depth(
    bam: pysam.AlignmentFile,
    contig: str,
    start: int,
    stop: int,
    *,
    min_base_quality: int = 20,
) -> list[int]:
    if not 0 <= start < stop <= bam.get_reference_length(contig):
        raise ValueError("coverage interval must be nonempty and within the contig")
    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:

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. Overlapping mates count separately. count_coverage does not use pileup BAQ/overlap correction; pass numeric span arguments because 0.24.1 does not size arrays from region=/end= aliases.

Convert low-depth bases into intervals

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:

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 sum(
                not item.is_del and not item.is_refskip
                and item.query_position is not None
                for item in column.pileups
            )
    return 0

get_num_aligned() includes D/N entries in 0.24.1, so this helper deliberately counts observed bases. BAQ is enabled when the reference is supplied. MAPQ255 is retained by the numeric threshold; apply an explicit unavailable-MAPQ policy if required. 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:

from collections import Counter


def snp_base_counts(
    bam: pysam.AlignmentFile,
    fasta: pysam.FastaFile,
    record: pysam.VariantRecord,
) -> Counter[str]:
    if (
        record.ref.upper() not in {"A", "C", "G", "T"}
        or not record.alts
        or any(alt.upper() not in {"A", "C", "G", "T"} for alt in record.alts)
        or any(alt.upper() == record.ref.upper() for alt in record.alts)
    ):
        raise ValueError("snp_base_counts only supports simple SNP records")

    # Validate REF without inferring a variant from a mismatched assembly.
    observed_ref = fasta.fetch(record.contig, record.start, record.start + 1)
    if observed_ref.upper() != record.ref.upper():
        raise ValueError("VCF REF does not match the selected FASTA")
    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:

import pysam


def annotate_depth(
    input_vcf: str,
    input_bam: str,
    reference_fasta: str,
    output_vcf: str,
) -> None:
    from pathlib import Path

    if Path(output_vcf).exists():
        raise FileExistsError(output_vcf)
    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

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:

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

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

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. Use --index --csi for BAM beyond BAI limits. Per-record filtering preserves RG/reference dictionaries and mate fields, but can remove one mate without fetching or repairing its partner.

Extract Strand-Aware BED Sequences

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 = pysam.reverse_complement(sequence)
        output.write(f">{gene.name}\n{sequence}\n")

This requires BED6 with a strand and extracts the full span, not spliced BED12 exons. Validate FASTA bounds and strand before extraction. 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:

with pysam.AlignmentFile("rna.bam", "rb") as bam:
    primary_reads = (
        read
        for read in bam.fetch("chr1")
        if not read.is_unmapped
        and not read.is_qcfail
        and 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:

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",
    "--multi-overlaps", ".",  # other ALT alleles become missing after splitting
    "-Oz",
    "-o", "normalized.vcf.gz",
    "input.vcf.gz",
    catch_stdout=False,
)
pysam.bcftools.index(
    "--csi",
    "normalized.vcf.gz",
    catch_stdout=False,
)

Splitting changes genotype/allele indexing: the explicit missing policy avoids turning another ALT into a reference call. It does not preserve the original multiallelic haplotype representation. Keep the source and validate GT/AD/PL semantics before interpreting normalized calls.

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:

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

Variant output:

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.