Variant Files: VCF and BCF

This reference targets pysam 0.24.1. Numeric pysam coordinates are 0-based, half-open even though VCF text uses a 1-based POS.

Open and Iterate

VariantFile auto-detects VCF, BGZF-compressed VCF, and BCF input:

import pysam

with pysam.VariantFile("cohort.vcf.gz", threads=4) as variants:
    for record in variants:
        print(record.contig, record.pos, record.ref, record.alts)

Common modes:

• r: VCF input
• rb: BCF input; input auto-detection usually makes an explicit mode unnecessary
• w: VCF output; a .vcf.gz suffix selects BGZF-compressed VCF
• wb: compressed BCF output
• wz: explicitly BGZF-compressed VCF output
• wb0: BCF in BGZF blocks with compression level zero

The advertised wbu mode is rejected as conflicting format specifiers in 0.24.1; use wb0. Do not mistake level-zero BGZF for a raw unframed BCF stream.

Writing requires a VariantHeader.

Useful constructor options:

• index_filename= for a nonstandard or remote index
• drop_samples=True to skip sample data
• threads= for compression/decompression
• ignore_truncation=True for missing BGZF EOF markers; not compatible with threads > 1

Coordinates and Region Queries

with pysam.VariantFile("cohort.vcf.gz") as variants:
    # Numeric coordinates: [999_999, 2_000_000)
    numeric = variants.fetch("chr1", 999_999, 2_000_000)

    # Region string: 1-based inclusive
    text = variants.fetch(region="chr1:1000000-2000000")

VariantFile.fetch() explicitly defines numeric start/stop as 0-based, half-open. Region strings follow samtools notation.

For each VariantRecord:

• contig / chrom: contig name
• pos: 1-based VCF position
• start: 0-based inclusive position
• stop: 0-based exclusive record end
• rlen: reference span

Do not subtract one from a numeric start passed to fetch(). Use record.start when integrating with BAM, BED, or FASTA APIs.

Random access requires:

• BGZF VCF: .tbi or .csi
• BCF: .csi

An unindexed VCF.gz or BCF can still be read sequentially with normal iteration. In 0.24.1, fetch() with no contig/region rewinds to the start of records and returns a sequential iterator without requiring an index. Normal for record in variants iteration starts at the current position. This differs from the broad index-driven wording in the upstream fetch docstring; verified against tagged source and native VCF/BCF tests.

Use the canonical contig= keyword. In 0.24.1, supplying only the legacy reference= alias takes the no-contig sequential branch and ignores the intended region bounds.

For simultaneous iterators, VariantFile.fetch() uses reopen=True:

first = variants.fetch("chr1", reopen=True)
second = variants.fetch("chr2", reopen=True)

Header Model

with pysam.VariantFile("cohort.vcf.gz") as variants:
    header = variants.header
    print(list(header.samples))

    for name, contig in header.contigs.items():
        print(name, contig.length)

    for name, field in header.info.items():
        print(name, field.number, field.type, field.description)

Important collections:

• header.contigs
• header.samples
• header.filters
• header.info
• header.formats
• header.records

INFO and FORMAT values can only be interpreted safely when their definitions exist in the header. Do not infer missing Number/Type metadata.

Record Fields

for record in variants:
    alleles = record.alleles       # (REF, ALT1, ALT2, ...)
    alt_alleles = record.alts      # tuple or None
    identifier = record.id         # string or None
    quality = record.qual          # float or None
    filters = tuple(record.filter.keys())

record.alleles_variant_types classifies alleles with values such as REF, SNP, MNP, INDEL, BND, OVERLAP, and OTHER.

Handle non-simple alleles:

• multiallelic records can have several ALT alleles
• symbolic alleles include <DEL>, <DUP>, <INS>, and others
• paired breakends use bracket notation; single breakends use A. or .A
• spanning deletion uses *
• gVCF records often use <NON_REF> or <*>

Do not apply single-nucleotide logic to every record.

INFO Fields

for record in variants:
    depth = record.info.get("DP")
    frequencies = record.info.get("AF")
    if frequencies is not None:
        for allele_index, frequency in enumerate(frequencies, start=1):
            alt = record.alleles[allele_index]
            print(alt, frequency)

Number semantics matter:

• 1: one value
• A: one value per ALT allele
• R: one value per REF+ALT allele
• G: one value per genotype
• .: variable number

For FORMAT fields, 0.24.1 also supports P (ploidy), local-allele LA, LR, LG, and modification M cardinalities. Preserve the source VCF version and required companion fields such as LAA; these are not interchangeable with ordinary A/R/G. Consult VCF 4.5 and test the actual caller's schema.

Flag fields are represented as booleans. Missing values may be None, a tuple containing None, or absent from the mapping depending on the field.

FILTER Semantics

filters = tuple(record.filter.keys())

if "PASS" in filters:
    status = "pass"
elif not filters:
    status = "unfiltered"  # VCF '.'
else:
    status = "failed"

PASS, an empty filter set (.), and a failed filter are distinct states. Do not treat . as equivalent to PASS without an explicit policy.

Sample and Genotype Data

for sample_name, call in record.samples.items():
    genotype = call.get("GT")
    depth = call.get("DP")
    genotype_quality = call.get("GQ")
    phased = call.phased
    allele_strings = call.alleles

Genotype allele integers index record.alleles:

• 0: REF
• 1: first ALT
• 2: second ALT
• None: missing allele

