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.Tuning — frozen dataclass of six explicit scheduling knobs (concurrent_chroms, reader_workers, overshard, dense_cap, merge_threads, sample_interval, all int | None) passed as tuning= to every SparseVar2.from_* writer; replaces the removed legacy environment-variable configuration (see "Tuning" under "SparseVar2 — quick reference", below, and "migrating from environment variables" in docs/source/svar.md) — no environment variable configures genoray
• genoray.exprs — polars filter expressions for .gvi indexes
• genoray.cosmic_signatures — fetch/cache COSMIC reference signatures
• genoray.fit_signatures — mutation-catalogue signature refit; strategy= selects Forward() (default, sparse forward selection) or Spa() (SigProfilerAssignment-faithful backward elimination)
• genoray.Forward / genoray.Spa — fit_signatures strategy dataclasses (frozen)
• genoray.Strategy / genoray.Criterion / genoray.Metric / genoray.ActivityScale — the associated type aliases; genoray.SPA_CONNECTED_GROUPS holds Spa's four default co-occurring-signature groups

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). INFO fields come back as top-level columns, not a nested INFO struct, named with the header's spelling of each ID. info= names are matched case-insensitively against the INFO IDs the header declares; None means all of them, [] means none, and a name the header does not declare raises ValueError (it used to be silently dropped along with every other INFO field, #139). An INFO ID that collides with a non-INFO column name (e.g. INFO/QUAL) also raises rather than flattening ambiguously.
• genoray/_pgen.py — PGEN class: constructor, read, chunk, read_ranges, chunk_ranges, mode constants near the top of the class; var_records(contig, starts=0, ends=POS_MAX) — the identity of the variants a read over the same ranges returns, as a pl.DataFrame of range (u32, which query range the row belongs to), CHROM, POS (1-based), REF, ALT, aligned 1:1 with the variant axis of read/read_ranges and honoring the reader's filter. ALT is always list[str] here regardless of how the .gvi stored it, and CHROM is the contig name as the file spells it (which may differ from the queried contig). Use this rather than var_idxs + the private _index when an artifact has to outlive the PGEN it came from
• 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/ — cosmic_signatures, fit_signatures, Forward, Spa, Strategy, Metric, ActivityScale, SPA_CONNECTED_GROUPS. fit_signatures(..., strategy=Forward()|Spa()) selects the refit algorithm; Forward is genoray's own greedy forward selection (the default), Spa reimplements SigProfilerAssignment's cosmic_fit (backward elimination from the saturated signature set, then add-remove refinement)
• 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.

Tuning (genoray.Tuning)

Every SparseVar2.from_* write method (from_vcf, from_pgen, from_vcf_list, from_svar1) accepts tuning: Tuning | None = None. Tuning is the explicit scheduling-knob object that replaced the old environment-variable configuration — no environment variable configures genoray; every knob is now a constructor argument (or a CLI flag, see "CLI" below).

from genoray import SparseVar2, Tuning

SparseVar2.from_vcf(
    "out.svar2", "file.vcf.gz", "ref.fa",
    tuning=Tuning(reader_workers=4, dense_cap=64),
)

Six fields, all int | None, keyword-only, on a frozen dataclass:

field	meaning	minimum
concurrent_chroms	contigs converted concurrently	1
reader_workers	independent indexed shard readers per concurrent contig	1
overshard	work units per reader (consulted only when a contig has no exact record count)	1
dense_cap	depth of the dense-chunk channel between reader and executor	1
merge_threads	gather threads for the per-contig var_key merge tail	1
sample_interval	monitor sampling cadence in seconds; 0 disables it	0

None (the default for every field) means "let the planner choose" — not "off" — the same value the planner would pick with no Tuning at all. A field you do set is honoured or refused, never silently shrunk: an explicit reader_workers that cannot fit max_mem raises InsufficientMemory rather than being downgraded to something that fits. This is the whole point of the object — the old environment variables could be silently ignored by a backend that never read them; Tuning cannot be. tuning= must be a Tuning instance or None — passing e.g. a dict raises TypeError naming the parameter, since the field names could otherwise look like plausible kwargs.

Not every backend can use every knob. Setting one it can't use raises ValueError naming the offending field(s) rather than being silently ignored:

field	from_vcf	from_pgen	from_vcf_list	from_svar1
concurrent_chroms	yes	yes	no	yes
reader_workers	yes	no	no	no
overshard	yes	no	no	no
dense_cap	yes	yes	yes	yes
merge_threads	yes	yes	yes	yes
sample_interval	yes	yes	yes	yes

Why:

• reader_workers/overshard are sharded-VCF only. from_pgen pins reader_workers=1 and never shards within a contig — pgenlib holds the GIL through decode, so sub-contig sharding there is pure overhead. Neither from_vcf_list (already one file descriptor per input file per contig) nor from_svar1 (no VCF/htslib reader at all) shards within a contig either.
• concurrent_chroms is unavailable on from_vcf_list because that pipeline walks contigs sequentially by design.

A field you do set is reported in the pipeline config log line tagged explicit; a field left None is reported tagged planner — so the log line itself tells you which knobs came from you versus the planner.

Tuning is a pure scheduling knob everywhere it applies: output is byte-identical at every value (gated by a store-hash oracle).

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, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, info_fields=None, format_fields=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None, 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.

• tuning.reader_workers (set via tuning=Tuning(reader_workers=N), see "Tuning" above) — independent indexed shard readers per concurrent contig, the knob that sets sub-contig read parallelism. Leaving it 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 at Tuning construction time, before conversion starts — 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 rejects reader_workers/overshard — see the applicability table above.) 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"/log_filter=None — 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. Accepted, case-insensitive: "off", "critical", "error", "warning", "info" (default), "debug", a logging module integer constant (logging.DEBUG and friends — an in-between int rounds UP to the more severe named level, matching logging's own gate), or that integer's decimal spelling ("10"), which is how the CLI reaches the int path. Canonical ranks: off=0, error=1, warning=2, info=3, debug=4. "critical" is accepted but is an alias for "error" — there is no distinct CRITICAL rank. "warn" is rejected (ValueError), on purpose: logging.warn() was removed in Python 3.13, so genoray does not accept that spelling either — use "warning". "off" disables everything, including the per-contig summaries and progress rendering (a pure no-op, zero overhead); "info" additionally includes thread-budget selection, per-contig start/finish, and contig-name resolution against the reference when it differs from the source's own spelling; "debug" additionally surfaces per-record detail (a record excluded for a REF/FASTA mismatch, each indel that gets left-aligned). See parse_log_level for the exact mapping.

  log_filter is an optional tracing-style EnvFilter directive string applied to genoray's own stderr diagnostic output, independent of log_level (which gates only the structured progress/summary events above). None (default) leaves stderr silent; an unparseable directive silences stderr rather than raising. Use it for module/target-scoped tracing that log_level can't express, e.g. log_filter="genoray::monitor=trace".

  No environment variable configures genoray. log_level/log_filter (and every Tuning field, above) replace all eight former environment variables that used to configure this pipeline — see "migrating from environment variables" in docs/source/svar.md for the full mapping, including the migration for the log-tracing variable's genoray::monitor=trace directive to log_filter="genoray::monitor=trace".

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

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, tuning=None, log_level="info", log_filter=None) -> 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. The per-contig chunk charge is 2 chunk-buffers, not the 10 the VCF path is charged: the 8-chunk reorder-backlog ceiling is a from_vcf mechanism that the PGEN pipeline (one reader per contig) never enforces, so budgets between those two brackets that used to raise ValueError now plan.

• 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"/log_filter=None — same as from_vcf (above). tuning= accepts concurrent_chroms, dense_cap, merge_threads, sample_interval; setting reader_workers/overshard raises (from_pgen pins a single reader per contig — see "Tuning" 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, tuning=None, log_level="info", log_filter=None) -> 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 entries, so a missing hap and a hom-ref hap both produce zero entries and cannot be told apart via decode or region_counts. The -1 missing sentinel is a dense genoray.VCF/genoray.PGEN convention and is not part of SVAR2's decode. (The distinction is real inside the merge, but it is discarded when genotypes are packed into the sparse carrier bit-grid — this matches from_vcf, so the two paths stay in parity.)

• The merge is join-on-atom — a variant is one shared row across files iff its normalized (pos, ref, alt) atom matches exactly, not merely its position.

• Each input file's records must already be position-sorted per contig (same assumption from_vcf makes for its single input) — an unsorted file raises ValueError naming the offending file and positions rather than silently corrupting the k-way merge.

