FASTA, FASTQ, and Tabix-Indexed Files

This reference targets pysam 0.24.0.

Indexed FASTA

FastaFile provides random access through a faidx index.

python

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 needs <name>.fai. 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:

python

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

python

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. Invalid regions raise ValueError or IndexError; do not silently clip unless that is the documented workflow.

Fetch variant context

python

def variant_context(
    fasta: pysam.FastaFile,
    contig: str,
    pos_1based: int,
    ref: str,
    flank: int = 20,
) -> tuple[str, bool]:
    start = pos_1based - 1
    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:

python

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


def reverse_complement(sequence: str) -> str:
    return sequence.translate(IUPAC_COMPLEMENT)[::-1]


sequence = fasta.fetch("chr1", start, stop)
if strand == "-":
    sequence = 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.

python

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:

python

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

python

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:

python

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

python

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

python

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. tabix_index() does not verify sort order.
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:

python

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:

python

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

Querying

python

with pysam.TabixFile(
    "regions.bed.gz",
    parser=pysam.asBed(),
    threads=4,
) 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/GFF-like fields and 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:

python

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

python

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

Data	Interface	Access
Reference FASTA	`FastaFile`	indexed random access
FASTA/FASTQ reads	`FastxFile`	sequential
BED/GFF/GTF/custom table	`TabixFile`	BGZF + tabix/CSI random access
VCF/BCF	`VariantFile`	structured records, sequential or indexed

Common Pitfalls

Mixing numeric 0-based coordinates with 1-based region strings
Expecting FastaFile to open without .fai
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