Do not assume diploidy:

def called_alleles(genotype):
    if genotype is None:
        return ()
    return tuple(allele for allele in genotype if allele is not None)


called = called_alleles(record.samples["sample_A"].get("GT"))
alternate_count = sum(allele > 0 for allele in called)

call.phased records separator semantics but the Python GT remains a tuple.

Efficient Sample Subsetting

Call subset_samples() before retrieving any record:

samples = ["sample_A", "sample_B"]

with pysam.VariantFile("cohort.bcf") as source:
    source.subset_samples(samples)
    output_header = source.header.copy()

    with pysam.VariantFile(
        "subset.bcf", "wb", header=output_header, threads=4
    ) as destination:
        for record in source:
            destination.write(record)

This reduces decoding and memory. Do not copy a header, clear its sample collection, and manually reconstruct records; that approach is error-prone. Validate that every requested sample is present before subsetting.

Writing Unchanged Records

When the destination uses the input header:

with pysam.VariantFile("input.vcf.gz") as source, pysam.VariantFile(
    "passing.vcf.gz",
    "w",
    header=source.header,
    threads=4,
) as destination:
    for record in source:
        if "PASS" in record.filter:
            destination.write(record)

Write to a new file. Reopen and validate it before replacing any source.

Safely Add Header Fields

Records are tied to their originating header. When the destination header is changed, copy and translate records:

with pysam.VariantFile("input.vcf.gz") as source:
    output_header = source.header.copy()
    output_header.info.add(
        "BAM_DP",
        number=1,
        type="Integer",
        description="Aligned base depth from the selected BAM and filters",
    )

    with pysam.VariantFile(
        "annotated.vcf.gz", "w", header=output_header
    ) as destination:
        for input_record in source:
            record = input_record.copy()
            record.translate(output_header)
            record.info["BAM_DP"] = 27
            destination.write(record)

Declare INFO, FORMAT, FILTER, and contig metadata before assigning values. Use distinct field names when the new value has semantics different from an existing field.

Construct New Records

header = pysam.VariantHeader()
header.contigs.add("chr1", length=248_956_422)
header.info.add(
    "DP",
    number=1,
    type="Integer",
    description="Total depth",
)
header.formats.add(
    "GT",
    number=1,
    type="String",
    description="Genotype",
)
header.add_sample("sample_A")

with pysam.VariantFile("new.vcf.gz", "w", header=header) as output:
    record = output.header.new_record(
        contig="chr1",
        start=99_999,
        stop=100_000,
        alleles=("A", "G"),
        id="example",
        qual=60,
    )
    record.filter.add("PASS")
    record.info["DP"] = 40
    record.samples["sample_A"]["GT"] = (0, 1)
    output.write(record)

VariantHeader() already starts with ##fileformat=VCFv4.2. Do not append a second fileformat line with add_meta; that creates a duplicate rather than changing the version. To use newer-version fields, copy a valid header read from that version and check header.version.

start/stop here are numeric Python coordinates. start=99_999 writes VCF POS=100000.

Filtering Patterns

Make missing-value policy explicit:

def passes(record, *, min_qual=30.0, min_dp=10) -> bool:
    if record.qual is None or record.qual < min_qual:
        return False
    depth = record.info.get("DP")
    if depth is None or depth < min_dp:
        return False
    return "PASS" in record.filter

For ALT-frequency filtering:

frequencies = record.info.get("AF")
keep = (
    frequencies is not None
    and any(value is not None and value >= 0.01 for value in frequencies)
)

For genotype predicates, ignore missing alleles and preserve ploidy:

gt = record.samples["sample_A"].get("GT")
has_alt = gt is not None and any(
    allele is not None and allele > 0 for allele in gt
)

Compression and Indexing

Plain VCF must be coordinate sorted before compression/indexing:

import pysam

pysam.tabix_compress("sorted.vcf", "sorted.vcf.gz")
pysam.tabix_index("sorted.vcf.gz", preset="vcf")

The last call creates TBI. To choose CSI instead, compress first and then use bcftools indexing:

import pysam.bcftools

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

Do not use ordinary gzip for a random-access VCF. BGZF permits blocked random access.

tabix_index() can automatically compress an uncompressed input and delete the original unless keep_original=True; the explicit two-step pattern above is safer. Do not use force=True unless overwrite is intended.

For BCF:

pysam.bcftools.index("--csi", "cohort.bcf", catch_stdout=False)

Use CSI when contigs may exceed the legacy TBI coordinate limit. Read coordinates_and_indexing.md for index selection.

Header Translation Across Inputs

Do not merge VCFs by grouping Python records and appending sample calls. Inputs can differ in contig order, INFO/FORMAT definitions, normalization, and allele representation. Prefer bcftools merge, bcftools concat, or a dedicated variant-merging tool after normalization and header reconciliation.

record.translate(destination_header) remaps header dictionaries but does not normalize alleles, split multiallelic records, or resolve conflicting metadata.

Validation Checklist

• Reopen the written file with VariantFile.
• Confirm contig order and sample order.
• Confirm each assigned INFO/FORMAT/FILTER field is declared with correct Number and Type.
• Check missing, haploid, polyploid, multiallelic, symbolic, and filtered records.
• Index the final coordinate-sorted output and fetch known edge intervals.
• Use bcftools validation/normalization commands when format-level guarantees are required.