• Every input file must use the same contig naming scheme (all chr1-style or all 1-style, not a mix) — the merge matches contigs by an exact per-file string, so a cohort mixing schemes raises ValueError up front (naming the conflicting files/spellings) instead of silently producing a store where half the cohort's samples decode as all-zeros on the "wrong-spelled" contigs.

• Opens all N input files concurrently (one file descriptor per file per contig) — at large N (roughly N > (soft RLIMIT_NOFILE - 64) / 2, often around N ≈ 480 at a default 1024 soft limit) this raises ValueError with the ulimit -n remedy instead of htslib's more confusing "is there a .tbi or .csi file?" error for some arbitrary file near the ceiling. There is no batched/hierarchical merge to fall back on for very large cohorts (future work) — raise the open-file limit instead.

• no_reference=True is supported, same as from_vcf/from_pgen: skips REF validation and left-alignment, reconstructing each atom's REF from the record's own REF bytes. The reference/no_reference exactly-one-of check and the signatures+no_reference incompatibility are otherwise identical to from_vcf.

  • Caveat specific to this entry point: because the merge is a per-contig k-way join keyed on each atom's normalized (pos, ref, alt), skipping left-alignment means a site shared across files only joins into one output row if every input already represents it identically — same anchor base, same padding (e.g. all files came from the same caller, or were all already run through bcftools norm against the same reference). Two files encoding the same indel with different normalization will not join under no_reference: they silently become two separate variants in the output store instead of one shared row. This is not a failure mode that raises — verify upstream normalization is consistent before relying on no_reference with from_vcf_list.

