Prebuilt wheels are available for supported macOS and Linux platforms. A
source build needs a C compiler and HTSlib build dependencies; read the
official installation guide linked from references/sources.md.
First Decide
Before writing code:
Identify the real format, compression, sort order, and available index.
Decide whether coordinates are numeric Python coordinates or a region
string. Do not mix them.
For CRAM, identify the exact reference assembly and FASTA.
Prefer indexed region access; use sequential iteration only when intended.
Preserve headers when writing and write to a new path by default.
State filtering semantics: mapping/base quality, flags, overlap handling,
duplicate handling, and pileup depth cap.
For unfamiliar files, start with the bundled read-only inspector:
All scripts refuse to overwrite existing outputs. Run each with --help for
coordinate, index, and privacy notes.
Coordinate Contract
Numeric coordinates accepted by pysam APIs are 0-based, half-open. This
includes numeric AlignmentFile.fetch(), VariantFile.fetch(),
FastaFile.fetch(), TabixFile.fetch(), and pileup() arguments.
Region strings are samtools-style: 1-based and inclusive.
# The same 100 bases:
bam.fetch("chr1", 99, 199) # [99, 199)
bam.fetch(region="chr1:100-199") # 1-based inclusive
VCF text uses 1-based POS, while record properties expose both systems:
Read references/coordinates_and_indexing.md for format conversions, overlap
semantics, index choices, and contig-name checks.
Alignment Files
Use context managers and explicit modes:
import pysam
with pysam.AlignmentFile("sample.bam", "rb", threads=4) as bam:
for read in bam.fetch("chr1", 1_000, 2_000):
if (
not read.is_unmapped
and not read.is_secondary
and not read.is_supplementary
and read.mapping_quality >= 30
):
print(read.query_name, read.reference_start, read.cigarstring)
Use fetch(until_eof=True) to stream every record in file order, including
unplaced unmapped reads, without requiring an index:
with pysam.AlignmentFile("sample.bam", "rb") as bam:
for read in bam.fetch(until_eof=True):
...
Important distinctions:
fetch() returns alignment records overlapping a region.
count() counts records and defaults to read_callback="nofilter".
count_coverage() returns A/C/G/T base counts and defaults to base quality
15 plus read_callback="all".
pileup() exposes per-column reads and has its own filtering, base-quality,
overlap, orphan, and max_depth=8000 defaults.
For exact-region pileups, set truncate=True and explicit filters:
with pysam.FastaFile("reference.fa") as fasta, pysam.AlignmentFile(
"sample.bam", "rb"
) as bam:
for column in 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,
):
print(column.reference_pos, column.get_num_aligned())
Read references/alignment_files.md for flags, CIGAR operations, tags,
modified bases, writing records, pileup details, and iterator lifetime.
Variant Files
Input format is auto-detected. Numeric fetch coordinates remain 0-based:
import pysam
with pysam.VariantFile("cohort.vcf.gz", threads=4) as variants:
for record in variants.fetch("chr1", 999_999, 2_000_000):
print(record.contig, record.pos, record.ref, record.alts)
for sample_name, call in record.samples.items():
print(sample_name, call.get("GT"))
Subset samples before retrieving records:
with pysam.VariantFile("cohort.bcf") as variants:
variants.subset_samples(["sample_A", "sample_B"])
for record in variants:
...
When changing a header, copy each record and translate it to the destination
header before assigning newly declared INFO/FORMAT/FILTER fields. Do not
manually clear and rebuild header.samples.
Read references/variant_files.md for safe headers, writing, sample
subsetting, missing genotypes, symbolic alleles, filtering, translation, and
indexing.
FASTA, FASTQ, and Tabix
Indexed FASTA uses numeric 0-based coordinates:
with pysam.FastaFile("reference.fa") as fasta:
sequence = fasta.fetch("chr1", 999, 1_099)
FastxFile is sequential. persist=False is faster but yielded records become
invalid after iteration advances:
with pysam.FastxFile("reads.fastq.gz", persist=False) as reads:
for read in reads:
qualities = read.get_quality_array()
...
Tabix input must be coordinate-sorted and BGZF-compressed, not ordinary gzip.
Use a non-destructive two-step workflow:
pysam.tabix_compress("regions.bed", "regions.bed.gz")
pysam.tabix_index("regions.bed.gz", preset="bed")
with pysam.TabixFile("regions.bed.gz", parser=pysam.asBed()) as tbx:
for interval in tbx.fetch("chr1", 1_000, 2_000):
print(interval.contig, interval.start, interval.end)
Read references/sequence_files.md for FASTA/FASTQ records and safe tabix
creation.
CRAM, Remote I/O, and Threads
pysam 0.24 changed inherited HTSlib behavior:
Newly written CRAM defaults to CRAM 3.1, not 3.0.
HTSlib no longer contacts the EBI reference server by default.
Prefer reference_filename="reference.fa" for deterministic local reads and
writes.
with pysam.AlignmentFile(
"sample.cram",
"rc",
reference_filename="reference.fa",
threads=4,
) as cram:
for read in cram.fetch("chr1", 1_000, 2_000):
...
Only configure REF_PATH/REF_CACHE when reference-by-MD5 lookup is
intentional. Do not assume a CRAM is self-contained. threads= accelerates
compression/decompression; it does not parallelize Python analysis.
Read references/cram_and_performance.md before CRAM conversion, remote access,
or concurrent iteration.
Wrapped samtools and bcftools
Import command modules explicitly. Pass each command-line token as a separate
string:
Dispatchers capture stdout by default. For large or binary output, use the
tool's -o option with catch_stdout=False, or save_stdout=..., rather than
returning the complete output in memory.
try:
pysam.samtools.quickcheck("-v", "sample.bam")
except pysam.SamtoolsError as error:
messages = pysam.samtools.quickcheck.get_messages()
raise RuntimeError(messages or str(error)) from error
Use the Python API for record-level logic and dispatchers for mature bulk
operations such as sort, index, merge, view, and normalization. Never compose
dispatcher arguments by splitting an untrusted shell command.
Writing Rules
Copy or construct a valid header before opening output.
Write to a new path; do not use force=True unless replacement is explicit.
Preserve sort order if the output will be indexed.
Set query_sequence before query_qualities.
Prefer pysam.CIGAR_OPS enum members; top-level constants such as
pysam.CMATCH are compatibility aliases slated for future removal.
Validate outputs with pysam.samtools.quickcheck() for alignments and reopen
variant/sequence outputs before downstream use.
Use CSI rather than BAI/TBI when references or coordinates exceed legacy
index limits.
