# Genoray API

> Use when writing or modifying Python code that imports `genoray` to read genotypes/dosages from VCF, PGEN, or SparseVar (`.svar`) files. Covers the public API surface, mode constants, range queries, chunking, filtering, and the SparseVar workflow. Skip for unrelated bioinformatics work.

- Skill: `d-laub/genoray-api` (Agent Skill)
- Install (CLI): `npx skillmds@latest add d-laub/genoray-api`
- Raw SKILL.md: https://api.skillmd.com/api/skills/d-laub/genoray-api/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Integrations & APIs
- Author: d-laub (https://skillmd.com/u/d-laub)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/d-laub/genoray-api

---


# genoray public API

`genoray` is a NumPy-first range-query layer over VCF/BCF (cyvcf2), PGEN
(pgenlib), and a sparse memmap format (`SparseVar` / `.svar`).

## Public surface

`import genoray` exposes exactly:

- `genoray.PGEN` — PLINK 2 PGEN reader
- `genoray.Reference` — indexed-FASTA reference genome reader
- `genoray.VCF` — VCF/BCF reader
- `genoray.Filter` — VCF filter value object bundling a cyvcf2 record predicate (`record`) with its matching `.gvi` polars expression (`expr`)
- `genoray.SparseVar` — sparse `.svar` reader/writer
- `genoray.SparseVar2` — next-gen sparse variant store (VCF/BCF → SVAR2 conversion via `from_vcf` (supports `regions=`/`samples=`/`merge_overlapping=`/`regions_overlap=`), PLINK2 PGEN → SVAR2 conversion via `from_pgen`, N single-sample VCFs/BCFs → one SVAR2 store via a native k-way merge in `from_vcf_list` (`reference`/`no_reference` supported like `from_vcf`, absent sites fill hom-ref; supports `regions=`/`merge_overlapping=`/`regions_overlap=` but **no `samples=`** — the cohort is the file set), SVAR1 (`SparseVar`) → SVAR2 native migration via `from_svar1` (reads no VCF/htslib; biallelic SVAR1 only; supports `regions=`/`samples=`/`merge_overlapping=`/`regions_overlap=` like `from_vcf`/`from_pgen`); range queries via `decode`/`region_counts`/`read_ranges`; mutational-signature support (SBS96/DBS78/ID83) via `annotate_mutations`/`mutation_matrix`/`assign_signatures`, or classify during the write with `from_vcf(signatures=True)`/`from_pgen(signatures=True)`/`from_svar1(signatures=True)`; scalar-numeric INFO/FORMAT field extraction during the write via `from_vcf(info_fields=, format_fields=)`/`from_vcf_list(info_fields=, format_fields=)` (`from_vcf_list` merges INFO first-carrier-wins, FORMAT per-sample); `from_pgen` instead stores per-sample **dosage** tracks as FORMAT fields via `dosages=Sequence[DosageField]` (from the hardcall `.pgen` itself via `source="self"`, or a separate `.pgen`) — it still has no `info_fields=`/`format_fields=` (PGEN has no VCF INFO/FORMAT); `from_svar1` carries SVAR1's existing fields through selectively via `fields=` (`None` default = all, `[]` = none, or a name subset) — read back opt-in via `fields=`/`with_fields`/`available_fields` and attached to `decode`'s result)
- `genoray.InfoField` / `genoray.FormatField` — frozen dataclasses (`name`, `dtype=None`, `default=None`) configuring a single INFO/FORMAT field for `SparseVar2.from_vcf`; a bare `str` name uses inferred defaults instead
- `genoray.DosageField` — frozen dataclass (`name="dosage"`, `source="self"|Path`, `dtype="f16"|"f32"="f32"`, `default=None`) configuring a PGEN dosage FORMAT field for `SparseVar2.from_pgen`
- `genoray.exprs` — polars filter expressions for `.gvi` indexes
- `genoray.cosmic_signatures` — fetch/cache COSMIC reference signatures
- `genoray.fit_signatures` — sparse forward-selection signature refit

Nothing else is public. Anything starting with `_` (e.g. `genoray._vcf`) is
internal — do not import it from user code.

## Where to look for details

Prefer reading these over guessing:

- `docs/source/index.md` — narrative tour with full examples (VCF, PGEN, filtering, chunking)
- `docs/source/svar.md` — SparseVar usage
- `genoray/__init__.py` — confirms the public surface
- `genoray/_vcf.py` — `VCF` class: constructor, `read`, `chunk`, mode constants near the top of the class; `get_record_info(contig=None, start=None, end=None, fields=None, info=None, lazy=False)` — non-FORMAT record-level fields (including INFO) for a range or the whole file, returns `pl.DataFrame` (or `pl.LazyFrame` when `lazy=True`)
- `genoray/_pgen.py` — `PGEN` class: constructor, `read`, `chunk`, `read_ranges`, `chunk_ranges`, mode constants near the top of the class
- `genoray/_svar.py` — `SparseVar`: `__init__`, `from_vcf`, `from_pgen`, `read_ranges`, `read_ranges_with_length(contig, starts=0, ends=POS_MAX, samples=None)` (length-guaranteed range read; returns the same type as `read_ranges` — a `Ragged` or fields-augmented record), `with_fields`, `annotate_mutations`, `mutation_matrix`, `assign_signatures`, `annotate_with_gtf(gtf, level_filter=1, write_back=True, *, strand_encoding=None, codon_null_token=None)` (GTF CDS annotation entry point, returns `pl.DataFrame` with `varID`/`gene_id`/`strand`/`codon_pos`), `cache_afs()` (computes and persists an `AF` column to the `.gvi` index; returns `None`)
- `genoray/_svar2.py` — `SparseVar2`: `__init__(path, *, fields=None)`, `with_fields(fields)` (new reader over the same store with those fields selected), `available_fields` (`dict[str, StoredField]`, set in `__init__`), `from_vcf` (VCF/BCF → SVAR2 conversion entry point, `signatures=` classifies during the write, `info_fields=`/`format_fields=` extract scalar-numeric fields during the write; supports `regions=`/`samples=`/`merge_overlapping=`/`regions_overlap=`), `from_pgen` (PLINK2 PGEN → SVAR2 conversion entry point; diploid-only, no `ploidy=`/`info_fields=`/`format_fields=`; `dosages=Sequence[DosageField]` stores per-sample dosage tracks as FORMAT fields, read from the hardcall `.pgen` itself (`source="self"`) or a separate `.pgen`; supports `regions=`/`samples=`/`merge_overlapping=`/`regions_overlap=` like `from_vcf`), `from_vcf_list` (N single-sample VCFs/BCFs → one SVAR2 store via a native k-way merge; `sources` accepts a `Sequence`/directory/manifest, resolved by module-level `_resolve_vcf_sources`; `reference`/`no_reference` supported (no_reference skips left-alignment, so cross-file joins require pre-normalized inputs); `info_fields=`/`format_fields=` supported — INFO merges first-carrier-wins, FORMAT stays per-sample; supports `regions=`/`merge_overlapping=`/`regions_overlap=` like `from_vcf`, but **no `samples=`** — the cohort is the file set), `from_svar1` (SVAR1 (`SparseVar`) → SVAR2 native migration entry point; reads no VCF/htslib, `ploidy` from SVAR1 metadata, biallelic SVAR1 only, no `info_fields=`/`format_fields=` (those are VCF-specific) — instead `fields=Sequence[str] | None` selects which SVAR1 fields carry through (`None` default = all, `[]` = none, a subset carries only those names, unknown name raises `ValueError`); `mutcat` is never selectable this way and is always dropped; supports `regions=`/`samples=`/`merge_overlapping=`/`regions_overlap=` like `from_vcf`/`from_pgen`, though regions filter per-record rather than narrowing a covering range up front); `n_samples`/`available_samples`/`contigs`/`ploidy` metadata. Read/query methods live in the mixins: `genoray/_svar2_decode.py` (`decode` — attaches one `Ragged` per selected field, `region_counts`), `genoray/_svar2_batch.py` (public `read_ranges`; internal gvl-only `_overlap_batch`/`_find_ranges`/`_gather_ranges`), and `genoray/_svar2_mutcat.py` (`annotate_mutations`, `mutation_matrix`, `assign_signatures` — COSMIC mutational-signature workflow, mirroring `SparseVar`'s but backed by a per-contig Rust sidecar instead of a `.gvi`-attached field)
- `genoray/_svar2_fields.py` — `InfoField`/`FormatField`/`DosageField` dataclasses + `FieldDtype` and the header/dtype validation used by `from_vcf(info_fields=, format_fields=)`; `_parse_cli_field_specs` (internal — parses bcftools-style `INFO/x`/`FORMAT/x`/`FMT/x` CLI field strings, used by the `genoray write vcf --fields` CLI); `StoredField` (frozen dataclass: `name`, `category`, `dtype`, `default`, `key`) is the read-side manifest entry type returned by `SparseVar2.available_fields` — not exported at top-level `genoray`, only reached via that dict
- `genoray/_cli/__main__.py` — the `genoray` CLI (`index`, `write vcf`/`write pgen`/`write svar1` (all → SVAR2), top-level `write-svar1` (legacy VCF/PGEN → SVAR1), `view` / `view svar1`, `concat`, `split`)
- `genoray/_signatures.py` — `cosmic_signatures`, `fit_signatures`
- `genoray/_reference.py` — `Reference`: `from_path`, `fetch`, `contig_array`
- `genoray/exprs.py` — the *complete* set of pre-built filter expressions (currently 7: `is_snp`, `is_indel`, `is_biallelic`, `is_symbolic`, `is_breakend`, `is_imprecise`, `ILEN`)

When a signature, kwarg, or shape is unclear, **read the docstring in the
source** rather than reasoning from first principles.

## Cross-cutting conventions

- Ranges are 0-based, half-open `[start, end)`.
- `max_mem` accepts strings like `"4g"`, `"512m"`, `"2GB"` — **except**
  `SparseVar2.from_vcf`, `SparseVar2.from_pgen`, and
  `SparseVar2.from_vcf_list`'s `max_mem`, all a **whole-process** planning
  budget, not a per-chunk cap; see their entries under "Conversion" below
  before assuming they mean the same thing as everywhere else this name
  appears (`VCF.chunk`/`chunk_ranges`, `PGEN.chunk`/`chunk_ranges`, etc.,
  where it caps one chunk directly).
- Contig names auto-normalize: `"chr1"` and `"1"` both work regardless of file convention (`ContigNormalizer`).
- Missing genotype = `-1` (int). Missing dosage = `np.nan` (float32).
- Ploidy is 2 by default; `SparseVar.from_vcf`/`from_pgen` (and `genoray write-svar1`) accept `haploid=True` / `--haploid`, which OR-collapses haplotypes into a single haploid call per sample and records `ploidy=1` in metadata (intended for unphased somatic data).
- All return arrays are NumPy; `mode` selects which arrays you get back.

## Sample accessors — canonical name + why the idioms diverge

`available_samples` (a `list[str]`) is the canonical "all samples in the
file" accessor — present on all four readers (`VCF`, `PGEN`, `SparseVar`,
`SparseVar2`).

`VCF` and `PGEN` additionally expose:
- `current_samples` — the currently-selected subset (read-only property).
- `set_samples(samples) -> Self` — a stateful call that mutates the reader
  in place to select a subset (or restore all samples with `None`), then
  returns `self`.

`SparseVar` and `SparseVar2` have **no** `current_samples`/`set_samples`.
Instead, every read method (`read_ranges`, `read_ranges_with_length`, etc.)
takes samples as a per-call `samples=` kwarg.

**Why the two idioms differ (performance):** subsetting samples on VCF/PGEN
is costly — it re-initializes the backend reader — so it's a deliberate,
stateful `set_samples()` call made once and reused across reads. On
`SparseVar`/`SparseVar2`, subsetting is ~free (it's just an index selection
over already-memory-mapped data), so it's exposed as a lightweight per-call
`samples=` kwarg instead of a persistent reader state. This is an
intentional divergence, not an inconsistency — don't "fix" one to match
the other.

## Mode constants — gotcha

Modes are **class attributes**, not top-level names:

```python
genoray.VCF.Genos8           # not genoray.Genos8
genoray.PGEN.GenosPhasingDosages
```

To discover the available modes for a class, read the class body in
`_vcf.py` / `_pgen.py` (search for `Genos` near the top).

When a mode bundles multiple arrays, the return tuple follows the order in
the constant name. `PGEN.GenosPhasingDosages` returns `(genos, phasing,
dosages)`; `VCF.Genos8Dosages` returns `(genos, dosages)`.

## VCF — quick reference

```python
vcf = genoray.VCF(
    "file.vcf.gz",
    phasing=True,             # constructor-time, not per-read
    dosage_field="DS",        # required to read dosages; FORMAT field with Number=A
    filter=genoray.Filter(
        record=lambda v: ...,                 # cyvcf2.Variant -> bool
        expr=~genoray.exprs.is_symbolic,       # matching .gvi index predicate
    ),
)

# Single range
arr = vcf.read("chr1", start=0, end=1_000_000, mode=genoray.VCF.Genos8)

# Chunked
for chunk in vcf.chunk("chr1", start=0, end=1_000_000,
                       max_mem="2g", mode=genoray.VCF.Genos8Dosages):
    ...
```

- Shape with `phasing=False`: `(samples, ploidy=2, variants)`.
- Shape with `phasing=True`: `(samples, ploidy+1=3, variants)` — the 3rd row along the ploidy axis is `0` (unphased) / `1` (phased), matching cyvcf2.
- Dosage arrays drop the ploidy axis: `(samples, variants)`, dtype `float32`.
- VCF intentionally has **no `read_ranges`** — benchmarking showed no throughput benefit.
- `read(out=...)` is **VCF-only** — pass a pre-allocated array to fill in place. PGEN random-access reads allocate fresh and have no `out=` buffer.

## PGEN — quick reference

```python
pgen = genoray.PGEN(
    "hardcalls.pgen",                # hardcalls live in the main path
    dosage_path="dosages.pgen",      # optional; defaults to the main path
    filter=genoray.exprs.is_snp & genoray.exprs.is_biallelic,
)
```

Important: when you have a dosage-only PGEN and a separate hardcalls PGEN,
**hardcalls go in the main path** and dosages go in `dosage_path`. If you
only pass one path, both hardcalls and dosages come from it (with the
hardcalls inferred from dosage threshold — see PLINK 2 docs).

A `.gvi` index file is created next to the PGEN on first construction.
Don't delete it.

```python
# Single range
genos = pgen.read("chr2", start=0, end=1000)

# Multiple ranges in one call (PGEN-only optimization)
data, offsets = pgen.read_ranges(
    "chr2",
    starts=[0, 1000, 2000],
    ends=[1000, 2000, 3000],
    mode=genoray.PGEN.GenosPhasingDosages,
)
# `data` matches the mode (tuple when mode bundles multiple arrays)
# `offsets` shape: (n_ranges + 1,). Slice range i with: arr[..., offsets[i]:offsets[i+1]]

# Chunked variants of both
for chunk in pgen.chunk("chr2", 0, 1000, max_mem="4g"): ...
for range_iter in pgen.chunk_ranges("chr2", starts, ends, max_mem="4g"):
    for chunk in range_iter: ...
```

Genotype dtype: `int32`. Dosage dtype: `float32`. Phasing is a separate
`bool` array of shape `(samples, variants)` — *not* an extra row in the
genotype array (unlike VCF with `phasing=True`).

## SparseVar (`.svar`) — quick reference

Build:

```python
# From a configured VCF reader
vcf = genoray.VCF("file.vcf.gz", dosage_field="DS")
genoray.SparseVar.from_vcf("out.svar", vcf, max_mem="4g",
                           with_dosages=True, overwrite=True)

# Or from a PGEN
genoray.SparseVar.from_pgen("out.svar", "file.pgen", max_mem="4g")

# Unphased somatic data: collapse to a single haploid call per sample (ploidy=1)
genoray.SparseVar.from_vcf("out.svar", vcf, max_mem="4g", haploid=True)
```

`SparseVar.from_vcf` / `from_pgen` inherit and apply the source's filter — filter the VCF/PGEN to filter the SVAR.

`SparseVar.from_vcf` / `from_pgen` accept `regions=`, `samples=`,
`merge_overlapping=`, `regions_overlap=` to subset by region and/or sample
during conversion (same semantics as `SparseVar.write_view`); a sample subset
drops MAC=0 variants from the output.

Read:

```python
# Plain ragged: data is just variant indices
svar = genoray.SparseVar("out.svar")
ragged = svar.read_ranges("chr1", starts=[0, 50_000], ends=[10_000, 60_000],
                          samples=["S1", "S2"])
# shape: (ranges, samples, ploidy, ~variants) — last axis is ragged

# With extra fields attached
svar = genoray.SparseVar("out.svar", fields={"dosages": np.float32})
# or, on an existing instance:
svar_with = svar.with_fields({"dosages": np.float32})
result = svar_with.read_ranges("chr1", [0], [10_000])
result.genos     # Ragged of variant indices (uint32)
result.dosages   # Ragged of dosages (float32)
```

`with_fields(False)` drops all extras and returns a plain
`Ragged[V_IDX_TYPE]` again from subsequent reads.

Each leaf value in the ragged result is a **variant index** — a row number
into `svar.index`, a polars `DataFrame` with at least `CHROM, POS, REF,
ALT (list[str]), ILEN`. To map indices back to chrom/pos/ref/alt, row-index
that DataFrame.

```python
v_idxs = ragged[0, 0, 0].to_numpy()
rows = svar.index[v_idxs.tolist()].select("CHROM", "POS", "REF", "ALT")
```

`svar.index.POS` is **1-based** (VCF convention), while query coordinates
are **0-based half-open**. Don't conflate them.

## SparseVar2 (`.svar2`) — quick reference

`SparseVar2` is the next-gen sparse variant store (VariantKey-style inline
encoding + per-variant dense/sparse cost model). Two halves: **conversion**
(`from_vcf`, below) writes a store; **range queries** (`decode` / `region_counts`
/ `read_ranges`, further below) read it back. All coordinates are 0-based
half-open `[start, end)`, as everywhere else in genoray.

### Conversion

```python
from genoray import SparseVar2

dropped = SparseVar2.from_vcf(
    "out.svar2", "file.vcf.gz", "ref.fa",   # reference: validates REF + left-aligns indels
    overwrite=True,
)

# Pre-normalized input (e.g. `bcftools norm`'d): skip REF validation/left-align
dropped = SparseVar2.from_vcf("out.svar2", "file.vcf.gz", no_reference=True)
```

Signature: `from_vcf(out, source, reference=None, *, regions=None, samples=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, ploidy=2, chunk_size=25_000, threads=None, reader_workers=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, info_fields=None, format_fields=None, check_ref="e", progress=False, log_level="info", max_mem=None) -> int`

- `source` — a bgzipped VCF (`.vcf.gz`, or the equivalent `.vcf.bgz` spelling)
  or BCF (`.bcf`). Auto-indexes (`.csi`) if no `.csi`/`.tbi` is found. For a PLINK2 PGEN source, use `from_pgen` instead
  (below).
- **`regions=`/`merge_overlapping=`/`regions_overlap=`** — restricts conversion
  to one or more indexed VCF fetch intervals. Region strings use the existing
  genoray convention (`"chrom:start-end"` is 1-based inclusive, converted to
  0-based half-open; tuple/BED/frame inputs are already 0-based half-open).
  Overlapping regions raise unless `merge_overlapping=True`. `regions_overlap`
  picks one of three modes, matching bcftools `--regions-overlap`: `"pos"`
  (default; POS inside `[start,end)`), `"record"` (POS in `[start,end+1)`, so
  an indel at the region's last base is kept), or `"variant"` (the
  anchor-trimmed variant extent overlaps the region). In `"variant"` mode a
  multiallelic record is kept whole if ANY of its alleles truly overlaps the
  region; individual non-overlapping alleles are not dropped. `"variant"`
  currently requires at most one region per contig; multiple regions per
  contig raise — use `"pos"`/`"record"`, or convert separately.
- **`samples=`** — selects and reorders VCF samples by name: preserves caller
  order, de-duplicates first occurrences, raises `ValueError` on an unknown
  name. `available_samples` and every decoded column match the caller's order
  exactly, regardless of each sample's original VCF header position.
- Exactly one of `reference` (a FASTA path, used to validate REF and left-align
  indels) or `no_reference=True` (trusts pre-normalized input, skips
  validation/left-align) is required — passing both or neither raises
  `ValueError`.
- The `reference=` FASTA may use a different contig naming scheme than the
  variant source (e.g. source `chr1`, FASTA `1`, or either side's mito contig
  spelled as `M`/`MT`/`chrM`/`chrMT`); genoray resolves the source's contig
  names against the FASTA's own naming (`chr`-prefix and mito aliases
  included) before validating REF/left-aligning. The output store keeps the
  source's contig spelling regardless of the FASTA's.
