skillZs
★ LIVE SKILL TAGS ★
>>> LIVE SKILLS INDEX <<<
* OPEN SOURCE *
NO LOGIN, NO TRACKING
※ REAL INSTALL DATA ※
← back to all skills
mcvickerlab/genvarloader129 installs

genvarloader

Use when writing or reading GenVarLoader (gvl) datasets — preparing VCF/PGEN/SVAR variant sources with bcftools/plink2, calling gvl.write, configuring gvl.Dataset for haplotype/reference/annotated/variants output modes, attaching BigWig or Table tracks, setting up spliced haplotypes from a GTF, choosing track insertion-fill strategies for indels, or filtering variants by allele frequency.

How do I install this agent skill?

npx skills add https://github.com/mcvickerlab/genvarloader --skill genvarloader
view source ↗

Is this agent skill safe to install?

  • Gen Agent Trust Hubpass

    This skill provides technical documentation and API usage guidelines for GenVarLoader, a genomic data processing library. It includes instructions for variant preprocessing, dataset creation, and integration with common bioinformatics tools like bcftools and plink2. No malicious patterns or security vulnerabilities were detected.

  • Socketpass

    No alerts

  • Snykpass

    Risk: LOW · No issues

What does this agent skill do?

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 the var-key window for each (region, sample, ploid) that holds a variant, under <dataset>/genotypes/svar2_ranges/ (sized to the dataset's selected samples, not the full .svar2 cohort) — as a sparse region-CSR table, 28 bytes per non-empty window, small at cohort scale (~504 MB for the All of Us chr22 grid; 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 bounds the genoray chunk stream gvl.write reads while producing the genotype range cache — ranges are produced in per-sample chunks sized to fit the budget rather than a whole contig at once. It does not bound the per-contig entry accumulator, which peaks at roughly 60 bytes per entry on the largest contig (see "Common gotchas" below), nor the permanent genotypes/svar2_ranges/ cache's on-disk size, which scales with variants observed inside regions (not regions x samples x ploidy) and is governed by disk space (see 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. Planning is streaming — memory scales with regions plus samples, not their product — so the resident cost is the merged offsets array (one int64 per merged (region, sample[, ploid]) slot) plus the inputs' own offsets arrays (together roughly twice that figure, since the inputs partition the merged grid), not the merge plan. 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:

kindReturnsUse when
"reference"Reference sequence (S1)Baseline / no personalization
"haplotypes"Personalized haplotypes with indels (S1)Standard variant-aware modeling
"annotated"AnnotatedHaps (haps + var_idxs + ref_coords)Need to map back to variants/ref coords
"variants"RaggedVariants (variants only, no seq)Variant-centric tasks
"variant-windows"FlatVariantWindows (per-allele window/allele token buffers; flat mode only)Tokenized model input around each variant
NoneNo sequencesTracks-only datasets

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

fmtContainer types returnedDefault?
"ragged"Ragged / RaggedVariants / RaggedAnnotatedHaps (all seqpro._core.Ragged-backed)Yes
"flat"Pure-numpy FlatRagged / FlatVariants / FlatAnnotatedHaps / FlatIntervalsNo

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).

StrategyBehavior
gvl.Repeat5p() (default)Repeat the value at variant POS across the insertion.
gvl.Repeat5pNormalized()Repeat track[POS] / (insertion_len + 1). Preserves sum.
gvl.Constant(value=nan)Constant value (default NaN) across the insertion.
gvl.FlankSample(flank_width=5)Resample with replacement from a 2·flank_width+1 window around POS.
gvl.Interpolate(order=1)Polynomial interp (order 1/2/3) between flanking reference values.

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 RaggedVariants). A string key (rv["start"]) returns a bare Ragged field, not a RaggedVariants. reshape takes the new shape either as unpacked ints — e.g. rv.reshape(1, 2, None) — or as a single tuple rv.reshape((1, 2, None)); the base Ragged signature accepts both. squeeze(axis=None) is a real axis-squeeze (base semantics) — it squeezes any size-1 axis, not a fixed "drop axis 0". An int index collapses the leading axis (numpy-consistent); slice/array indexing preserves it. Indexing an allele field (rv.alt[b][h], rv.ref[b][h]) yields the alleles of that haplotype per variant, aligned one-for-one with rv.start[b][h] / rv.ilen[b][h], so zip(rv.start[b][h], rv.alt[b][h]) pairs each position with its own allele. This requires seqpro >= 0.22; on older seqpro the same expression returned every allele of the group concatenated into a single bytes, which silently misaligned against the scalar fields for any group containing a multi-byte (indel) allele. Named properties (.alt, .ref, .start, .ilen, .end) are the primary access point; extra fields (e.g. AF, custom FORMAT fields) are also accessible via rv["field"] or rv.field (via __getattr__). RaggedVariants itself does not define __eq__ (wrapper-level == is Python object-identity, not element-wise); to compare contents, compare individual fields — e.g. rv["alt"] == other_alt or rv.start == other_start — which use seqpro.rag.Ragged's ufunc-based (element-wise) comparison. Domain methods retained on RaggedVariants: .rc_(), .pad(), .to_nested_tensor_batch(); derived read-only properties: .ilen, .end; fields: .alt, .ref, .start, .dosage.
  • gvl.FlatRagged — flat analog of Ragged: .data (flat numpy array), .offsets (int64), .shape. Methods: .to_ragged(), .to_fixed(length), .to_padded(pad_value), .reshape(shape), .squeeze(axis). Source: python/genvarloader/_flat.py.
  • gvl.FlatIntervals — flat-buffer interval container returned by with_tracks(kind="intervals") + with_output_format("flat"). Fields .starts/.ends/.values are FlatRagged; .to_ragged() → RaggedIntervals; .reshape(...), .squeeze(...), .shape. Source: python/genvarloader/_ragged.py.
  • gvl.FlatAnnotatedHaps — flat analog of RaggedAnnotatedHaps: fields .haps, .var_idxs, .ref_coords (each a FlatRagged). Methods: .to_ragged(), .to_fixed(length), .to_padded(), .reshape(shape), .squeeze(axis). Source: python/genvarloader/_flat.py.
  • gvl.FlatVariants — flat analog of RaggedVariants: .fields dict mapping field names to FlatRagged (scalar fields: start/ilen/dosage/info) or FlatAlleles (alt/ref). .shape delegates to fields["start"].shape. .flank_tokens: optional ride-along FlatRagged of shape (b, p, ~v, 2L) (set when with_settings(flank_length=L, ...) is configured; None otherwise). Methods: .to_ragged(), .reshape(shape), .squeeze(axis). Source: python/genvarloader/_dataset/_flat_variants.py.
  • gvl.FlatAlleles — two-level flat bytestring for allele fields: .byte_data (uint8), .seq_offsets (per-variant byte offsets, int64), .var_offsets (per-(batch×ploidy)-row variant offsets, int64), .shape. Methods: .to_ragged(), .reshape(shape), .squeeze(axis). Source: python/genvarloader/_dataset/_flat_variants.py.
  • gvl.FlatVariantWindows — returned by with_seqs("variant-windows", VarWindowOpt(...)) in flat mode. .fields: dict of scalar FlatRagged (start/ilen/dosage/info; raw byte alleles are dropped). Per-allele token buffers — exactly one of .ref_window (flanked ref window, "window" mode) or .ref (bare ref allele tokens, "allele" mode) is set; same for .alt_window / .alt. Each non-None buffer is a two-level token buffer (internal _FlatWindow, not the public FlatRagged) of shape (b, p, ~v, ~len) with its own .to_ragged(). The container's .shape delegates to fields["start"].shape. Methods: .to_ragged() (returns dict of ragged parts), .reshape(shape), .squeeze(axis). Source: python/genvarloader/_dataset/_flat_variants.py.
  • gvl.VarWindowOpt — frozen config dataclass for with_seqs("variant-windows", ...). Fields: flank_length (int), token_alphabet (str | bytes | sp.NucleotideAlphabet, e.g. sp.alphabets.DNA; normalized to bytes on construction), unknown_token (int), ref ∈ {"window","allele"}, alt ∈ {"window","allele"}. ref and alt are chosen independently. "window" = flanked + tokenized reference read (ref) or flank·alt·flank assembly (alt); "allele" = bare tokenized allele with no flanks. Source: python/genvarloader/_dataset/_flat_variants.py.
  • gvl.DummyVariant — frozen dataclass used with with_settings(dummy_variant=...). Fields and defaults: start: int = -1, ilen: int = 0, dosage: float = 0.0, ref: bytes = b"N", alt: bytes = b"N", info: dict = {}. Unspecified info keys default to 0 for integer columns and NaN for float columns. Source: python/genvarloader/_dataset/_flat_variants.py.
  • gvl.migrate(path) — upgrade a pre-2.0 (array-of-structs) dataset to format 2.0 (struct-of-arrays) in place. Streaming, idempotent, crash-safe; converts intervals/<track>/ and annot_intervals/<track>/ interval storage and bumps metadata.json. A no-op (with leftover-AoS cleanup) on an already-2.0 dataset. Source: python/genvarloader/_dataset/_migrate.py. (Distinct from gvl.migrate_svar_link, which upgrades legacy SVAR symlink layouts.)
  • gvl.to_nested_tensor(ragged) — convert to a PyTorch nested tensor (requires torch).
  • gvl.get_dummy_dataset() — small in-memory dataset for examples/tests.
  • gvl.RefDataset — reference-only dataset (no genotypes).
  • gvl.Table — core interval-track source backed by a Rust COITrees overlap engine. Zero-based half-open coordinates; positive-width intervals assumed. Usable directly as a tracks= or annot_tracks= source in gvl.write and gvl.update. No extra install required — ships as part of the core package. CI-covered via a brute-force numpy oracle + property tests.
  • gvl.data_registry.fetch(name) — download public test/demo datasets.

