/pysam-genomic-files

name: pysam-genomic-files
description: "Read/write SAM/BAM/CRAM, VCF/BCF, FASTA/FASTQ. Region queries, pileup, variant filtering, read groups. Python htslib wrapper exposing samtools/bcftools CLI. Use STAR/BWA for alignment; GATK/DeepVariant for variant calling."
license: MIT

Pysam — Genomic File Toolkit

Overview

Pysam provides a Pythonic interface to htslib for reading, manipulating, and writing genomic data files. It handles SAM/BAM/CRAM alignments, VCF/BCF variants, and FASTA/FASTQ sequences with efficient region-based random access. Also exposes samtools and bcftools as callable Python functions.

When to Use

• Reading and querying BAM/CRAM alignment files (region extraction, read filtering)
• Analyzing VCF/BCF variant files (genotype access, variant filtering, annotation)
• Extracting reference sequences from indexed FASTA files
• Calculating per-base coverage and pileup statistics
• Building custom bioinformatics pipelines that combine alignment + variant + sequence data
• Quality control of NGS data (mapping quality, flag filtering, coverage)
• For **alignment from FASTQ** (read mapping), use STAR, BWA, or minimap2 instead
• For **variant calling from BAM**, use GATK or DeepVariant instead

Prerequisites

pip install pysam

**Note**: Requires htslib C library (bundled with pip install on most platforms). On some Linux systems, may need `libhts-dev` or equivalent. Index files (`.bai`, `.tbi`, `.fai`) required for random access — create with `pysam.index()`, `pysam.tabix_index()`, or `pysam.faidx()`.

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: mappingQualityFloor
    kind: required
    source: user
    ask: "Below what mapping confidence should reads be ignored when piling up or counting?"
    default: "0 - no filter"

  - id: D2
    param: baseQualityFloor
    kind: required
    source: user
    ask: "How confident must an individual base call be to be counted at a position?"
    default: 15

  - id: D3
    param: regionTruncation
    kind: required
    source: user
    ask: "Should a pileup report only positions inside the requested region, or every position touched by an overlapping read?"
    default: "not truncated - overlapping reads extend the reported range"

  - id: D4
    param: indexRequirement
    kind: derived
    source: data
    ask: "Is the file indexed, or must records be streamed sequentially?"
    default: "indexed random access"

  - id: D5
    param: iteratorSafety
    kind: never_ask
    source: user
    reason: "Whether simultaneous iterators are allowed is an API-usage concern with a small overhead cost, not a result-changing choice"
    default: false

D3 is the one that surprises people: an untruncated pileup returns positions outside the region that was asked for, so per-region statistics computed straight off the result cover a wider span than intended.

Quick Start

import pysam

# Read BAM file, fetch reads in a region
with pysam.AlignmentFile("sample.bam", "rb") as bam:
    for read in bam.fetch("chr1", 1000, 2000):
        print(f"{read.query_name}: pos={read.reference_start}, mapq={read.mapping_quality}")
    print(f"Total reads in region: {bam.count('chr1', 1000, 2000)}")

Core API

1. Alignment Files (SAM/BAM/CRAM)

Read, query, and write aligned sequencing reads.

import pysam

# Open BAM (binary) or SAM (text) file
bam = pysam.AlignmentFile("sample.bam", "rb")  # rb=read BAM, r=read SAM, rc=read CRAM

# Fetch reads overlapping a region (requires .bai index)
for read in bam.fetch("chr1", 10000, 20000):
    print(f"Name: {read.query_name}")
    print(f"  Position: {read.reference_start}-{read.reference_end}")
    print(f"  MAPQ: {read.mapping_quality}")
    print(f"  CIGAR: {read.cigarstring}")
    print(f"  Sequence: {read.query_sequence[:30]}...")
    break

# Count reads in region (fast, no iteration needed)
n_reads = bam.count("chr1", 10000, 20000)
print(f"Reads in region: {n_reads}")

# Filter reads by quality and flags
for read in bam.fetch("chr1", 10000, 20000):
    if read.mapping_quality >= 30 and not read.is_unmapped and not read.is_duplicate:
        pass  # Process high-quality, mapped, non-duplicate reads

bam.close()

# Write filtered reads to a new BAM file
with pysam.AlignmentFile("input.bam", "rb") as inbam:
    with pysam.AlignmentFile("filtered.bam", "wb", header=inbam.header) as outbam:
        for read in inbam.fetch("chr1", 10000, 20000):
            if read.mapping_quality >= 30:
                outbam.write(read)

# Index the output
pysam.index("filtered.bam")
print("Created filtered.bam + filtered.bam.bai")

2. Coverage and Pileup Analysis

Calculate per-base coverage statistics.

import pysam
import numpy as np

bam = pysam.AlignmentFile("sample.bam", "rb")

# Pileup: per-base coverage with read-level detail
for pileup_col in bam.pileup("chr1", 10000, 10100, min_mapping_quality=30):
    bases = [p.alignment.query_sequence[p.query_position]
             for p in pileup_col.pileups if not p.is_del and p.query_position is not None]
    print(f"Pos {pileup_col.reference_pos}: depth={pileup_col.nsegments}, bases={''.join(bases[:5])}")

# Quick coverage count per region (faster than pileup)
coverage = bam.count_coverage("chr1", 10000, 10100, quality_threshold=20)
# Returns tuple of 4 arrays (A, C, G, T counts per position)
total_cov = np.array(coverage).sum(axis=0)
print(f"Mean coverage: {total_cov.mean():.1f}x")

bam.close()

3. Variant Files (VCF/BCF)

Read, query, and filter genetic variants.

import pysam

# Open VCF/BCF file
vcf = pysam.VariantFile("variants.vcf.gz")

# Iterate all variants
for record in vcf.fetch("chr1", 10000, 50000):
    print(f"{record.chrom}:{record.pos} {record.ref}>{','.join(record.alts or [])}")
    print(f"  QUAL={record.qual}, FILTER={list(record.filter)}")
    print(f"  INFO: {dict(record.inf