FASTA, FASTQ, and Tabix-Indexed Files

This reference targets pysam 0.24.1.

Indexed FASTA

FastaFile provides random access through a faidx index.

import pysam

pysam.faidx("reference.fa")

with pysam.FastaFile("reference.fa") as fasta:
    print(fasta.references)
    print(fasta.lengths)
    print(fasta.get_reference_length("chr1"))

An uncompressed FASTA uses <name>.fai. With a default index path, FastaFile can build a missing index on open. Pass an existing filepath_index to prevent implicit index creation (as the bundled inspector does). A BGZF-compressed FASTA also needs a .gzi compressed-offset index. Ordinary gzip is not suitable for indexed random access.

The constructor can use nonstandard index paths:

fasta = pysam.FastaFile(
    "reference.fa.gz",
    filepath_index="indexes/reference.fa.gz.fai",
    filepath_index_compressed="indexes/reference.fa.gz.gzi",
)

Numeric and string coordinates

with pysam.FastaFile("reference.fa") as fasta:
    # Numeric: 0-based, half-open
    sequence = fasta.fetch("chr1", 999, 1_099)

    # Region string: 1-based, inclusive
    same_sequence = fasta.fetch(region="chr1:1000-1099")

If start or end is omitted, pysam uses the sequence boundary. Negative or inverted coordinates raise, but an end beyond the contig clips and a start at/beyond the end returns an empty string in 0.24.1. Check 0 <= start <= stop <= fasta.get_reference_length(contig) when exact length is required; absence of an exception is not a bounds check.

Fetch variant context

def variant_context(
    fasta: pysam.FastaFile,
    contig: str,
    pos_1based: int,
    ref: str,
    flank: int = 20,
) -> tuple[str, bool]:
    if pos_1based < 1 or flank < 0 or not ref:
        raise ValueError("positive POS, nonempty REF, and nonnegative flank required")
    start = pos_1based - 1
    if start + len(ref) > fasta.get_reference_length(contig):
        raise ValueError("REF extends beyond the contig")
    context_start = max(0, start - flank)
    context_stop = min(
        fasta.get_reference_length(contig),
        start + len(ref) + flank,
    )
    context = fasta.fetch(contig, context_start, context_stop)
    observed_ref = fasta.fetch(contig, start, start + len(ref))
    return context, observed_ref.upper() == ref.upper()

Validate the entire REF allele, not only its first base. A mismatch often means the VCF and FASTA use different assemblies, contig aliases, or normalization.

Strand-aware extraction

Coordinates do not encode strand. Reverse-complement after fetching:

sequence = fasta.fetch("chr1", start, stop)
if strand == "-":
    sequence = pysam.reverse_complement(sequence)

Confirm annotation coordinates before conversion: BED is normally 0-based half-open; GFF/GTF text is normally 1-based inclusive.

Sequential FASTA and FASTQ

FastxFile streams FASTA, FASTQ, or mixed FASTX records. It does not implement random access.

with pysam.FastxFile("reads.fastq.gz") as reads:
    for record in reads:
        print(record.name)
        print(record.sequence)
        print(record.comment)
        print(record.quality)

The current class is FastxFile; FastqFile is an old compatibility name.

Record persistence

persist=True is the default and copies each record so it remains valid. For high-throughput scans, persist=False avoids the copy:

with pysam.FastxFile("reads.fastq.gz", persist=False) as reads:
    for record in reads:
        process(record.name, record.sequence)
        # Do not save `record` for use after the iterator advances.

With persist=False, records are read-only proxies and cease to be valid after iteration advances. Copy needed strings or use the default when retaining records.

Quality scores

with pysam.FastxFile("reads.fastq") as reads:
    for record in reads:
        if record.quality is None:  # FASTA record
            continue
        qualities = record.get_quality_array()
        mean_quality = (
            sum(qualities) / len(qualities) if qualities else None
        )

Pysam converts the FASTQ quality string to numeric Phred scores. Confirm the source encoding for legacy FASTQ; modern data is normally Phred+33.

Safe writing

str(record) preserves whether the record is FASTA or FASTQ and includes the comment and quality line correctly:

with pysam.FastxFile("reads.fastq.gz") as source, open(
    "filtered.fastq", "x", encoding="utf-8"
) as destination:
    for record in source:
        qualities = record.get_quality_array()
        if qualities and sum(qualities) / len(qualities) >= 20:
            destination.write(str(record) + "\n")