Full list lives in python/genvarloader/__init__.py __all__.

On-disk layout (quick reference)

ds.gvl/
├── metadata.json              # version, samples, contigs, ploidy, max_jitter, svar_link?
├── input_regions.arrow        # BED + region index map
├── genotypes/                 # variant_idxs.npy, dosages.npy, variants.arrow
│                              # (absent when sourced from .svar; see svar_link)
│                              # svar2_ranges/ present iff sourced from .svar2 (see svar2_link)
├── intervals/<track>/         # per-sample track data (BigWigs / Table)
└── annot_intervals/<track>/   # sample-independent annotation track data

In format 2.0, each intervals/<track>/ (and annot_intervals/<track>/) directory stores its intervals as struct-of-arrays — three contiguous files starts.npy (int32), ends.npy (int32), values.npy (float32), sharing one offsets.npy (int64) — replacing the format 1.x single intervals.npy record array. This lets the contiguous memmaps cross the Python→Rust boundary zero-copy. Upgrade a 1.x dataset with gvl.migrate(path) (see the format version gate above).

See docs/source/format.md for the full schema, versioning, and SVAR-link details.

Where to look next

For…Read…
End-to-end RNA-seq exampledocs/source/geuvadis.ipynb
Splicing tutorialdocs/source/splicing.ipynb
Deep-learning eval pipelinedocs/source/basenji2_eval.ipynb
BED / BigWig / bcftools recipesdocs/source/write.md
Dataset shapes, ragged/var/fixeddocs/source/dataset.md
On-disk format + SVAR resolutiondocs/source/format.md
FAQ (with_* design, typing)docs/source/faq.md
Auto-generated referencedocs/source/api.md → https://genvarloader.readthedocs.io
gvl.write / gvl.update internalspython/genvarloader/_dataset/_write.py
gvl.concat internals + preconditionspython/genvarloader/_dataset/_concat.py, _concat_validate.py
Sharded-build / merge workflowdocs/source/write.md ("Merging datasets"), docs/source/faq.md
Track re-alignment internalspython/genvarloader/_dataset/_tracks.py, _reconstruct.py
Insertion fill internalspython/genvarloader/_dataset/_insertion_fill.py
SVAR back-reference / migrationpython/genvarloader/_dataset/_svar_link.py
.svar2 back-reference / read-bound wiringpython/genvarloader/_dataset/_svar2_link.py, _svar2_haps.py
.svar2 sparse range-cache layout (write/read/concat)python/genvarloader/_dataset/_svar2_ranges.py
Format 1.x → 2.0 migration internalspython/genvarloader/_dataset/_migrate.py
Flat-buffer ragged containerspython/genvarloader/_flat.py
Flat variants + alleles typespython/genvarloader/_dataset/_flat_variants.py
Flank fetch + tokenization + windowspython/genvarloader/_dataset/_flat_flanks.py