- `skip_out_of_scope=False` — when `True`, drops out-of-scope (symbolic
  `<DEL>`/`<INS>`/… and breakend) ALTs instead of erroring; the strict default
  errors on the first one. The two classes are **not** distinguishable at this
  layer — there's no separate "symbolic only" vs. "breakend only" toggle.
- Returns the number of dropped out-of-scope ALTs as an `int` (always `0`
  unless `skip_out_of_scope=True`).
- `check_ref: Literal["e", "x"] = "e"` — policy for a record whose REF
  disagrees with the reference FASTA (ignored when `no_reference=True`).
  `"e"` (default) raises and aborts the build, matching `bcftools norm
  --check-ref e`. `"x"` drops the offending record (including a REF that
  runs past the contig end) and continues, logging a per-contig count.
  Comparison is case-insensitive (soft-masked lowercase reference bases
  match). Any other value raises `ValueError` before conversion starts.
- No `dosages=` kwarg here (unlike `from_pgen`, below) — VCF dosage-like data
  goes through `format_fields=` instead (e.g. a `DS` FORMAT field). No
  `haploid=` OR-collapse, no `max_mem`-based chunking (use `chunk_size`
  instead) — those two remain `SparseVar` (SVAR 1.0)-only for now.
- `threads=None` — total thread budget (autodetected if `None`). Drives contig
  concurrency and, through the planner, the per-contig reader count.
