Coordinates and Indexing

Coordinate mistakes and stale indexes are the most common causes of plausible but wrong genomic results. This reference targets pysam 0.24.0.

The Pysam Rule

For Python API numeric arguments and properties, pysam uses 0-based, half-open intervals:

[start, stop)

The first base is 0; start is included and stop is excluded. Interval length is stop - start.

The main exception is a textual samtools-style region string, which is 1-based, inclusive:

chr1:100-199

These refer to the same 100 bases:

file.fetch("chr1", 99, 199)
file.fetch(region="chr1:100-199")

This rule applies to:

• AlignmentFile.fetch(), count(), count_coverage(), and pileup()
• VariantFile.fetch()
• FastaFile.fetch()
• TabixFile.fetch()

Do not treat numeric VariantFile.fetch() arguments as VCF text coordinates.

Format Conversion Table

For VCF structural variants, symbolic alleles, breakends, and records with INFO/END, use VariantRecord.start and VariantRecord.stop rather than reconstructing the interval from len(REF).

Single Positions

A 1-based position p becomes the one-base Python interval:

start = p - 1
stop = p

For a VCF record:

assert record.start == record.pos - 1
base = fasta.fetch(record.contig, record.start, record.start + 1)

Overlap Versus Containment

Region fetches are overlap queries. An alignment or variant can begin before the requested interval and still overlap it.

To require complete containment:

def fully_contained(read, start: int, stop: int) -> bool:
    return (
        read.reference_start is not None
        and read.reference_end is not None
        and read.reference_start >= start
        and read.reference_end <= stop
    )

For point-based logic, define exactly what "overlap" means for deletions, reference skips, symbolic alleles, and breakends.

pileup() has an additional trap: without truncate=True, it can emit columns outside the requested interval because reads overlap the interval.

Parser Coordinates

Pysam parser objects normalize coordinates:

• asBed().start / .end: 0-based, half-open
• asGTF().start / .end: exposed in Python coordinate convention
• asVCF().pos: parser-specific lightweight field; use VariantFile for full VCF record semantics

When creating a custom tabix index:

• seq_col, start_col, and end_col are 0-based column indices
• file coordinates default to 1-based unless zerobased=True
• later TabixFile.fetch() numeric query coordinates are still 0-based

These are three distinct concepts: Python column index, coordinate encoding in the stored table, and coordinate encoding in the query.

Contig Identity

Coordinate conversion does not solve contig mismatches. Check:

• chr1 versus 1
• mitochondrial names (chrM, MT, M)
• alternate loci and decoys
• assembly version (for example GRCh37 versus GRCh38)
• contig order and lengths

alignment_contigs = dict(zip(bam.references, bam.lengths))
fasta_contigs = dict(zip(fasta.references, fasta.lengths))

shared = alignment_contigs.keys() & fasta_contigs.keys()
length_mismatches = {
    name: (alignment_contigs[name], fasta_contigs[name])
    for name in shared
    if alignment_contigs[name] != fasta_contigs[name]
}

Do not silently strip or add chr across arbitrary assemblies. Use an explicit reviewed mapping.

Index Matrix

An index is a view of a specific file. If the data file changes, rebuild the index. A stale index may fail loudly or return wrong/incomplete regions.

BAI/TBI Versus CSI

Legacy BAI and standard TBI indexes have a maximum coordinate near 2^29 (512 Mi bases). This is insufficient for some plant, animal, and synthetic references. CSI is parameterized and supports larger coordinates.

Create a BAM CSI:

import pysam

pysam.index("-c", "large-reference.bam", catch_stdout=False)

Create a tabix CSI:

pysam.tabix_index(
    "large-reference.bed.gz",
    preset="bed",
    csi=True,
    min_shift=14,
)

Create a VCF/BCF CSI:

import pysam.bcftools

pysam.bcftools.index(
    "--csi",
    "variants.vcf.gz",
    catch_stdout=False,
)

Prefer CSI when reference sizes are unknown or potentially large. Confirm downstream tools support it.

Sort Before Indexing

Indexing does not sort records.

BAM:

import pysam.samtools

pysam.samtools.sort(
    "-@", "4",
    "-o", "sorted.bam",
    "input.bam",
    catch_stdout=False,
)
pysam.samtools.index(
    "-@", "4",
    "sorted.bam",
    catch_stdout=False,
)

VCF:

import pysam.bcftools

pysam.bcftools.sort(
    "-Oz",
    "-o", "sorted.vcf.gz",
    "input.vcf",
    catch_stdout=False,
)
pysam.bcftools.index(
    "--csi",
    "sorted.vcf.gz",
    catch_stdout=False,
)

Tabix tables must be sorted before tabix_index(). The Python function does not verify sort order.

Safe Tabix Creation

Prefer separate compression and indexing:

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

Calling tabix_index("regions.bed") can automatically create regions.bed.gz and remove the original. Use keep_original=True if relying on that one-step path.

Do not pass force=True by default. Existing output should trigger review, not silent replacement.

Nonstandard and Remote Index Locations

Pass an explicit index:

with pysam.AlignmentFile(
    "sample.bam",
    "rb",
    index_filename="indexes/sample.csi",
) as bam:
    ...

with pysam.VariantFile(
    "cohort.vcf.gz",
    index_filename="indexes/cohort.vcf.gz.csi",
) as variants:
    ...

Remote random access additionally depends on:

• an HTSlib build with the relevant network/plugin support
• a reachable index
• server range requests
• stable data and index URLs

Pass index_filename explicitly when automatic URL derivation is unreliable. Read cram_and_performance.md before remote or CRAM access.

Index Checks

Alignment:

with pysam.AlignmentFile("sample.bam", "rb") as bam:
    if not bam.has_index():
        raise ValueError("random access requires a BAM/CRAM index")
    bam.check_index()

Variant and tabix constructors open a discovered index automatically; a region fetch fails when none is available. Reopen output and test known regions rather than checking only that an index filename exists.

Boundary Tests

For important pipelines, test:

• first base of a contig
• exact interval start and stop
• a record spanning the query boundary
• a zero-length or invalid interval
• contig end
• a high coordinate beyond 512 Mi bases when CSI is expected
• missing and aliased contigs
• region string and numeric equivalents

One useful invariant:

numeric = list(file.fetch(contig, start, stop))
region = f"{contig}:{start + 1}-{stop}"
textual = list(file.fetch(region=region))

For the same indexed file and valid nonempty interval, these queries should select equivalent records.