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:

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

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

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.

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

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

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

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

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:

  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

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.

  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

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)