# 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. - Skill: `mcvickerlab/genvarloader` (Agent Skill) - Install (CLI): `npx skillmds@latest add mcvickerlab/genvarloader` - Raw SKILL.md: https://api.skillmd.com/api/skills/mcvickerlab/genvarloader/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: AI & ML - Author: mcvickerlab (https://skillmd.com/u/mcvickerlab) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/mcvickerlab/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 ```python 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 (``, ``, …) 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. ```bash # 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`: ```python 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 `/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=)` 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 ```python 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 `/annot_intervals//`. 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 ```python 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 `/intervals//`. - `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 `/annot_intervals//`. - `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.replace`d into place), so a reader can never see a half-written track. Source: `python/genvarloader/_dataset/_write.py`. ## `gvl.concat` — merge on-disk datasets ```python 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 ```python 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 `/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/"` / `"FORMAT/"`. 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[]` 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: | `kind` | Returns | Use 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 | | `None` | No sequences | Tracks-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: | `fmt` | Container types returned | Default? | |------------|----------------------------------------------------------------------------|----------| | `"ragged"` | `Ragged` / `RaggedVariants` / `RaggedAnnotatedHaps` (all `seqpro._core.Ragged`-backed) | Yes | | `"flat"` | Pure-numpy `FlatRagged` / `FlatVariants` / `FlatAnnotatedHaps` / `FlatIntervals` | No | 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`). ```python 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:** ```python 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). | Strategy | Behavior | |--------------------------------|-----------------------------------------------------------------------| | `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. ```python 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`. ```python 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. ```python 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)