Genvarloader

GenVarLoader Public API

GenVarLoader (gvl) reconstructs personalized haplotypes and re-aligns functional genomic tracks on the fly from a reference + variants + BigWig/Table tracks, without writing personalized genomes to disk. Variable-length output is the norm (indels make lengths region- and sample-dependent).

This skill is a pointer-dense overview. Symbol names link to where to find the authoritative docstring or source.

End-to-end shape

import genvarloader as gvl

# 1. Preprocess variants outside Python (see "Variant preprocessing")
# 2. Write the dataset
gvl.write(
    path="ds.gvl",
    bed="rois.bed",
    variants="normed.bcf",  # or .pgen, or .svar directory
    tracks=[gvl.BigWigs.from_table("signal", "bw_table.tsv")],
    max_jitter=128,
)

# 3. Open and configure (chainable fluent API)
ds = (
    gvl.Dataset
    .open("ds.gvl", reference="ref.fa")
    .with_seqs("haplotypes")
    .with_tracks(["signal"])
    .with_insertion_fill(gvl.Repeat5pNormalized())
    .with_len(2048)  # or "ragged" / "variable"
    .with_settings(jitter=32, deterministic=False)
)

# 4. Eager indexing: dataset[region_idx, sample_idx]
batch = ds[0:8, :]  # shape depends on with_* state — see "Output shapes"

Variant preprocessing requirements

Variants passed to gvl.write must be left-aligned, bi-allelic, atomized (no MNPs or compound MNP-indels), and free of symbolic (<DEL>, <INS>, …) and breakend ALT alleles. gvl expands every ALT into literal haplotype sequence and cannot reconstruct symbolic or breakend records — gvl.write raises ValueError (with per-class counts of multi-allelic, symbolic, and breakend variants) if any are present, for VCF, PGEN, and SVAR inputs alike. VCFs must be indexed.

# VCF/BCF
bcftools norm -f ref.fa \
    -a --atom-overlaps . \
    -m -any --multi-overlaps . \
    -O b -o normed.bcf in.vcf.gz
bcftools index normed.bcf

# PGEN
plink2 --make-bpgen --pfile in --out tmp
plink2 --make-pgen --normalize --ref-from-fa --fa ref.fa --bpfile tmp --out normed

See docs/source/write.md for the canonical recipe and BED/BigWig table layouts.

When to use SVAR vs BCF/PGEN

.svar is a sparse columnar variant archive (from genoray). Pass it to gvl.write(variants="x.svar") exactly like a BCF or PGEN — the resulting dataset stores a back-reference instead of duplicating per-variant arrays.

Use SVAR when:

• You need allele-frequency filtering at read time (Dataset.open(min_af=..., max_af=...) is supported for .svar only — a .svar2-backed dataset raises NotImplementedError).
• Many datasets share the same variant source — SVAR avoids duplicating variant_idxs.npy/dosages.npy/variants.arrow into each .gvl directory.
• You're working at population scale and want compact on-disk variant storage.

Use BCF/PGEN directly when you have a one-off dataset and don't need AF filtering.

Create an SVAR from a normalized VCF/PGEN with genoray:

from genoray import VCF, SparseVar

SparseVar.from_vcf("normed.svar", VCF("normed.bcf"), max_mem="4g")  # writes a .svar/ directory

SVARs are resolved at Dataset.open time via metadata.json → caller svar= arg → recorded relative path → recorded absolute path → sibling *.svar. See docs/source/format.md ("SVAR resolution at open time") and _dataset/_svar_link.py. Legacy symlink-based SVAR layouts: run gvl.migrate_svar_link(path) once to upgrade.

.svar2 — the read-bound sparse variant format

.svar2 is genoray's newer sparse columnar variant store. Pass it to gvl.write exactly like a .svar, BCF, or PGEN — gvl.write(path, bed, variants="cohort.svar2") or variants=SparseVar2("cohort.svar2"). Like .svar, the dataset stores a back-reference (metadata.json → svar2_link) instead of duplicating per-variant arrays, so the .svar2 store must remain accessible at read time.