- `reader_workers=None` — independent indexed shard readers per concurrent
  contig, the knob that sets sub-contig read parallelism. `None` derives it
  from the core budget: a quarter of usable cores is reserved for the merge
  tail, contig concurrency is chosen preferring depth (~8 readers per contig),
  and the rest goes to readers. An explicit value must be `None` or an
  integer `>= 1`; anything below 1 raises `ValueError` before conversion
  starts (checked in Python, ahead of the planner — a value that instead
  fits `max_mem` but cannot otherwise be honoured raises `InsufficientMemory`
  rather than being silently reduced). Output is byte-identical at every
  valid value. (`from_vcf_list`, the N-single-sample-VCF merge path, does not
  shard within a contig and does not accept this argument.) See "Parallel
  conversion" in `docs/source/svar.md` for scaling numbers.
- `signatures=False` — when `True`, classifies every SNP/indel into its
  SBS96/ID83 mutation-type code during the write and stores a `mutcat`
  sidecar per contig (factored into the write's dense/var_key cost model).
  Requires a reference (`reference=`); raises `ValueError` if combined with
  `no_reference=True`. There is no public read-side API for the SVAR2
  `mutcat` sidecar yet (unlike `SparseVar.annotate_mutations`/
  `mutation_matrix`, below) — this flag only controls whether the sidecar is
  written.
- `info_fields=`/`format_fields=` — `Sequence[str | InfoField]` /
  `Sequence[str | FormatField]`, `None` by default. Extracts **scalar-numeric**
  INFO/FORMAT fields into the store during the write:

  ```python
  from genoray import SparseVar2, InfoField, FormatField

  SparseVar2.from_vcf(
      "out.svar2", "file.vcf.gz", "ref.fa",
      info_fields=["AC", InfoField("AF", dtype="f16")],
      format_fields=[FormatField("DS", default=0.0)],
  )
  ```

  - **Scope: scalar-numeric only.** Header `Type` must be `Integer`, `Float`,
    or `Flag`; `Number` must be `1`, biallelic-split `A`, or `0` (`Flag`,
    INFO-only). Anything else (`Number=R`/`G`/`.`, `String`/`Character`
    fields) raises `ValueError` at config time, before conversion starts. A
    bare `str` name uses inferred defaults (`dtype=None`, no `default`); pass
    an `InfoField`/`FormatField` to override.
  - **`dtype`** (`FieldDtype = Literal["bool","i8","u8","i16","u16","i32","u32","f16","f32"]`):
    `None` (default) auto-resolves — `Integer`/`Flag` are **losslessly
    auto-narrowed** to the smallest width fitting the observed global range
    (plus a reserved missing sentinel); `Float` always resolves to `f32`
    (never silently downcast). An explicit `dtype` is validated at
    conversion time against both the header type (e.g. `Float` cannot target
    an int width) and the observed range — overflow, or `f16`'s ~65504
    range, raises `ValueError`. `f16` is the only lossy option and must be
    requested explicitly.
  - **`default`** — the value written for VCF-missing entries; otherwise a
    reserved sentinel at the extreme of the chosen width (`INT*_MIN` for
    signed widths, `u*::MAX` for unsigned widths — auto-narrowing prefers
    unsigned when the observed range is non-negative — and `NaN` for float
    widths). `Flag` fields are never missing (absent ⇒ `false`/`0`).
  - **FORMAT is genotype-aligned, not independently lossless**: a FORMAT
    value is stored only where the genotype has a call — one value per
    carrier call in var_key-routed variants, or a full dense per-sample
    column (non-carrier slots filled with `default`/sentinel) in
    dense-routed variants. Non-carrier FORMAT values (e.g. an imputed
    dosage at a ref/ref genotype) are **dropped by design** in this version;
    an independent lossless FORMAT stream is deferred to a future spec.
  - **Read path:** see "Reading INFO/FORMAT fields (SVAR2)" below —
    `SparseVar2(path, fields=…)` / `.with_fields(…)` / `.available_fields`
    opt into decoding these back out via `decode()`.
- **`max_mem: int | str | None = None`** — byte budget for the **concurrency
  planner**: how many contigs convert at once, chosen so cohort-baseline
  memory plus each concurrent contig's in-flight chunk buffers fit inside it
  (in addition to the existing core-count bound). Same string forms as the
  module-level `max_mem` convention above (`"4g"`, `"512m"`, `"2GB"`, parsed
  by `parse_memory`), and **the same whole-process meaning as
  `from_vcf_list`'s `max_mem`** (below) — `from_vcf_list` just has no fitted
  concurrency planner to spend it on (its contigs run strictly sequentially),
  so it derives its own per-chunk `chunk_size` from this budget instead.
  **`None` (the default) means a DETECTED budget — 80% of the cgroup memory
  limit (or `/proc/meminfo` total outside a cgroup) — NOT unbounded.** This
  is a deliberate default behavior change from the pre-`max_mem` planner. If
  detection itself fails (no cgroup limit and no readable `/proc/meminfo` —
  always true on macOS), genoray warns and falls back to the old
  core-bound-only planning rather than raising. Pass an explicit value to
  raise or lower the budget, or a very large value to approximate unbounded
  planning. **Practical floor:** the planner's RAM law has a fixed
  cohort-baseline term plus a per-concurrent-contig term, so any budget that
  can't cover baseline plus one concurrent contig is rejected with
  `ValueError`, even for a tiny cohort. The floor is **backend-specific**: the
  VCF law's raw LP coefficients are ~457 MB baseline plus ~111 MB per
  concurrent contig, but the per-contig bracket's `kappa` term dominates
  those two numbers completely, so the real floor for even a tiny cohort —
  evaluated at the `cc=1, w=1` point the planner actually lands on when the
  budget is tight — is **roughly 1.0 GB** (1,016 MB at S=4,000; see the table
  below), not ~600 MB — anything much below that is rejected in practice.
  `from_pgen`'s floor is roughly 2.7 GB plus ~210 MB per concurrent contig,
  putting its floor nearer ~3 GB.

  **The 2026-08-11 envelope refit roughly quadrupled `from_vcf`'s real-world
  floor versus the pre-refit law**, evaluated at a fixed `chunk_size=25_000,
  reader_workers=3, cc=1` illustrative point: the minimum `max_mem` for one
  concurrent contig went ~1.15 GB → ~1.38 GB at S=4,000, ~7.8 GB → ~26.4 GB
  at S=128,000, and ~27.9 GB → ~101.5 GB at S=500,000. Task 3 then changed
  the per-`w` charge from `kappa * (2w - 1)` to `w * (kappa + 2) + 8` per
  chunk-MB, and `from_vcf` no longer pins `reader_workers` at a fixed
  default — `None` derives it, and `plan_sharded` scans `w` downward from
  `w_max` to `1` before giving up a contig, so a tight budget now yields a
  smaller `w` instead of a refusal. The floor the shipped planner actually
  enforces at its new default is the `cc=1, w=1` point, well below either
  `w=3` figure below. The "old law" column is the 2026-08-11 refit number
  quoted just above; the "current law" column is what an explicit
  `reader_workers=3` actually costs today under the `w * (kappa + 2) + 8`
  charge — the two are not the same number, so don't read them as one:

  | cohort | old law, `w=3` (2026-08-11 refit) | current law, `w=3` | actual floor now (`w=1`) |
  |---|---|---|---|
  | S=4,000 | 1,380 MB | 1,421 MB | **1,016 MB** |
  | S=128,000 | 26,400 MB | 27,833 MB | **14,864 MB** |
  | S=500,000 | 101,480 MB | 107,069 MB | **56,408 MB** |

  An explicit `reader_workers` raises the floor above this `w=1` minimum
  (a larger `w` costs more per the `w * (kappa + 2) + 8` charge above),
  because an explicit value is honoured or refused rather than silently
  degraded to whatever `w` fits — passing `reader_workers=3` demands the
  current-law `w=3` figure above (e.g. 107,069 MB at S=500,000), or raises
  `PlanError::InsufficientMemory`, never a silent downgrade to `w=1`. The
  direction is still safe (a larger requirement means more over-allocation
  or an outright refusal to plan, never an OOM). At S=500,000, a **64 GB
  host** (`max_mem` defaults to 80% of detected RAM, i.e. `52,429 MB`) is
  still below the `w=1` floor of `56,408 MB`, so it still raises
  `PlanError::InsufficientMemory` — naming
  both remedies in its message, "raise `max_mem` or lower `chunk_size`" —
  though the margin is now narrow (52.4 vs 56.4 GB) rather than the old
  law's enormous gap. A **128 GB host** (`104,858 MB`) is no longer a near
  miss: it plans `cc=1, w=2`, which needs `81,739 MB` -- about 23 GB of
  headroom. (It stops at `w=2` because `w=3` would need `107,069 MB`, just
  over the budget.)
- **`progress=False`/`log_level="info"`** — write-time progress/logging,
  shared by `from_vcf`/`from_pgen`/`from_vcf_list`/`from_svar1`/`write_view`.
  `progress=True` renders live progress: in a terminal or Jupyter, a `rich`
  bar (one row per in-flight contig); elsewhere, compact heartbeat lines
  throttled to roughly one per 5s per contig (`"chr1 42% (12,345/29,000)
  ..."`). Regardless of `progress`, a one-line `"[svar2] chrom done: N kept,
  M excluded (Ts)"` summary prints per contig once it finishes, unless
  `log_level="off"`. `log_level` is the minimum severity for structured
  write-time log lines — `"off"` (disables everything, including the
  per-contig summaries and progress rendering — a pure no-op, zero
  overhead), `"warning"`, `"info"` (default; also includes thread-budget
  selection, per-contig start/finish, and contig-name resolution against the
  reference when it differs from the source's own spelling), or `"debug"`
  (additionally surfaces per-record detail: a record excluded for a
  REF/FASTA mismatch, and each indel that gets left-aligned). The
  `GENORAY_LOG` environment variable overrides the `log_level` argument when
  set to one of the same four values (e.g. `GENORAY_LOG=debug`), without
  touching call sites.
  Structured log lines render their fields inline as ` key=value` pairs
  after the message (e.g. `pipeline config concurrent_chroms=8
  reader_workers=4 exact_counts=true planned_units=32`), matching what
  `GENORAY_LOG`'s stderr layer emits.

### Conversion from PGEN

```python
from genoray import SparseVar2

dropped = SparseVar2.from_pgen(
    "out.svar2", "file.pgen", "ref.fa",   # reference: validates REF + left-aligns indels
    overwrite=True,
)
```

Signature: `from_pgen(out, source, reference=None, *, regions=None, samples=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, chunk_size=None, max_mem=None, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, dosages=None, check_ref="e", progress=False, log_level="info") -> int`

- `source` — a `.pgen` file. Variant metadata is read from the sibling
  `.pvar`/`.pvar.zst`, sample names from the sibling `.psam`.
  `reference`/`no_reference`, `skip_out_of_scope`, `overwrite`,
  `long_allele_capacity`, `signatures`, and `check_ref` all mean the same as
  `from_vcf` (above), and return the same `int` (dropped out-of-scope ALTs).
- Unlike `from_vcf`, PGEN sub-contig **sharding is disabled** (single reader
  per contig) and `threads` never changes a single output byte. Reason
  (measured, chr21c ~1M variants x 3202 samples): single-reader conversion is
  already fast (~33s) and bound by the shared executor/writer + reference I/O,
  not by pgenlib decode -- so sharding cannot beat that floor and measured
  as *slower* (44.9s at threads=24 vs 32.6s serial). Bumping `pgenlib` to a
  GIL-releasing build (>=0.94.x, which parallelizes decode via `prange`) does
  not help either: the conversion is flat at ~33s across `OMP_NUM_THREADS`
  1..32, so decode parallelism buys nothing. The sharding machinery exists and
  is byte-identical (validated to 1M variants) for re-enablement only if a
  future change shifts the bottleneck onto decode.
- **Diploid only** — no `ploidy=` kwarg (`from_vcf`'s default `ploidy=2` is
  implicit and fixed here).
- `chunk_size=None` — unlike `from_vcf`'s fixed `25_000` default, `None` here
  derives a variant-count budget from sample count (a packed dense chunk costs
  `chunk_size * n_samples * 2 / 8` bytes), so a fixed constant that's fine at
  200 samples doesn't blow memory at 500k. Pass an explicit `int` to override.
  Warns if the derived value falls below 256 variants — see
  `from_vcf_list`'s `chunk_size` entry below for the details. `dosages`
  counts as `n_format_fields` here.
- **`max_mem: int | str | None = None`** — byte budget for the **concurrency
  planner**: how many contigs convert at once, chosen so cohort-baseline
  memory plus each concurrent contig's in-flight chunk buffers fit inside it
  (in addition to the existing core-count bound, also capped at 8 concurrent
  contigs regardless of budget). Same string forms as the module-level
  `max_mem` convention above (`"4g"`, `"512m"`, `"2GB"`, parsed by
  `parse_memory`), and **the same whole-process meaning as `from_vcf`'s
  `max_mem`** (above) — both pipelines have a fitted concurrency planner and
  spend the budget on concurrency the same way, just with
  separately-fitted RAM-law coefficients (a PGEN chunk decodes both
  haplotypes at once, so its per-variant cost is higher). `from_vcf_list`'s
  `max_mem` means the same whole-process budget too, but that path has no
  concurrency planner to spend it on (its contigs run strictly
  sequentially), so it derives its own per-chunk `chunk_size` from the
  budget instead — see its entry below. **`None` (the default) means a
  DETECTED budget — 80% of the cgroup memory limit (or `/proc/meminfo` total
  outside a cgroup) — NOT unbounded.** This is a deliberate default behavior
  change from the pre-`max_mem` planner. If detection itself fails (no
  cgroup limit and no readable `/proc/meminfo` — always true on macOS),
  genoray warns and falls back to the old core-bound-only planning rather
  than raising. Pass an explicit value to raise or lower the budget, or a
  very large value to approximate unbounded planning. **Practical floor:**
  the planner's RAM law has a fixed cohort-baseline term of roughly 2.7 GB
  (PGEN's own fitted coefficients, higher than `from_vcf`'s ~457 MB), so any
  budget that can't cover baseline plus one concurrent contig's chunk
  buffers is rejected with `ValueError`, even for a tiny cohort. That
  baseline scales with cohort size (`~0.0158 MB/sample`), so this isn't just
  a small-cohort concern: at ~500k samples it alone predicts ~10.6 GB, so a
  *detected* budget on a smaller host will reject the conversion — pass an
  explicit `max_mem` sized to the host in that case.
- **`regions=`/`merge_overlapping=`/`regions_overlap=`** — same convention,
  semantics, and three overlap modes (`"pos"`/`"record"`/`"variant"`) as
  `from_vcf`, restricting conversion to one or more `.pvar` variant-index
  ranges. As with `from_vcf`, `"variant"` mode keeps a multiallelic record
  whole if ANY of its alleles truly overlaps the region.
- **`samples=`** — selects and reorders `.psam` samples by name (same
  convention as `from_vcf`): preserves caller order, de-duplicates first
  occurrences, raises `ValueError` on an unknown name. `available_samples`
  and every decoded column match the caller's order exactly, regardless of
  each sample's original `.psam` position.
- **No `info_fields=`/`format_fields=`** — PGEN carries no FORMAT, and `.pvar`
  INFO extraction is not implemented.
- **`dosages=Sequence[DosageField]`** — stores per-sample dosage tracks as
  FORMAT fields. Each `DosageField(name="dosage", source="self"|Path,
  dtype="f16"|"f32"="f32", default=None)`: `source="self"` reads dosages from
  the hardcall `.pgen` (`source` above) itself; a `Path` reads from a
  *separate* `.pgen` (e.g. a VAF/CCF file kept apart because pgenlib derives
  hardcalls from dosage when both live in one file) — it must share the
  hardcall `.psam`'s samples and align 1:1 on the hardcall `.pvar`'s variants.
  Stored genotype-aligned like any FORMAT field: under var_key routing a
  non-carrier's dosage is dropped (harmless for VAF/CCF-style fields, ~0 for
  non-carriers). Read back the same way as other FORMAT fields — see "Reading
  INFO/FORMAT fields (SVAR2)" below.

  ```python
  from genoray import SparseVar2, DosageField

  SparseVar2.from_pgen("out.svar2", "cohort.pgen", "ref.fa",
                       dosages=[DosageField(name="DS", source="self")])
  # separate dosage file (e.g. VAF stored as dosage):
  SparseVar2.from_pgen("out.svar2", "hardcalls.pgen", "ref.fa",
                       dosages=[DosageField(name="VAF", source="vaf.pgen")])
  ```
- Unphased heterozygotes resolve haplotypes in the allele-code order
  `pgenlib` returns — the same caveat `from_vcf` carries for unphased `GT`.
- **`progress=False`/`log_level="info"`** — same as `from_vcf` (above).

### Conversion from a list of single-sample VCFs

```python
from genoray import SparseVar2

# Explicit list
dropped = SparseVar2.from_vcf_list("out.svar2", ["s1.vcf.gz", "s2.bcf"], "ref.fa")

# A directory of single-sample files
# (non-recursive: all *.vcf.gz/*.vcf.bgz, then all *.bcf)
dropped = SparseVar2.from_vcf_list("out.svar2", "vcfs/", "ref.fa")

# A manifest file (one path per line; blank/`#`-comment lines skipped;
# relative entries resolved against the manifest's directory)
dropped = SparseVar2.from_vcf_list("out.svar2", "manifest.txt", "ref.fa")
```

Signature: `from_vcf_list(out, sources, reference=None, *, regions=None, merge_overlapping=False, regions_overlap="pos", no_reference=False, skip_out_of_scope=False, ploidy=2, chunk_size=None, max_mem=None, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, info_fields=None, format_fields=None, check_ref="e", progress=False, log_level="info") -> int`

Builds **one** SVAR2 store from **N single-sample** VCFs/BCFs with different
site lists, via a native k-way merge — no `bcftools merge`, no intermediate
multi-sample VCF.

- **`regions=`/`merge_overlapping=`/`regions_overlap=`** — same convention,
  semantics, and three overlap modes (`"pos"`/`"record"`/`"variant"`) as
  `from_vcf`, applied identically to every input file in the merge. As with
  `from_vcf`, `"variant"` mode keeps a multiallelic record whole if ANY of its
  alleles truly overlaps the region.
- **No `samples=` parameter** — unlike `from_vcf`/`from_pgen`/`from_svar1`,
  `from_vcf_list` has no cohort to subset by name: each input file is already
  single-sample, and the cohort is exactly the file set passed via `sources`.

- **Each input file must be single-sample** — exactly one sample column;
  `ValueError` if any file has zero or more than one. That sample's VCF
  header name becomes its sample name in the store; duplicate sample names
  across input files raise `ValueError`.
- `sources` — one of three forms, resolved by module-level
  `_resolve_vcf_sources`:
  - a `Sequence[str | Path]` — explicit files, in the given order.
  - a single directory `Path` — every bgzipped VCF (`*.vcf.gz`/`*.vcf.bgz`)
    then every `*.bcf` directly inside it (non-recursive), each group
    `natsort`-ordered.
  - a single file `Path` — `.vcf.gz`/`.vcf.bgz`/`.bcf` is taken as one file; anything
    else is a manifest (one path per line, blank/`#`-comment lines skipped,
    relative entries resolved against the manifest's parent directory).
  - Resolving to zero files raises `ValueError`.
- **Absent site → hom-ref `0`.** A site called in file A but not present at
  all in file B fills `0` (hom-ref) for B's sample at that site.
- **A within-file `./.` is not observable after the merge.** SVAR2's sparse
  layout stores only ALT-carrying e

…(truncated)