Use exclusive creation ("x") when replacement was not requested. Manual four-line FASTQ formatting can lose comments or create sequence/quality length mismatches.

Streaming statistics

from __future__ import annotations


def fastx_stats(path: str) -> dict[str, float | int | None]:
    record_count = 0
    base_count = 0
    quality_sum = 0
    quality_count = 0

    with pysam.FastxFile(path, persist=False) as records:
        for record in records:
            record_count += 1
            base_count += len(record.sequence)
            if record.quality is not None:
                qualities = record.get_quality_array()
                quality_sum += sum(qualities)
                quality_count += len(qualities)

    return {
        "records": record_count,
        "bases": base_count,
        "mean_length": (
            base_count / record_count if record_count else None
        ),
        "mean_quality": (
            quality_sum / quality_count if quality_count else None
        ),
    }

This is a full scan but has constant memory.

BGZF and Tabix

TabixFile provides random access to coordinate-sorted, BGZF-compressed tabular files.

Non-destructive compression and indexing

import pysam

pysam.tabix_compress("regions.bed", "regions.bed.gz")
pysam.tabix_index("regions.bed.gz", preset="bed")

This leaves regions.bed intact. By contrast, calling tabix_index() directly on an uncompressed file may automatically compress it and remove the original unless keep_original=True.

Do not use force=True unless replacing existing compressed data or indexes is explicitly intended.

Input requirements

• Sort by contig and coordinate first. The wrapper does not sort or perform a full preflight; HTSlib rejects detected ordering errors during indexing.
• Use BGZF, not ordinary gzip.
• Select the correct preset: commonly bed, gff, sam, or vcf.
• Presets define columns and coordinate conventions.

For a custom table:

pysam.tabix_index(
    "custom.tsv.gz",
    seq_col=0,
    start_col=1,
    end_col=2,
    zerobased=True,
)

Python column indices are 0-based. File coordinates are assumed 1-based unless zerobased=True. This is separate from query coordinates, which are always numeric 0-based in the Python API.

Use CSI for references beyond legacy TBI limits, then pass its path to TabixFile(..., index="regions.bed.gz.csi") (the default only looks for TBI):

pysam.tabix_index(
    "regions.bed.gz",
    preset="bed",
    csi=True,
    min_shift=14,
)

Querying

with pysam.TabixFile(
    "regions.bed.gz",
    parser=pysam.asBed(),
) as regions:
    for row in regions.fetch("chr1", 1_000, 2_000):
        print(row.contig, row.start, row.end, row.name)

Numeric query coordinates are 0-based, half-open. Region strings are samtools-style 1-based inclusive.

Useful parsers:

• pysam.asTuple(): tuple-like fields
• pysam.asBed(): BED fields with 0-based start/end
• pysam.asGTF(): GTF fields and quoted attributes
• pysam.asGFF3(): GFF3 fields and key=value attributes
• pysam.asVCF(): lightweight tabix VCF parser

For complete VCF semantics, use VariantFile, not TabixFile(asVCF()).

Without a parser, each result is the raw tab-delimited string:

with pysam.TabixFile("annotations.tsv.gz") as table:
    for line in table.fetch("chr1", 1_000, 2_000):
        fields = line.split("\t")

Headers and multiple iterators

with pysam.TabixFile("annotations.gff.gz") as table:
    header_lines = list(table.header)
    first = table.fetch("chr1", multiple_iterators=True)
    second = table.fetch("chr2", multiple_iterators=True)

Header lines are yielded without trailing newlines. Each multiple_iterators=True iterator reopens the file and adds overhead.

Choosing the Right Interface

Common Pitfalls

• Mixing numeric 0-based coordinates with 1-based region strings
• Unexpected .fai creation when opening FastaFile without an explicit index
• Using ordinary gzip for indexed FASTA or tabix data
• Retaining a persist=False FASTX proxy
• Assuming every FASTX record has qualities
• Manually formatting FASTQ and losing comments or quality alignment
• Letting tabix_index() remove the uncompressed source unexpectedly
• Indexing an unsorted file
• Forgetting zerobased=True for custom 0-based table coordinates
• Using a TBI index for coordinates beyond its legacy range