← Về thư mục
📄 / / profiles / r-and-d / raw_evidence / scientific-agent-skills / skills / pysam / references / sequence_files.md

FASTA, FASTQ, and Tabix-Indexed Files

This reference targets pysam 0.24.0.

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 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:

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

Fetch variant context

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:

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.

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

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

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:

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

Querying

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:

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

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