Unlike .svar (whose read path builds an interval search tree + a per-read dense-union over the queried window), a .svar2-backed dataset reconstructs via a read-bound path: gvl.write caches per-(region, sample, ploid) variant-key ranges under <dataset>/genotypes/svar2_ranges/ (sized to the dataset's selected samples, not the full .svar2 cohort) — not small at cohort scale, see the "Common gotchas" bullet below — and at read time gvl gathers directly off that cache and calls all-Rust kernels — no interval-search-tree build and no dense-union rebuild per read. .svar2 stores are also typically smaller on disk than .svar, especially for large cohorts. See docs/source/faq.md.

.svar2 is resolved at Dataset.open time in the same order as .svar: caller svar2= arg → recorded relative path → recorded absolute path → sibling *.svar2. Dataset.open(path, svar2=<override>) mirrors svar=. See docs/source/format.md (".svar2 resolution at open time").

Phase-1 scope — unsupported combinations raise NotImplementedError. .svar2-backed datasets support all four output modes (haplotypes, variants, variant-windows, and haplotype-realigned tracks) byte-identical to the .svar/union-oracle backend (except pure-deletion ALT bytes — see below), and with_seqs("variant-windows") (ref="window", alt ∈ {"window", "allele"}), unphased_union (for both "variants" and "variant-windows" output), and var_fields-selected store INFO/FORMAT fields (also for both "variants" and "variant-windows"; see var_fields under Dataset.open below) are all fully wired for .svar2. Plain haplotype output also supports splicing, var_filter="exonic", and automatic reverse-complementation of negative-strand regions. The following are still not yet wired and raise a clear error instead of silently mis-computing:

• Splicing with variants, variant-windows, or haplotype-realigned tracks.
• var_filter="exonic" with non-haplotype output or haplotype-realigned tracks.
• min_af / max_af filtering.
• annotated haplotypes (with_seqs("annotated")).
• VarWindowOpt(ref="allele") (bare-allele REF mode; blocked upstream of .svar2 too — REF alleles aren't stored).
• Reverse-complement with haplotype-realigned tracks.
• Fixed-length (integer output_length) haplotype-realigned track output (plain haplotype output at a fixed length is fine — only the track kernel is guarded).
• variants / variant-windows output on a dataset written with max_jitter>0 or read with jitter>0 (the read-bound decode does not right-clip to the post-jitter window; haplotypes and tracks are unaffected and support jitter fully).
• gvl.write(..., extend_to_length=False) for a .svar2 variant source (write-time; raises NotImplementedError — .svar2 sources must use the default extend_to_length=True).

(FlankSample insertion-fill for tracks spanning multiple contigs in one query is now supported and byte-identical to the .svar backend — issue #267.)

variants/variant-windows ALT bytes differ from .svar for pure deletions (format convention, not a bug). For a pure deletion (e.g. VCF GTA>G), with_seqs("variants") on a .svar dataset yields the VCF anchor base as ALT (b"G"), while a .svar2 dataset yields the atomized empty ALT (b"") — this is how genoray's .svar2 format represents pure deletions. The same convention carries into with_seqs("variant-windows"): ref_window is byte-identical between .svar/.svar2, but alt/alt_window differ for pure-deletion records for the same reason. Reconstructed haplotypes are byte-identical between the two backends (both consume the ALT identically when building sequence); only the raw allele/window bytes differ for pure-deletion records. See docs/source/faq.md.

Symbolic/breakend variants are rejected the same as .svar, but for .svar2 the rejection happens upstream, at .svar2 conversion time (the store format cannot represent them) — a .svar2 must be built from an already-filtered source; gvl cannot re-filter a materialized .svar2 any more than it can a materialized .svar.

gvl.write — key arguments

gvl.write(
    path,
    bed,
    variants=None,
    tracks=None,
    annot_tracks=None,
    samples=None,
    max_jitter=None,
    overwrite=False,
    max_mem="4g",
    extend_to_length=True,
)

Notable:

• bed: path or polars DataFrame with chrom, chromStart, chromEnd (0-based). Optional strand (+/-/.) controls reverse-complement on read. Extra columns are preserved on Dataset.regions.
• tracks: a gvl.BigWigs (or a list of them), or a gvl.Table. Each must have a unique .name. BigWigs need a sample→path mapping (dict or table with sample, path columns; see BigWigs.from_table). gvl.Table is a core interval-track source backed by a Rust COITrees overlap engine (zero-based half-open coordinates); pass it directly as a tracks= or annot_tracks= source in gvl.write.
• annot_tracks: dict[str, str | Path | pl.DataFrame | pl.LazyFrame] | None — sample-independent annotation tracks, written to <path>/annot_intervals/<name>/. Each value is either a path to an interval table/bigWig file, or a polars DataFrame/LazyFrame with BED-like columns (chrom, chromStart, chromEnd, score). Annotation tracks are sample-independent and can be read without a per-sample variant source.
• max_jitter: max read-time jitter; pads stored data on both sides of every region by this many bases so Dataset.with_settings(jitter=j) works for any j <= max_jitter.
• extend_to_length=True keeps reading past the BED end until every haplotype is ≥ the region length (matters when deletions would shorten output); set False for faster writes if shorter haps are acceptable. Not supported for a .svar2 variant source — extend_to_length=False raises NotImplementedError there; only BCF/PGEN/.svar sources may disable it.
• samples: which samples to include; None (default) takes every sample available across variants and all tracks. Either way the dataset's sample order is the lexicographic sort of that selection, which is not the numeric order a phenotype table usually carries ("1000" sorts before "999"). Align external tables to Dataset.samples by name, never by position.
• Inner-joins samples across variants and all tracks.

Parallelism: gvl.write now parallelizes over write categories. Variants are processed first (serially). Then per-sample tracks and annot_tracks run concurrently (joblib loky backend). The max_mem budget is divided across the concurrently-running categories.

max_mem and .svar2: for a .svar2 variant source, max_mem also bounds the genotype range-cache write — ranges are produced in per-sample chunks sized to fit the budget rather than a whole contig at once. It does not bound the permanent genotypes/svar2_ranges/ cache's on-disk size; that scales with regions x samples x ploidy and is governed by disk space (see "Common gotchas" below and format.md).

Source: python/genvarloader/_dataset/_write.py.

Atomic creation: gvl.write builds into a private sibling temp directory and publishes via an atomic os.replace. A best-effort filelock avoids redundant rebuilds when parallel jobs share the same destination, but correctness relies on the rename — the lock is advisory only. Datasets do not auto-rebuild; if the on-disk artifact is missing or corrupt, re-run gvl.write.

Out-of-scope: genoray .gvi index files and pysam .fai/.gzi index files are created by those libraries and are not covered by gvl's atomic/locked creation. Concurrent jobs that trigger index creation for those files depend on the upstream libraries' behavior.

gvl.update — add tracks to an existing dataset

gvl.update(
    dataset,          # str | Path | Dataset
    tracks=None,      # BigWigs | Table | Sequence[BigWigs | Table] | None
    annot_tracks=None,  # dict[str, str | Path | pl.DataFrame | pl.LazyFrame] | None
    *,
    overwrite=False,
    max_mem="4g",
) -> None

Adds tracks to an existing on-disk GVL dataset without rewriting it from scratch.

• dataset: path to a dataset directory, or an opened Dataset (its .path is used). A live dataset can be read during update; it will not observe the new track until reopened.
• tracks: per-sample BigWigs or Table sources. The track's sample set must match the dataset's exactly (no missing, no extra); samples are reordered to dataset order automatically. Written to <path>/intervals/<track>/.
• annot_tracks: sample-independent sources, identical to gvl.write's annot_tracks (path to interval table, path to bigWig, or polars DataFrame/LazyFrame with BED-like columns). Written to <path>/annot_intervals/<name>/.
• overwrite=True: replace a same-named existing track; False (default) raises FileExistsError if the name already exists.
• max_mem: approximate memory budget, divided across concurrently-running categories.

Each track subdirectory is published atomically (built into a temp sibling, then os.replaced into place), so a reader can never see a half-written track.

Source: python/genvarloader/_dataset/_write.py.

gvl.concat — merge on-disk datasets

gvl.concat(
    path,              # str | Path — destination directory
    datasets,          # Sequence[str | Path | Dataset] — 2+ inputs
    axis,               # "regions" | "samples"
    *,
    overwrite=False,
    max_mem="4g",       # advisory; accepted for symmetry with gvl.write
) -> None

Merges datasets that all share one variant source — the same PGEN/VCF variant table (checked via a content fingerprint, not just backend type) or the same .svar/.svar2 store; merging datasets built from different variant sources raises ValueError. Requires at least two inputs.

• axis="regions": inputs must have identical samples in identical order; their regions are concatenated (sorted by contig/position) into the output.
• axis="samples": inputs must have identical regions and disjoint sample sets; their samples are merged into sorted order.

Both axes merge per-sample tracks and annotation tracks alongside the genotypes. Supported backends: PGEN/VCF, .svar, .svar2. For a PGEN/VCF source, variants.arrow is hardlinked from the first input and a fingerprint is recorded on the merged metadata.json; Dataset.open re-verifies it, so an out-of-band rewrite of the variant index raises instead of silently misresolving.

Cost: gvl.concat streams bytes rather than re-deriving anything, so its cost tracks I/O, not compute — it moves roughly the full size of the merged dataset. Only worth it against the alternative of re-extracting genotypes from scratch, not as a routine step. See "Merging datasets" in docs/source/write.md for the full cost discussion and the genoray.SparseVar2.concat store-level alternative for contig-sharded .svar2 workflows.

Source: python/genvarloader/_dataset/_concat.py, _concat_validate.py.

Dataset.open — key arguments

gvl.Dataset.open(
    path, reference=None, jitter=0, rng=None,
    deterministic=True, rc_neg=True,
    min_af=None, max_af=None,           # .svar only — raises NotImplementedError on .svar2
    region_names=None,
    splice_info=None,                    # see "Spliced haplotypes"
    var_filter=None,                     # None | "exonic"
    var_fields=None,                     # list[str] | None — see below
    *, svar=None, svar2=None,
)

Without reference=, a genotypes-only dataset opens with the "variants" view by default (yielding RaggedVariants) — Dataset.open(path) just works, no with_seqs needed. The "haplotypes", "annotated", and "reference" views all require a reference; requesting one via with_seqs on a reference-less dataset raises a clear ValueError. svar= overrides the recorded SVAR location; svar2= mirrors it for a .svar2-backed dataset.

with_settings(dummy_variant=...) — inserts a gvl.DummyVariant into every empty (region, sample, ploid) variant group so that every group has at least one variant. Only fills groups that are empty; non-empty groups are unchanged. Valid for both "variants" and "variant-windows" output kinds; indexing raises ValueError if dummy_variant is set and the output kind is any other kind ("haplotypes", "annotated", "reference", or no seqs) — the check is order-independent with with_seqs. False disables dummy padding; None (default) leaves the current setting unchanged. Setting dummy_variant=False when the output is an unsupported kind is a harmless no-op.

For token outputs (ride-along FlatVariants.flank_tokens and "variant-windows" fields ref_window/alt_window/bare ref/alt), each empty group's dummy entry is filled entirely with the configured unknown_token (the user-supplied integer that out-of-alphabet bytes map to — alphabet-agnostic, so N for DNA/RNA, X for amino acids tokenize to it). The DummyVariant.ref/.alt bytes only determine the dummy allele's byte-length, not the token value. Fill lengths per empty (region, sample, ploid) group: ride-along flank_tokens → 2·flank_length unknown tokens; ref_window/alt_window → 2·flank_length + len(dummy allele) unknown tokens; bare ref/alt → len(dummy allele) unknown tokens. (Default b"N" allele → length 1.)

Scalar fields (start/ilen/dosage/info[...]) are still filled from DummyVariant values as before. The dummy fill applies before reverse-complementing on the "variants"-output alleles only, so a non-b"N" dummy allele is reverse-complemented on negative-strand regions like any real allele — the default b"N" is rc-invariant. Flank tokens and the variant-window token buffers are reference-oriented and are NOT RC'd.

with_settings(unphased_union=...) — fold the stored diploid haplotypes onto a single haploid sequence: the union of called ALTs per (region, sample). When True, ds.ploidy reports 1 (instead of the stored 2); n_variants(...) reports a single ploidy slot (shape (..., 1)), with counts equal to the naive per-haplotype sum (a hom call appears twice — once per haplotype — with no dedup). "variants" and "variant-windows" output decode at ploidy 1; ALT occurrences are concatenated across haplotypes with no sort and no dedup. Phase is discarded — intended for haploid somatic modeling of unphased somatic calls. Requires a dataset with genotypes (raises ValueError on reference-only datasets). Incompatible with "haplotypes" / "annotated" output — with_seqs("haplotypes") or with_seqs("annotated") (or setting this flag while one of those is the active output kind) raises ValueError. See issue #222.

Format validation: Dataset.open validates the dataset's format_version and structural integrity (file presence + sizes). A corrupt dataset raises a ValueError instructing regeneration with gvl.write. Datasets do not auto-rebuild.

Format version gate (2.0): the current on-disk format is 2.0.0. Opening a dataset written by genvarloader < 2.0 (or any unversioned dataset) raises a ValueError whose message points at gvl.migrate(path); a dataset written by a newer major raises a ValueError telling you to upgrade genvarloader. Run gvl.migrate(path) once to upgrade a pre-2.0 dataset in place — it is streaming (peak extra disk is one track's interval store), idempotent, and crash-safe (metadata is bumped only after every track's struct-of-arrays files are durable, then the old array-of-structs files are deleted). It converts the track-interval storage only; genotypes, regions, and reference are untouched.

• var_fields: list[str] | None — Variant fields to include on RaggedVariants output. Defaults to the minimum useful set ["alt", "ilen", "start"]. Pass additional names (e.g. "ref", "dosage", or any numeric info column in the source variants table) to load them eagerly at open time. Must be a subset of Dataset.available_var_fields. Can be reconfigured later via Dataset.with_settings(var_fields=...), which lazily loads any newly-requested columns. "dosage" must be requested explicitly — it is not added automatically even when dosages.npy exists on disk. Beyond the built-ins (alt, start, ref, ilen, dosage) and per-variant INFO columns, a genoray .svar may register arbitrary per-call (Number=G) FORMAT fields in <svar>/metadata.json["fields"]; these appear in Dataset.available_var_fields and can be requested via Dataset.open(..., var_fields=[...]) or with_settings(var_fields=[...]). Each surfaces in variants, variant-windows, and flat outputs as a per-call ragged field aligned with the genotypes. A FORMAT field shadows a same-named INFO column.

  On .svar2, available_var_fields is ["alt", "ilen", "start"] plus whatever store fields the .svar2 was written with — "ref" and "dosage" are not available fields for a .svar2 source and requesting either raises (they are only valid for BCF/PGEN/.svar). Beyond that, var_fields additionally exposes the store's own scalar-numeric INFO/FORMAT fields — whichever ones the .svar2 was written with via genoray.SparseVar2.from_vcf(info_fields=[...], format_fields=[...]) (bare str names also work there). Only scalar-numeric fields can exist in a .svar2 store at all — INFO/FORMAT Type=Integer/Float with Number=1 or Number=A, plus INFO Type=Flag (stored as bool); anything else is rejected by genoray at write time and never reaches gvl. gvl only reads whatever the store already carries — it cannot add fields, and re-requesting a field the store doesn't have raises (it isn't in available_var_fields). Dataset.available_var_fields advertises each store field's key, sourced from genoray.SparseVar2.available_fields: the bare field name when it's unique across the store's INFO/FORMAT namespace, else "INFO/<name>" / "FORMAT/<name>". A builtin name (alt/start/ref/ilen/dosage) always wins — a store field that happens to be named e.g. alt is never advertised and cannot shadow the builtin. Values keep the store's dtype exactly, with no widening (an i32 field decodes int32, an f32 field float32), and a VCF-missing entry carries the store's stored default verbatim (NaN for a float field declared with no default). Both entry points route to the same svar2 reconstructor: gvl.Dataset.open(path, reference=..., var_fields=[...]) and ds.with_seqs("variants")/.with_settings(var_fields=[...]). Supported on both output modes: "variants" (the field appears on the returned RaggedVariants, e.g. rv["AF"], sharing alt/start/ilen's variant offsets) and "variant-windows" (the field appears in win.fields["AF"] alongside start/ilen). A FORMAT field's value is the value for the sample that row belongs to (not sample 0). Empty (region, sample, ploid) groups fill each store field via the same DummyVariant.info[<key>] mechanism as any other scalar field (see with_settings(dummy_variant=...) above): the user-supplied value if given, else NaN for a float column or 0 for an integer column.

Output modes — with_seqs × with_tracks

with_seqs(kind) selects the sequence output channel:

"variant-windows" requires a VarWindowOpt second argument (with_seqs("variant-windows", gvl.VarWindowOpt(...))), with_output_format("flat"), and a reference genome. It does not inherit flank settings from with_settings.

variant-windows may be combined with tracks when with_settings(realign_tracks=False) is set; the returned tracks/intervals are reference-coordinate (as-is). Float tracks come back as FlatRagged, interval tracks as FlatIntervals.

with_tracks(tracks=..., kind=...) selects tracks:

• tracks: None (default), False (disable), a single name, or a list of names.
• kind: "tracks" (re-aligned numeric values) or "intervals" (raw interval representation).

Track re-alignment to haplotype coordinates is controlled by with_settings(realign_tracks=True) (default). Set realign_tracks=False for reference-coordinate ("as-is") tracks. realign_tracks=False is required for kind="intervals" with any variant-aware seq mode, and for variant-windows + tracks. with_insertion_fill requires realign_tracks=True.

with_settings(parallel=...) — parallelism policy for this dataset's reads. True forces the Rust kernels multithreaded, False forces them serial, "auto" (default) decides per batch from the output size. An explicit True/False overrides the GVL_FORCE_PARALLEL environment variable — precedence is explicit setting > environment > size gate — so a script's parallelism can be read off the script rather than depending on ambient state. "auto" defers to the environment, so datasets that never set it behave as before. The setting is per-dataset and travels with it into dataloader worker processes. It governs whether to parallelize, not thread count: the worker count is fixed at import from GVL_NUM_THREADS because rayon reads it at global-pool init, so it cannot vary per dataset. See issue #352.

with_len(L) controls output shape:

• "ragged" (default): returns gvl.Ragged (variable length per item).
• "variable": NumPy array right-padded to the batch's longest item (N for seqs, 0 for tracks).
• integer L: fixed length; jitter/random shift/truncate/pad-with-more-personalized-data combine to meet L. Must satisfy L + 2·jitter ≤ min(region_length) + 2·max_jitter.

Returns either a RaggedDataset or ArrayDataset (frozen dataclass views) based on with_len. See docs/source/dataset.md for diagrams.

with_output_format(fmt) selects the container type returned by eager indexing:

In "flat" mode the hot path is zero-awkward; the returned containers carry .data (flat numpy array) and .offsets (int64). Every flat type has .to_ragged() back to its _core.Ragged-backed form. Densification escape hatches vary by type: FlatRagged has .to_fixed(length) and .to_padded(pad_value); FlatAnnotatedHaps has .to_fixed(length) and .to_padded() (no arg — uses per-field pad defaults); FlatVariants/FlatAlleles expose only .to_ragged() (plus .reshape/.squeeze).

ds_flat = ds.with_output_format("flat")
result = ds_flat[0:8, :]   # FlatRagged or FlatAnnotatedHaps or FlatVariants
# direct tensorization — no awkward round-trip
import torch
t = torch.from_numpy(result.data)
# or convert back
ragged = result.to_ragged()

with_output_format is orthogonal to and composes with with_len and subset_to.

Scope note: In "flat" mode, kind="intervals" tracks return FlatIntervals (.to_ragged() → RaggedIntervals; fields .starts/.ends/.values are each a FlatRagged). Float (numeric) tracks return FlatRagged in flat mode. Seqs/haplotypes/annotated-haps/reference and variants outputs are also flattened as before.

Flat variants extras — ride-along flank tokens and variant windows:

import seqpro as sp
import genvarloader as gvl

# Both paths need tracks disabled — the flat variants/windows channel is not
# produced when tracks are active (see gotchas).

# (a) ride-along flank tokens on the "variants" output
fv = (ds.with_tracks(False).with_seqs("variants").with_output_format("flat")
        .with_settings(flank_length=128, token_alphabet=sp.DNA.alphabet,
                       unknown_token=len(sp.DNA)))[0:8]
fv.flank_tokens          # FlatRagged, shape (b, p, ~v, 2*128), or None if not configured

# (b) per-allele windows: ref as a flanked window, alt as a bare tokenized allele
fw = (ds.with_tracks(False).with_output_format("flat")
        .with_seqs("variant-windows",
                   gvl.VarWindowOpt(flank_length=128, token_alphabet=sp.DNA.alphabet,
                                    unknown_token=len(sp.DNA), ref="window", alt="allele")))[0:8]
fw.ref_window            # flanked ref window tokens (two-level token buffer)
fw.alt                   # bare alt allele tokens (no flanks); fw.alt_window is None
fw.ref_window.shape      # the window buffer's own shape: (b, p, ~v, ~len)

Ride-along FlatVariants.flank_tokens (with_seqs("variants") + with_settings(flank_length=L, token_alphabet=..., unknown_token=...)): appends a FlatRagged of shape (b, p, ~v, 2L) to the returned FlatVariants. Per variant the buffer holds [flank5 | flank3] reference-context tokens (each L long). Coordinate rule: flank5 = [start-L, start), flank3 = [end, end+L) where end = start - min(ilen, 0) + 1. token_alphabet (str, bytes, or seqpro.NucleotideAlphabet — e.g. sp.alphabets.DNA / sp.DNA.alphabet; normalized to bytes at the with_settings/build_token_lut boundary) and unknown_token (int) together build a 256-entry byte→token LUT (seqpro-style): each alphabet byte → its 0-based index; every other byte (including N and out-of-bounds padding) → unknown_token. flank_length=0/None disables; both token_alphabet and unknown_token must be set together. Token dtype is uint8 when max token id ≤ 255, else int32; offsets are int64. When with_settings(dummy_variant=...) is set, each empty (region, sample, ploid) group's flank_tokens row is a 2L-long run of unknown_token.

VarWindowOpt / FlatVariantWindows (with_seqs("variant-windows", opt)): each variant gets a fixed-length token buffer in two modes selected independently for ref and alt via VarWindowOpt.ref / VarWindowOpt.alt ∈ {"window", "allele"}:

• "window": flanked + tokenized — ref-window = tokenized [start-L, end+L) reference read; alt-window = tokenized flank5 · alt-allele · flank3 assembly.
• "allele": the bare tokenized allele (ref or alt bases) with no flanks.

FlatVariantWindows sets exactly one of .ref_window / .ref (the other is None) and one of .alt_window / .alt (the other is None). .fields is a dict of scalar FlatRagged (start/ilen/dosage/info; raw byte alleles are dropped). Flanks and windows are reference-oriented — NOT reverse-complemented even when rc_neg=True. Splicing is not supported with "variant-windows". When with_settings(dummy_variant=...) is set, each empty (region, sample, ploid) group is padded with one all-unknown_token entry: length 2·flank_length + len(dummy allele) tokens for ref_window/alt_window ("window" mode), or len(dummy allele) tokens for bare ref/alt ("allele" mode).

Track insertion fill (only when haps + tracks together)

Indels make track length differ from reference length. Dataset.with_insertion_fill(fill) controls what gets written into inserted positions. Only valid when the dataset returns both haplotypes and tracks — pure-ref and pure-hap datasets ignore it (raises if attempted).

Pass a single strategy (applies to every track) or a dict[track_name, strategy] (missing tracks fall back to Repeat5p). Source: python/genvarloader/_dataset/_insertion_fill.py.

Spliced haplotypes

Splicing is opt-in at Dataset.open (or via with_settings). It groups the BED rows for one transcript and concatenates exon-level sequences/tracks per sample.

splice_bed = gvl.get_splice_bed("annotation.gtf", transcript_support_level="1")
gvl.write(path="splice.gvl", bed=splice_bed, variants="normed.svar")

sds = gvl.Dataset.open(
    "splice.gvl",
    reference="ref.fa",
    splice_info=("transcript_id", "exon_number"),  # tuple = (group_col, order_col)
    var_filter="exonic",  # optional: drop intronic variants
)

splice_info accepts:

• a column name string (single grouping column, order inferred from BED row order), or
• a (group_col, order_col) tuple (explicit ordering, e.g. exon number).

get_splice_bed does GTF→BED with TSL filtering and an optional "CDS length multiple of 3" filter. To roll your own splice BED, just include transcript_id (or any grouping column) and exon_number columns on the BED. See docs/source/splicing.ipynb.

RefDataset splicing

gvl.RefDataset accepts the same splice_info argument as Dataset.open. Pass either a transcript-ID column name (rows already in splice order) or a (group_col, sort_col) tuple to reorder exons. with_settings(splice_info=False) disables splicing on an existing RefDataset; pass a new value to re-enable. Splicing requires output_length in {"ragged", "variable"}, jitter=0, and deterministic=True. subset_to(transcript_ids) works the same as for Dataset.

ref = gvl.Reference.from_path("hg38.fa.bgz")
bed = gvl.get_splice_bed("annotations.gtf")
ref_ds = gvl.RefDataset(ref, bed, splice_info="transcript_id")
seqs = ref_ds[:]  # Ragged[S1], one row per transcript

Site-only variants (e.g. ClinVar)

Use gvl.sites_vcf_to_table(vcf) → pl.DataFrame (bi-allelic SNPs only), then wrap an ArrayDataset[AnnotatedHaps, ...] with gvl.DatasetWithSites(ds, sites, max_variants_per_region=1). Returns (wt_haps, mut_haps, flags[, tracks]); flags encode applied / deleted-overlap / already-existing. See _variants/_sitesonly.py.

Prefetching dataloader (mode=... on to_dataloader)

Dataset.to_dataloader() accepts an optional mode to coarsen fetching: gvl's throughput scales with fetch size (internal multithreading amortizes overhead), so one big dataset[r_idx, s_idx] call sliced into mini-batches outperforms many small per-batch calls.

loader = ds.to_dataloader(
    batch_size=32,
    mode="double_buffered",   # or "buffered", or None
    buffer_bytes=2 * 1024**3, # total RAM budget; split across slots in double mode
    copy=True,                # zero-copy opt-out (default True = safe)
    heartbeat_seconds=60.0,   # double_buffered: max wait per chunk before liveness check
)

Modes:

• None (default) — plain torch.utils.data.DataLoader; existing behavior.
• "buffered" — main process fetches one chunk per call (sized to buffer_bytes), slices into mini-batches. Refill latency visible but amortized over many batches.
• "double_buffered" — subprocess producer fills one shm slot while consumer drains the other; refill latency hidden. Requires a file-backed Dataset.open(path).

Flat output composes with buffered modes. ds.with_output_format("flat").to_dataloader(mode="buffered" | "double_buffered") yields Flat* mini-batches (FlatRagged / FlatVariants / FlatAnnotatedHaps / FlatVariantWindows) instead of Ragged / RaggedVariants, with zero awkward on the transport path — the double_buffered producer writes flat buffers and the consumer reads them back without re-wrapping. Densify each batch with .to_fixed(length) / .to_padded(pad), or wrap batch.data / batch.offsets with torch.from_numpy (offsets are int64; cast to int32 for a torch Nested tensor). Re-wrap to the _core.Ragged-backed types with .to_ragged() (element-identical to ragged mode). The mode's preconditions are unchanged: double_buffered still requires a file-backed Dataset.open(path), still rejects spliced datasets and non-default insertion_fill, and haplotype/annotated output still needs deterministic=True. "variant-windows" output (with_seqs("variant-windows", VarWindowOpt(...))) and flat "variants" output carrying ride-along flank tokens (with_settings(flank_length=...)) are supported over both mode="buffered" and mode="double_buffered", byte-identical to mode=None — including with dummy_variant set (the producer subprocess replays it, and the empty-group dummy fill is counted in the double_buffered slot sizing).

Preconditions (all raise ValueError at construction):

• with_seqs in {"haplotypes", "annotated"} requires deterministic=True (set via with_settings(deterministic=True)). reference and variants modes have no determinism requirement.
• Spliced datasets are not supported.
• num_workers > 0 is rejected — the new loader IS the concurrency strategy.
• A single mini-batch whose exact footprint exceeds the per-slot capacity raises with the offending size and remediation knobs (batch_size↓, buffer_bytes↑).

Footprint is computed exactly via Dataset._output_bytes_per_instance(...) (uses haplotype_lengths, n_variants, and allele offset tables) — no Zipf-style worst-case slack.

Other public surface (one-liners)

• gvl.Reference.from_path(fasta, contigs=None, in_memory=True) — wrap a FASTA (path to a .fa/.fa.bgz, or a .gvlfa cache dir). Builds/reuses a sibling .gvlfa cache directory (self-describing, fingerprint-validated; legacy .fa.gvl caches auto-migrate). The cache is built atomically (temp + os.replace) under a best-effort lock, so concurrent builders sharing one reference are safe; the cache auto-rebuilds from its source when stale or missing. in_memory=False reads on-demand from a memory map (lower RAM) but keeps FASTA contig order, so contigs must be None or exactly the full FASTA order — reordering or subsetting contigs requires in_memory=True (otherwise raises ValueError).
• gvl.read_bedlike(path) / gvl.with_length(bed, L) — BED helpers (re-exported from seqpro).
• gvl.Ragged, gvl.RaggedAnnotatedHaps, gvl.RaggedVariants, gvl.RaggedIntervals — ragged return containers. All are backed by seqpro.rag.Ragged (_core.Ragged Rust backend); not awkward. RaggedVariants is a subclass of seqpro.rag.Ragged (class RaggedVariants(seqpro.rag.Ragged)), so isinstance(rv, Ragged) is True. Structural methods — indexing, reshape, squeeze, to_packed — are inherited from the base and preserve the RaggedVariants type (positional/structural operations return `Ra

…(truncated)