• info_fields=/format_fields= — same declaration API as from_vcf (resolved against the FIRST file in sources's header). Merge semantics differ from a single-file conversion because there are now N source columns per site:

  • INFO fields merge first-carrier-wins. When a site is shared across files, the stored INFO value comes from the lowest-numbered (earliest in sources order) file that carries the atom — not the last file, and not an aggregate (e.g. max/sum) of the carriers' values.
  • FORMAT fields stay per-sample, exactly as in from_vcf: each sample gets its own file's value; a sample that doesn't carry the atom at all gets the field's default (reserved sentinel/NaN, or an explicit default=).

• chunk_size=None — unlike from_vcf's fixed 25_000 default, None here derives a budget-based chunk size from the cohort size (_auto_chunk_size), so one packed dense chunk stays within a max_mem-derived per-chunk target (see max_mem below), up to a ~256 MiB ceiling; this is the same default from_pgen/from_svar1 already use. The budget accounts for both the packed genotype grid AND any staged format_fields (n_format_fields * n_samples * 4 bytes/variant — this term can dominate the grid by 32 * F / ploidy, e.g. 112x at F=7, ploidy=2), so requesting FORMAT fields on a large cohort shrinks the auto chunk size accordingly. Scope: it bounds only the dense-chunk term, which is a small fraction of peak RAM at typical cohort sizes — a large-cohort guardrail, not a fix for overall RAM scaling in the number of inputs. Pass an int to override with a fixed count. If the derived chunk_size falls below 256 variants, _auto_chunk_size warns (naming the budget, the derived chunk_size, and n_samples/ploidy/n_format_fields) and keeps the small value rather than raising — e.g. at ploidy=2, n_format_fields=7 this fires above n_samples≈37,118 against the default ~256 MiB budget. Raise max_mem or request fewer format_fields to clear it. Shared by all three converters that derive chunk_size via this helper (from_pgen, from_svar1, from_vcf_list).

• max_mem: int | str | None = None — byte budget the whole process may use (same string forms as the module-level max_mem convention, e.g. "4g", parsed by parse_memory) — the same meaning as from_vcf's max_mem (above), not a per-chunk cap. This budget only sizes the dense-chunk term — it does NOT cap the reader's per-input-file overhead (from_vcf_list opens all N files at once per contig), which scales with the number of input files and is not bounded by max_mem at all; don't read "whole process" as a hard total-RSS guarantee.

  BREAKING CHANGE: earlier versions treated this max_mem as a direct cap on the bytes of one in-flight dense chunk. It is now a whole-process budget. A call like max_mem="512MiB" that used to mean "let one dense chunk use up to 512 MiB" now means "the dense-chunk term should use at most 512 MiB", which derives a MUCH smaller chunk_size. Re-tune any carried-over value.

  from_vcf_list has no fitted concurrency planner the way the sharded from_vcf reader does — its contigs convert strictly sequentially — so instead of planning concurrency it derives a per-chunk byte target: chunk_target = min(_DENSE_CHUNK_TARGET_BYTES, max_mem // (concurrent_jobs * in_flight_chunks_per_job)). concurrent_jobs and in_flight_chunks_per_job are not hardcoded Python literals — they are read live from the Rust extension module (_core.VCF_LIST_CONCURRENT_CHROMS, currently 1, and _core.VCF_LIST_DENSE_CHANNEL_CAP + 2, currently 8: the reader/executor channel's capacity plus one chunk each of the single reader and single executor thread may hold outside it), so this can't silently drift out of sync with the orchestrator's actual architecture the way a duplicated constant could — not a fitted memory law either way (none exists for this pipeline, unlike the sharded path's RAM law). _DENSE_CHUNK_TARGET_BYTES stays a ceiling, so a large budget can't grow chunks past today's size — only a tight budget shrinks them below it. On top of that, the resulting target is still a worst-case ceiling: the estimate assumes every variant routes dense, so cohorts whose variants route sparse (e.g. private somatic calls) use considerably less. Ignored when chunk_size is passed explicitly — including skipping memory-budget DETECTION entirely (no attempt, no warning) when chunk_size is explicit. None (the default) means a DETECTED budget, same detection as from_vcf; if detection fails, this degrades to the historical fixed ~256 MiB dense-chunk target (with a warning) rather than raising.

• ploidy, skip_out_of_scope, threads, overwrite, long_allele_capacity, signatures, check_ref all mean the same as from_vcf, and the return value is the same int (dropped out-of-scope ALTs). check_ref is applied per input file during the merge (ignored when no_reference=True): under "x", a bad record is excluded from its own file only (not the whole merged site), and the per-contig log reports the total excluded across every input file.

• progress=False/log_level="info"/log_filter=None — same as from_vcf (above). tuning= accepts only dense_cap, merge_threads, sample_interval — setting concurrent_chroms, reader_workers, or overshard raises (this pipeline walks contigs sequentially and never shards within a contig — see "Tuning" above).

Conversion from SVAR1

from genoray import SparseVar2

dropped = SparseVar2.from_svar1(
    "out.svar2", "old.svar", "ref.fa",   # reference: validates REF + left-aligns indels
    overwrite=True,
)

Signature: from_svar1(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, threads=None, overwrite=False, long_allele_capacity=8*1024*1024, signatures=False, fields=None, check_ref="e", progress=False, tuning=None, log_level="info", log_filter=None) -> int

Migrates an existing SVAR 1.0 (SparseVar) store to SVAR2 natively — reads no VCF and no htslib; SVAR1 is already sparse, so this reconstructs variant records from SVAR1's arrays and reuses the same conversion spine as from_vcf.

• source — a SparseVar store directory (SVAR1). 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).
• ploidy is read from SVAR1's metadata — no ploidy= kwarg.
• chunk_size=None derives a variant-count budget from cohort size the same way as from_pgen/from_vcf_list (_auto_chunk_size) and warns under the same below-256-variant condition — see from_vcf_list's chunk_size entry above for the details. The budget counts the FORMAT fields fields= carries (all of them by default), so a wide SVAR1 store derives a smaller chunk, the same as from_pgen does for dosages=.
• Biallelic SVAR1 only — raises ValueError if the source store has multiallelic variants (SVAR1's geno==1 model); re-create the SVAR1 store biallelically first.
• regions=/merge_overlapping=/regions_overlap= — same convention, semantics, and three overlap modes ("pos"/"record"/"variant") as from_vcf/from_pgen. "variant" mode keeps a record whole if ANY of its alleles truly overlaps the region (though SVAR1 is itself biallelic-only, so this only ever judges a single ALT). Unlike from_pgen, SVAR1 has no on-disk covering-range index to narrow against up front — a selected contig's local variants are still scanned in full; the per-record filter is what actually restricts the output, so this costs a full-contig scan rather than a range-restricted one.
• samples= — selects and reorders SVAR1 samples by name (same convention as from_vcf/from_pgen): 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 SVAR1 position.
• fields=Sequence[str] | None — selects which SVAR1 FORMAT fields (e.g. dosages) carry through, keyed by their SVAR1 name: None (default) carries all of them (the prior, lossless-by-default behavior), [] carries none, and a subset of names carries only those — an unknown name raises ValueError listing the available fields. mutcat is never selectable this way and is always dropped — pass signatures=True to recompute signatures from the reference instead of carrying SVAR1's.
• Field parity caveat: because SVAR1 never stored non-carrier FORMAT values, field output is byte-identical to from_vcf only for var_key (carrier-only) routed variants — for dense-routed variants, non-carrier cells are filled with the field's default/missing sentinel rather than the source VCF's true value. Genotype streams themselves (not fields) are byte-identical to from_vcf under matching normalization regardless of routing.
• No info_fields=/format_fields= kwargs — those names are VCF-specific; use fields= (above) instead to select which SVAR1 fields carry.
• progress=False/log_level="info"/log_filter=None — same as from_vcf (above). tuning= accepts concurrent_chroms, dense_cap, merge_threads, sample_interval; setting reader_workers/overshard raises (from_svar1 reads no VCF/htslib and never shards within a contig — see "Tuning" above).

Range queries

Open a finished store, then query per contig. Construction reads meta.json and opens one native reader per contig, exposing .available_samples (list; the canonical sample-name accessor shared with VCF/PGEN/SparseVar), .n_samples, .contigs, .ploidy, .format_version.

from genoray import SparseVar2

sv = SparseVar2("out.svar2")
regions = [(0, 40), (1_000, 2_000)]   # 0-based half-open [start, end)

# Analysis path — decode to a seqpro Ragged record (one call per contig)
rag = sv.decode("chr1", regions)      # fields pos (i32), ilen (i32), allele (ALT bytes)
                                      # + one per selected field (see "Reading
                                      # INFO/FORMAT fields" below); shape
                                      # (R, S, P, None); pure-DEL ALT is empty

# Decode-free per-(region, sample, ploid) variant count — replaces SVAR 1.0's var_ranges
counts = sv.region_counts("chr1", regions)   # np.ndarray, shape (R, S, P)

• decode(contig, regions) returns a seqpro.rag.Ragged whose layout is byte-identical to gvl's RaggedVariants (pos/ilen numeric, allele opaque-string ALT, one shared variant-axis offsets object). ALT is empty for a pure deletion (the reference base is not re-emitted). Requires seqpro.
• region_counts(contig, regions) is the decode-free count (offset diffs + dense-mask popcount) — the simplified stand-in for SparseVar.var_ranges (SVAR2 has no unified variant table, so variant indices no longer exist).
• Queries are per contig — cross-contig batching is the caller's job. Regions are an iterable of (start, end) pairs.
• The contig argument to decode/region_counts/read_ranges accepts alternate naming schemes — chr-prefixed vs unprefixed (chr1 ↔ 1) and the mitochondrial aliases {M, MT, chrM, chrMT} — resolved via ContigNormalizer to the store's own spelling. An unresolvable contig raises ValueError.

The user-facing SVAR2 query API is decode / region_counts / read_ranges (above). read_ranges(contig, starts, ends, samples=None) is a fused search+gather; starts/ends are parallel 1D arrays (mirrors SparseVar.read_ranges), samples selects/reorders a subset by name. It returns the raw two-channel BatchResult → numpy dict, a TypedDict with a fixed field set: vk_pos/vk_key/vk_off, dense_pos/dense_key/ dense_range/dense_present/dense_present_off, lut_bytes/lut_off, and scalars n_regions/n_samples/ploidy.

SparseVar2 also has _overlap_batch/_find_ranges/_gather_ranges (underscore-prefixed) — an internal, gvl-only numpy-dict wire contract for the search/gather split used by a write-time overlap cache. They are not part of the public API, are not covered by semver, and may change or disappear without notice; don't call them from user code.

Reading INFO/FORMAT fields (SVAR2)

Fields written by from_vcf(info_fields=…, format_fields=…) (above) are read back by opting in — they are not decoded by default (each one costs extra I/O).

sv = SparseVar2("out.svar2")
sv.available_fields                    # {"AF": StoredField(...), "DS": StoredField(...)}

sv = sv.with_fields(["AF", "DS"])       # or SparseVar2("out.svar2", fields=["AF", "DS"])
rag = sv.decode("chr1", [(0, 10_000)])
rag["AF"]                              # Ragged, sharing offsets with pos/ilen/allele

• SparseVar2(path, *, fields=None) / .with_fields(fields) — fields is a Sequence[str] of canonical keys (see available_fields below). with_fields returns a new SparseVar2 over the same store; it does not mutate the original in place. fields=None (the constructor default) selects nothing — fields are opt-in.
• available_fields -> dict[str, StoredField] — every field declared in the store's meta.json, keyed canonically: the bare field name when it is unique across INFO and FORMAT, else bcftools-style INFO/DP / FORMAT/DP when a name is used by both categories. StoredField (defined in genoray._svar2_fields, not exported at top-level genoray) is a frozen dataclass: name, category ("info"/"format"), dtype (np.dtype), default (float | None), key.
• decode(contig, regions) attaches one Ragged per selected field to the returned record Ragged, alongside pos/ilen/allele — every one sharing a single variant-axis offsets object, shape (R, S, P, None). Access a field's data via rag["KEY"] (Ragged.__getitem__), not rag.fields["KEY"] — Ragged.fields is just the list[str] of field names on the record.
• Dtype is preserved as stored. SVAR2 losslessly auto-narrows integer fields at write time, so e.g. an AC field may come back as int8; nothing is widened on read.
• Missing values are the field's default if one was set at write time, else a reserved sentinel (NaN for floats, iinfo.min/iinfo.max for ints) — returned as-is, never translated.
• FORMAT fields are genotype-aligned (see from_vcf's format_fields= above) — decode() only ever emits carrier records, so the "non-carrier values aren't stored" caveat from the write path is invisible on this read surface.

Mutational signatures (SBS96 / DBS78 / ID83)

Same COSMIC workflow as SparseVar (see "Mutation catalogues" above), backed by a Rust per-contig sidecar instead of a .gvi-attached field. Annotation is required before mutation_matrix — either post-hoc, or by passing signatures=True to from_vcf (above):

sv = SparseVar2("out.svar2")
ref = genoray.Reference.from_path("hg38.fa")

sv.annotate_mutations(ref)                   # post-hoc; writes the mutcat sidecar
sv.annotate_mutations(ref, contigs=["chr1"])  # restrict to a subset of contigs

df = sv.mutation_matrix("SBS96")                     # count="allele" (default)
df = sv.mutation_matrix("DBS78", count="sample")
act = sv.assign_signatures("SBS96")                  # mutation_matrix + fit_signatures

• annotate_mutations(reference, *, gtf=None, contigs=None) -> None — reference is a genoray.Reference or a FASTA path; contigs=None (default) annotates every contig. contigs= accepts alternate naming — chr-prefixed vs unprefixed and the mitochondrial aliases {M, MT, chrM, chrMT} — resolved via ContigNormalizer to the store's own spelling; raises ValueError if every requested contig fails to resolve. Unlike SparseVar.annotate_mutations, there is no write_back= toggle — SVAR2 always persists the sidecar to disk. gtf= optionally supplies a GTF/GFF gene model path; when given, each SNV is additionally classified by transcriptional-strand class (from feature == "gene" footprints) and persisted to a strand.bin sidecar, which unlocks the "SBS192"/"SBS384" catalogs below. Each sidecar and meta.json is written atomically (same-directory temp + rename), so a crash never leaves a truncated file; re-running is unconditional and therefore also the recovery path for a missing or suspect sidecar. Re-running without a gtf= removes a previously written strand.bin (stale strand classes are never kept alongside a fresh gtf-less annotation).
• mutation_matrix(kind, *, count="allele"|"sample", contigs=None) -> pl.DataFrame — a MutationType column (fixed COSMIC codebook order) plus one column per sample. kind ∈ {"SBS96", "DBS78", "ID83", "SBS192", "SBS384"}. count="allele" counts every non-ref allele copy; count="sample" counts each category at most once per sample, OR-combined across contigs. contigs= restricts the sum to a subset (alternate naming accepted, duplicates collapsed; a name absent from the store raises ValueError) and requires only those contigs to be annotated; None (default) uses every contig. Raises ValueError if a contig in scope is not annotated (no on-disk sidecar) — annotate first, either via annotate_mutations or from_vcf(..., signatures=True). "SBS192"/"SBS384" additionally require strand annotation (annotate_mutations(..., gtf=...)) and raise ValueError if the scoped contigs lack it. A sidecar that is truncated or stale (length does not match the store's variant counts) raises OSError telling you to re-run annotate_mutations. assign_signatures does not accept "SBS192"/"SBS384" — see below.
• assign_signatures(kind, *, reference=None, count="allele", contigs=None, strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame — mutation_matrix(kind, count=..., contigs=...) then genoray.fit_signatures(...). reference accepts a pl.DataFrame, a TSV path, or None (defaults to genoray.cosmic_signatures(kind)). contigs= (contig subset, alternate naming accepted) is passed through to mutation_matrix. strategy=/max_delta=/min_activity=/criterion= are forwarded to fit_signatures exactly like the top-level function — pass strategy=Spa() for SPA's algorithm, or leave it None to use the max_delta/min_activity/criterion forward-selection shorthand. Calls fit_signatures once, so passing strategy= together with any of the three shorthand kwargs raises ValueError, exactly like fit_signatures.
• Same classification rules as v1 (shared Rust classifier): DBS78 arises only from isolated adjacent same-haplotype SNV pairs — runs of ≥3 adjacent SNVs stay as individual SBS96 entries, native MNVs > 2bp are atomized into SBS96, and each isolated doublet is counted once (not once per constituent SNV).
• No public read-side access to the raw per-genotype mutcat codes for SVAR2 (unlike v1's fields=["mutcat"]) — only the aggregated mutation_matrix output is exposed.

Strand-resolved catalogs (SBS192 / SBS384)

SparseVar2 also supports the transcriptional-strand-bias catalogs, which require a gene model (GTF) at annotation time:

sv2.annotate_mutations(reference, gtf="gencode.v45.annotation.gtf.gz")
sbs384 = sv2.mutation_matrix("SBS384")   # 384 rows: [T, U, N, B] x 96
sbs192 = sv2.mutation_matrix("SBS192")   # 192 rows: the {T, U} sub-view = SBS384[:192]

• SBS384 = 96 trinucleotide channels x 4 strand categories, SigProfiler order [T, U, N, B]: Transcribed, Untranscribed, Nontranscribed (intergenic), Bidirectional (position covered by genes on both strands).
• SBS192 is the {T, U} sub-view (SBS384[:192]).
• Strand rule (pyrimidine-folded): a genic SNV is Untranscribed iff the pyrimidine of its ref/alt pair sits on the gene's coding strand, else Transcribed. Gene footprints come from feature == "gene" rows (full gene body); pre-filter the GTF to restrict biotypes.
• Without a gtf=, mutation_matrix("SBS192"/"SBS384") raises. Write-time from_vcf(..., signatures=True) stays strand-free; obtain strand catalogs via a post-hoc annotate_mutations(reference, gtf=...).
• assign_signatures("SBS192"/"SBS384") raises NotImplementedError: COSMIC publishes no strand-resolved reference set. Use mutation_matrix for strand-bias analysis.

Mutation clusters (SigProfilerClusters)

annotate_clusters ports SigProfilerClusters' per-mutation subclassification (doublet / MBS / omikli / kataegis / other) onto a finished store as the cluster_class FORMAT field. It is post-hoc only (no write-time flag) and needs no reference — the classifier uses no sequence context.

sv = SparseVar2("out.svar2")

sv.annotate_clusters(imd_cutoff=1000.0)                  # no-VAF mode
sv.annotate_clusters(imd_cutoff=1000.0, vaf_field="VAF") # VAF mode
sv.annotate_clusters(imd_cutoff={"S1": 500.0, "S2": 1_000.0}, vaf_field="VAF")
sv.annotate_clusters(imd_cutoff=1000.0, vaf_field="VAF", contigs=["chr1"])

# The original instance is stale for the new field (`available_fields` is read
# at construction); re-open the store and opt in before decoding.
sv = SparseVar2("out.svar2").with_fields(["cluster_class"])
rag = sv.decode("chr1", [(0, 1_000_000)])
labels = rag["cluster_class"]   # 1:1 with rag["pos"], the per-call flat order

• Signature: annotate_clusters(*, imd_cutoff: float | Mapping[str, float], vaf_field: str | None = None, vaf_cut: float = 0.1, contigs: Sequence[str] | None = None) -> None.
• imd_cutoff — inter-mutational-distance cutoff(s), in bases: a mutation is clustered iff its minimum neighbour distance is <= its sample's cutoff. Either one float applied to every sample, or a mapping from sample name to cutoff that must name every available_samples entry and no others (ValueError otherwise); cutoffs must be > 0. Cutoffs come from the caller, not the store: upstream derives them per sample by simulation and persists them in imds.pickle — feed the same values and every event upstream classifies gets the same label (SigProfilerClusters v1.2.2).
• vaf_field — optional FORMAT field (canonical key or bare name, as in available_fields) holding per-call VAF/CCF. When given, VAF consistency participates in the decision tree and failed events are greedily re-split. Must resolve to a FORMAT field with a 2- or 4-byte float dtype (f16/f32); anything else raises ValueError. Missing values are treated as upstream's -1.5 sentinel. Omit (default None) for the no-VAF path.
• vaf_cut — maximum |delta VAF| for adjacent mutations to stay in one event (upstream's default 0.1; 0.25 is its CCF setting). Must be > 0.
• contigs — restrict annotation to a subset (alternate naming resolved as in annotate_mutations: chr-prefixed vs unprefixed and the mito aliases); a name absent from the store raises ValueError. Every out-of-scope contig is 255-filled so the field stays decodable there. None (default) annotates every contig. Re-running is unconditional.
• Stamps meta.json with cluster_version (currently 1), cluster_contigs, cluster_cutoff, cluster_vaf_field, cluster_vaf_cut, and adds the cluster_class FORMAT entry (u8).

Label codebook — one label per SNP carrier call (a 1/2 genotype gets a label per recorded alt; indels are never subclassified):

code	label	upstream
0	nonclustered	not in the clustered file
1	doublet	ClassIA
2	mbs	ClassIB
3	omikli	ClassIC
4	kataegis	ClassII
5	other	ClassIII
255	not_annotated	out-of-scope contig / non-SNP call

Read back through the existing field machinery — with_fields(["cluster_class"]) then decode(contig, ranges). The field rides the same variant axis as pos/ilen/allele, so rag["cluster_class"] is in the store's per-call flat order, 1:1 with rag["pos"]. On disk, var_key_* streams store one label per call; dense_* streams store one per (dense_row, sample) in dense_row * n_samples + sample order, 255 for non-carriers; indel streams are 255 throughout. Labels are total over in-scope SNP calls: where upstream drops events (clustered singletons, no-VAF failures) genoray writes other, and non-clustered SNP calls are nonclustered.

Merge and split by contig

SVAR2 contigs are fully independent on disk, so recombining or subsetting whole contigs is a cheap metadata-rewrite + file-copy operation — unlike write_view (see the CLI section below), none of these methods re-run conversion or the var_key/dense cost model.

from genoray import SparseVar2

sv = SparseVar2("out.svar2")
sv.subset_contigs("chr1.svar2", "chr1")                # single contig
sv.subset_contigs("subset.svar2", ["chr1", "chr2"])    # multiple, source order preserved
paths = sv.split_by_contig("by_contig/")               # one store per contig, out_dir/{contig}.svar2

SparseVar2.concat("merged.svar2", ["chr1.svar2", "chr2.svar2"])  # disjoint-contig merge

• subset_contigs(output, contigs, *, mode="copy", overwrite=False) -> None — write a new store containing only contigs (a single contig name or a sequence of names). contigs accepts alternate naming — chr-prefixed vs unprefixed and the mitochondrial aliases {M, MT, chrM, chrMT} — resolved via ContigNormalizer to the store's own spelling. Pure metadata rewrite + file copy of the kept contig directories, preserving the source store's contig order. Raises ValueError if any name is unresolvable against self.contigs, or if output resolves to this store's own path (in-place subsetting is rejected, mirroring write_view's in-place guard). Raises FileExistsError if output exists and overwrite=False.
• split_by_contig(out_dir, *, mode="copy", overwrite=False) -> list[Path] — explode into one single-contig store per contig at out_dir/{contig}.svar2; returns the output paths in self.contigs order. Implemented as one subset_contigs call per contig.
• SparseVar2.concat(output, sources, *, mode="copy", overwrite=False) -> None (classmethod) — concatenate stores with disjoint contig sets into one. sources is a sequence of paths (or SparseVar2 instances); all sources must agree on samples, ploidy, format_version, and fields — disagreement on any of those, or a contig name appearing in more than one source, raises ValueError. The merged contig list is natsorted, independent of the order sources were passed in.
• mode (all three methods) is the shared Mode literal — "copy"|"hardlink"|"symlink"|"move" — controlling how each contig directory is transplanted into the output store.

Errors

genoray raises standard Python builtins, by category:

• ValueError — bad input content: contig/sample not found, REF disagrees with the reference FASTA, or a symbolic/breakend ALT with skip_out_of_scope=False.
• FileNotFoundError — a required input file is missing.
• OSError — a corrupt/truncated store sidecar or an underlying disk I/O failure.
• RuntimeError — an internal genoray bug (a worker thread panicked); please report it.

CLI

genoray write has three subcommands — write vcf, write pgen, write svar1 — and all three target SVAR2. There is no bare auto-detecting genoray write SOURCE OUT anymore; you must name the source kind. The previous SVAR 1.0 (SparseVar) write path lives at the top-level genoray write-svar1 command (hyphenated, not a write subcommand) — it takes a VCF or PGEN source, same as before. genoray view still defaults to SVAR2, with the previous SVAR 1.0 behavior under view svar1. genoray concat/genoray split are SVAR2-only (no SVAR1 equivalent).

genoray write vcf / genoray write pgen / genoray write svar1

# write vcf — VCF/BCF (or a directory/manifest of single-sample VCFs/BCFs) → SVAR2
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa
genoray write vcf file.vcf.gz out.svar2 --no-reference
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa --skip-symbolics-and-breakends --threads 4
genoray write vcf file.vcf.gz out.svar2 --reference ref.fa --fields INFO/AF --fields FORMAT/DP
genoray write vcf vcf_dir/ out.svar2 --no-reference --regions chr1:1-1000   # vcf-list form

# write pgen — PLINK2 PGEN → SVAR2 (no --ploidy; PGEN is diploid-only)
genoray write pgen file.pgen out.svar2 --reference ref.fa
genoray write pgen file.pgen out.svar2 --no-reference --regions chr1:1-1000 --samples A,B
genoray write pgen file.pgen out.svar2 --reference ref.fa --dosages DS=self
genoray write pgen hardcalls.pgen out.svar2 --reference ref.fa --dosages VAF=vaf.pgen

# write svar1 — SVAR1 (SparseVar) → SVAR2
genoray write svar1 store.svar out.svar2 --no-reference --samples A,B
genoray write svar1 store.svar out.svar2 --no-reference --fields dosages
genoray write svar1 store.svar out.svar2 --no-reference --empty-fields

# write-svar1 (legacy, top-level) — VCF or PGEN → SVAR 1.0, dosages, --haploid, --max-mem
genoray write-svar1 file.vcf.gz out.svar --max-mem 4g --haploid

All three write subcommands share --regions/-r, --regions-file/-R, --samples/-s, --samples-file/-S, --merge-overlapping, --regions-overlap (pos/record/variant), --reference XOR --no-reference (required), --chunk-size, --threads/-@, --overwrite, --long-allele-capacity (advanced), a single --skip-symbolics-and-breakends flag (maps to skip_out_of_scope=; the SVAR2 core can't expand either symbolic ALTs (<DEL>, <INS>, …) or breakends into nucleotides, so they're dropped together and print a Dropped {n} out-of-scope (symbolic/breakend) ALT alleles. line when set), --check-ref {e,x} (default e, ignored with --no-reference; e aborts on the first REF/FASTA disagreement, x drops the offending record and continues — mirrors bcftools norm --check-ref), --progress/--no-progress + --log-level LEVEL (map to progress=/log_level=; see "Conversion" above for behavior — default --no-progress --log-level info. LEVEL is one of off, critical, error, warning, info, debug, case-insensitive — every name parse_log_level accepts. warn is rejected: logging.warn() was removed in Python 3.13, so it is a deprecated spelling, not a convention. The parameter is a validated str rather than a fixed choice list precisely so it cannot drift from parse_log_level, which is the single source of truth; the integer logging constants that function also accepts are reachable from this flag too, spelled out: --log-level 10 == --log-level debug), and --log-filter DIRECTIVE (maps to log_filter=; a tracing-style EnvFilter string, e.g. --log-filter "genoray::monitor=trace", independent of --log-level; unset/None by default). write vcf's vcf-list form forwards these to from_vcf_list; its single-file form forwards them to from_vcf.

Every write subcommand also exposes Tuning's always-applicable knobs as --dense-cap N, --merge-threads N, and --sample-interval N (map to tuning=Tuning(dense_cap=, merge_threads=, sample_interval=), see "Tuning" above). --concurrent-chroms N is additionally available on write vcf, write pgen, and write svar1 — not on the vcf-list form (from_vcf_list walks contigs sequentially by design; see "Tuning" above). As with the Python tuning= kwarg, an explicit flag value is honoured or refused — never silently reduced — and passing a flag a backend can't use raises rather than being ignored.

write vcf additionally accepts --reader-workers N and --overshard N (single-file input only; passing either with a directory/manifest raises — the vcf-list form never shards within a contig). The single-file form forwards both straight through to tuning=Tuning(reader_workers=N, overshard=N), so N < 1 raises the same ValueError there — the CLI has no separate check for the lower bound.

• genoray write vcf (SparseVar2.from_vcf/from_vcf_list): source is a single .vcf.gz/.vcf.bgz/.bcf → from_vcf; anything else (a directory, or a file that isn't .vcf.gz/.vcf.bgz/.bcf) → the vcf-list form (a directory of single-sample VCFs/BCFs, or a manifest listing them) → from_vcf_list — a .svar (SVAR1) source belongs under write svar1 instead, not here. --samples/--samples-file work only for the single-file form — they raise for the vcf-list form (each input file already contributes exactly one sample, so there's no cohort to subset). --fields (-f, repeatable) takes bcftools-style INFO/x/FORMAT/x/FMT/x specs, parsed by _parse_cli_field_specs and forwarded as info_fields=/format_fields=; defaults to unset (no fields carried, genotypes only). --chunk-size defaults to 25000. --ploidy (default 2) is accepted here.
• genoray write pgen (SparseVar2.from_pgen): source is a .pgen. No --ploidy (PGEN is diploid-only). --dosages (repeatable) takes NAME=self (read dosage from source itself) or NAME=/path/to/vaf.pgen (read from a separate PGEN), each becoming a DosageField(name=NAME, source=...) passed as dosages=. --chunk-size defaults to a memory-derived value (None). --max-mem (default None = a DETECTED budget, not unbounded) is the same whole-process concurrency-planner budget as from_pgen(max_mem=) (above).
• genoray write svar1 (SparseVar2.from_svar1): source is a *.svar (SVAR1) directory. --fields (repeatable) selects which SVAR1 FORMAT fields carry through (default: all); --empty-fields overrides --fields to carry none. --chunk-size defaults to a memory-derived value (None).
• genoray write-svar1 (top-level, legacy SparseVar.from_vcf/from_pgen → SVAR 1.0): unchanged prior behavior — VCF or PGEN source (auto-detected), --dosages (a FORMAT field name for VCF, or a dosage .pgen path for PGEN), --max-mem (default "1g"), --haploid, --no-symbolic/ --no-breakend (independent flags here, unlike the SVAR2 write subcommands' single --skip-symbolics-and-breakends), --threads/-@, --overwrite. No --regions/--samples/--fields/--reference/ --check-ref/--progress/--log-level/--log-filter, and none of the Tuning-backed flags (--dense-cap/--merge-threads/--sample-interval/ --concurrent-chroms/--reader-workers/--overshard) — those are all SVAR2-write-only (the legacy SparseVar.from_vcf/from_pgen backends don't accept progress=/log_level=/log_filter=/tuning=).

genoray view

# SVAR2 (default) — thin CLI over SparseVar2.write_view
genoray view in.svar2 out.svar2 -r chr1:1-1000 -s A,B
genoray view in.svar2 out.svar2 -r chr1:1-1000          # all samples
genoray view in.svar2 out.svar2 -s A,B                  # all variants (one region per contig)
genoray view in.svar2 out.svar2 -r chr1:1-1000 --no-reroute   # representation-preserving, low-memory view
genoray view in.svar2 out.svar2 -r chr1:1-1000 --reroute      # force the size-optimal re-route

# SVAR 1.0 (previous default) — unchanged SparseVar.write_view CLI
genoray view svar1 in.svar out.svar -r chr1:1-1000 -s A,B --progress

Both subcommands share the same -r/--regions, -R/--regions-file, -s/--samples, -S/--samples-file, -f/--fields, --merge-overlapping, --regions-overlap, --overwrite, -@/--threads, --progress/ --no-progress options and the same no-op guard (at least one of regions/samples is required) and mutex checks (--regions/--regions-file and --samples/--samples-file are each mutually exclusive). genoray view (SVAR2) additionally has --log-level LEVEL (the same six names as write, above) and --log-filter DIRECTIVE (same semantics as write, above); genoray view svar1 does not — its SparseVar.write_view backend has no log_level=/log_filter= kwarg. Neither view subcommand has any Tuning-backed flag — write_view isn't a conversion pipeline, so none of Tuning's knobs apply.

• genoray view (SVAR2, thin wrapper over SparseVar2.write_view): when --regions/--regions-file is omitted, "all variants" defaults to one region per contig (SparseVar2.contigs, since SVAR2 has no contig-length metadata) spanning [0, 2**31 - 1) — every real POS is smaller. --fields defaults to None, meaning no fields are carried through (genotypes only) — this always succeeds, even on a store that has INFO/FORMAT fields. Both --reroute and --no-reroute go through the same slicer backend and carry --fields/--reference identically — there is no longer a fields-carrying vs. genotypes-only split between them:

  • --reroute reruns the var_key/dense routing cost model over the subset — size-optimal (each variant re-routed to whichever representation is smaller for the subset's sample/carrier counts).
  • --no-reroute (reroute=False) slices each variant's existing on-disk representation directly (no cost model, byte-level slice) — representation-preserving regardless of the subset's sample/carrier counts. Recommended for somatic/all-rare cohorts (nearly every variant is already var_key-routed) or memory-constrained runs.
  • Omitting both flags (the default) is "auto": resolves to --no-reroute's behavior when any FORMAT field is carried, to --reroute's otherwise. WHY: a dense→var_key flip stores one value per carrier call and has no slot for a non-carrier sample's FORMAT value, so re-routing a source-dense variant under a FORMAT-carrying view would silently drop it — "auto" prefers fidelity whenever FORMAT is in play and takes the size-optimal re-route otherwise (genotype-only / INFO-only views have no per-sample slot to lose).

  Both --reference (recomputes mutcat from scratch on the subset) and -@/--threads (caps contigs sliced concurrently; autodetected when omitted) are real on both --reroute and --no-reroute — there is no longer an "accepted but ignored/unused" caveat on either path. --progress, --log-level, and --log-filter are all real here — see "write_view progress bar" below for the coarse, one-line-per-contig rendering and log-level semantics. write_view's underlying reroute= kwarg only accepts "auto", True, or False — any other value (e.g. reroute=1) raises ValueError rather than silently falling through to the reroute=False slicer.

• genoray view svar1: unchanged SVAR 1.0 behavior — "all variants" defaults from SparseVar's _contig_stats ([0, pos_max + 1) per contig); --fields defaults to all available fields (use an explicit empty selection to carry none); no --reference/--reroute/--log-level options; --progress shows a real phase-level bar (see below).

genoray concat / genoray split

genoray concat merged.svar2 part1.svar2 part2.svar2       # disjoint-contig merge
genoray split in.svar2 out_dir/                           # explode into out_dir/{contig}.svar2
genoray split in.svar2 subset.svar2 --contigs chr1,chr2    # subset into one store

Both accept --mode (Literal["copy", "hardlink", "symlink", "move"], default "copy" — see SparseVar2.concat/split_by_contig/subset_contigs docstrings) and --overwrite.

Filtering

VCF: pass a genoray.Filter(record=, expr=) value object to filter=. record is a Callable[[cyvcf2.Variant], bool] applied during the genotype scan; expr is the matching polars pl.Expr applied to the .gvi index — VCF requires both halves, bundled together so they can never diverge.

To change a VCF's filter after construction, assign a Filter (or None to clear it) to the vcf.filter setter; the in-memory index is invalidated. The getter returns the Filter | None currently in effect, so vcf.filter = vcf.filter round-trips.

from genoray import VCF, Filter

vcf = VCF("file.vcf", filter=Filter(
    record=lambda v: not v.INFO.get("SVTYPE"),   # cyvcf2 record predicate
    expr=~genoray.exprs.is_symbolic,               # matching .gvi index predicate
))
vcf.filter = None                                  # clear
f = vcf.filter                                      # -> Filter | None

The former two-argument constructor (a separate polars-expression keyword argument alongside filter=) and its tuple-valued vcf.filter getter/setter are removed in 3.0.0 — migrate any code passing the record predicate and polars expression separately to the single Filter(record=, expr=) object shown above.

PGEN: pass a polars pl.Expr returning a boolean mask, operating on the .gvi index columns. Built-in expressions in genoray.exprs (the complete list):

• is_snp (True if all ALT alleles have ILEN == 0; rows with any null ILEN → False)
• is_indel (True if all ALT alleles have ILEN != 0; rows with any null ILEN → False)
• is_biallelic
• is_symbolic (True if any ALT is a VCF 4.x symbolic allele, i.e. starts with <)
• is_breakend (True if any ALT is a VCF 4.x breakend in mate-pair / single-breakend notation, e.g. G[chr2:321[, ]chr2:321]G, .TGCA, TGCA.. A distinct ALT class from symbolic alleles — is_symbolic does not flag breakends)
• is_imprecise (True if any ALT's ILEN is null — an un-sizable symbolic allele or a breakend)
• ILEN (a List[Int32] expression — one value per ALT allele, not a boolean)

ILEN semantics for symbolic SVs. For precise <DEL>/<INS>/<DUP>, ILEN is computed at index-build time from INFO fields: -|SVLEN| for <DEL>, +|SVLEN| for <INS>/<DUP> (falls back to |END - POS| when SVLEN is absent). For VCF, INFO fields are read from header-declared columns (via oxbow); for PGEN, they are parsed from the PVAR INFO string. Non-symbolic ALTs use the literal len(ALT) - len(REF).

Un-sizable symbolic alleles carry null ILEN. An allele is un-sizable when: the IMPRECISE INFO flag is set, SVLEN/END are both missing, the symbolic type is unsupported (<BND>, <CNV>, <INV>, <*>/<NON_REF>), or the ALT is a breakend in mate-pair / single-breakend notation (e.g. G[chr2:321[). At NumPy materialization, null ILEN is coerced to 0 (treated as a point variant).

Filtering guidance (use filter= — a bare pl.Expr for PGEN, a genoray.Filter for VCF):

• ~genoray.exprs.is_symbolic — drops all symbolic alleles (precise or not). Required for haplotype consumers (e.g. genvarloader) that cannot expand any symbolic ALT into literal sequence:

  # PGEN
  pgen = genoray.PGEN("file.pgen", filter=~genoray.exprs.is_symbolic)
  # VCF (both halves required, bundled in a Filter)
  vcf = genoray.VCF(
      "file.vcf.gz",
      filter=genoray.Filter(
          record=lambda rec: not any(a.startswith("<") for a in rec.ALT),
          expr=~genoray.exprs.is_symbolic,
      ),
  )
  

• ~genoray.exprs.is_imprecise — keeps precise symbolic SVs (correctly sized/spanned) and drops only the un-sizable ones (including breakends, which are always un-sizable). Suitable for range/overlap queries where precise SVs are queryable:

  pgen = genoray.PGEN("file.pgen", filter=~genoray.exprs.is_imprecise)
  

• For haplotype consumers, drop all un-expandable ALTs (symbolic and breakends) — breakends are not caught by ~is_symbolic:

  hap_safe = ~genoray.exprs.is_symbolic & ~genoray.exprs.is_breakend
  pgen = genoray.PGEN("file.pgen", filter=hap_safe)
  

For anything else, write pl.col(...) against the .gvi schema — read genoray/exprs.py for the available columns. Combining two exprs expressions with & / | works without importing polars; you only need import polars as pl to build custom predicates.

Reference — quick reference

genoray.Reference is a pysam-backed indexed-FASTA reader used to supply flanking context for mutation-catalogue classification.

ref = genoray.Reference.from_path("hg38.fa")          # auto-creates .fai if absent
ref = genoray.Reference.from_path("hg38.fa", contigs=["chr1", "chr2"])

seq: np.ndarray = ref.fetch("chr1", start=1_000_000, end=1_000_010)
# returns uint8 NDArray, 0-based half-open [start, end)
# bytes(seq) gives the ASCII sequence

Key properties:

• from_path(fasta, contigs=None) — fasta is a str | Path; auto-calls pysam.faidx if the .fai index is missing. contigs filters which contigs the caller cares about (defaults to all in the FASTA).
• fetch(contig, start, end) — 0-based half-open [start, end). Positions outside the contig are N-padded. Returns NDArray[np.uint8].
• contig_array(contig) — the full contig sequence as a cached NDArray[np.uint8]. Shares the one-contig-in-memory cache with fetch. Accepts chr-prefixed or unprefixed names.
• Contig-name agnostic: "chr1" and "1" both resolve correctly (ContigNormalizer under the hood).
• One contig is cached in memory at a time; sequential per-contig access is efficient.

Mutation catalogues (SBS-96 / DBS-78 / ID-83)

write_view progress bar

SparseVar.write_view(..., progress=False) accepts an opt-in progress keyword. When True, a phase-level rich progress bar is shown while the view is written (one tick per major step: counting, genotypes, each carried field, the index build, and mutation annotation when reference= is given). It defaults to False — no bar and no overhead — so library and pipeline callers are unaffected. The genoray view svar1 CLI exposes the same option as --progress (also default off):

genoray view svar1 in.svar out.svar -r chr1:1-1000 -s A,B --progress

The bar is cosmetic: output bytes, schema, and dtypes are identical whether or not it is enabled.

SparseVar2.write_view(..., progress=False, log_level="info", log_filter=None) renders live write progress the same way the from_* writers do (see "Conversion" above), with one difference: unlike the from_* writers, write_view has no per-record stream to sample from, so its progress is COARSE — one line per contig, no within-contig bar movement. In a terminal or Jupyter, progress=True shows a live-updating list of in-flight/finished contigs; elsewhere, a compact "chrom done" line prints as each contig finishes. Regardless of progress, a one-line "[svar2] chrom done: N kept, 0 excluded (Ts)" summary prints per contig once it finishes, unless log_level="off" (slicing never excludes variants, so excluded is always 0). log_level (including the "critical" alias and "warn" rejection) and log_filter behave identically to the from_* writers, above; note write_view has no tuning= parameter — it isn't a conversion pipeline, so none of Tuning's knobs apply. The genoray view CLI exposes both as --progress/--no-progress, --log-level, and --log-filter (see the genoray view CLI section above); genoray view svar1 (SVAR 1.0) exposes only --progress — its SparseVar.write_view backend (above) has no log_level/log_filter kwarg.

Atomic crash-safe writes

Writes are crash-safe and atomic. from_vcf, from_pgen, and write_view build the .svar directory in a hidden sibling staging directory (.<name>.tmp… next to the output) and atomically rename it into place only after the write fully succeeds; .gvi index files are written the same way. A crash mid-write never leaves a partial or corrupt output, and overwriting an existing output preserves it until the replacement is complete. Output bytes are unchanged — this is a durability guarantee only.

Overview

SparseVar supports COSMIC-style mutation catalogues. The workflow is:

1. Call svar.annotate_mutations(reference) once to classify every variant and write mutcat.npy to the .svar directory.
2. Call svar.mutation_matrix(kind) to get a per-sample count matrix.

SparseVar.annotate_mutations

svar = genoray.SparseVar("out.svar")
ref  = genoray.Reference.from_path("hg38.fa")

svar.annotate_mutations(ref)                   # write_back=True (default)
svar.annotate_mutations(ref, write_back=False) # in-memory only; not persisted
svar.annotate_mutations("hg38.fa")             # path accepted directly

Signature: annotate_mutations(reference, *, contigs=None, write_back=True) -> None

• reference — a genoray.Reference instance or a path to a FASTA file (auto-wraps via Reference.from_path).
• contigs=None — if given (a list of contig names), only variants on those contigs are classified; entries on all other contigs are marked NOT_ANNOTATED (sentinel -4) and their contigs are never fetched from the reference. Names match via the ContigNormalizer (chr1/1 both work). Requested contigs absent from the .svar index are skipped with a warning; a listed contig present in the index but absent from the reference raises (omit it from the list to exclude it cleanly). None (default) classifies all contigs. When write_back=True, the normalized scope is recorded in metadata.json as mutcat_contigs (None = all).
• write_back=True — persists mutcat.npy and updates metadata.json so that subsequent SparseVar(dir, fields=["mutcat"]) opens will see the field. Note: write_view never copies mutcat positionally to the output (see below); pass reference= to write_view to recompute it on the subset, or call annotate_mutations on the output view yourself. (This is SparseVar/v1 behavior; on SparseVar2.write_view, reference= recomputes mutcat from scratch on the subset on both reroute=True and reroute=False.)
• write_back=False — the mutcat field lives only in memory (svar.fields["mutcat"]); reopening the file will NOT find it.
• After the call, svar.fields["mutcat"] is populated regardless of write_back.

What it classifies:

Variant type	Channel
Isolated SNV	SBS-96 (trinucleotide context)
Adjacent SNV pair on the same haplotype	DBS-78 (5' entry = DBS code, 3' entry = DBS_PARTNER sentinel)
Runs of ≥ 3 adjacent SNVs	SBS (each stays independent; no DBS collapse)
Native 2 bp MNV in the VCF	DBS-78
MNV > 2 bp, symbolic, non-ACGT	UNCLASSIFIED
Insertion / deletion	ID-83 (size, repeat-context bucketing)
Variant on a contig outside contigs=	NOT_ANNOTATED (excluded from all matrices)

SparseVar.mutation_matrix

svar = genoray.SparseVar("out.svar", fields=["mutcat"])  # pre-load field
df = svar.mutation_matrix("SBS96")                        # default count="allele"
df = svar.mutation_matrix("DBS78", count="sample")
df = svar.mutation_matrix("ID83",  count="allele")

Signature: mutation_matrix(kind, *, count="allele") -> pl.DataFrame

• kind — one of "SBS96", "DBS78", "ID83".
• count="allele" — counts every non-ref allele copy (diploid homozygous = 2).
• count="sample" — counts each category at most once per sample (presence/absence).
• Returns a Polars DataFrame with a MutationType string column followed by one Int64 column per sample. Rows are in fixed COSMIC codebook order (96 / 78 / 83 rows respectively).
• Requires the mutcat field to be available: either loaded at open time with fields=["mutcat"], or already in memory from a prior annotate_mutations call, or present on disk from a prior annotate_mutations(write_back=True). Raises ValueError if none of those hold.

The mutcat field

mutcat is an int16 field stored per genotype entry (same ragged layout as genos). The int16 code space is:

Range	Channel
[0, 96)	SBS-96
[96, 174)	DBS-78
[174, 257)	ID-83
-1	DBS_PARTNER — 3' half of an adjacent SNV pair; never counted
-2	UNCLASSIFIED — symbolic / complex / MNV > 2 bp / non-ACGT
-3	MISSING — reserved sentinel (defined in the code space but not emitted by annotate_mutations v1; SparseVar stores only ALT-carrying entries, so no-call slots do not appear in the ragged field)
-4	NOT_ANNOTATED — entry on a contig outside the contigs= annotation scope; never counted

To read a previously annotated file:

svar = genoray.SparseVar("out.svar", fields=["mutcat"])
# svar.fields["mutcat"] is a Ragged[int16] mirroring svar.genos

v1 scope limits (no strand-bias; calibrated against PCAWG/SigProfiler rules)

• No strand-bias separation (SBS-192 / SBS-384) — v1 SparseVar only. Use SparseVar2.annotate_mutations(reference, gtf=...) for transcriptional strand-resolved catalogs.
• DBS collapse applies only to isolated adjacent pairs on the same haplotype. Runs of ≥ 3 adjacent SNVs stay as individual SBS entries.
• Indel channel (ID-83) bucketing follows PCAWG/SigProfiler published rules and is pinned by the unit tests in tests/test_mutcat.py. Cross-validation against SigProfilerMatrixGenerator is deferred (it is not a declared dependency).

Signature refitting (COSMIC)

Decompose a catalogue into per-sample COSMIC signature activities.

import genoray
from genoray import Forward, Spa

ref = genoray.cosmic_signatures("SBS96")        # pooch-fetched + cached
cat = svar.mutation_matrix("SBS96")              # MutationType + sample cols
act = genoray.fit_signatures(cat, ref)           # activities + cosine_similarity

# strategy=: pick the refit algorithm and its parameters explicitly
act = genoray.fit_signatures(cat, ref, strategy=Forward(max_delta=0.02))
act = genoray.fit_signatures(cat, ref, strategy=Spa())  # SPA's cosmic_fit, faithfully

# convenience: mutation_matrix -> fit_signatures in one call
act = svar.assign_signatures("SBS96")                       # default COSMIC ref
act = svar.assign_signatures("SBS96", reference=ref, min_activity=0.01)
act = svar.assign_signatures("SBS96", reference="my_sigs.txt")  # TSV path
act = svar.assign_signatures("SBS96", strategy=Spa())       # same strategy= as fit_signatures

Signatures:

• cosmic_signatures(kind, *, version="3.4", genome="GRCh38") -> pl.DataFrame — fetches/caches the COSMIC reference set for kind ∈ {"SBS96","DBS78","ID83"}. Returns a MutationType column (canonical codebook order) + one column per signature. genome is ignored for ID83.

• fit_signatures(catalogue, reference, *, strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame — refits catalogue against reference. Aligns rows by joining on MutationType (raises ValueError if the catalogue has a type missing from the reference). Returns one row per sample: Sample, one Float column per signature (activities; 0.0 if unselected), and cosine_similarity. n_jobs=1 (default) is serial; n_jobs=-1 uses all cores. Results are identical regardless of n_jobs/backend.

  strategy= selects the algorithm and cannot be combined with max_delta/min_activity/criterion (raises ValueError if both are given, even at their default values) — those three are permanent shorthand for strategy=Forward(...), not deprecated:

  • Forward(max_delta=0.01, min_activity=0.005, criterion="cosine") (the default when strategy is omitted) — genoray's own sparse greedy forward selection: NNLS + criterion-guided add + min-activity prune. criterion="cosine" stops on cosine-similarity plateau (scale-invariant, so blind to mutation burden); criterion="bic" stops on a Poisson-BIC test (burden-aware, more consistent above ~1,000 mutations, but over-selects at very low burden). Returned activities are raw NNLS weights — they need not sum to anything in particular.
  • Spa(metric="l2", initial_remove_penalty=0.05, add_penalty=0.05, remove_penalty=0.01, background_sigs=("SBS1","SBS5"), connected_sigs=True, activity_scale="burden") — a from-scratch reimplementation of SigProfilerAssignment's cosmic_fit (solver="nnls", pcawg_rule=False), pure numpy/scipy/polars with no SigProfiler dependency: saturate over every reference signature, prune backward on relative L2 error, then refine with add-remove layers. background_sigs names are force-added to the working set at the start of every refinement layer, but they are not absolutely protected from the removal sweeps: SPA's own protection decays partway through a sweep and genoray reproduces that, so a background signature the sample does not need is still dropped (names absent from reference are ignored, which is what makes the SBS1/SBS5 default inert on DBS78/ID83 references); connected_sigs=True force-adds co-occurring SBS partners (e.g. SBS2/SBS13) once any member is selected, False disables it, or pass your own Sequence[Sequence[str]] of signature-name groups. Under the default activity_scale="burden", activities are rescaled and integer-rounded so each sample's signature columns sum to its mutation burden; activity_scale="raw" returns unscaled NNLS weights instead, like Forward does. Deterministic (touches no RNG), unlike upstream SPA. Recovers more true signatures than Forward at every burden measured but, like Forward's "cosine" criterion, is scale-invariant and so not a consistent estimator either. Roughly an order of magnitude slower per sample than Forward (measured ~7-27x on real COSMIC SBS96, 86 signatures, depending on the catalogue).

• SparseVar.assign_signatures(kind, *, reference=None, count="allele", strategy=None, max_delta=0.01, min_activity=0.005, criterion="cosine", n_jobs=1, backend="loky") -> pl.DataFrame — mutation_matrix(kind, count=...) then fit_signatures(...). reference accepts a pl.DataFrame, a TSV path, or None (defaults to cosmic_signatures(kind)). Forwards n_jobs/backend to fit_signatures for per-sample parallelism (n_jobs=1 (default) is serial; n_jobs=-1 uses all cores). strategy= and criterion= are forwarded exactly like fit_signatures itself — pass strategy=Spa() for SigProfilerAssignment's algorithm, or leave strategy None (default) and use max_delta/min_activity/criterion to configure the Forward shorthand. Passing strategy= together with any of max_delta/min_activity/criterion raises ValueError, exactly like fit_signatures (omitted kwargs take their defaults). SparseVar2.assign_signatures has the identical signature, except that it additionally accepts contigs= to restrict the catalogue to a subset of the store's contigs.

Refit strategies

fit_signatures(strategy=...) selects the algorithm. Both strategies are frozen dataclasses carrying only their own parameters.

• Forward(max_delta=0.01, min_activity=0.005, criterion="cosine") — greedy forward selection from the empty set. The default. Activities are raw NNLS weights and need not sum to the sample's mutation burden.
• Spa(metric="l2", initial_remove_penalty=0.05, add_penalty=0.05, remove_penalty=0.01, background_sigs=("SBS1","SBS5"), connected_sigs=True, activity_scale="burden") — SigProfilerAssignment's cosmic_fit. Every default is SPA's own. Under the default activity_scale="burden", activities are integer-valued floats summing exactly to the sample's burden; activity_scale="raw" gives unscaled NNLS weights instead.

background_sigs names absent from the reference are ignored, so the default is inert for DBS78 and ID83. connected_sigs=True uses SPA's four SBS groups (genoray.SPA_CONNECTED_GROUPS: SBS2/13, SBS7a-d, SBS10a/b, SBS17a/b); pass your own groups as a sequence of sequences, or False to disable.

Neither Forward(criterion="cosine") nor Spa is burden-aware. For whole-genome catalogues where consistency matters, use Forward(criterion="bic").

The output schema is unchanged by the strategy: Sample, one Float column per reference signature, and a trailing cosine_similarity.

Out of scope (v1): de novo extraction, opportunity normalization, bootstrap CIs, plotting.

Common mistakes

Mistake	Fix
genoray.Genos8	genoray.VCF.Genos8 (class attribute)
vcf.read(..., phasing=True)	Set phasing=True on the VCF() constructor
Reading dosages from a VCF without dosage_field=	Pass dosage_field="DS" (or appropriate Number=A field) on the constructor
Putting a dosage-only PGEN in the main path when you also have hardcalls	Hardcalls in main path, dosages in dosage_path=
Importing from genoray._vcf import VCF	Use from genoray import VCF
Expecting VCF to have read_ranges	VCF doesn't; loop over single-range read calls, or use PGEN/SparseVar
Treating svar.index["POS"] as 0-based	It's 1-based; subtract 1 to compare with query coords
Calling read_ranges and assuming a flat array	PGEN returns (data, offsets); SparseVar returns a Ragged (or awkward record with fields)
Calling mutation_matrix without a mutcat field	Run annotate_mutations first, or open with fields=["mutcat"]
Expecting mutation_matrix to auto-run annotation	It does not; call annotate_mutations separately
Re-opening SparseVar and losing the mutcat field	Use write_back=True (default) in annotate_mutations; then open with SparseVar(dir, fields=["mutcat"])
Calling write_view and expecting mutcat to be in the output	write_view never copies mutcat positionally (subsetting invalidates DBS adjacency codes). Pass reference= to write_view to recompute it on the subset, or call annotate_mutations on the output view yourself. Explicitly including "mutcat" in fields= without a reference= raises ValueError. On SparseVar2.write_view, reference= recomputes mutcat from scratch on both reroute=True and reroute=False.
Passing a FORMAT field in fields= to SparseVar2.write_view(..., reroute=True) and expecting it dropped/rejected	Both reroute=True and reroute=False carry fields through now (previously reroute=True raised ValueError). Watch the reroute="auto" default instead: it resolves to reroute=False whenever any FORMAT field is carried, because a dense→var_key flip has no slot for a non-carrier sample's FORMAT value.
Passing the source dataset directory as output to write_view (even with overwrite=True)	Raises ValueError — writing in place would delete the source before the view is written. Pass a different output path.
Passing a FASTA path directly to annotate_mutations	Supported — it auto-wraps via Reference.from_path
rag.fields["AF"] on a SparseVar2.decode() result	Ragged.fields is a list[str] of names, not a mapping; index the field itself with rag["AF"]
Expecting SparseVar2.decode() to include INFO/FORMAT fields	Fields are opt-in — pass fields=[...] to SparseVar2(...) or call .with_fields([...]) first

When this skill needs updating

Any PR that adds, removes, renames, or changes the semantics of a public name (anything reachable from import genoray without underscores) must update this skill alongside the code change. See the project CLAUDE.md.