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 fromvcf (supports regions=/samples=/mergeoverlapping=/regionsoverlap=), PLINK2 PGEN → SVAR2 conversion via frompgen, N single-sample VCFs/BCFs → one SVAR2 store via a native k-way merge in fromvcflist (reference/noreference supported like fromvcf, absent sites fill hom-ref; supports regions=/mergeoverlapping=/regionsoverlap= but no samples= — the cohort is the file set), SVAR1 (SparseVar) → SVAR2 native migration via fromsvar1 (reads no VCF/htslib; biallelic SVAR1 only; supports regions=/samples=/mergeoverlapping=/regionsoverlap= like fromvcf/frompgen); range queries via decode/regioncounts/readranges; mutational-signature support (SBS96/DBS78/ID83) via annotatemutations/mutationmatrix/assignsignatures, or classify during the write with fromvcf(signatures=True)/frompgen(signatures=True)/fromsvar1(signatures=True); scalar-numeric INFO/FORMAT field extraction during the write via fromvcf(infofields=, formatfields=)/fromvcflist(infofields=, formatfields=) (fromvcflist merges INFO first-carrier-wins, FORMAT per-sample); frompgen 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 infofields=/formatfields= (PGEN has no VCF INFO/FORMAT); fromsvar1 carries SVAR1's existing fields through selectively via fields= (None default = all, [] = none, or a name subset) — read back opt-in via fields=/withfields/availablefields and attached to decode's result)
• genoray.InfoField / genoray.FormatField — frozen dataclasses (name, dtype=None, default=None) configuring a single INFO/FORMAT field for SparseVar2.from_vcf; a bare str name uses inferred defaults instead
• genoray.DosageField — frozen dataclass (name="dosage", source="self"|Path, dtype="f16"|"f32"="f32", default=None) configuring a PGEN dosage FORMAT field for SparseVar2.from_pgen
• genoray.exprs — polars filter expressions for .gvi indexes
• genoray.cosmic_signatures — fetch/cache COSMIC reference signatures
• genoray.fit_signatures — sparse forward-selection signature refit

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

Where to look for details

Prefer reading these over guessing:

• docs/source/index.md — narrative tour with full examples (VCF, PGEN, filtering, chunking)
• docs/source/svar.md — SparseVar usage
• genoray/init.py — confirms the public surface
• genoray/vcf.py — VCF class: constructor, read, chunk, mode constants near the top of the class; getrecord_info(contig=None, start=None, end=None, fields=None, info=None, lazy=False) — non-FORMAT record-level fields (including INFO) for a range or the whole file, returns pl.DataFrame (or pl.LazyFrame when lazy=True)
• genoray/pgen.py — PGEN class: constructor, read, chunk, readranges, chunk_ranges, mode constants near the top of the class
• genoray/svar.py — SparseVar: init, fromvcf, frompgen, readranges, readrangeswithlength(contig, starts=0, ends=POSMAX, samples=None) (length-guaranteed range read; returns the same type as readranges — a Ragged or fields-augmented record), withfields, annotatemutations, mutationmatrix, assignsignatures, annotatewithgtf(gtf, levelfilter=1, writeback=True, *, strandencoding=None, codonnulltoken=None) (GTF CDS annotation entry point, returns pl.DataFrame with varID/geneid/strand/codonpos), cache_afs() (computes and persists an AF column to the .gvi index; returns None)
• genoray/svar2.py — SparseVar2: init(path, *, fields=None), withfields(fields) (new reader over the same store with those fields selected), availablefields (dict[str, StoredField], set in init), fromvcf (VCF/BCF → SVAR2 conversion entry point, signatures= classifies during the write, infofields=/formatfields= extract scalar-numeric fields during the write; supports regions=/samples=/mergeoverlapping=/regionsoverlap=), frompgen (PLINK2 PGEN → SVAR2 conversion entry point; diploid-only, no ploidy=/infofields=/formatfields=; 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=/mergeoverlapping=/regionsoverlap= like fromvcf), fromvcflist (N single-sample VCFs/BCFs → one SVAR2 store via a native k-way merge; sources accepts a Sequence/directory/manifest, resolved by module-level resolvevcfsources; reference/noreference supported (noreference skips left-alignment, so cross-file joins require pre-normalized inputs); infofields=/formatfields= supported — INFO merges first-carrier-wins, FORMAT stays per-sample; supports regions=/mergeoverlapping=/regionsoverlap= like fromvcf, but no samples= — the cohort is the file set), fromsvar1 (SVAR1 (SparseVar) → SVAR2 native migration entry point; reads no VCF/htslib, ploidy from SVAR1 metadata, biallelic SVAR1 only, no infofields=/formatfields= (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=/mergeoverlapping=/regionsoverlap= like fromvcf/frompgen, though regions filter per-record rather than narrowing a covering range up front); nsamples/availablesamples/contigs/ploidy metadata. Read/query methods live in the mixins: genoray/svar2decode.py (decode — attaches one Ragged per selected field, regioncounts), genoray/svar2batch.py (public readranges; internal gvl-only overlapbatch/findranges/gatherranges), and genoray/svar2mutcat.py (annotatemutations, mutationmatrix, assignsignatures — COSMIC mutational-signature workflow, mirroring SparseVar's but backed by a per-contig Rust sidecar instead of a .gvi-attached field)
• genoray/svar2fields.py — InfoField/FormatField/DosageField dataclasses + FieldDtype and the header/dtype validation used by fromvcf(infofields=, formatfields=); parseclifieldspecs (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.availablefields — not exported at top-level genoray, only reached via that dict
• genoray/_cli/main.py — the genoray CLI (index, write vcf/write pgen/write svar1 (all → SVAR2), top-level write-svar1 (legacy VCF/PGEN → SVAR1), view / view svar1, concat, split)
• genoray/signatures.py — cosmicsignatures, fit_signatures
• genoray/reference.py — Reference: frompath, fetch, contig_array
• genoray/exprs.py — the complete set of pre-built filter expressions (currently 7: issnp, isindel, isbiallelic, issymbolic, isbreakend, isimprecise, 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.fromvcf, SparseVar2.frompgen, and SparseVar2.fromvcflist's maxmem, 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/chunkranges, 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.fromvcf/frompgen (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 currentsamples/setsamples. Instead, every read method (readranges, readrangeswithlength, 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.fromvcf / frompgen inherit and apply the source's filter — filter the VCF/PGEN to filter the SVAR.

SparseVar.fromvcf / frompgen accept regions=, samples=, mergeoverlapping=, regionsoverlap= 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)

withfields(False) drops all extras and returns a plain Ragged[VIDX_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 (fromvcf, below) writes a store; range queries (decode / regioncounts / read_ranges, further below) read it back. All coordinates are 0-based half-open [start, end), as everywhere else in genoray.

Conversion

from genoray import SparseVar2

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

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

Signature: fromvcf(out, source, reference=None, , regions=None, samples=None, mergeoverlapping=False, regionsoverlap="pos", noreference=False, skipoutofscope=False, ploidy=2, chunksize=25000, threads=None, overwrite=False, longallelecapacity=81024*1024, signatures=False, infofields=None, formatfields=None, checkref="e", progress=False, loglevel="info", maxmem=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=/mergeoverlapping=/regionsoverlap= — 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 mergeoverlapping=True. regionsoverlap 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.

• skipoutof_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 skipoutof_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 formatfields= instead (e.g. a DS FORMAT field). No haploid= OR-collapse, no maxmem-based chunking (use chunk_size instead) — those two remain SparseVar (SVAR 1.0)-only for now.

• threads=None — thread budget (autodetected if None). For single-file

input, conversion shards within a contig once the budget clears HTSlib's decode-thread allocation (roughly 15+ threads); below that it runs one un-sharded reader. Sub-contig sharding is driven entirely by this existing threads value — no separate knob — and output is byte-identical to serial conversion at every thread count (fromvcflist, the N-single-sample-VCF merge path, does not shard within a contig). 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/varkey cost model). Requires a reference (reference=); raises ValueError if combined with noreference=True. There is no public read-side API for the SVAR2 mutcat sidecar yet (unlike SparseVar.annotatemutations/ mutationmatrix, below) — this flag only controls whether the sidecar is written.

• infofields=/formatfields= — Sequence[str | InfoField] /

Sequence[str | FormatField], None by default. Extracts scalar-numeric INFO/FORMAT fields into the store during the write:

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

SparseVar2.fromvcf( "out.svar2", "file.vcf.gz", "ref.fa", infofields=["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 (INTMIN 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 varkey-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=…) / .withfields(…) / .availablefields 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 maxmem convention above ("4g", "512m", "2GB", parsed by parsememory), and the same whole-process meaning as fromvcflist's maxmem (below) — fromvcflist just has no fitted concurrency planner to spend it on (its contigs run strictly sequentially), so it derives its own per-chunk chunksize 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-maxmem 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 at fromvcf's own defaults (chunksize=25000, readerworkers=3) the per-contig bracket's kappa term dominates those two numbers completely, so the real floor for even a tiny cohort is roughly 1.38 GB, not ~600 MB (see the S=4,000 figure below) — anything much below that is rejected in practice. frompgen'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 quadruples fromvcf's real-world floor at its own defaults (chunksize=25000, readerworkers=3, cc=1) versus the pre-refit law: the minimum maxmem for one concurrent contig goes ~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. The direction is safe (a larger requirement means more over-allocation or an outright refusal to plan, never an OOM), but the size of the jump is large enough to flip outcomes at production scale. At S=500,000 with fromvcf's own defaults (chunksize=25000, readerworkers=3), the floor for concurrentchroms=1 is ~101,480 MB: a 64 GB host (maxmem defaults to 80% of detected RAM, i.e. 51,200 MB) now raises PlanError::InsufficientMemory — naming both remedies in its message, "raise maxmem or lower chunksize" — where the pre-refit law planned concurrentchroms=2 and ran; a 128 GB host** (102,400 MB) still clears and plans concurrent_chroms=1, but by under 1% of headroom, where the pre-refit law planned 4.

• progress=False/log_level="info" — write-time progress/logging,

shared by fromvcf/frompgen/fromvcflist/fromsvar1/writeview. 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 loglevel="off". loglevel is the minimum severity for structured write-time log lines — "off" (disables everything, including the per-contig summaries and progress rendering — a pure no-op, zero overhead), "warning", "info" (default; also includes thread-budget selection, per-contig start/finish, and contig-name resolution against the reference when it differs from the source's own spelling), or "debug" (additionally surfaces per-record detail: a record excluded for a REF/FASTA mismatch, and each indel that gets left-aligned). The GENORAYLOG environment variable overrides the loglevel argument when set to one of the same four values (e.g. GENORAYLOG=debug), without touching call sites. Structured log lines render their fields inline as key=value pairs after the message (e.g. pipeline config concurrentchroms=8 readerworkers=4), matching what GENORAYLOG's stderr layer emits.

Conversion from PGEN

from genoray import SparseVar2

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

Signature: frompgen(out, source, reference=None, , regions=None, samples=None, mergeoverlapping=False, regionsoverlap="pos", noreference=False, skipoutofscope=False, chunksize=None, maxmem=None, threads=None, overwrite=False, longallelecapacity=81024*1024, signatures=False, dosages=None, checkref="e", progress=False, log_level="info") -> int

• source — a .pgen file. Variant metadata is read from the sibling

.pvar/.pvar.zst, sample names from the sibling .psam. reference/noreference, skipoutofscope, overwrite, longallelecapacity, signatures, and checkref all mean the same as fromvcf (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 OMPNUMTHREADS 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).

• chunksize=None — unlike fromvcf's fixed 25_000 default, None here

derives a variant-count budget from sample count (a packed dense chunk costs chunksize nsamples 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 fromvcflist's chunksize entry below for the details. dosages counts as nformat_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 maxmem convention above ("4g", "512m", "2GB", parsed by parsememory), and the same whole-process meaning as fromvcf's maxmem (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). fromvcflist's maxmem 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 chunksize 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-maxmem 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 fromvcf's ~457 MB), so any budget that can't cover baseline plus one concurrent contig's chunk buffers is rejected with ValueError, even for a tiny cohort. That baseline scales with cohort size (~0.0158 MB/sample), so this isn't just a small-cohort concern: at ~500k samples it alone predicts ~10.6 GB, so a detected budget on a smaller host will reject the conversion — pass an explicit max_mem sized to the host in that case.

• regions=/mergeoverlapping=/regionsoverlap= — same convention,

semantics, and three overlap modes ("pos"/"record"/"variant") as fromvcf, restricting conversion to one or more .pvar variant-index ranges. As with fromvcf, "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 fromvcf): preserves caller order, de-duplicates first occurrences, raises ValueError on an unknown name. availablesamples and every decoded column match the caller's order exactly, regardless of each sample's original .psam position.

• No infofields=/formatfields= — PGEN carries no FORMAT, and .pvar

INFO extraction is not implemented.

• dosages=Sequence[DosageField] — stores per-sample dosage tracks as

FORMAT fields. Each DosageField(name="dosage", source="self"|Path, dtype="f16"|"f32"="f32", default=None): source="self" reads dosages from the hardcall .pgen (source above) itself; a Path reads from a separate .pgen (e.g. a VAF/CCF file kept apart because pgenlib derives hardcalls from dosage when both live in one file) — it must share the hardcall .psam's samples and align 1:1 on the hardcall .pvar's variants. Stored genotype-aligned like any FORMAT field: under var_key routing a non-carrier's dosage is dropped (harmless for VAF/CCF-style fields, ~0 for non-carriers). Read back the same way as other FORMAT fields — see "Reading INFO/FORMAT fields (SVAR2)" below.

```python from genoray import SparseVar2, DosageField

SparseVar2.frompgen("out.svar2", "cohort.pgen", "ref.fa", dosages=[DosageField(name="DS", source="self")]) # separate dosage file (e.g. VAF stored as dosage): SparseVar2.frompgen("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/loglevel="info" — same as fromvcf (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: fromvcflist(out, sources, reference=None, , regions=None, mergeoverlapping=False, regionsoverlap="pos", noreference=False, skipoutofscope=False, ploidy=2, chunksize=None, maxmem=None, threads=None, overwrite=False, longallelecapacity=81024*1024, signatures=False, infofields=None, formatfields=None, checkref="e", progress=False, loglevel="info") -> int

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

• regions=/mergeoverlapping=/regionsoverlap= — same convention,

semantics, and three overlap modes ("pos"/"record"/"variant") as fromvcf, applied identically to every input file in the merge. As with fromvcf, "variant" mode keeps a multiallelic record whole if ANY of its alleles truly overlaps the region.

• No samples= parameter — unlike fromvcf/frompgen/from_svar1,

fromvcflist 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

resolvevcf_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 regioncounts. 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 fromvcf, 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.

• noreference=True is supported, same as fromvcf/from_pgen: skips

REF validation and left-alignment, reconstructing each atom's REF from the record's own REF bytes. The reference/noreference exactly-one-of check and the signatures+noreference incompatibility are otherwise identical to fromvcf. - 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 noreference: 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 noreference with fromvcf_list.

• infofields=/formatfields= — 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=).

• chunksize=None — unlike fromvcf's fixed 25_000 default, None here

derives a budget-based chunk size from the cohort size (autochunksize), so one packed dense chunk stays within a maxmem-derived per-chunk target (see maxmem below), up to a ~256 MiB ceiling; this is the same default frompgen/fromsvar1 already use. The budget accounts for both the packed genotype grid AND any staged formatfields (nformatfields nsamples 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 chunksize falls below 256 variants, autochunksize warns (naming the budget, the derived chunksize, and nsamples/ploidy/nformatfields) and keeps the small value rather than raising — e.g. at ploidy=2, nformatfields=7 this fires above nsamples≈37,118 against the default ~256 MiB budget. Raise maxmem or request fewer formatfields to clear it. Shared by all three converters that derive chunksize via this helper (frompgen, fromsvar1, fromvcf_list).

• max_mem: int | str | None = None — byte budget the whole process may

use (same string forms as the module-level maxmem convention, e.g. "4g", parsed by parsememory) — the same meaning as fromvcf's maxmem (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 (fromvcflist 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 maxmem as a direct cap on the bytes of one in-flight dense chunk. It is now a whole-process budget. A call like maxmem="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.

fromvcflist has no fitted concurrency planner the way the sharded fromvcf reader does — its contigs convert strictly sequentially — so instead of planning concurrency it derives a per-chunk byte target: chunktarget = min(DENSECHUNKTARGETBYTES, maxmem // (concurrentjobs * inflightchunksperjob)). concurrentjobs and inflightchunksperjob are not hardcoded Python literals — they are read live from the Rust extension module (core.VCFLISTCONCURRENTCHROMS, currently 1, and core.VCFLISTDENSECHANNELCAP + 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). DENSECHUNKTARGETBYTES 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 chunksize is passed explicitly — including skipping memory-budget DETECTION entirely (no attempt, no warning) when chunksize 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, skipoutof_scope, threads, overwrite,

longallelecapacity, signatures, checkref all mean the same as fromvcf, and the return value is the same int (dropped out-of-scope ALTs). checkref is applied per input file during the merge (ignored when noreference=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/loglevel="info" — same as fromvcf (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: fromsvar1(out, source, reference=None, , regions=None, samples=None, mergeoverlapping=False, regionsoverlap="pos", noreference=False, skipoutofscope=False, chunksize=None, threads=None, overwrite=False, longallelecapacity=81024*1024, signatures=False, fields=None, checkref="e", progress=False, loglevel="info") -> 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,

skipoutofscope, overwrite, longallelecapacity, signatures, and checkref 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 frompgen/fromvcflist (autochunksize) and warns under the same below-256-variant condition — see fromvcflist's chunksize entry above for the details. Known gap: this call site always passes nformatfields=0, even though fields= (below) selects SVAR1 FORMAT fields and defaults to carrying all of them — so unlike frompgen, the budget here does not account for staged FORMAT bytes and can under-size the chunk when fields= carries a wide FORMAT set. Tracked in #157.

• 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=/mergeoverlapping=/regionsoverlap= — same convention,

semantics, and three overlap modes ("pos"/"record"/"variant") as fromvcf/frompgen. "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 fromvcf/frompgen): 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 fromvcf only for varkey (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 infofields=/formatfields= kwargs — those names are VCF-specific;

use fields= (above) instead to select which SVAR1 fields carry.

• progress=False/loglevel="info" — same as fromvcf (above).

Range queries

Open a finished store, then query per contig. Construction reads meta.json and opens one native reader per contig, exposing .availablesamples (list; the canonical sample-name accessor shared with VCF/PGEN/SparseVar), .nsamples, .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/regioncounts/readranges 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 / regioncounts / readranges (above). readranges(contig, starts, ends, samples=None) is a fused search+gather; starts/ends are parallel 1D arrays (mirrors SparseVar.readranges), samples selects/reorders a subset by name. It returns the raw two-channel BatchResult → numpy dict, a TypedDict with a fixed field set: vkpos/vkkey/vkoff, densepos/densekey/ denserange/densepresent/densepresentoff, lutbytes/lutoff, and scalars nregions/n_samples/ploidy.

SparseVar2 also has overlapbatch/findranges/gatherranges (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 fromvcf(infofields=…, 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 availablefields below). withfields 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.svar2fields, 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 fromvcf's formatfields=

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 mutationmatrix — either post-hoc, or by passing signatures=True to fromvcf (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.annotatemutations, there is no writeback= 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.

• mutation_matrix(kind, *, count="allele"|"sample") -> 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. Raises ValueError if called before the store is annotated (no on-disk sidecar for every contig) — annotate first, either via annotatemutations or fromvcf(..., signatures=True). "SBS192"/"SBS384" additionally require strand annotation (annotatemutations(..., gtf=...)) and raise ValueError if the store lacks it. assignsignatures does not accept "SBS192"/"SBS384" — see below.

• assignsignatures(kind, *, reference=None, count="allele", maxdelta=0.01, minactivity=0.005, njobs=1, backend="loky") -> pl.DataFrame

— mutationmatrix(kind, count=...) then genoray.fitsignatures(...). reference accepts a pl.DataFrame, a TSV path, or None (defaults to genoray.cosmic_signatures(kind)).

• 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

fromvcf(..., signatures=True) stays strand-free; obtain strand catalogs via a post-hoc annotatemutations(reference, gtf=...).

• assign_signatures("SBS192"/"SBS384") raises NotImplementedError: COSMIC

publishes no strand-resolved reference set. Use mutation_matrix for strand-bias analysis.

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 writeview (see the CLI section below), none of these methods re-run conversion or the varkey/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.

• splitbycontig(out_dir, *, mode="copy", overwrite=False) -> list[Path] —

explode into one single-contig store per contig at outdir/{contig}.svar2; returns the output paths in self.contigs order. Implemented as one subsetcontigs 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 skipoutof_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 skipoutofscope=; 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), and --progress/--no-progress + --log-level {off,warning,info,debug} (map to progress=/loglevel=; see "Conversion" above for behavior — default --no-progress --log-level info). write vcf's vcf-list form forwards both to fromvcflist; its single-file form forwards both to from_vcf.

• genoray write vcf (SparseVar2.fromvcf/fromvcf_list): source is a

single .vcf.gz/.vcf.bgz/.bcf → fromvcf; 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) → fromvcflist — 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 parseclifieldspecs and forwarded as infofields=/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 frompgen(maxmem=) (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.fromvcf/frompgen

→ 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 — those are SVAR2-write-only (the legacy SparseVar.fromvcf/frompgen backends don't accept progress=/log_level=).

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 {off,warning,info,debug}; genoray view svar1 does not — its SparseVar.writeview backend has no loglevel= kwarg.

• 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 varkey/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 varkey-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 and --log-level are both real here — see "writeview progress bar" below for the coarse, one-line-per-contig rendering and log-level semantics. writeview'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 contigstats ([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/splitbycontig/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 <)
• isbreakend (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 — issymbolic 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:

``python # PGEN pgen = genoray.PGEN("file.pgen", filter=~genoray.exprs.issymbolic) # 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.issymbolic, ), ) ``

• ~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:

``python 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:

``python hapsafe = ~genoray.exprs.issymbolic & ~genoray.exprs.isbreakend pgen = genoray.PGEN("file.pgen", filter=hapsafe) ``

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.writeview(..., progress=False, loglevel="info") renders live write progress the same way the from writers do (see "Conversion" above), with one difference: unlike the from writers, writeview 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 loglevel="off" (slicing never excludes variants, so excluded is always 0). loglevel and the GENORAYLOG env override behave identically to the from* writers. The genoray view CLI exposes both as --progress/--no-progress and --log-level (see the genoray view CLI section above); genoray view svar1 (SVAR 1.0) exposes only --progress — its SparseVar.writeview backend (above) has no log_level kwarg.

Atomic crash-safe writes

Writes are crash-safe and atomic. fromvcf, frompgen, 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.

1. 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: annotatemutations(reference, *, contigs=None, writeback=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 NOTANNOTATED (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 writeback=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: writeview never copies mutcat positionally to the output (see below); pass reference= to writeview to recompute it on the subset, or call annotatemutations on the output view yourself. (This is SparseVar/v1 behavior; on SparseVar2.writeview, 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:

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 annotatemutations call, or present on disk from a prior annotatemutations(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:

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

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

# 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

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.

• fitsignatures(catalogue, reference, *, maxdelta=0.01, minactivity=0.005, njobs=1, backend="loky") -> pl.DataFrame

— sparse forward-selection refit (NNLS + cosine-guided add + min-activity prune). 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 (counts; 0.0 if unselected), and cosinesimilarity. njobs=1 (default) is serial; njobs=-1 uses all cores. Results are identical regardless of njobs/backend.

• SparseVar.assignsignatures(kind, *, reference=None, count="allele", maxdelta=0.01, minactivity=0.005, njobs=1, backend="loky") -> pl.DataFrame

— mutationmatrix(kind, count=...) then fitsignatures(...). reference accepts a pl.DataFrame, a TSV path, or None (defaults to cosmicsignatures(kind)). Forwards njobs/backend to fitsignatures for per-sample parallelism (njobs=1 (default) is serial; n_jobs=-1 uses all cores).

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

Common mistakes

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.