Common gotchas

  • Pre-2.0 datasets must be migrated once before opening. Dataset.open rejects any dataset written by genvarloader < 2.0 (or unversioned) with a ValueError pointing at gvl.migrate(path). Run it once (in place, idempotent, crash-safe). A dataset written by a newer major raises a different ValueError asking you to upgrade genvarloader. Note gvl.migrate (format upgrade) is not the same as gvl.migrate_svar_link (SVAR symlink-layout upgrade).
  • gvl.update does not hot-reload open datasets. A Dataset instance that was opened before gvl.update ran will not see the new track; reopen the dataset to pick it up. The update itself is safe to run while readers are active — each track is published atomically so a reader never sees a half-written track.
  • Dataset.write_annot_tracks has been removed. Use gvl.update(dataset, annot_tracks={"name": source}) instead, or pass annot_tracks= to gvl.write at creation time.
  • gvl.Table is a core public API. No extra install required. It uses a Rust COITrees overlap engine and is CI-covered. Import it as gvl.Table (re-exported from the top-level package).
  • Symbolic / breakend variants are rejected, not skipped. Remove them before gvl.write — e.g. bcftools view -e 'ALT~"<" || ALT~"\["' (drop SVs and breakends), or construct the genoray reader with filter=genoray.exprs.is_biallelic & ~genoray.exprs.is_symbolic & ~genoray.exprs.is_breakend. SVAR inputs must be built from an already-filtered source, since gvl validates but cannot re-filter a materialized .svar. For .svar2 the same variants are rejected upstream at .svar2 conversion time (genoray), not at gvl.write time — the store format cannot represent them at all.
  • .svar2 has a Phase-1 unsupported-combination matrix. min_af/max_af, annotated haplotypes, VarWindowOpt(ref="allele"), splicing or var_filter="exonic" with non-haplotype/track output, reverse-complemented or fixed-length haplotype-realigned tracks, variants/variant-windows output with jitter (max_jitter>0 at write or jitter>0 at read), and extend_to_length=False at write time all raise NotImplementedError on a .svar2-backed dataset instead of silently mis-computing. Plain haplotypes support splicing, exonic filtering, and negative-strand reverse-complementation. with_seqs("variant-windows"), unphased_union, var_fields-selected store INFO/FORMAT fields (on both "variants" and "variant-windows"), and multi-contig FlankSample track fills (issue #267) are also supported for .svar2. See ".svar2 — the read-bound sparse variant format" above.
  • .svar2 variants/variant-windows ALT bytes differ from .svar for pure deletions. .svar keeps the VCF anchor base (b"G" for GTA>G); .svar2 decodes the atomized empty ALT (b""). Reconstructed haplotypes are byte-identical either way; ref_window is also byte-identical — only raw ALT/alt_window bytes differ for pure-deletion records.
  • Opening a genotypes-only dataset without a reference= defaults to the "variants" view (RaggedVariants), not "haplotypes". Only "variants" is available without a reference; with_seqs("haplotypes" | "annotated" | "reference") raises ValueError if no reference was provided.
  • with_insertion_fill raises unless the dataset has both haplotypes AND tracks active.
  • min_af / max_af are supported for .svar only. A non-SVAR-backed dataset (plain BCF/PGEN) raises RuntimeError; a .svar2-backed dataset raises NotImplementedError specifically.
  • with_len(L) requires L + 2·jitter ≤ min(region_length) + 2·max_jitter — set max_jitter accordingly at write time.
  • Tracks must have unique .name; the on-disk layout is intervals/<name>/.
  • BED strand of . is treated as +. Reverse-complement happens automatically when rc_neg=True (default) and strand == "-".
  • Splicing is a read-time setting on a flat BED of exons — do not pre-concatenate exons before gvl.write.
  • extend_to_length=False at write time will produce haplotypes shorter than the BED region when deletions are present; downstream code must tolerate < region length. Not an option for a .svar2 variant source — gvl.write(..., variants=<.svar2>, extend_to_length=False) raises NotImplementedError.
  • Missing a dosage field on a RaggedVariants output you expected? Check var_fields — dosage must be requested explicitly even if dosages.npy exists on disk.
  • Allele fields index per-variant, and need seqpro>=0.22. rv.alt[b][h] gives one allele per variant, matching rv.start[b][h]. On seqpro <0.22 it returned the group's alleles concatenated as one bytes — a mismatch that is invisible for SNP-only data (every allele is 1 byte, so the counts coincide) and only appears once an indel is in the window. gvl requires seqpro>=0.22; if you see a concatenated blob, check the installed seqpro version first.
  • FlatRagged / FlatVariants offsets are int64. PyTorch nested tensors require int32 offsets — cast with .astype(np.int32) or tensor.to(torch.int32) before passing to torch.nested.narrow.
  • kind="intervals" cannot be re-aligned: combining it with a variant-aware seq mode (haplotypes/annotated/variants/variant-windows) raises unless with_settings(realign_tracks=False). (Breaking change: haplotypes+intervals previously returned un-realigned intervals silently under the default.)
  • with_insertion_fill raises when realign_tracks=False.
  • with_seqs("variant-windows") requires a VarWindowOpt second argument, with_output_format("flat"), and a reference genome. Querying in "ragged" mode raises. It does not inherit flank settings from with_settings.
  • Flat ride-along flank tokens need tracks disabled. FlatVariants.flank_tokens (ride-along on "variants" output) is only produced on the pure flat variants channel — with active tracks the variants output falls back to ragged RaggedVariants (no flank_tokens). Call with_tracks(False) for the ride-along. "variant-windows" may be combined with tracks via with_settings(realign_tracks=False) (tracks are returned reference-coordinate alongside the windows output).
  • Flank tokens (FlatVariants.flank_tokens) and variant windows (FlatVariantWindows) are reference-oriented — they are NOT reverse-complemented even when rc_neg=True. Only the alt/ref allele fields of FlatVariants / RaggedVariants are RC'd (not their scalar fields, not flank tokens, not windows).
  • Token dtype is uint8 when max token id ≤ 255, else int32; offsets are int64.
  • VarWindowOpt(ref="allele") (bare allele mode) requires REF alleles on disk. flank_length ride-along on FlatVariants requires token_alphabet and unknown_token to be set together in the same or a prior with_settings call.
  • dummy_variant padding applies to both "variants" and "variant-windows" outputs. Setting dummy_variant=<DummyVariant> and then indexing with any other kind ("haplotypes", "annotated", "reference", or no seqs) raises ValueError. For token fields (flank_tokens, ref_window/alt_window, bare ref/alt), the dummy fill is all-unknown_token — the DummyVariant.ref/.alt bytes only set the dummy allele's byte-length, not the token value. dummy_variant=False with an unsupported output kind is silently ignored.
  • A non-b"N" DummyVariant.alt (or .ref) is reverse-complemented on negative-strand regions, exactly like a real variant allele. The default b"N" is rc-invariant; use it if you want a strand-neutral sentinel.
  • unphased_union=True + with_seqs("haplotypes") / with_seqs("annotated") raises — unphased_union only applies to "variants" / "variant-windows" output.
  • The .svar2 range-cache write-time accumulator is not bounded by max_mem. gvl.write with a .svar2 source builds the permanent genotypes/svar2_ranges/ cache as a sparse region-CSR table — 28 bytes per non-empty (region, sample, ploid) window, ~504 MB for the All of Us chr22 grid, small at cohort scale (see format.md). The remaining memory surprise is at write time: max_mem bounds the genoray chunk stream but not the per-contig entry accumulator, which peaks at roughly 60 bytes per entry on the largest contig — the per-chunk (region int32, cell int32, 24-byte entry) blocks plus the merged (cell, entry) output live at the same time (~1.1 GB at chr22, low single-digit GB on the largest contigs).
  • Dataset.samples is sorted lexicographically, not numerically. Cohorts with integer-like IDs of mixed digit counts come out in an order that differs from the numeric sort a phenotype table typically has ("1000" < "999" as strings). A positional join between the two is silently wrong; join on the sample name.
  • gvl.concat requires one shared variant source and at least two inputs. Mismatched variant sources (checked by fingerprint, not just backend type), axis="regions" with differing sample sets/order, or axis="samples" with differing regions or overlapping samples all raise ValueError before any bytes move. It moves roughly the full merged-dataset size at sequential-IO speed — not a routine step.

Maintaining this skill

Whenever a PR changes the public API (anything in python/genvarloader/__init__.py __all__, or the docstring/signature of gvl.write, Dataset.open, or any Dataset.with_* method), the author must also update this SKILL.md. New public symbols, removed symbols, renamed args, changed defaults, and new output modes are all in scope. CLAUDE.md enforces this as part of the contribution checklist.

Add the canonical catalog link to the repository README so users can inspect current installs and available audits. The publishing guide covers the complete discovery path.

<a href="https://skillzs.dev/skills/mcvickerlab/genvarloader/genvarloader">View genvarloader on skillZs</a>