Repertoire signatures

A signature is a fixed-order, name-addressed feature vector for one repertoire, on a scale a downstream model can consume without fitting a scaler of its own. vdjtools emits the statistics half, vsig; the geometry half, rsig, comes from mirpy and shares this contract exactly.

It is built in three stages, and a corpus — a large published reference collection of repertoires — fixes the last two:

raw features  ->  winsorization bounds  ->  rotation + per-PC scaling  ->  vsig
(this sample                (all three from one corpus artifact)
 alone)                                                              +  channels, untouched

In plain terms, reading left to right:

1. Raw features are the ordinary measurements of one repertoire — how diverse, how clonal, which V and J genes, what junction lengths, which isotypes. They come from this sample and nothing else.

2. Winsorization clips each measurement to a range measured on the corpus, so that one pathological sample cannot dominate everything downstream. Clipping edits a real number, so how much of the row was clipped is always reported — see winsor_frac below.

3. The rotation re-expresses those correlated measurements as uncorrelated axes (principal component s) ordered by how much variation each carries across the corpus, and each axis is put on a robust z-score scale: how far this sample sits from the corpus’s typical one, in robust deviations. That is what makes your matrix and a collaborator’s directly comparable without either of you fitting a scaler.

Channel s sit outside all of this. They are carried in their own units, never rotated and never clipped, because a number whose job is to tell you whether to trust a row must not be mixed into the row.

If you want the underlying measurements back rather than the rotated axes, ask for them: --named emits the raw groups the rotation is fitted on under their own names. See Getting the named statistics back.

Quickstart

# nine corpora are published; name one and it is fetched on first use
vdjtools signature --corpus blood samples/*.tsv.gz -o sig.tsv

# what exactly will I get?
vdjtools signature --corpus blood --components 32 --describe

# pre-warm the cache instead of fetching lazily -- the install-time step
vdjtools corpus --fetch all

# or build your own synthetic corpus, which needs no cohort at all
vdjtools corpus --corpus synthetic-tissue -o /tmp/st.npz
vdjtools signature --corpus /tmp/st.npz samples/*.tsv.gz -o sig.tsv

from vdjtools.signature import vsig, vsig_cohort
from vdjtools.signature.corpus import Corpus, bundled_path

corpus = Corpus.load(bundled_path("synthetic-blood"))
row = vsig({"TRB": trb, "IGH": igh}, corpus)
frame = vsig_cohort({"S1": sample1, "S2": sample2}, corpus, n_jobs=0)

Pick synthetic-blood or synthetic-tissue by which compartment your samples come from; both are described below, with what they reproduce and what they do not.

Choosing a corpus

This is the first decision and the one that matters most, so start here. Match the corpus to how your samples were produced, not to what you are testing for.

Your samples are	Use	Why
Bulk RNA-seq of blood or PBMC	blood	Fitted on 11,117 real blood samples across 947 study groups. The default choice for peripheral-blood work.
Bulk RNA-seq of a solid tissue or tumour	tissue	Fitted on 21,131 real non-blood samples across 1,934 study groups. Tissue repertoires are shallower and more skewed than blood, and a blood corpus clamps them harder.
Targeted or amplicon deep TCR sequencing	deep-tcr	Fitted on 3,936 amplicon samples. TRA and TRB only – it has no other locus to offer.
Any of the above, and you need the per-study cap’s effect	blood-uncapped, tissue-uncapped	The same populations with no cap at 30 samples per study, so the difference is measurable rather than assumed.
Something none of these describe	synthetic-blood, synthetic-tissue	Drawn across the measured depth and clone-size range of that compartment, using no cohort sample at all – so they rebuild bit-for-bit anywhere and carry no population you have to accept.
A reference point rather than a cohort	naive, memory	The pure regimes: what the recombination model emits, and the same under a Zipf rank-abundance clone-size law.

Two practical consequences:

A real corpus fits real data better, and by a measurable margin. On four TCGA tissue samples at --components 32, the fraction of the row that the bounds clamped was 0.199 through tissue against 0.595 through synthetic-tissue – roughly three times less clamping. deep-tcr on the same samples clamps 0.483 with a median absolute z of 2.1, correctly far out: amplicon TCR is the wrong reference for RNA-seq, and the numbers say so rather than quietly absorbing it.

Check ``winsor_frac`` before you trust a matrix. vsig:qc:-:winsor_frac is the fraction of the row the corpus’s bounds clamped. A high value means your samples sit outside the corpus’s range, so the standardisation is extrapolating – pick a closer corpus rather than proceeding.

A corpus is required

There is no default. A signature is comparable to another one only if both were rotated through the same corpus, so the artifact has to be named, and the emitted row records which one it was. Two matrices whose column names match but whose corpora differ are not comparable, and nothing about the numbers says so — which is why the choice is not allowed to be implicit.

All nine corpora are published. Four are synthetic – every receptor drawn from the bundled recombination models, so anyone can rebuild the artifact bit-for-bit and no cohort sample is needed to use one – and five are real, fitted on repertoires. Both kinds exist on purpose: a synthetic corpus spans a measured range by construction, while a real one carries the joint structure the generative model does not produce (selection-shaped V/J usage, isotype and SHM structure, and the cross-locus covariance of libraries prepared together).

corpus	kind	fitted on	what it is
naive	synthetic	10,000 drawn	every clone size 1 – what the recombination model emits
memory	synthetic	10,000 drawn	Zipf rank-abundance, sampled by multinomial, zeros dropped
synthetic-blood	synthetic	10,000 drawn	a naive/memory mixture across measured blood ladders (TRB richness 74-3,162, f1 0.215-0.949)
synthetic-tissue	synthetic	10,000 drawn	the same across measured tissue ladders (TRB richness 30-2,977, IGH 40-10,352)
blood	real	11,117 samples	public bulk RNA-seq blood, 947 study groups, capped at 30 per study
blood-uncapped	real	22,441 samples	the same population with no cap, so the cap’s effect is measurable
tissue	real	21,131 samples	public bulk RNA-seq non-blood, 1,934 study groups, capped at 30 per study
tissue-uncapped	real	33,874 samples	the same population with no cap
deep-tcr	real	3,936 samples	targeted/amplicon deep TCR, TRA and TRB only

Artifacts are fetched, not bundled

At 256 components a vsig artifact is about 10 MB of float32 rotation, and the nine corpora across both halves are roughly 110 MB. Making every pip install carry all of that in order to use one of them is the wrong trade, so the wheel ships a few KB of index – corpora.json, naming every corpus with its size and SHA-256 – and the artifact itself is a GitHub release asset fetched on first use into $VDJTOOLS_CORPUS_DIR (default ~/.cache/vdjtools/signature) and verified against that digest.

The resolution order is a local path, then the wheel, then the cache, then the release. A local path always wins, so a caller who fitted their own corpus and passes its path is never served a download of the same name. A download whose digest does not match is deleted rather than cached: a half-written rotation that loads is worse than one that is missing.

vdjtools corpus --fetch all (and mir corpus --fetch all) pre-warms the cache – the install-time step, since anything not pre-warmed is fetched lazily. The two halves are published by the library that owns each and share one cache directory, because a vsig/rsig pair belongs in one place.

Which one? Use blood or tissue when your samples are bulk RNA-seq of that compartment and deep-tcr for amplicon TCR: a real corpus is the closest reference to real data. Use synthetic-blood / synthetic-tissue when you need a reference that is reproducible from the library alone, with no cohort behind it, or when your depths sit outside what the real corpora cover. naive and memory are the two pure regimes and neither is a bulk sample – naive has no clone-size structure at all, and memory at a fixed nominal size has a read count of exactly size × reads_per_clone in every sample – so they are reference points for what selection does to a rotation, not descriptions of a cohort.

What comes out

Column names are four colon-separated parts throughout, so one parse() reads every kind:

vsig:pc:TRB:PC01          a rotated coordinate
vsig:cov:TRB:cstar        a channel, carried through untouched
vsig:qc:-:winsor_frac     a channel that is not per-locus

Channels are never rotated and never clamped. A provenance number that has been mixed with the measurements it was supposed to qualify is no longer provenance. Read these three before trusting any rotated column of a row:

channel	what it tells you
vsig:cov:<locus>:cstar	The coverage this sample actually attained, in [0, 1]. Always emitted, including when the coverage-standardised diversity features turn out not to be estimable — which is the case it matters most in.
vsig:mask:<locus>:present / :estimable	Why a column is a hole. “This locus is absent” and “this sample is too shallow to reach the target” are different facts that would otherwise both render as one nan.
vsig:qc:-:winsor_frac	What fraction of this row the corpus’s bounds clamped. See Winsorization: by percentile, and one-sided where the metric is.

vsig:qc:<locus>:v_fallback_frac and :j_fallback_frac are the other two worth reading: an unrecognised V or J call raises nowhere, because the germline lookup falls back to the maximum observed distance, so a cohort on Adaptive nomenclature or an older IMGT release yields a fully populated, entirely plausible, systematically wrong result. These fractions are the only thing that says so.

A hole is nan, never 0.

Winsorization: by percentile, and one-sided where the metric is

Bounds are estimated once, on the corpus, at a percentile — not as a multiple of a robust standard deviation. A robust-SD bound is not a fixed amount of probability: it depends on the distribution’s shape, so an 8-SD bound trims nothing at all on a right-skewed count and a third of the corpus on a near-degenerate column.

The side that gets trimmed comes from each feature’s declared support, and never from its transform:

support	trimmed	because
nonneg	top only	[0, inf) — cannot run away downward. Counts, richness, Hill numbers, standard deviations.
nonpos	bottom only	(-inf, 0] — a log-probability is bounded above by 0.
real	both	(-inf, inf) — log-ratios, clr and logit coordinates, PC scores.
unit	neither	[a, b] — already bounded on both sides. Raw proportions, cstar.

Three raw features declare transform="none" and have three different supports. Deriving the side from the transform silently trims the wrong end of two of them, which is why support is a declared field rather than an inference.

Two modes, and the mode fixes how the rotation was fitted

--winsorize	at emit time	the corpus PCA was fitted by
features	clamp raw features to the corpus bounds, then rotate	winsorized corpus, median/MAD, plain SVD. Deterministic and bit-reproducible; the robustness sits in a step you can inspect.
pcs	rotate, then clamp the PC scores	RobustScaler and a full-solver PCA — there the corpus still has its tails at fit time, so the estimator has to resist them itself.
none	clamp nothing	whichever of the two the artifact records

Both bound sets and both percentiles (0.01 and 0.05) are stored, so --winsorize and --winsor-p are apply-time choices that need no refit.

The clamp is never silent. Clamping is many-to-one: it overwrites a real measurement of a real sample, and the sample cannot tell afterwards. In the system this replaces, a reference whose centre was wrong by 57 robust deviations went unnoticed because every affected sample was clamped to the bound and looked like an ordinary tail case. So every row carries vsig:qc:-:winsor_frac, and a value near 1.0 means this sample does not belong to this corpus — not that it is unusual.

How many components

--components / n_components= takes either kind of answer:

• an integer — a fixed count, identical at every locus. At fit time it is a cap: a block narrower than the request (the five cross-locus ratios against a request of 256) gets its full rotation, recorded per locus in the artifact.

• a fraction in (0, 1) — cumulative variance explained. The count then differs per locus and per corpus, because some corpora are genuinely easier to explain. That is biology, not a defect.

The artifact stores the rotation only up to the count it was fitted at, plus the full eigenvalue spectrum. So truncating downward at apply time is free and exact — the components are ordered, so the first n of a longer rotation are the same vectors. A count above the stored one is capped per locus, exactly as the fit capped it; a fraction the stored spectrum cannot reach raises, and the error quotes what it does reach. Neither ever pads: no column appears under a name the rotation does not hold.

The coverage target is a runtime argument

Hill numbers are read at a coverage level, and that level is not in the artifact. A per-sample quantity has no business in a corpus artifact: in the system this replaces it lived there, and one shipped reference ended up standardising tissue samples to a level measured on blood, identical to 17 significant digits across all seven loci while its location and scale had genuinely been refitted.

--cstar-target	meaning
min (default)	This cohort’s own per-locus minimum attained coverage. Comparable across its samples, and a property of this cohort rather than of somebody’s corpus. Costs a cheap extra pass over the count column — and resolves deferred samples twice, so pass a number for a one-pass run.
a number	That level for every locus. One pass.
own	Each sample at its own attained coverage: the diversity features are then observed rather than standardised, and are not comparable across samples.

Whatever the target, vsig:cov:<locus>:cstar reports what the sample reached and vsig:mask:<locus>:estimable reports whether the target was reachable without extrapolating.

Raw feature groups

Every group at a locus is concatenated into one feature vector and rotated together, so the PCA can find cross-group structure. One rotation per (half, locus), plus one for the cross-locus group — not one joint rotation over all loci, because under a joint rotation one absent locus silently moves every component. Per-locus blocking is what keeps a hole a hole.

group	width	what
div	6	Coverage-standardised Hill numbers 0D_c/1D_c/2D_c, clonality, 0D_chao, d50
depth	3	reads, richness, S_unseen (Chao’s estimate of what was never drawn)
clon	3	Clone-size composition f1/f2 in clr coordinates, plus the top clone’s share
len	3	Weighted junction-length mean, sd, skew
vus / jus	germline	V and J gene usage over the arda germline vocabulary, clr
spec	germline x 26	Spectratype: junction length resolved per V gene, clr
kmer	400	Junction 2-mer composition, arcsine
aa	20	Single-residue composition, arcsine
pchem	30	Physicochemistry over two junction regions
iso / shm	5 / 1	IGH only: isotype composition in clr, and mean V identity
pair	5	Cross-locus: five log read-count ratios, depth divided out

Widths per locus run from 685 (TRD) to 3,744 (IGH), 13,483 raw features in total.

Note

spec is by far the widest group and the sparsest. 121 IGH V genes x 26 lengths is 3,146 cells against a median 464 IGH clonotypes per real sample — about 0.15 observations per cell, mostly structural zeros. That sparsity is real and is not hidden: a corpus of 10,000 repertoires puts ~4.6M observations behind the same cells, so the rotation’s directions are well determined even where one sample’s scores along them are noisy. How noisy is readable per sample from vsig:cov:<locus>:cstar.

columns= selects output, it does not skip work

Because a locus is rotated jointly over all its groups, every group must be computed before any of that locus’s components exist. columns= therefore restricts the emitted frame; it does not avoid the computation. Declining a whole locus does skip it, and is the way to make a run cheaper.

Building a corpus

vdjtools corpus --corpus synthetic-blood -o synthetic-blood.npz       # 10,000 repertoires
vdjtools corpus --corpus naive  --out naive.npz
vdjtools corpus --corpus memory --size n_eff --components 0.95 -o m.npz
vdjtools corpus --smoke -o /tmp/smoke.npz                             # minutes, for tests

The synthetic corpora use no samples from anybody’s cohort: every receptor is drawn from the bundled recombination models.

The build runs across worker processes, one contiguous block of samples each; --jobs/-j sets how many and defaults to every core the process may use (-j 1 stays in-process). A sample is a pure function of its index – its generator is seeded from (seed, locus, index) – so each worker draws its own repertoires from memory-mapped pools and returns one row of numbers, and the artifact is identical at every --jobs. Measured on Aldan-3, seven loci at --samples 10000 --size 10000 on 72 cores: 157 s to generate the receptor pools (4,850 s in one process) and about 50 samples/s to featurise them.

Note

Two settings, two layers. --jobs is worker processes; POLARS_MAX_THREADS and OMP_NUM_THREADS are kernel threads. The builder sets its workers to one kernel thread each, because jobs x cores threads is how a 15 s stage was once measured at 296 s.

The build is a deterministic function of (corpus, loci, samples, size, seed, source) and those models, and is required to be byte-identical across processes, worker counts and thread counts:

vdjtools corpus --corpus naive --loci TRB,TRG --samples 40 --size 150 -o /tmp/a.npz
OMP_NUM_THREADS=1 POLARS_MAX_THREADS=1 \
  vdjtools corpus --corpus naive --loci TRB,TRG --samples 40 --size 150 -o /tmp/b.npz
vdjtools corpus --corpus naive --loci TRB,TRG --samples 40 --size 150 -j 4 -o /tmp/c.npz
cmp /tmp/a.npz /tmp/b.npz && cmp /tmp/a.npz /tmp/c.npz     # all three must be identical

Important

Each repertoire’s depth is drawn, log-uniformly over the measured real p05–p95 spread for its locus, not fixed. A corpus at one fixed depth has exactly zero variance in depth:reads, depth:richness and all five pair: ratios, so those contribute nothing to the rotation and the cross-locus block cannot be fitted at all — and it fails quietly, as an artifact that simply omits them. Real repertoires span 4.1x (TRB) to 11.0x (IGH) between their 5th and 95th percentiles.

The artifact verifies itself on load: the raw column order is re-derived from the stored vocabulary through the current layout and compared, and the recorded bundled-model version is checked against the installed one. Both raise rather than producing numbers in a different coordinate system that look perfectly reasonable — retraining the bundled models moves the synthetic corpora by 0.1–1.6% per locus.

Parallelism

--jobs/n_jobs is worker processes. There is exactly one concurrency setting, and its help text says which layer it reaches; a flag named --threads that started processes is a mistake this subsystem has made before. Workers are spawned, not forked — polars cannot be combined with fork — and a pool that cannot start raises rather than falling back to one process, because a correctness-preserving fallback is what hid a 20x slowdown here once.

Getting the named statistics back

A signature is a rotation, and a rotated coordinate has no name a domain reader can use. But the rotation is fitted on features that do: diversity, depth, clone-size fractions, junction length, isotype composition, SHM, cross-locus yield. Those are computed for every sample either way.

named= returns them:

Python

from vdjtools.signature import vsig, vsig_cohort
from vdjtools.signature.corpus import Corpus, bundled_path

corpus = Corpus.load(bundled_path("blood"))

vsig(sample, corpus, n_components=32, named=True)
vsig(sample, corpus, n_components=32, named=("div", "depth"))   # or pick blocks
vsig_cohort(samples, corpus, n_components=32, named=True)

Command line

vdjtools signature --corpus blood --components 32 --named all -o sig.tsv
vdjtools signature --corpus blood --components 32 --named div,depth -o sig.tsv

named=(), the default, emits exactly the rotated columns and the channels – the output is byte-identical to a run without the argument. The reportable blocks are:

half	block	what it holds
vsig	div	Coverage-standardised Hill numbers 1D_c/0D_c/2D_c, 0D_chao, clonality, d50
vsig	depth	reads, richness, S_unseen
vsig	clon	Singleton, doubleton and top-clone fractions
vsig	len	Junction-length mean, sd, skew
vsig	iso / shm	IGH isotype composition, and mean V identity
vsig	pair	Cross-locus yield log-ratios
rsig	depth / band / band_igh	n_eff and observed mass; clone-size and isotype shares of Phi

Warning

These carry their declared transform, not their natural scale. A log10 diversity comes back as log10 and a clr composition as a log-ratio with no unique inverse. Do not read vsig:div:TRB:1D_c = 2.386 as a clone count – it is log10 of one. channel_table() reports the transform for every feature, and --describe prints it per emitted column.

The reason this exists rather than being a convenience: the diversity floor of any study using this library is made of exactly these blocks. Without a documented route to them there is no floor, and a signature that beats nothing gets reported as if it beat something.

Fitting a corpus on your own cohort

The published corpora are one route; fitting the rotation on your own training half is the other, and it is the first thing anybody with a cohort asks for.

from vdjtools.signature.corpus import Corpus, fit_cohort

corpus = fit_cohort(train_samples, sig="vsig", name="my-cohort",
                    loci=("TRA", "TRB", "IGH"), n_components=64, n_jobs=0)
corpus.save("my-cohort.npz")

# and it is the same artifact type the shipped corpora are
mine = Corpus.load("my-cohort.npz")
vsig(held_out_sample, mine, n_components=64)

Samples are {sample_id: {locus: frame}}, or an iterable of frames, or picklable zero-argument callables that defer the read into the worker. mir.signature must be imported before sig="rsig" will resolve; nothing in vdjtools imports mir, so the featuriser reaches the corpus through a registry rather than an import.

Warning

A cross-validated score computed on a rotation fitted inside the same cohort is not evidence that the rotation generalises.

Measured on raw repertoire embeddings, 5,376 columns over seven loci, one logistic head, with a balanced train/test split of one cohort and a second cohort held out entirely. At matched widths:

read-out	in-cohort leads at	median, in-cohort	median, corpus
training out-of-fold (what selects)	4 of 5 widths	0.5039	0.4968
held-out half	2 of 5	0.5222	0.5635
external cohort	0 of 5	0.4757	0.5270

Fitting the rotation in-cohort improves the number the configuration is selected on, and neither held-out read-out. A rotation fitted on 612 samples of one trial learns that trial’s covariance, which is what a within-cohort cross-validation rewards and what does not travel.

Both routes are legitimate and they answer different questions. Fit here when you want a basis for this cohort; name a shipped corpus when the number has to travel.

Further reading

How signatures were fitted is the measurement record: how each corpus’s population was selected, what the quantile ladders are, how far a corpus carries across sequencing depths, and which alternatives were tried and rejected. Channels documents the columns that are never rotated and never clamped.

API

Repertoire signatures: raw features, a corpus-fitted rotation, and the emitted vector.

Four modules, one per stage:

• layout – the column contract. Which raw features exist, what each one’s support is, which columns pass through as channels. Computes nothing, imports nothing heavy.

• transform – the variance-stabilising transforms, applied where a raw feature is computed so its denominator is still in scope.

• features – the raw features and channels for one sample. A pure function of that sample plus the germline vocabulary.

• corpus – the only module that needs a corpus: build one, winsorize it, fit the rotation and the per-PC scaling, write and read the artifact.

vsig() and vsig_cohort() in signature put the four together. The geometry half (rsig) lives in mir.signature and registers its groups into this same layout registry; nothing here imports mir.

class vdjtools.signature.Corpus(sig, name, vocab, fits, meta=<factory>)

A fitted corpus: one LocusFit per locus, plus the manifest.

A load verifies rather than trusts: the raw column order is re-derived from the stored vocabulary through the current layout and compared, so a layout change that would silently re-index the rotation raises instead. The bundled-model version is checked the same way, because the synthetic corpora are drawn from those models and retraining them moves the result by 0.1-1.6% per locus.

Parameters:

• sig (str)

• name (str)

• vocab (dict[str, dict[str, list[str]]])

• fits (dict[str, LocusFit])

• meta (dict)

property k: dict[str, int]

{locus: components}. A locus absent here emits no columns – how a corpus that could not fit a locus declares it, rather than shipping a rotation of zeros.

columns(n_components=None, named=())

The signature columns this corpus emits: rotated, then channels, then named blocks.

Parameters:

• n_components (int | float | None)

• named (bool | Sequence[str])

Return type:

list[str]

resolve_k(n_components=None)

{locus: k} after apply-time truncation.

Parameters:

n_components (int | float | None) – None keeps what was fitted. An int truncates every locus to that count. A float in (0, 1) truncates each locus to the fewest components reaching that cumulative variance.

Raises:

ValueError – If an int exceeds what a locus was fitted with, or a float target is not reached by the stored spectrum. Neither ever pads: a count is capped at each locus’s stored k (the same cap the fit applied, so a 5-feature cross-locus block yields 5 whether you asked for 5 or 128), and a variance fraction the stored spectrum cannot reach raises rather than silently returning the whole rotation under the requested name.

Return type:

dict[str, int]

save(path)

Write <path>.npz (arrays) and <path>.json (manifest and vocabulary).

Parameters:

path (str | Path)

Return type:

Path

classmethod load(path, *, verify=True)

Read an artifact, and check it is the one the caller thinks it is.

Parameters:

• path (str | Path) – The .npz path; the .json sidecar must sit beside it.

• verify (bool) – Re-derive the raw column order from the stored vocabulary through the current layout and compare, and check the recorded bundled-model version against the installed one. Turn it off only to inspect a deliberately stale artifact.

Raises:

ValueError – On a column-order or model-version mismatch. Both would otherwise produce numbers in a different coordinate system that look perfectly reasonable.

Return type:

Corpus

verify()

Raise unless the stored column order and model version match this installation.

Return type:

None

vdjtools.signature.fit(rows, vocab, *, sig, **kw)

Fit a corpus from raw feature rows – one dict per repertoire.

Convenient for a small corpus and for tests. For a full-size build use locus_matrix() / fill_row() / fit_matrices(), which never holds more than one sample’s dict at a time. To fit on your own cohort of clonotype frames rather than on feature rows you assembled yourself, use fit_cohort().

Parameters:

• rows (list[dict[str, float]])

• vocab (dict[str, dict[str, list[str]]])

• sig (str)

Return type:

Corpus

vdjtools.signature.raw_and_channels(frames, vocab, *, cstar_target=None, weight='log2p1', prefiltered=False, on_duplicate='error')

Every raw feature and every channel for one sample.

Parameters:

• frames (dict[str, DataFrame]) – {locus: clonotype frame}, raw (this sanitises).

• vocab (dict[str, dict[str, list[str]]]) – {locus: {"vus": [...], "jus": [...], "spec": [...]}} from the corpus artifact, or from gene_vocab() at fit time.

• cstar_target (float | dict[str, float] | None) – Coverage level the Hill numbers are read at. A float applies to every locus, a dict is per-locus, and None uses each locus’s own attained coverage – which makes the diversity features observed rather than standardised, and is the only default that does not smuggle in a constant from somewhere the caller cannot see.

• weight (str) – Clone-size weight, a key of WEIGHTS.

• prefiltered (bool) – When the caller already removed non-productive rows, report nonstd_aa_frac as nan rather than a confident floor of zero.

• on_duplicate (str) – Forwarded to sanitise().

Returns:

(raw, channels) – both {column: value}, holes as nan.

Return type:

tuple[dict[str, float], dict[str, float]]

vdjtools.signature.synthesize(corpus_name, *, loci=('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL'), n_samples=10000, size=None, seed=20260927, n_components=128, mode='features', winsor_p=0.01, source='olga', fit_corpus=True, progress=None, n_jobs=1, depth_spread=None)

Build a synthetic corpus, and fit it.

Parameters:

• corpus_name (str) – One of SYNTHETIC – "naive" and "memory" are the pure regimes at a nominal size; "synthetic-blood" and "synthetic-tissue" are the naive/memory mixture drawn across the named cohort’s measured per-locus richness, read-depth and singleton-fraction bands (COHORT), and are the two that describe a real bulk cohort rather than a regime.

• loci (tuple[str, ...]) – Loci to build. All seven by default.

• n_samples (int) – Repertoires in the corpus.

• size (int | str | None) – Receptors per repertoire. None, the default, is the corpus’s own – 10,000 for the pure regimes, the geometric centre of the cohort’s measured richness band for a synthetic-* one. Also takes "n_eff" for the per-locus real medians in N_EFF, "p05" / "p95" for the sweep endpoints, or a cohort name.

• seed (int) – Base seed; every draw is a recorded offset from it.

• n_components (int | float) – Components per locus, or a variance fraction.

• mode (str) – Winsorization mode to fit at.

• winsor_p (float) – Percentile for the fitted bounds. All of WINSOR_PS are stored.

• source (str) – Bundled model set – "olga", "learned" or "arda".

• fit_corpus (bool) – False returns the raw samples instead of fitting, for scoring a held-out draw against a corpus fitted on another.

• progress – Optional callable(locus, done, total).

• n_jobs (int) – Worker processes (not kernel threads). 0 means every available core; 1, the default, runs in-process. The result is bit-identical at every value.

• depth_spread (float | str | None) – Multiplicative depth range each repertoire’s size is drawn log-uniformly across, around size. None is the corpus’s own – DEPTH_SPREAD for a pure regime (2.4x-11.0x), the cohort’s measured richness band for a synthetic-* one (43x on blood TRB, 259x on tissue IGH). Pass a number to override it: 1000 with size=3162 spans 100 to 100,000 receptors per locus. The bounds, centre and per-PC scaling are all estimated from the draw, so a corpus describes only the depths it was drawn across.

Returns:

(corpus, mats) when fitting – mats being {locus: (matrix, columns)}, the corpus matrix itself – else a list of {locus: frame} samples.

The whole build is a deterministic function of (corpus_name, loci, n_samples, size, seed, source) and the bundled models, so two machines at different thread counts must produce byte-identical artifacts. That is an acceptance criterion, not a hope: generate was once irreproducible across processes while being deterministic within one, because a marginal-table aggregation left its group order unspecified.

vdjtools.signature.vsig(sample, corpus, *, mode=None, winsor_p=None, n_components=None, cstar_target=None, weight='log2p1', prefiltered=False, on_duplicate='error', named=(), columns=None)

The statistics half of the signature for one sample.

Parameters:

• sample – {locus: frame}, one frame with a locus column, or a zero-argument callable returning either.

• corpus (Corpus) – A fitted Corpus for sig="vsig".

• mode (str | None) – Winsorization mode override – "features", "pcs" or "none". Defaults to whichever the corpus was fitted with.

• winsor_p (float | None) – Which stored percentile to clamp at. Defaults to the fitted one.

• n_components (int | float | None) – Truncate to this many components, or this cumulative variance fraction.

• cstar_target (float | dict[str, float] | None) – Coverage level the Hill numbers are read at. None uses each locus’s own attained coverage, which makes them observed rather than standardised – honest for one sample, and not comparable across samples. Use vsig_cohort() for a cohort, which resolves a shared level.

• weight (str) – Clone-size weight; a key of WEIGHTS.

• prefiltered (bool) – Report qc:*:nonstd_aa_frac as nan because the caller already filtered.

• on_duplicate (str) – "error" or "sum", for a frame repeating an amino-acid clonotype key.

• named (bool | Sequence[str]) –

  Also return the reportable raw blocks – True for all of them, or a sequence such as ("div", "depth", "clon"). (), the default, emits exactly the rotated columns and channels. Values carry their declared transform (log10, clr, logit, …), which channel_table() reports.

  These are the diversity, depth, clone-size, junction-length, isotype, SHM and cross-locus yield numbers. They are computed either way, because the rotation is fitted on them; without this argument there is no supported route to reading them back, and a study that cannot report its own diversity floor has no floor. See channel_table() for the full list and its units.

• columns (list[str] | None) –

  Restrict the output to these columns, in layout order.

  Unlike the system this replaces, this does not skip work: the rotation for a locus is fitted jointly over every raw group at that locus, so every group must be computed before any of that locus’s components exist. Declining a locus entirely does skip it.

Returns:

{column: value} in layout order, holes as nan.

Return type:

dict[str, float]

vdjtools.signature.vsig_cohort(samples, corpus, *, n_jobs=1, cstar_target='min', columns=None, **kw)

One row per sample, sample_id first.

Parameters:

• samples – {sample_id: sample} or an iterable of pairs. A sample may be a zero-argument picklable callable, which defers the read into the worker and keeps peak memory at O(n_jobs) samples rather than the whole cohort.

• corpus (Corpus) – A fitted corpus.

• n_jobs (int) – Worker processes. 1 stays in-process, 0 uses every available core.

• cstar_target (float | dict[str, float] | str | None) – "min" (default) reads every sample’s attained coverage first and standardises the whole cohort to the per-locus minimum, so the diversity columns are comparable across samples and the level is a property of this cohort rather than of somebody’s corpus. That is a cheap extra pass over the count column – but it does resolve deferred samples twice, so pass an explicit float or dict for a one-pass run. None reads each sample at its own coverage: one pass, not comparable across samples.

• columns (list[str] | None) – Restrict the output columns.

• **kw – Forwarded to vsig() – including named=, which adds the reportable raw blocks in their own units.

Returns:

A frame whose columns are sample_id then the corpus’s signature columns.

Return type:

DataFrame

class vdjtools.signature.Channel(sig, name, features, loci=None)

A named family carried through untouched – never winsorized, never rotated.

Parameters:

• sig (str) – Owning signature.

• name (str) – Channel name.

• features (dict[str, str]) – {name: support}. The support is recorded for documentation and for a range assertion; no trimming is applied either way.

• loci (tuple[str, ...] | None) – As RawGroup.

class vdjtools.signature.RawGroup(sig, name, features=<factory>, loci=None, dynamic=None, named=False)

One named family of raw features – an input to the rotation, never an output.

Parameters:

• sig (str) – Owning signature, "vsig" or "rsig".

• name (str) – Group name, unique within a sig.

• features (dict[str, tuple[str, str]]) – {name: (transform, support)}, most easily built with feats(). Empty for a group whose column names are not knowable without the germline; such a group declares dynamic instead.

• loci (tuple[str, ...] | None) – Loci the group is emitted for. None means all of LOCI; an empty tuple means the group is not per-locus and uses NO_LOCUS.

• dynamic (tuple[str, str] | None) – For a group whose width depends on the germline vocabulary (vus, jus, spec): one (transform, support) pair shared by every column, with the names resolved at fit time and then recorded in the corpus artifact. Keeping the vocabulary out of here is what lets import vdjtools stay free of arda.

• named (bool) – This group’s features are interpretable domain scalars – a diversity index, a read count, a singleton fraction – that a caller may legitimately want back in their own units rather than only as a rotation input. Declared here rather than inferred, because “is this number reportable” is a property of the feature and not of its width: pchem is 30 static columns and is not reportable, shm is one and is. Groups marked here are what named=True selects; see named_groups().

property emitted_loci: tuple[str, ...]

Loci this group emits for; (NO_LOCUS,) when it is not per-locus.

spec(feature)

(transform, support) for one of this group’s features.

Parameters:

feature (str)

Return type:

tuple[str, str]

columns(locus, names=None)

Raw column names for one locus.

Parameters:

• locus (str) – A member of emitted_loci.

• names (list[str] | None) – Required for a dynamic group – the germline-resolved column names, in the order the artifact recorded them.

Return type:

list[str]

vdjtools.signature.arcsine(x, m)

Anscombe’s variance-stabilising arcsine transform, asin(sqrt((x·m + 3/8)/(m + 3/4))).

The right transform for a sparse composition — residue and gene-usage profiles, where most cells are structurally zero at shallow depth. Unlike a CLR it is defined at zero without any replacement step, and unlike a raw proportion its variance does not collapse near the boundary. Bounded in [0, π/2], so it cannot produce the heavy tail a log-ratio would.

Parameters:

• x – Proportion(s) in [0, 1].

• m – Denominator(s) the proportion was observed on.

vdjtools.signature.channel_columns(sig)

Every pass-through channel column for one sig, in emitted order.

Parameters:

sig (str)

Return type:

list[str]

vdjtools.signature.channels(sig=None)

Registered pass-through channels, in declaration order.

Parameters:

sig (str | None)

Return type:

list[Channel]

vdjtools.signature.clr(parts, m=None, *, keys=None)

Centred log-ratio of a composition, with multiplicative zero replacement.

A CLR is the natural coordinate for a composition whose ratios carry the meaning: it is the log of each part over the geometric mean of all of them, so it is invariant to the total and a difference between two coordinates is a log-ratio of two parts.

Zeros are replaced multiplicatively, not additively: each zero part is set to delta = 0.5/m and the non-zero parts are scaled down by 1 − n_zero·delta so the composition still closes. Adding a constant to every part instead — the common shortcut — distorts the ratios among the parts that were observed, which are the only ratios the coordinate system is about.

Compute the CLR over the whole composition, then select coordinates. A CLR of a sub-composition is a different number from the corresponding coordinate of the full one, so a tier that ships four of six isotype parts must still divide by the six-part geometric mean. Doing it the other way would make the narrower tier stop being a slice of the wider one, which is the contract the layout exists to guarantee.

Parameters:

• parts – Non-negative part values as a mapping {name: value} or an array. They need not sum to 1; they are closed here.

• m – Total count the composition was observed on, setting the replacement scale. Defaults to the sum of parts when they are counts.

• keys – Part order when parts is an array. Ignored for a mapping.

Returns:

{name: clr} when parts is a mapping (or keys is given), else an array. The coordinates sum to zero by construction, which is why the layout ships all but one of them: the last is exactly determined by the rest and would make any unregularised design matrix singular.

Raises:

ValueError – If fewer than two parts are given, or any part is negative.

vdjtools.signature.feats(transform, support, *names)

{name: (transform, support)} for a run of features sharing both.

Groups are heterogeneous – a diversity group carries log10 Hill numbers beside a logit evenness – so both properties belong to the feature, not the group. This keeps the homogeneous runs terse anyway.

Parameters:

• transform (str)

• support (str)

• names (str)

Return type:

dict[str, tuple[str, str]]

vdjtools.signature.gene_vocab(locus, organism='human')

Column names for the three germline-width raw groups at one locus.

arda’s germline is the single source of truth for V/J vocabulary across this ecosystem, so the usage and spectratype widths are a property of the germline release rather than of any cohort. The corpus artifact records what came back here, and apply time reads it from the artifact – never from arda again, because a germline update would silently re-index the rotation.

Returns:

{"vus": [...], "jus": [...], "spec": [...]}. spec is <V gene>_<length> over SPECTRATYPE_LENGTHS, V-major.

Parameters:

• locus (str)

• organism (str)

Return type:

dict[str, list[str]]

vdjtools.signature.log10(x, floor=1.0)

log10 of a positive quantity, floored so an empty locus maps to 0 rather than -inf.

Used for counts and for Hill numbers. The floor is 1 because both are counts of things: one clonotype, one read, one effective species. Zero of them is the same as the floor for every downstream purpose, and the presence mask already records that the locus was empty.

Parameters:

floor (float)

vdjtools.signature.log1p(x)

log(1+x) for a non-negative quantity that genuinely reaches zero.

Distinct from log10() in intent: this is for norms, dispersions and hit counts, where zero is a real, attainable value rather than an empty measurement.

vdjtools.signature.logit(x, m)

Haldane–Anscombe logit of a proportion observed on a denominator m.

log((x·m + 1/2) / ((1−x)·m + 1/2)). Adding half an observation to each side is the standard remedy for an empty cell; its side effect is exactly the behaviour wanted here — the transform of 0 depends on how many chances there were to see something:

logit(0, m=3)   ->  -1.95      a fifth of the repertoire could hide here
logit(0, m=500) ->  -6.91      it is really absent

Parameters:

• x – Proportion(s) in [0, 1].

• m – Denominator(s) the proportion was observed on. Broadcasts against x.

Returns:

The transformed value, finite for every input including exactly 0 and exactly 1.

vdjtools.signature.parse(column)

Split a column name into (sig, block, locus, feature).

Also the definition of “is this a signature column”, used as an allow-list wherever a frame may carry a caller’s own joined columns – so an age or n_reads column can never be silently winsorized or rotated.

Raises:

ValueError – If the name is not four colon-separated parts.

Parameters:

column (str)

Return type:

tuple[str, str, str, str]

vdjtools.signature.pc_columns(sig, k)

Rotated column names, given the component count per locus.

Parameters:

• sig (str) – "vsig" or "rsig".

• k (dict[str, int]) – {locus: n_components}, from the corpus artifact. A locus absent from k emits nothing, which is how a corpus that could not fit a locus declares it.

Return type:

list[str]

vdjtools.signature.raw_columns(sig, locus, vocab=None)

Every raw feature name for one (sig, locus), in rotation-input order.

This is the column order of one row of the corpus matrix M_L, so it is also the row order of the rotation. It must be reproduced exactly at apply time, which is why the corpus artifact stores it verbatim rather than recomputing it.

Parameters:

• sig (str) – "vsig" or "rsig".

• locus (str) – A locus, or NO_LOCUS for the cross-locus group.

• vocab (dict[str, list[str]] | None) – {group_name: [column names]} for every dynamic group emitted at this locus. Required unless no dynamic group applies.

Return type:

list[str]

vdjtools.signature.raw_groups(sig=None)

Registered raw groups, in declaration order.

Parameters:

sig (str | None)

Return type:

list[RawGroup]

vdjtools.signature.register_channel(*channels)

Add pass-through channels to the registry.

Parameters:

channels (Channel)

Return type:

None

vdjtools.signature.register_raw(*groups)

Add raw groups to the registry. Used by mir.signature for the geometry half.

Parameters:

groups (RawGroup)

Return type:

None

vdjtools.signature.sanitise(df, *, strict=True, on_duplicate='error')

Drop non-productive clonotypes; return the frame and the dropped weight fraction.

Dropped by weight, not by row: losing one dominant clone matters more than losing fifty singletons, and a row fraction would hide that.

Parameters:

• df (DataFrame) – A clonotype frame.

• strict (bool) – Raise on unparseable characters (see assert_parseable()). False drops them, for a corpus known to carry ambiguity codes.

• on_duplicate (str) – What to do when the frame has no junction_nt and repeats an amino-acid clonotype key – "error" (default) or "sum". The default refuses because the two readings give different richness and clonality and the frame cannot say which is meant.

Returns:

(kept_frame, dropped_weight_fraction).

Return type:

tuple[DataFrame, float]

vdjtools.signature.signature_columns(sig, k)

The emitted signature for one half: rotated columns then channels, in that order.

Parameters:

• sig (str)

• k (dict[str, int])

Return type:

list[str]

vdjtools.signature.support_of(column, vocab=None)

The declared support of a raw feature or channel column.

Parameters:

• column (str) – A raw or channel column name.

• vocab (dict[str, dict[str, list[str]]] | None) – Unused for lookup – a dynamic group’s support is declared on the group, not per column – and accepted so callers can pass their artifact’s vocab uniformly.

Raises:

ValueError – If no registered group or channel declares the column.

Return type:

str

vdjtools.signature.work_frame(df, weight='log2p1')

Overwrite frequency with the normalised clone weight, so every profiler agrees.

Order matters. Call this after filtering, and never call filter_functional, downsample or select_top afterwards – each recomputes frequency from the counts and would silently restore read weighting.

Parameters:

• df (DataFrame)

• weight (str)

Return type:

DataFrame

The signature column contract: raw feature groups, supports, pass-through channels.

A signature is a fixed-order, name-addressed feature vector for one sample. It is built in three steps, and this module owns the vocabulary of all three while computing none of them:

raw features  ->  rotation (per locus, from a corpus artifact)  ->  PC columns
                                                                +  pass-through channels

Names are four colon-separated parts throughout, so one parse() reads every kind:

vsig:div:TRB:1D_c        a raw feature      (an input to the rotation)
vsig:pc:TRB:PC01         a rotated column   (an output)
vsig:cov:TRB:cstar       a channel          (carried through untouched)
vsig:qc:-:n_loci_present a channel that is not per-locus

Two things here exist because getting them wrong produced silent wrong answers before.

``support`` is declared, never inferred. Winsorization trims the tail a metric can run away into, and which tail that is depends on the metric’s support, not on its transform. Three raw features declare transform="none" and have three different supports – a log-probability is bounded above by 0, a standard deviation is bounded below by 0, and a log-ratio is bounded neither way. Deriving the trimming side from the transform silently trims the wrong end of two of them, so the side comes from SUPPORTS and nothing else.

Channels are not rotated. A provenance number that went through a rotation is no longer a provenance number: it has been mixed with the measurements it was supposed to qualify. So the fallback fractions, the coverage the sample actually attained, and the masks pass through in their own units, and the rotation never sees them.

vsig (statistics) is declared here. rsig (geometry) is computed in mir.signature and registers its groups into this same registry with register_raw() / register_channel(); nothing here imports mir.

vdjtools.signature.layout.LOCI: tuple[str, ...] = ('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL')

The seven human receptor loci, in canonical signature order.

vdjtools.signature.layout.NO_LOCUS = '-'

Placeholder in the locus slot of a column that is not per-locus.

vdjtools.signature.layout.PC_BLOCK = 'pc'

Block name reserved for rotated output. <sig>:pc:<locus>:PC01.

vdjtools.signature.layout.TRANSFORMS: tuple[str, ...] = ('none', 'log10', 'log1p', 'logit', 'clr', 'arcsine')

Variance-stabilising transforms a raw feature may declare, applied where the feature is computed, while its denominator is still in scope. See vdjtools.signature.transform.

vdjtools.signature.layout.SUPPORTS: dict[str, tuple[bool, bool]] = {'nonneg': (False, True), 'nonpos': (True, False), 'real': (True, True), 'unit': (False, False)}

Supports a raw feature may declare, and the winsorization side each implies. This is the closed vocabulary the trimming side is read from.

support	trimmed	because
nonneg	top only	[0, inf) – cannot run away downward
nonpos	bottom only	(-inf, 0] – a log-probability, bounded by 0 above
real	both	(-inf, inf) – log-ratios, clr, logit, PC scores
unit	neither	[a, b] – already bounded on both sides

vdjtools.signature.layout.AMINO_ACIDS: str = 'ACDEFGHIKLMNPQRSTVWY'

Amino acids, anchored order. Also the k-mer alphabet.

vdjtools.signature.layout.SPECTRATYPE_LENGTHS: tuple[int, ...] = (6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31)

Junction lengths the spectratype resolves, in residues. Fixed rather than data-driven so the raw column set is the same for every sample and every corpus; anything outside falls in the two saturating end bins.

vdjtools.signature.layout.KMER_K: int = 2

k for the k-mer group. 2, giving 400 columns per locus over the plain 20-letter alphabet. k=3 is 8,000 cells against roughly 1,200 junction tokens at the corpus median depth – over 85% structural zeros, whose log-ratio coordinates read as a depth measurement wearing a motif label. k=2 is ~1,200 observations over 400 parts, which is thin but genuinely estimated.

vdjtools.signature.layout.feats(transform, support, *names)

{name: (transform, support)} for a run of features sharing both.

Groups are heterogeneous – a diversity group carries log10 Hill numbers beside a logit evenness – so both properties belong to the feature, not the group. This keeps the homogeneous runs terse anyway.

Parameters:

• transform (str)

• support (str)

• names (str)

Return type:

dict[str, tuple[str, str]]

class vdjtools.signature.layout.RawGroup(sig, name, features=<factory>, loci=None, dynamic=None, named=False)

One named family of raw features – an input to the rotation, never an output.

Parameters:

• sig (str) – Owning signature, "vsig" or "rsig".

• name (str) – Group name, unique within a sig.

• features (dict[str, tuple[str, str]]) – {name: (transform, support)}, most easily built with feats(). Empty for a group whose column names are not knowable without the germline; such a group declares dynamic instead.

• loci (tuple[str, ...] | None) – Loci the group is emitted for. None means all of LOCI; an empty tuple means the group is not per-locus and uses NO_LOCUS.

• dynamic (tuple[str, str] | None) – For a group whose width depends on the germline vocabulary (vus, jus, spec): one (transform, support) pair shared by every column, with the names resolved at fit time and then recorded in the corpus artifact. Keeping the vocabulary out of here is what lets import vdjtools stay free of arda.

• named (bool) – This group’s features are interpretable domain scalars – a diversity index, a read count, a singleton fraction – that a caller may legitimately want back in their own units rather than only as a rotation input. Declared here rather than inferred, because “is this number reportable” is a property of the feature and not of its width: pchem is 30 static columns and is not reportable, shm is one and is. Groups marked here are what named=True selects; see named_groups().

property emitted_loci: tuple[str, ...]

Loci this group emits for; (NO_LOCUS,) when it is not per-locus.

spec(feature)

(transform, support) for one of this group’s features.

Parameters:

feature (str)

Return type:

tuple[str, str]

columns(locus, names=None)

Raw column names for one locus.

Parameters:

• locus (str) – A member of emitted_loci.

• names (list[str] | None) – Required for a dynamic group – the germline-resolved column names, in the order the artifact recorded them.

Return type:

list[str]

class vdjtools.signature.layout.Channel(sig, name, features, loci=None)

A named family carried through untouched – never winsorized, never rotated.

Parameters:

• sig (str) – Owning signature.

• name (str) – Channel name.

• features (dict[str, str]) – {name: support}. The support is recorded for documentation and for a range assertion; no trimming is applied either way.

• loci (tuple[str, ...] | None) – As RawGroup.

vdjtools.signature.layout.register_raw(*groups)

Add raw groups to the registry. Used by mir.signature for the geometry half.

Parameters:

groups (RawGroup)

Return type:

None

vdjtools.signature.layout.register_channel(*channels)

Add pass-through channels to the registry.

Parameters:

channels (Channel)

Return type:

None

vdjtools.signature.layout.raw_groups(sig=None)

Registered raw groups, in declaration order.

Parameters:

sig (str | None)

Return type:

list[RawGroup]

vdjtools.signature.layout.channels(sig=None)

Registered pass-through channels, in declaration order.

Parameters:

sig (str | None)

Return type:

list[Channel]

vdjtools.signature.layout.raw_columns(sig, locus, vocab=None)

Every raw feature name for one (sig, locus), in rotation-input order.

This is the column order of one row of the corpus matrix M_L, so it is also the row order of the rotation. It must be reproduced exactly at apply time, which is why the corpus artifact stores it verbatim rather than recomputing it.

Parameters:

• sig (str) – "vsig" or "rsig".

• locus (str) – A locus, or NO_LOCUS for the cross-locus group.

• vocab (dict[str, list[str]] | None) – {group_name: [column names]} for every dynamic group emitted at this locus. Required unless no dynamic group applies.

Return type:

list[str]

vdjtools.signature.layout.channel_columns(sig)

Every pass-through channel column for one sig, in emitted order.

Parameters:

sig (str)

Return type:

list[str]

vdjtools.signature.layout.named_groups(sig=None)

Raw groups whose features are reportable in their own units (RawGroup.named).

Parameters:

sig (str | None)

Return type:

list[RawGroup]

vdjtools.signature.layout.resolve_named(sig, named)

Normalise a named= argument to a tuple of group names.

True means every group declared named; False and () mean none. A sequence is taken literally and checked, so a typo raises here rather than silently emitting nothing – a caller who asks for "diversity" and gets no diversity columns has no way to tell that from a corpus that could not compute them.

Parameters:

• sig (str)

• named (bool | Sequence[str])

Return type:

tuple[str, …]

vdjtools.signature.layout.named_columns(sig, named=True, vocab=None)

Raw column names carried through in natural units, in layout order.

Parameters:

• sig (str) – "vsig" or "rsig".

• named (bool | Sequence[str]) – True for every reportable group, or an explicit sequence of group names.

• vocab (dict[str, dict[str, list[str]]] | None) – {locus: {group: [names]}} from the corpus. When given it also restricts the loci: a raw feature is only computed for a locus the corpus models, so emitting vsig:div:TRA:* from a TRB-only corpus would hand back a column of holes that never had a chance of being anything else. NO_LOCUS is always kept, because the cross-locus group is computed whatever the corpus models.

Return type:

list[str]

vdjtools.signature.layout.channel_table(sig=None)

Every feature this half can emit, with what it is and what units it is in.

One row per (block, feature) with a kind saying how it reaches the output:

rotated

A rotation input. You do not get this column; you get <sig>:pc:<locus>:PCnn.

named

A rotation input that is also reportable on its own, via named= on vsig / rsig. Same number, natural units.

channel

Carried through untouched – never winsorized, never rotated.

The point of the kind column is that “what can I get, and in what units” was previously three lookups and a reading of corpus.py.

Returns:

A list of dicts, ready for polars.DataFrame(...). Plain dicts so that importing layout stays free of polars.

Parameters:

sig (str | None)

Return type:

list[dict[str, str]]

vdjtools.signature.layout.pc_columns(sig, k)

Rotated column names, given the component count per locus.

Parameters:

• sig (str) – "vsig" or "rsig".

• k (dict[str, int]) – {locus: n_components}, from the corpus artifact. A locus absent from k emits nothing, which is how a corpus that could not fit a locus declares it.

Return type:

list[str]

vdjtools.signature.layout.signature_columns(sig, k)

The emitted signature for one half: rotated columns then channels, in that order.

Parameters:

• sig (str)

• k (dict[str, int])

Return type:

list[str]

vdjtools.signature.layout.parse(column)

Split a column name into (sig, block, locus, feature).

Also the definition of “is this a signature column”, used as an allow-list wherever a frame may carry a caller’s own joined columns – so an age or n_reads column can never be silently winsorized or rotated.

Raises:

ValueError – If the name is not four colon-separated parts.

Parameters:

column (str)

Return type:

tuple[str, str, str, str]

vdjtools.signature.layout.support_of(column, vocab=None)

The declared support of a raw feature or channel column.

Parameters:

• column (str) – A raw or channel column name.

• vocab (dict[str, dict[str, list[str]]] | None) – Unused for lookup – a dynamic group’s support is declared on the group, not per column – and accepted so callers can pass their artifact’s vocab uniformly.

Raises:

ValueError – If no registered group or channel declares the column.

Return type:

str

Raw vsig features and pass-through channels, for one sample.

Everything here is a pure function of one sample plus frozen vocabulary. Nothing is fitted, no corpus is consulted, and no value is standardised, winsorized or rotated – those need a corpus and live in vdjtools.signature.corpus. That split is the point: two people who never share a cohort compute the same raw features, and only the artifact they rotate through can differ.

Two axes this module keeps apart, because collapsing them is how wrong answers got shipped:

• Raw features feed the rotation. Each carries its declared transform, applied here while its denominator is still in scope, because a proportion separated from the count it was observed on cannot be stabilised afterwards.

• Channels never reach the rotation. A provenance number that has been mixed with the measurements it was meant to qualify is no longer provenance, so the fallback fractions, the masks, and the coverage the sample actually attained pass through in their own units.

The coverage target is a runtime argument, not a frozen constant. Carrying a per-sample quantity in a corpus artifact is how one shipped reference ended up standardising tissue samples to a level measured on blood, digit for digit across all seven loci.

vdjtools.signature.features.VALID_AA = '^[ACDEFGHIKLMNPQRSTVWY]+$'

The 20 proteinogenic amino acids, anchored. A junction containing anything else – a stop codon, an ambiguity code, the legacy out-of-frame marker – is dropped before anything is computed. Not a crash guard: most downstream code accepts a malformed junction silently and returns a finite, meaningless number, so this filter is the only thing between a contaminated feature and a plausible-looking wrong answer. The dropped fraction is reported as qc:*:nonstd_aa_frac.

vdjtools.signature.features.PARSEABLE_AA = '^[ACDEFGHIKLMNPQRSTVWY*_]+$'

The 20 amino acids plus the two characters a WELL-FORMED AIRR table can legitimately carry: * for a stop codon, _ for the legacy out-of-frame marker. Anything outside this is a broken file, not a kind of receptor – see sanitise().

vdjtools.signature.features.ISOTYPES: dict[str, tuple[str, ...]] = {'IgA': ('IGHA1', 'IGHA2'), 'IgD': ('IGHD',), 'IgE': ('IGHE',), 'IgG': ('IGHG1', 'IGHG2', 'IGHG3', 'IGHG4'), 'IgM': ('IGHM',)}

Isotype classes and the constant-gene names that map onto them. IGHGP is a pseudogene and IGHC is ambiguous, so neither is called.

vdjtools.signature.features.WEIGHTS = {'anscombe': <function <lambda>>, 'distinct': <function <lambda>>, 'duplicate_count': <function <lambda>>, 'log1p': <ufunc 'log1p'>, 'log2p1': <function <lambda>>}

Clone-size weights g, the one definition both halves of the signature use. Concave by default: raw read weighting lets a single dominant clone be the entire profile, presence weighting throws the expansion signal away, and log2(1+a) sits between.

It lives in the shared dependency because vsig and rsig must weight the same repertoire identically or they are not describing one measure – and a second copy that drifted would not raise anywhere, it would quietly make the two halves disagree.

vdjtools.signature.features.PAIRS: tuple[tuple[str, str], ...] = (('TRA', 'TRB'), ('TRG', 'TRB'), ('TRD', 'TRB'), ('IGK', 'IGL'), ('IGH', 'TRB'))

T-vs-B balance, the gamma-delta share, the light-chain balance. A ratio rather than two counts, because the sequencing depth that drives both cancels.

Type:

Locus pairs whose read-count ratio is a compartment read-out

vdjtools.signature.features.assert_parseable(df)

Raise if any junction_aa is not a parseable amino-acid string.

A data-integrity gate, not a filter, and the first of three axes to keep apart: parseable (can we read the string at all – here), productive (does the rearrangement encode a chain – AIRR), functional (is the germline gene real – IMGT F/ORF/P).

Measured across 6,047,716 rows of a clinical AIRR store: zero violations. The gate therefore costs nothing on well-formed data and is only ever reached by a real problem.

Raises:

ValueError – If any junction carries an unparseable character.

Parameters:

df (DataFrame)

Return type:

None

vdjtools.signature.features.sanitise(df, *, strict=True, on_duplicate='error')

Drop non-productive clonotypes; return the frame and the dropped weight fraction.

Dropped by weight, not by row: losing one dominant clone matters more than losing fifty singletons, and a row fraction would hide that.

Parameters:

• df (DataFrame) – A clonotype frame.

• strict (bool) – Raise on unparseable characters (see assert_parseable()). False drops them, for a corpus known to carry ambiguity codes.

• on_duplicate (str) – What to do when the frame has no junction_nt and repeats an amino-acid clonotype key – "error" (default) or "sum". The default refuses because the two readings give different richness and clonality and the frame cannot say which is meant.

Returns:

(kept_frame, dropped_weight_fraction).

Return type:

tuple[DataFrame, float]

vdjtools.signature.features.work_frame(df, weight='log2p1')

Overwrite frequency with the normalised clone weight, so every profiler agrees.

Order matters. Call this after filtering, and never call filter_functional, downsample or select_top afterwards – each recomputes frequency from the counts and would silently restore read weighting.

Parameters:

• df (DataFrame)

• weight (str)

Return type:

DataFrame

vdjtools.signature.features.gene_vocab(locus, organism='human')

Column names for the three germline-width raw groups at one locus.

arda’s germline is the single source of truth for V/J vocabulary across this ecosystem, so the usage and spectratype widths are a property of the germline release rather than of any cohort. The corpus artifact records what came back here, and apply time reads it from the artifact – never from arda again, because a germline update would silently re-index the rotation.

Returns:

{"vus": [...], "jus": [...], "spec": [...]}. spec is <V gene>_<length> over SPECTRATYPE_LENGTHS, V-major.

Parameters:

• locus (str)

• organism (str)

Return type:

dict[str, list[str]]

vdjtools.signature.features.div_group(df, level)

Hill numbers standardised to a coverage level, plus distribution shape.

Coverage standardisation is what makes a hundred clonotypes and a hundred thousand comparable: both are read at the same completeness rather than at their own depth. It only works where the sample actually reaches level; beyond that the estimator extrapolates, and extrapolation is not a mild approximation here. Measured on one real repertoire subsampled across a 200x depth range: at a level every subsample attained, Shannon diversity agreed within 4%; at a level the shallow ones extrapolated to, the same statistic inflated roughly tenfold.

Whether this sample reached it is reported by mask:*:estimable, and the level it actually attained by cov:*:cstar. When it did not reach the target the features are nan – a hole a downstream model can see, rather than a confident wrong number.

Parameters:

• df (DataFrame)

• level (float)

Return type:

dict[str, float]

vdjtools.signature.features.estimable(df, level, est=None)

Whether a coverage-standardised diversity is a measurement rather than an extrapolation.

Two necessary conditions: the estimator must not have extrapolated, and the target depth must be within twice the observed one. The second catches the case where it interpolates nominally but is leaning on almost no data.

Parameters:

• df (DataFrame)

• level (float)

• est (DataFrame | None)

Return type:

bool

vdjtools.signature.features.depth_group(df, n_reads, richness)

How much was seen, and how much was not.

S_unseen is Chao’s estimate of the clonotypes that exist but were never drawn. Carried explicitly rather than folded into a diversity estimate because it is the honest statement of what the sample could not resolve.

Parameters:

• df (DataFrame)

• n_reads (float)

• richness (float)

Return type:

dict[str, float]

vdjtools.signature.features.clon_group(df, n_reads, richness)

Clone-size distribution shape as a composition, plus the top clone’s share.

f1/f2/f3plus – seen once, twice, more – is a three-part composition read in log-ratio coordinates. Two of the three ship because the third is exactly determined by them.

Parameters:

• df (DataFrame)

• n_reads (float)

• richness (float)

Return type:

dict[str, float]

vdjtools.signature.features.len_group(df)

Weighted moments of the junction length distribution.

Moments rather than a pooled histogram: the per-V-gene histogram is the spec group’s job, and a locus-pooled 26-bin histogram is neither of the two useful things.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.aa_group(df)

Weighted single-residue composition of the junctions, arcsine-stabilised.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.kmer_group(df)

Weighted k-mer composition of the junctions, arcsine-stabilised.

k=2 over the plain 20-letter alphabet: 400 parts against roughly 1,200 junction tokens at the corpus median depth. Thin, but genuinely estimated – k=3 would be 8,000 cells, over 85% of them structural zeros, whose coordinates read as a depth measurement wearing a motif label.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.pchem_group(df, regions=('all', 'center'))

Weighted mean physicochemistry of the junction, over two regions.

Two regions rather than five, and means rather than four quantiles apiece: at a hundred clonotypes the quantiles of a per-clonotype property are noise, while the weighted mean is a well-behaved average over every residue seen.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.usage_group(df, col, genes)

Weighted gene-usage composition over a fixed germline vocabulary, in clr coordinates.

The vocabulary is fixed by the artifact, not by the sample, so a gene this sample never used is a structural zero of a known composition rather than a missing column. An unrecognised call is different and is not silently folded in here – it lands in the closing residual part and is reported by qc:*:{v,j}_fallback_frac.

The clr is taken over the vocabulary plus that residual, then the residual is dropped: a clr of a sub-composition is a different number from the corresponding coordinate of the full one.

Parameters:

• df (DataFrame)

• col (str)

• genes (list[str])

Return type:

dict[str, float]

vdjtools.signature.features.spec_group(df, names, v_genes)

Spectratype: the junction-length distribution resolved per V gene, in clr coordinates.

A (V gene x length) table, closed by one residual cell for calls or lengths outside the vocabulary. This is the widest group by a wide margin – 121 IGH V genes x 26 lengths is 3,146 cells against a median 464 IGH clonotypes per real sample, so a single sample’s table is ~0.15 observations per cell and the great majority of it is structural zeros.

That sparsity is real and is not hidden: a corpus of 10,000 repertoires puts ~4.6M observations behind the same cells, so the rotation’s directions are well determined even where one sample’s scores along them are noisy. How noisy is readable per sample from cov:*:cstar, which is why that channel is always emitted.

Parameters:

• df (DataFrame)

• names (list[str])

• v_genes (list[str])

Return type:

dict[str, float]

vdjtools.signature.features.iso_group(df)

Isotype composition of an IGH repertoire, in log-ratio coordinates.

The uncalled share is a real part of the composition, not a rounding error – roughly two fifths of IGH reads carry no constant-gene call – so it closes the composition rather than being silently dropped the way a usage profile would drop it.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.shm_group(df)

Mean somatic hypermutation load, via v_identity.

Absent from the canonical schema – vdjtools’ readers narrow to eight columns and v_identity is not one of them – so this masks out unless the caller kept it explicitly.

Parameters:

df (DataFrame)

Return type:

dict[str, float]

vdjtools.signature.features.pair_group(reads)

Log read-count ratios between loci – compartment balance, with depth divided out.

Parameters:

reads (dict[str, float])

Return type:

dict[str, float]

vdjtools.signature.features.coverage_of(df)

The sample’s own attained Chao coverage at its observed depth, in [0, 1].

Parameters:

df (DataFrame)

Return type:

float

vdjtools.signature.features.qc_channel(raw, clean, locus, nonstd_frac)

Whether this sample’s gene calls are in a vocabulary we recognise.

An unrecognised V or J call raises nowhere downstream: the germline distance lookup falls back to the maximum observed distance, so a cohort on Adaptive nomenclature or an older IMGT release yields a fully populated, entirely plausible, systematically wrong result. These fractions are the one number telling a collaborator their vector is not comparable to ours, which is why they are columns and not a warning nobody reads.

Reported as raw fractions in [0, 1], not logit-transformed: a channel is read by a human deciding whether to trust the row, and a logit is not.

Parameters:

• raw (DataFrame)

• clean (DataFrame)

• locus (str)

• nonstd_frac (float)

Return type:

dict[str, float]

vdjtools.signature.features.raw_and_channels(frames, vocab, *, cstar_target=None, weight='log2p1', prefiltered=False, on_duplicate='error')

Every raw feature and every channel for one sample.

Parameters:

• frames (dict[str, DataFrame]) – {locus: clonotype frame}, raw (this sanitises).

• vocab (dict[str, dict[str, list[str]]]) – {locus: {"vus": [...], "jus": [...], "spec": [...]}} from the corpus artifact, or from gene_vocab() at fit time.

• cstar_target (float | dict[str, float] | None) – Coverage level the Hill numbers are read at. A float applies to every locus, a dict is per-locus, and None uses each locus’s own attained coverage – which makes the diversity features observed rather than standardised, and is the only default that does not smuggle in a constant from somewhere the caller cannot see.

• weight (str) – Clone-size weight, a key of WEIGHTS.

• prefiltered (bool) – When the caller already removed non-productive rows, report nonstd_aa_frac as nan rather than a confident floor of zero.

• on_duplicate (str) – Forwarded to sanitise().

Returns:

(raw, channels) – both {column: value}, holes as nan.

Return type:

tuple[dict[str, float], dict[str, float]]

The corpus: build one, winsorize it, fit the rotation and the scaling, apply it to a sample.

This is the only module here that needs more than one sample, and everything it fits comes out of one pass over one matrix of repertoires. That is the whole design, and it is a correction: the system this replaces fitted its rotation on 10,000 individual clonotypes from a prototype panel while every column that rotation produced was a repertoire statistic, and its centre and scale arrived from a separate corpus of real samples. Two independent fits, stitched – which is how one shipped reference came to pair a centre of exactly 0.0 with a scale plainly fitted from data, putting a corpus-typical sample 81 robust deviations out.

So, per locus L:

M_L (N, p_L)  one row per repertoire, columns in layout order
  -> bounds       per column, BY PERCENTILE, side from the declared support
  -> centre/scale median and 1.4826*MAD, from OBSERVED entries only
  -> rotation     R_L (p_L, k_L)
  -> pc centre/scale/quantiles   from the corpus's own PC scores

Nothing per-sample enters the artifact. The coverage a sample attained is a per-sample statistic and ships as a channel; the coverage target is a runtime argument. Carrying a per-sample quantity in a corpus artifact is how one shipped reference ended up standardising tissue to a level measured on blood, digit for digit across all seven loci.

Winsorization is by percentile and one-sided wherever the metric is bounded on that side. A robust-SD bound is not a fixed amount of probability: it depends on the distribution’s shape, so an 8-SD bound trims nothing at all on a right-skewed count and a third of the corpus on a near-degenerate column. And the side comes from the declared support, never from the transform – transform="none" spans three different supports.

vdjtools.signature.corpus.MAD_TO_SD = 1.4826

1.4826 * MAD estimates the standard deviation of a normal distribution.

vdjtools.signature.corpus.WINSOR_PS: tuple[float, ...] = (0.01, 0.05)

Winsorization percentiles the artifact stores, so a caller picks one at apply time without a refit. 0.01 trims the 1st/99th percentile, 0.05 the 5th/95th.

vdjtools.signature.corpus.MIN_SPREAD = 1e-12

it is centred, not scaled, and contributes nothing to the rotation. Not a tolerance on a measurement – a structural zero of a composition really is constant, and dividing by its MAD would turn float noise into the corpus’s dominant direction.

Type:

A column whose corpus spread is at or below this is treated as never having varied

vdjtools.signature.corpus.DEFAULT_COMPONENTS: int = 128

Default components per locus. A count rather than a variance fraction, because the artifact stores the rotation only up to the count it was fitted at and a fraction on a wide sparse group resolves to enough components to make the artifact hundreds of megabytes. The full eigenvalue spectrum is stored either way, so the variance a given count reaches is always readable.

vdjtools.signature.corpus.MODES: tuple[str, ...] = ('features', 'pcs', 'none')

Modes for apply(). features clamps raw features to the corpus bounds before rotating; pcs rotates first and clamps the PC scores; none clamps nothing.

vdjtools.signature.corpus.bounds_for(x, supports, p)

Per-column winsorization bounds at percentile p, side chosen by declared support.

Parameters:

• x (ndarray) – (N, p_L) corpus matrix; non-finite entries are ignored.

• supports (list[str]) – One support name per column, from support_of().

• p (float) – Tail fraction, e.g. 0.01 for the 1st/99th percentile.

Returns:

(lo, hi), each (p_L,). A side that is not trimmed is -inf / +inf, so the clamp is a no-op there rather than a bound that happens not to bite.

Return type:

tuple[ndarray, ndarray]

vdjtools.signature.corpus.clamp(x, lo, hi)

Clamp to [lo, hi], and report which entries moved.

The second return value is what keeps the bound from being silent. Clamping is many-to-one: it overwrites a real measurement of a real sample, and the sample cannot tell afterwards. A reference whose centre was wrong by 57 robust deviations went unnoticed for months precisely because every affected sample was clamped to the bound and looked like an ordinary tail case.

Parameters:

• x (ndarray)

• lo (ndarray)

• hi (ndarray)

Return type:

tuple[ndarray, ndarray]

vdjtools.signature.corpus.robust_loc_scale(x)

Per-column (median, 1.4826*MAD) from observed entries only.

Observed entries only, before any imputation. Filling holes first and measuring afterwards deflates the scale in proportion to how sparse a column is, so the least-observed locus ends up with the largest apparent values and dominates every distance and every principal component.

Parameters:

x (ndarray)

Return type:

tuple[ndarray, ndarray]

class vdjtools.signature.corpus.LocusFit(columns, supports, loc, scale, rotation, eigenvalues, pc_loc, pc_scale, bounds, pc_bounds, n_obs, n_rows)

Everything fitted for one locus. Arrays are column-aligned to columns.

Parameters:

• columns (list[str])

• supports (list[str])

• loc (ndarray)

• scale (ndarray)

• rotation (ndarray)

• eigenvalues (ndarray)

• pc_loc (ndarray)

• pc_scale (ndarray)

• bounds (dict[float, tuple[ndarray, ndarray]])

• pc_bounds (dict[float, tuple[ndarray, ndarray]])

• n_obs (ndarray)

• n_rows (int)

variance_at(k)

Cumulative variance fraction the first k components reach.

Parameters:

k (int)

Return type:

float

class vdjtools.signature.corpus.Corpus(sig, name, vocab, fits, meta=<factory>)

A fitted corpus: one LocusFit per locus, plus the manifest.

A load verifies rather than trusts: the raw column order is re-derived from the stored vocabulary through the current layout and compared, so a layout change that would silently re-index the rotation raises instead. The bundled-model version is checked the same way, because the synthetic corpora are drawn from those models and retraining them moves the result by 0.1-1.6% per locus.

Parameters:

• sig (str)

• name (str)

• vocab (dict[str, dict[str, list[str]]])

• fits (dict[str, LocusFit])

• meta (dict)

property k: dict[str, int]

{locus: components}. A locus absent here emits no columns – how a corpus that could not fit a locus declares it, rather than shipping a rotation of zeros.

columns(n_components=None, named=())

The signature columns this corpus emits: rotated, then channels, then named blocks.

Parameters:

• n_components (int | float | None)

• named (bool | Sequence[str])

Return type:

list[str]

resolve_k(n_components=None)

{locus: k} after apply-time truncation.

Parameters:

n_components (int | float | None) – None keeps what was fitted. An int truncates every locus to that count. A float in (0, 1) truncates each locus to the fewest components reaching that cumulative variance.

Raises:

ValueError – If an int exceeds what a locus was fitted with, or a float target is not reached by the stored spectrum. Neither ever pads: a count is capped at each locus’s stored k (the same cap the fit applied, so a 5-feature cross-locus block yields 5 whether you asked for 5 or 128), and a variance fraction the stored spectrum cannot reach raises rather than silently returning the whole rotation under the requested name.

Return type:

dict[str, int]

save(path)

Write <path>.npz (arrays) and <path>.json (manifest and vocabulary).

Parameters:

path (str | Path)

Return type:

Path

classmethod load(path, *, verify=True)

Read an artifact, and check it is the one the caller thinks it is.

Parameters:

• path (str | Path) – The .npz path; the .json sidecar must sit beside it.

• verify (bool) – Re-derive the raw column order from the stored vocabulary through the current layout and compare, and check the recorded bundled-model version against the installed one. Turn it off only to inspect a deliberately stale artifact.

Raises:

ValueError – On a column-order or model-version mismatch. Both would otherwise produce numbers in a different coordinate system that look perfectly reasonable.

Return type:

Corpus

verify()

Raise unless the stored column order and model version match this installation.

Return type:

None

vdjtools.signature.corpus.fit_locus(x, columns, *, mode='features', n_components=128, winsor_p=0.01)

Fit bounds, centre, scale, rotation and PC scaling for one locus, in one pass.

Parameters:

• x (ndarray) – (N, p) raw feature matrix. Rows that observed nothing at this locus must already be dropped; individual holes are fine and are handled without imputing before measuring.

• columns (list[str]) – Raw column names, defining the rotation’s row order.

• mode (str) – "features" winsorizes the corpus and then takes a plain SVD – the covariance the SVD sees has no tails left to be dragged by, so the robustness sits in a step you can inspect, and the result is deterministic and bit-reproducible. "pcs" leaves the corpus untrimmed and uses a robust scaler with sklearn’s PCA, because there the tails are still present at fit time and the estimator has to resist them itself.

• n_components (int | float) – Count, or a variance fraction in (0, 1).

• winsor_p (float) – The percentile whose bounds are used when mode="features". All of WINSOR_PS are stored regardless.

Returns:

The fit, or None when the locus cannot support one (fewer rows than 2, or no column that ever varied). None is how a corpus declares a locus it could not fit; it is not an error, and it emits no columns.

Return type:

LocusFit | None

vdjtools.signature.corpus.locus_matrix(vocab, sig, locus, n)

A preallocated (n, p_L) buffer and its column order, or None if the locus is empty.

Preallocated and filled in place rather than built from a list of per-sample dicts. At the shipped corpus size a dict of 13,483 Python floats per sample is 4-5 GB of interpreter objects for a matrix that is 1.08 GB as float64 – and the dicts have to coexist with the array while it is assembled. See fill_row().

Parameters:

• vocab (dict[str, dict[str, list[str]]])

• sig (str)

• locus (str)

• n (int)

Return type:

tuple[ndarray, list[str]] | None

vdjtools.signature.corpus.fill_row(buf, cols, i, raw)

Write one sample’s values for one locus into row i of a preallocated buffer.

Parameters:

• buf (ndarray)

• cols (list[str])

• i (int)

• raw (dict[str, float])

Return type:

None

vdjtools.signature.corpus.fit_matrices(mats, vocab, *, sig, name, mode='features', n_components=128, winsor_p=0.01, meta=None)

Fit a corpus from per-locus matrices – the memory-bounded path.

Each locus is fitted only on the rows that observed it, so a corpus where half the samples have no TRD does not learn TRD from imputed values.

Parameters:

• mats (dict[str, tuple[ndarray, list[str]]])

• vocab (dict[str, dict[str, list[str]]])

• sig (str)

• name (str)

• mode (str)

• n_components (int | float)

• winsor_p (float)

• meta (dict | None)

Return type:

Corpus

vdjtools.signature.corpus.register_featuriser(sig, fn)

Declare how one half turns {locus: frame} into (raw, channels).

Called once at import by each half. fit_cohort() is the only consumer: it is what lets a caller hand in clonotype frames rather than pre-computed feature rows.

Parameters:

sig (str)

Return type:

None

vdjtools.signature.corpus.featuriser(sig)

The registered featuriser for sig, or a pointed error naming what to import.

Parameters:

sig (str)

vdjtools.signature.corpus.fit(rows, vocab, *, sig, **kw)

Fit a corpus from raw feature rows – one dict per repertoire.

Convenient for a small corpus and for tests. For a full-size build use locus_matrix() / fill_row() / fit_matrices(), which never holds more than one sample’s dict at a time. To fit on your own cohort of clonotype frames rather than on feature rows you assembled yourself, use fit_cohort().

Parameters:

• rows (list[dict[str, float]])

• vocab (dict[str, dict[str, list[str]]])

• sig (str)

Return type:

Corpus

vdjtools.signature.corpus.fit_cohort(samples, *, sig='vsig', name='custom', vocab=None, loci=None, organism='human', n_components=128, weight='log2p1', n_jobs=0, **kw)

Fit a corpus on your own cohort of clonotype frames.

The shipped corpora are one route and this is the other. fit takes feature rows, so a caller holding repertoires had to run the featuriser per sample and assemble the dicts themselves; this does that, in processes, and hands back the same Corpus type the published artifacts are. Corpus.save() writes it, Corpus.load() reads it back, and a collaborator can score against it exactly as against a shipped one.

Parameters:

• samples – {sample_id: sample}, an iterable of pairs, or an iterable of samples. A sample is {locus: frame}, one frame with a locus column, or a picklable zero-argument callable returning either.

• sig (str) – Which half to fit. "rsig" requires mir.signature to have been imported.

• name (str) – Recorded in the manifest, and reported in the preflight line wherever this corpus is used. Name it after the cohort, not after the analysis.

• vocab (dict[str, dict[str, list[str]]] | None) – {locus: {group: [names]}}. Defaults to the germline vocabulary for loci.

• loci (tuple[str, ...] | None) – Which loci to model. Defaults to every locus any sample carries.

• organism (str) – Germline organism, when vocab is derived rather than given.

• n_components (int | float) – Components per locus, as an int or a cumulative variance fraction.

• weight (str) – Clone-size weight, forwarded to the featuriser.

• n_jobs (int) – Worker processes; 0 uses every available core.

• **kw – Forwarded to fit_matrices() (mode, winsor_p, …).

Returns:

A fitted Corpus.

Return type:

Corpus

Warning

A cross-validated score computed on a rotation fitted inside the same cohort is not evidence that the rotation generalises. Measured on 5,376 raw columns over seven loci with one logistic head: a rotation fitted on 612 training repertoires beat the shipped blood artifact on the training out-of-fold number at 4 of 5 matched widths (median 0.5039 against 0.4968) – and on the external cohort it led at 0 of 5 (median 0.4757 against 0.5270). Fitting in-cohort improves the number the configuration is selected on and neither held-out read-out, because a rotation fitted on one trial learns that trial’s covariance.

Both routes are legitimate and they answer different questions. Fit here when you want a basis for this cohort; use a shipped corpus when the number has to travel.

vdjtools.signature.corpus.apply(raw, chan, corpus, *, mode=None, winsor_p=None, n_components=None, named=())

Rotate one sample’s raw features through a corpus, and carry its channels through.

Parameters:

• raw (dict[str, float]) – Raw features for one sample, as raw_and_channels() returns them.

• chan (dict[str, float]) – That sample’s channels. Copied out untouched apart from qc:-:winsor_frac.

• corpus (Corpus) – A fitted Corpus.

• mode (str | None) – Override the artifact’s winsorization mode. "pcs" on an artifact fitted with "features" is allowed and is a legitimate thing to want, but the bounds it applies were then measured on an already-trimmed corpus; the manifest records which was fitted.

• winsor_p (float | None) – Which stored percentile to clamp at. Defaults to the fitted one.

• n_components (int | float | None) – Truncate to this many components (or this variance fraction). Exact – the components are ordered, so the first n of a longer rotation are the same vectors.

• named (bool | Sequence[str]) –

  Also carry through the reportable raw blocks, untouched by the winsorization and by the rotation – True for every group declared reportable, or an explicit sequence of group names. (), the default, reproduces the rotation-plus-channels output exactly. These are the same numbers the rotation is fitted on, not a second computation of them: raw already holds every one.

  They carry each feature’s declared transform, not its natural scale: a log10 diversity comes back as log10, and a clr composition as a log-ratio that has no unique inverse. channel_table() gives the transform per feature. Do not read one of these as a clone count.

Returns:

{column: value} in signature_columns() order, followed by the requested named blocks in layout order.

Return type:

dict[str, float]

vdjtools.signature.corpus.model_fingerprint(loci, source='olga')

{locus: "<model source>@<model version>"} for the models a synthetic pool is drawn from.

The version is a hash of the germline the generator will actually draw from (every allele’s CDR3-region cut segment of the collapsed model), not the manifest’s declared version string. A declared version does not move when a germline is repaired: correcting TRBV4-3*02’s anchor changed what every TRB pool contains while leaving olga:human_T_beta@2.0.0 identical, so a corpus fitted before the repair would have loaded silently against models that no longer produce it. The hash cannot miss that.

This, not the library version, is what a synthetic corpus’s numbers rest on: every pool comes out of these models, and retraining one moves that locus by 0.1-1.6%. Gating on the library version instead was both too strict – a patch release invalidated every corpus for no reason – and wrong: vdjtools.__version__ is read from installed distribution metadata, so a build driven by PYTHONPATH against a different installed version records that version rather than the code that ran. The first four cluster artifacts recorded 3.6.0 for exactly that reason, from a 4.0.0 tree.

A real corpus is fitted on real repertoires and no model enters it, so its manifest carries no fingerprint and Corpus.verify() has nothing to check.

Parameters:

source (str)

Return type:

dict[str, str]

vdjtools.signature.corpus.CORPUS_RELEASE = 'corpora-v1'

GitHub release the corpus artifacts are published under, in the repo that owns each half. A tag of its own rather than the library version: an artifact changes far less often than the code, and keying the download on the library version would invalidate every cached corpus on a patch release – the same mistake the germline fingerprint replaced for the load gate.

vdjtools.signature.corpus.CORPORA_INDEX = 'corpora.json'

The shipped index of downloadable artifacts. A few KB of JSON rather than ~110 MB of rotations, so a plain install stays small while –corpus can still NAME every corpus, report its size before fetching anything, and verify what it downloaded.

vdjtools.signature.corpus.corpus_cache_dir()

Where downloaded corpus artifacts are cached.

VDJTOOLS_CORPUS_DIR overrides it; otherwise the usual per-user cache. Shared between the two halves on purpose – a vsig/rsig pair belongs in one place, and a cluster job that pre-warms one warms both.

Return type:

Path

vdjtools.signature.corpus.corpora_index(res_dir)

{name: {"npz": {"sha256", "bytes"}, "json": {...}}} from the shipped index, or {}.

Parameters:

res_dir (Path)

Return type:

dict

vdjtools.signature.corpus.fetch_artifact(name, *, sig, res_dir, repo, quiet=False, base_url=None)

Download one corpus artifact and its manifest into the cache, verified. Returns the .npz.

Both files, because Corpus.load() reads the .json beside the .npz – fetching only the arrays would leave a corpus that cannot say what it was fitted on.

base_url overrides where the assets come from – an internal mirror, or a file:// URL, which is what lets the resolution order be tested without a network.

Parameters:

• name (str)

• sig (str)

• res_dir (Path)

• repo (str)

• quiet (bool)

• base_url (str | None)

Return type:

Path

vdjtools.signature.corpus.resolve_artifact(name, *, sig, res_dir, repo)

A corpus artifact from a path, the wheel, the cache, or the release – in that order.

The order is what makes a local file always win: a caller who fitted their own corpus and passes its path must never be silently served a downloaded one of the same name.

Parameters:

• name (str | Path)

• sig (str)

• res_dir (Path)

• repo (str)

Return type:

Path | None

vdjtools.signature.corpus.bundled_names()

Every corpus name this version knows – installed, cached, or downloadable.

Return type:

list[str]

vdjtools.signature.corpus.bundled_path(name)

Resolve a corpus name or a filesystem path to a vsig artifact, else None.

Downloads on first use if the name is in the shipped index and not yet cached; see resolve_artifact() for the resolution order.

Parameters:

name (str | Path)

Return type:

Path | None

vdjtools.signature.corpus.SEED: int = 20260927

Seed for the shipped synthetic corpora. Every draw is SEED + offset, recorded in the manifest, so a rebuild is reproducible by anyone who installs the library.

vdjtools.signature.corpus.ZIPF_A: float = 1.5

clone i of n gets frequency proportional to i ** -ZIPF_A. Recorded in the artifact rather than assumed, because it is one of the two free parameters of the selected-repertoire regime.

NOTE this is the Zipf law over ranks, not numpy.random.Generator.zipf, which samples integers from a Zipf distribution. At a = 1.5 that distribution has infinite mean, so normalising a draw of it gives one clone almost all the mass: measured here, repertoires collapsed to as few as 1 surviving clonotype, which then made every coverage-standardised diversity feature a hole. The rank form is a well-behaved rank-abundance curve.

Type:

Rank-abundance exponent for the memory corpus

vdjtools.signature.corpus.READS_PER_CLONE: int = 20

Reads per clone in the memory multinomial. Sets how much of the rank-abundance tail survives sampling, which is the realistic part: at 20 reads per clone a 500-clone repertoire keeps 336 of them with 129 singletons, which is the shape a real library has. The zeros are dropped, so a memory repertoire has fewer distinct clones than a naive one at the same nominal size – exactly as a sampled selected repertoire does.

vdjtools.signature.corpus.DEFAULT_SIZE: int = 10000

Receptors per synthetic repertoire, and the pool each is drawn from. A pool much larger than a sample is what makes two synthetic repertoires nearly disjoint, the way two donors are.

vdjtools.signature.corpus.DEPTH_SPREAD: dict[str, float] = {'IGH': 10.981981981981981, 'IGK': 10.407407407407407, 'IGL': 10.027777777777779, 'TRA': 2.9182389937106916, 'TRB': 4.051282051282051, 'TRD': 2.4285714285714284, 'TRG': 4.2727272727272725}

Multiplicative p05-to-p95 spread of repertoire depth, per locus, from N_EFF. Each synthetic repertoire’s size is drawn log-uniformly across this range around the nominal size.

A corpus at one fixed depth is not usable, and it fails silently. Measured here: with every naive repertoire at exactly 10,000 receptors, depth:reads and depth:richness are identical in all N samples and all five pair: log-ratios are constant to the last bit – so their corpus spread is 0, they contribute nothing to the rotation, the cross-locus block cannot be fitted at all, and a real sample’s depth features get standardised against a degenerate reference. Real repertoires span 4.1x (TRB) to 11.0x (IGH) between their 5th and 95th percentiles, so drawing the depth is not an embellishment: it is what makes the corpus describe the axis it is going to be asked about.

vdjtools.signature.corpus.N_EFF: dict[str, tuple[int, int, int]] = {'IGH': (111, 376, 1219), 'IGK': (54, 203, 562), 'IGL': (36, 123, 361), 'TRA': (159, 276, 464), 'TRB': (195, 401, 790), 'TRD': (14, 34, 34), 'TRG': (22, 49, 94)}

Per-locus median effective clone count of real bulk blood repertoires, and the p05/p95 of the same, measured on 1,168 samples as n_eff = 1/sum(w^2) with w = log2(1+count)/sum. Used by size="n_eff" and by the depth sweep, so the artifact can report how its bounds and its scaling move with repertoire depth instead of leaving that unmeasured.

vdjtools.signature.corpus.COHORT: dict[str, dict[str, tuple]] = {'blood': {'IGH': (33245, (84, 271, 697, 1711, 6200), (2.11, 2.42, 2.88, 3.93, 12.83), (0.255, 0.528, 0.697, 0.819, 0.927), (0.204, -0.132, -0.717)), 'IGK': (43678, (80, 233, 549, 1172, 3432), (2.29, 2.83, 3.67, 5.41, 15.97), (0.143, 0.384, 0.555, 0.704, 0.846), (0.399, -0.301, -0.759)), 'IGL': (37601, (61, 160, 361, 790, 2468), (2.33, 2.86, 3.63, 5.37, 15.89), (0.122, 0.346, 0.514, 0.671, 0.817), (0.349, -0.259, -0.714)), 'TRA': (32732, (67, 148, 276, 554, 1724), (2.08, 2.38, 2.75, 3.5, 6.98), (0.197, 0.566, 0.745, 0.851, 0.942), (-0.045, 0.142, -0.547)), 'TRB': (34365, (74, 208, 457, 985, 3162), (2.1, 2.4, 2.83, 3.74, 7.91), (0.215, 0.58, 0.759, 0.86, 0.949), (-0.036, 0.128, -0.503)), 'TRD': (5620, (19, 48, 80, 121, 280), (2.25, 2.96, 4.17, 6.75, 17.36), (0.167, 0.418, 0.591, 0.714, 0.856), (-0.547, 0.484, -0.482)), 'TRG': (8688, (21, 48, 72, 108, 220), (2.64, 3.51, 4.74, 7.0, 15.38), (0.167, 0.438, 0.583, 0.686, 0.806), (-0.421, 0.463, -0.394))}, 'tissue': {'IGH': (47031, (40, 164, 602, 2458, 10352), (2.39, 3.3, 5.27, 10.0, 40.36), (0.055, 0.395, 0.561, 0.69, 0.855), (0.017, 0.079, -0.509)), 'IGK': (58706, (34, 150, 519, 1813, 6592), (2.66, 4.09, 7.05, 14.9, 69.34), (0.147, 0.29, 0.41, 0.546, 0.75), (0.134, -0.147, -0.582)), 'IGL': (51010, (29, 113, 343, 1188, 4152), (2.64, 3.88, 6.49, 13.09, 54.7), (0.125, 0.26, 0.377, 0.51, 0.717), (0.075, -0.094, -0.565)), 'TRA': (16799, (18, 84, 148, 417, 1839), (2.1, 2.42, 2.95, 4.68, 120.24), (0.111, 0.476, 0.668, 0.806, 0.93), (-0.32, 0.34, -0.688)), 'TRB': (22298, (30, 89, 158, 498, 2977), (2.11, 2.48, 3.0, 4.6, 94.74), (0.194, 0.517, 0.688, 0.817, 0.928), (-0.233, 0.265, -0.598)), 'TRD': (1180, (6, 21, 54, 96, 253), (2.55, 4.05, 9.29, 53.4, 373.45), (0.12, 0.4, 0.594, 0.728, 0.876), (-0.661, 0.373, -0.219)), 'TRG': (2235, (7, 31, 57, 93, 172), (2.85, 4.42, 6.82, 13.29, 102.0), (0.077, 0.309, 0.5, 0.661, 0.82), (-0.54, 0.462, -0.377))}}

What a REAL repertoire looks like per locus, in the two compartments the mixture corpora are named after: {cohort: {locus: (samples, richness, reads per expanded clone, singleton fraction, rank correlations)}}, each of the three quantities being its measured p05/p25/p50/p75/p95 ladder.

Measured 2026-09-27 in one streaming pass over the harmonized AIRR store, every sample with at least 100 reads in the locus, grouped by sample_id. Richness is distinct clonotypes; f1 is the fraction of them seen exactly once; reads per expanded clone is (reads - singletons) / (richness - singletons). The sample count leads each row, so a reader can see which loci the table is thin on (tissue TRD, 1,180 samples) and which it is not (tissue IGK, 58,706).

These three ladders are the entire parameterisation of a mixture corpus, and the middle one is deliberately not the obvious quantity. Reads per clonotype cannot be drawn independently of the singleton fraction: a repertoire with f1 singletons whose other clones all carry at least 2 reads has at least 2 - f1 reads per clonotype, so a pair drawn from those two marginals lands in a forbidden region about half the time on blood TRB and has to be clipped back out of it – which silently rewrites the marginal it was drawn from. Drawing the read count itself is worse still: measured here, a copula over (richness, reads, f1) asks for a repertoire below that floor for 21.1% of 20,000 blood TRB draws. Reads per expanded clone is >= 2 whatever f1 is, so the trio has no forbidden region at all.

Five quantiles, not two, and it is the resolution that makes the corpus land. The read count is a product of all three drawn quantities, so it is the sharpest check on the draw: against blood TRB’s measured median of 763 reads (n = 34,365 samples), a 40,000-draw simulation gives 762 at five quantiles and 846 at three; blood IGH’s 1,298 comes out 1,276 against 1,490 (n = 33,245). Tissue is the harder case and is stated rather than smoothed over: tissue TRB’s median 308 comes out 425 at five quantiles, against 756 at three (n = 22,298). A product of three heavy-tailed factors has a median above the product of their medians unless their joint tails are matched exactly, and no marginals-plus-copula draw does that.

The last field is the Spearman rank correlations among the three, in the order (richness vs expanded count, richness vs f1, expanded count vs f1), drawn through a Gaussian copula. That is for the rotation, not for the read count: a corpus’s rotation IS its covariance structure, and f1 against expansion size is -0.503 on blood TRB and -0.846 on blood IGK (n = 43,678) – fewer singletons, bigger expansions. A corpus that drew the two independently would hand the PCA a correlation the cohort does not have, which no amount of correct marginals repairs.

Against all this, DEPTH_SPREAD – 2.4x on TRD to 11.0x on IGH, from 1,168 deep blood samples, and what the naive/memory corpora draw across – covers about a tenth of the real axis: blood TRB richness spans 43x (74 to 3,162 clonotypes) and tissue IGH 259x (40 to 10,352). A corpus estimates its bounds, its centre and its per-PC scaling from its own draw, and not one of the three extrapolates beyond it.

vdjtools.signature.corpus.COHORT_QS: tuple[float, ...] = (0.05, 0.25, 0.5, 0.75, 0.95)

The quantiles every COHORT ladder is recorded at.

vdjtools.signature.corpus.SYNTHETIC: dict[str, tuple[str, str | None]] = {'memory': ('memory', None), 'naive': ('naive', None), 'synthetic-blood': ('mixed', 'blood'), 'synthetic-tissue': ('mixed', 'tissue')}

{name: (regime, cohort)}. There are four, and the name is the only thing a caller has to know.

naive and memory are the two pure regimes at one nominal size. Neither is representative of a real bulk sample and neither can be: naive has no clone-size structure at all, and memory has no naive background and, at a fixed size, no read-depth variance either – its read count is exactly size * READS_PER_CLONE in every sample. The two synthetic-* corpora are the MIXTURE, drawn across the measured per-locus richness, read-depth and singleton-fraction ranges of the cohort they are named after, so they belong in the same set as the real blood / tissue / deep-tcr corpora and mean what their names say.

Type:

The named synthetic corpora

vdjtools.signature.corpus.corpus_plan(name, *, size=None, depth_spread=None)

(regime, cohort, size, depth_spread) for a named corpus – the ONE place a name becomes a draw, so the two halves of a corpus cannot resolve the same name differently.

A cohort name stands in for both size and depth_spread, because both come out of the same measured band: the nominal size is its geometric centre and the spread is its width. An explicit value for either wins, which is what makes a deliberately wider or narrower variant of a named corpus a one-flag change rather than a new table.

Parameters:

• name (str)

• size (int | str | None)

• depth_spread (float | str | None)

Return type:

tuple

vdjtools.signature.corpus.cohort_bands(cohort, locus)

One locus of one cohort: the richness, expanded-count and f1 ladders, and their copula.

Parameters:

• cohort (str)

• locus (str)

Return type:

tuple[tuple, tuple, tuple, tuple]

vdjtools.signature.corpus.ladder_draw(q, u, log=True)

Inverse-CDF draw from a measured COHORT_QS ladder, at probabilities u.

Linear interpolation of the quantile function – in log space for a scale quantity – so the drawn marginal carries all five of the cohort’s recorded quantiles and spans exactly its p05 to p95.

That is the difference from drawing uniformly across the band, and it is not cosmetic: real singleton fractions sit near the top of theirs (blood TRB median 0.759 of a 0.215-0.949 band, n = 34,365 samples), so a uniform draw would centre the corpus at 0.582 and leave a cohort-median sample 0.65 robust SD off centre on every singleton-sensitive feature.

Parameters:

• q (tuple)

• u (ndarray)

• log (bool)

Return type:

ndarray

vdjtools.signature.corpus.copula_uniforms(rank_corr, n, rng)

(n, 3) probabilities carrying the measured rank correlations – a Gaussian copula.

Spearman’s rho converts to the Gaussian correlation it implies by 2 * sin(pi * rho / 6), the three of them are made a matrix, and correlated normals are pushed through the normal CDF. The result is uniform in every column whatever the correlations are, so ladder_draw() still reproduces each measured marginal exactly – a copula only decides how they move together, which is what the rotation is fitted on.

Measured over all 14 (cohort, locus) entries of COHORT: the realised rank correlations match their targets to within 0.015 on 20,000 draws.

The correlation matrix is factorised in closed form, not by an eigendecomposition, and that is a reproducibility requirement rather than a micro-optimisation. numpy.linalg.eigh returns eigenvectors whose SIGN is a LAPACK convention, not mathematics: negating a column leaves the covariance – and so every marginal and every correlation – untouched, but it changes the realised sample. A corpus built on a machine whose LAPACK signs a column differently would then differ byte-for-byte from the same build here while being statistically identical, and a corpus whose values depend on the builder’s linear-algebra library cannot be compared with one built anywhere else. The Cholesky factor of a positive-definite matrix is unique, and at 3x3 it is six arithmetic operations, so there is nothing to defer to.

Unit diagonal makes the factor’s rows unit-norm by construction, so the normals need no rescaling. An infeasible triple – three correlations that cannot come from one joint distribution – makes the last radicand negative and raises, rather than being clipped into something plausible.

Parameters:

• rank_corr (tuple)

• n (int)

• rng (Generator)

Return type:

ndarray

vdjtools.signature.corpus.depth_spread_of(locus, depth_spread=None)

The multiplicative depth range for one locus: measured, named, or a caller’s override.

None means DEPTH_SPREAD, the p05-p95 spread of that locus among 1,168 deep blood samples, which is what naive and memory draw across – 2.4x on TRD to 11.0x on IGH, so under one decade everywhere. A cohort name from COHORT means that cohort’s own measured richness band, which is 43x on blood TRB and 259x on tissue IGH. A number widens or narrows it deliberately: a corpus meant to describe samples across decades of sequencing depth has to be drawn across those decades, because the bounds, the centre and the per-PC scaling are all estimated from the draw and none of them extrapolates.

Parameters:

• locus (str)

• depth_spread (float | str | None)

Return type:

float

vdjtools.signature.corpus.draw_pool(locus, n, *, seed, source='olga')

n productive rearrangements from the bundled model for one locus.

Measured 2026-09-27, generate(load_bundled("IGH"), 50_000, productive_only=True) on one core of a 16-core M-series laptop: 4,200 sequences/s, and identical at 1 and 16 polars threads – the sampler is a per-sequence Python/numpy loop, so it does not thread. Pool generation is therefore the dominant serial cost of a synthetic build (about 15M sequences at size=10,000 across seven loci) and is why build_pools() fans out over processes.

NOTE an earlier docstring here claimed 20,372 seq/s on IGH. That figure is 4.8x the measured rate and it is what made a 77 s stage look like a 16 s one when the build was planned.

Parameters:

• locus (str)

• n (int)

• seed (int)

• source (str)

Return type:

DataFrame

vdjtools.signature.corpus.draw_sizes(nominal, n, spread, rng)

n repertoire sizes, log-uniform over a spread-fold range around nominal.

Log-uniform rather than normal because depth is a scale, and the geometric centre is the nominal size so the corpus is centred where the caller asked for it. See DEPTH_SPREAD for why a fixed depth is not an option.

Parameters:

• nominal (int)

• n (int)

• spread (float)

• rng (Generator)

Return type:

ndarray

vdjtools.signature.corpus.draw_sample(pool, size, regime, rng, frac=None, mexp=None)

One synthetic repertoire: size receptors drawn from pool, with clone sizes.

naive gives every clone a count of 1, which is what an unselected repertoire looks like and what the generator produces. memory draws Zipf frequencies, takes a multinomial sample at the same total, and drops the zeros – so a memory repertoire has fewer distinct clones than a naive one at the same nominal size, exactly as a selected repertoire does.

mixed is what a real bulk sample is – a singleton naive background plus an expanded memory component – and it needs frac (the singleton fraction) and mexp (reads per expanded clone), both drawn per sample by draw_plan() from the cohort’s measured ladders. Neither pure regime can stand in for it: naive has no clone-size structure at all, and memory at a fixed size has a read count of exactly size * READS_PER_CLONE in every sample, so its depth does not vary either.

Parameters:

• pool (DataFrame)

• size (int)

• regime (str)

• rng (Generator)

• frac (float | None)

• mexp (float | None)

Return type:

DataFrame

vdjtools.signature.corpus.resolved_size(locus, size)

The per-locus receptor count a size asks for – an integer, or a measured depth name.

Separate from sample_stream() because the manifest has to record what each locus was actually drawn at: size="n_eff" means seven different depths, and a manifest saying only "n_eff" cannot tell a reader which.

Parameters:

• locus (str)

• size (int | str)

Return type:

int

vdjtools.signature.corpus.draw_plan(loci, *, n_samples, size, seed, depth_spread=None, cohort=None)

Every sample’s drawn size, singleton fraction and expanded count: {locus: {name: array}}.

Drawn in the parent, once, from one generator per locus seeded at seed + 1000 + i, so every worker is handed the same plan rather than reproducing it – and so the plan itself can be inspected, which is how the realised corpus gets checked against the cohort it is named after.

cohort=None leaves the two mixture knobs unset, which is correct for naive and memory and keeps their draw byte-identical to what the shipped artifacts were built from.

Parameters:

• n_samples (int)

• size (int | str)

• seed (int)

• depth_spread (float | str | None)

• cohort (str | None)

Return type:

dict

vdjtools.signature.corpus.draw_one(pools, plan, regime, j, seed)

Sample j of a synthetic corpus – a pure function of j.

Each sample gets its own generator, seeded from (seed, locus index, j), so any process can draw any sample in any order and get the same repertoire. That is what lets the build run across processes while staying bit-identical at every n_jobs: the drawing is the only part of the pipeline carrying RNG state, and here it carries none between samples.

Parameters:

• pools (dict)

• plan (dict)

• regime (str)

• j (int)

• seed (int)

Return type:

dict

vdjtools.signature.corpus.pool_target(locus, size, n_samples, depth_spread=None, cohort=None)

How many receptors a locus’s pool needs: the largest repertoire, times POOL_FACTOR.

The pool has to cover the LARGEST depth the draw can ask for, since a sample is drawn from it without replacement. For a cohort corpus that is the top of the measured richness ladder, scaled the same way draw_plan() scales it; for a pure regime it is the nominal size times the half-spread.

Parameters:

• locus (str)

• size (int | str)

• n_samples (int)

• depth_spread (float | str | None)

• cohort (str | None)

Return type:

int

vdjtools.signature.corpus.build_pools(loci, *, size, n_samples, seed, source, tmp, n_jobs=1, progress=None, depth_spread=None, cohort=None)

Generate every locus’s pool across processes; return {locus: [chunk path, ...]}.

Generation is the single most expensive stage of a synthetic build and it is embarrassingly parallel – measured on Aldan-3, seven pools at size=10,000 took 4,850 s in one process, which is 81 minutes of a build whose featurisation is the part anyone cares about.

Written as uncompressed IPC and read back memory-mapped, so N workers share one physical copy of a multi-gigabyte pool through the page cache instead of each pickling its own.

Parameters:

• n_samples (int)

• seed (int)

• source (str)

• tmp (Path)

• n_jobs (int)

• depth_spread (float | str | None)

• cohort (str | None)

Return type:

dict

vdjtools.signature.corpus.read_pools(paths)

Memory-map the pool chunks and present one frame per locus, without copying them.

rechunk=False is load-bearing: a rechunk would materialise the whole pool in this process, which is the per-worker copy this design exists to avoid.

Parameters:

paths (dict)

Return type:

dict

vdjtools.signature.corpus.POOL_CHUNKS: int = 64

Chunks each locus’s pool is generated in. Fixed, so the pool is a function of the build and not of the machine: a chunk count taken from the core count would make an artifact built on 64 cores differ from the same build on 16. Sixty-four because ONE CHUNK IS THE FLOOR on wall time – it cannot be split further, however many cores the machine has. At 16 chunks IGH’s 3.3M-sequence pool is 207k per chunk, which is 49 s at the measured 4,200 seq/s; at 64 it is 12 s. Changing this value changes the pools and so the artifact, which is why it is a recorded constant.

vdjtools.signature.corpus.build_matrices(regime, *, sig, featurise, vocab, loci=('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL'), n_samples=10000, size=10000, seed=20260927, source='olga', n_jobs=1, progress=None, tmp=None, depth_spread=None, cohort=None)

The corpus matrix, {locus: (array, columns)}, built across n_jobs processes.

Shared by both halves of the signature: vsig and rsig differ only in featurise and sig, and the samples they see are identical – same pools, same per-sample seeds, same depths – which is what makes the two artifacts joinable on sample_id.

Parallel across samples, never inside one. A sample’s featurisation is a pure function of that sample, so a worker needs no coordination; and because every worker draws its own samples from memory-mapped pools, the only thing that crosses a process boundary is one row of numbers per sample. Shipping drawn repertoires instead would be ~38 GB of pickling at the shipped size.

Parameters:

• featurise – A picklable callable mapping {locus: frame} to (raw, channels) – a module-level function or a functools.partial over one, never a lambda.

• n_jobs (int) – Worker processes; 0 means every core available. 1 runs in-process and is the default, so that a doctest or a 12-sample test does not spawn a pool.

• tmp (Path | None) – Scratch directory for the pool files. A private temporary directory by default, removed when the build finishes.

• regime (str)

• sig (str)

• vocab (dict)

• loci (tuple[str, ...])

• n_samples (int)

• size (int | str)

• seed (int)

• source (str)

• depth_spread (float | str | None)

• cohort (str | None)

Returns:

{locus: (matrix, columns)}, the matrix being (n_samples, p_locus) float64.

Return type:

dict

vdjtools.signature.corpus.single_threaded_children()

Make spawned workers take ONE kernel thread each, by setting the env they inherit.

n_jobs processes each starting a kernel sized off the core count is cores x cores threads. Measured here: seven pool workers on a 16-core box, each with a 16-thread polars, turned 15 s of receptor generation into 296 s. The variables have to be set in the parent, before any child exists – a spawned child reads them while importing polars, which is strictly earlier than any initializer of ours can run. The parent’s own kernel is already up, so this does not throttle it.

Return type:

None

vdjtools.signature.corpus.sample_stream(regime, *, loci=('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL'), n_samples=10000, size=10000, seed=20260927, source='olga', progress=None, tmp=None, n_jobs=1, depth_spread=None, cohort=None)

Set up the pools in this process and return (ordered_loci, draw); draw(j) is sample j.

The in-process path, for callers that want the repertoires themselves rather than a corpus matrix. draw may be called in any order – see draw_one().

Parameters:

• regime (str)

• loci (tuple[str, ...])

• n_samples (int)

• size (int | str)

• seed (int)

• source (str)

• tmp (Path | None)

• n_jobs (int)

• depth_spread (float | str | None)

• cohort (str | None)

vdjtools.signature.corpus.synthesize(corpus_name, *, loci=('TRA', 'TRB', 'TRG', 'TRD', 'IGH', 'IGK', 'IGL'), n_samples=10000, size=None, seed=20260927, n_components=128, mode='features', winsor_p=0.01, source='olga', fit_corpus=True, progress=None, n_jobs=1, depth_spread=None)

Build a synthetic corpus, and fit it.

Parameters:

• corpus_name (str) – One of SYNTHETIC – "naive" and "memory" are the pure regimes at a nominal size; "synthetic-blood" and "synthetic-tissue" are the naive/memory mixture drawn across the named cohort’s measured per-locus richness, read-depth and singleton-fraction bands (COHORT), and are the two that describe a real bulk cohort rather than a regime.

• loci (tuple[str, ...]) – Loci to build. All seven by default.

• n_samples (int) – Repertoires in the corpus.

• size (int | str | None) – Receptors per repertoire. None, the default, is the corpus’s own – 10,000 for the pure regimes, the geometric centre of the cohort’s measured richness band for a synthetic-* one. Also takes "n_eff" for the per-locus real medians in N_EFF, "p05" / "p95" for the sweep endpoints, or a cohort name.

• seed (int) – Base seed; every draw is a recorded offset from it.

• n_components (int | float) – Components per locus, or a variance fraction.

• mode (str) – Winsorization mode to fit at.

• winsor_p (float) – Percentile for the fitted bounds. All of WINSOR_PS are stored.

• source (str) – Bundled model set – "olga", "learned" or "arda".

• fit_corpus (bool) – False returns the raw samples instead of fitting, for scoring a held-out draw against a corpus fitted on another.

• progress – Optional callable(locus, done, total).

• n_jobs (int) – Worker processes (not kernel threads). 0 means every available core; 1, the default, runs in-process. The result is bit-identical at every value.

• depth_spread (float | str | None) – Multiplicative depth range each repertoire’s size is drawn log-uniformly across, around size. None is the corpus’s own – DEPTH_SPREAD for a pure regime (2.4x-11.0x), the cohort’s measured richness band for a synthetic-* one (43x on blood TRB, 259x on tissue IGH). Pass a number to override it: 1000 with size=3162 spans 100 to 100,000 receptors per locus. The bounds, centre and per-PC scaling are all estimated from the draw, so a corpus describes only the depths it was drawn across.

Returns:

(corpus, mats) when fitting – mats being {locus: (matrix, columns)}, the corpus matrix itself – else a list of {locus: frame} samples.

The whole build is a deterministic function of (corpus_name, loci, n_samples, size, seed, source) and the bundled models, so two machines at different thread counts must produce byte-identical artifacts. That is an acceptance criterion, not a hope: generate was once irreproducible across processes while being deterministic within one, because a marginal-table aggregation left its group order unspecified.

vdjtools.signature.corpus.corpus_meta(regime, cohort, loci, *, n_samples, size, seed, source, depth_spread, **extra)

The manifest of a synthetic build – everything needed to reproduce the draw, per locus.

Shared by both halves, because a manifest that disagrees between vsig_<name> and rsig_<name> about the depths, the singleton fractions or the germline they were drawn from is a joinability claim nobody can check.

Parameters:

• regime (str)

• cohort (str | None)

• n_samples (int)

• size (int | str)

• seed (int)

• source (str)

Return type:

dict

vsig and vsig_cohort: raw features, a corpus rotation, and the emitted vector.

Three steps, each owned elsewhere and joined here:

features.raw_and_channels(sample, vocab)   ->  raw features + channels
corpus.apply(raw, chan, corpus)            ->  rotated columns + channels carried through

A corpus is required. There is no default, deliberately: a silently chosen rotation is the mixed-units failure this rewrite exists to end, and two matrices rotated through different corpora are not comparable no matter how alike their column names look.

vdjtools.signature.signature.vsig(sample, corpus, *, mode=None, winsor_p=None, n_components=None, cstar_target=None, weight='log2p1', prefiltered=False, on_duplicate='error', named=(), columns=None)

The statistics half of the signature for one sample.

Parameters:

• sample – {locus: frame}, one frame with a locus column, or a zero-argument callable returning either.

• corpus (Corpus) – A fitted Corpus for sig="vsig".

• mode (str | None) – Winsorization mode override – "features", "pcs" or "none". Defaults to whichever the corpus was fitted with.

• winsor_p (float | None) – Which stored percentile to clamp at. Defaults to the fitted one.

• n_components (int | float | None) – Truncate to this many components, or this cumulative variance fraction.

• cstar_target (float | dict[str, float] | None) – Coverage level the Hill numbers are read at. None uses each locus’s own attained coverage, which makes them observed rather than standardised – honest for one sample, and not comparable across samples. Use vsig_cohort() for a cohort, which resolves a shared level.

• weight (str) – Clone-size weight; a key of WEIGHTS.

• prefiltered (bool) – Report qc:*:nonstd_aa_frac as nan because the caller already filtered.

• on_duplicate (str) – "error" or "sum", for a frame repeating an amino-acid clonotype key.

• named (bool | Sequence[str]) –

  Also return the reportable raw blocks – True for all of them, or a sequence such as ("div", "depth", "clon"). (), the default, emits exactly the rotated columns and channels. Values carry their declared transform (log10, clr, logit, …), which channel_table() reports.

  These are the diversity, depth, clone-size, junction-length, isotype, SHM and cross-locus yield numbers. They are computed either way, because the rotation is fitted on them; without this argument there is no supported route to reading them back, and a study that cannot report its own diversity floor has no floor. See channel_table() for the full list and its units.

• columns (list[str] | None) –

  Restrict the output to these columns, in layout order.

  Unlike the system this replaces, this does not skip work: the rotation for a locus is fitted jointly over every raw group at that locus, so every group must be computed before any of that locus’s components exist. Declining a locus entirely does skip it.

Returns:

{column: value} in layout order, holes as nan.

Return type:

dict[str, float]

vdjtools.signature.signature.attained_coverage(sample, *, on_duplicate='error')

Per-locus attained Chao coverage for one sample – the cheap first pass.

Only the count column is touched, so this is a small fraction of the cost of the full feature set and is what lets a cohort agree on one coverage level without anybody choosing a constant.

Parameters:

on_duplicate (str)

Return type:

dict[str, float]

vdjtools.signature.signature.vsig_cohort(samples, corpus, *, n_jobs=1, cstar_target='min', columns=None, **kw)

One row per sample, sample_id first.

Parameters:

• samples – {sample_id: sample} or an iterable of pairs. A sample may be a zero-argument picklable callable, which defers the read into the worker and keeps peak memory at O(n_jobs) samples rather than the whole cohort.

• corpus (Corpus) – A fitted corpus.

• n_jobs (int) – Worker processes. 1 stays in-process, 0 uses every available core.

• cstar_target (float | dict[str, float] | str | None) – "min" (default) reads every sample’s attained coverage first and standardises the whole cohort to the per-locus minimum, so the diversity columns are comparable across samples and the level is a property of this cohort rather than of somebody’s corpus. That is a cheap extra pass over the count column – but it does resolve deferred samples twice, so pass an explicit float or dict for a one-pass run. None reads each sample at its own coverage: one pass, not comparable across samples.

• columns (list[str] | None) – Restrict the output columns.

• **kw – Forwarded to vsig() – including named=, which adds the reportable raw blocks in their own units.

Returns:

A frame whose columns are sample_id then the corpus’s signature columns.

Return type:

DataFrame

Variance-stabilising transforms.

Every signature column is a different kind of number — a read count, a Hill number, a share of a composition, a coordinate of an embedding — and a downstream model should not have to know which. Each raw feature therefore declares one transform (vdjtools.signature.layout), applied where the feature is computed, so that what reaches the rotation is dimensionless and roughly symmetric.

Standardisation, winsorization and rotation are not here – they need a corpus, and they live in vdjtools.signature.corpus. This module is a pure function of one sample.

Small denominators are the whole problem. At the depths this signature has to work at — a median of order a hundred clonotypes, a quarter of samples below ten — a “fraction” is very often 0/3. A transform that sees only the ratio cannot tell 0/3 from 0/500, and maps both to the same place, asserting a precision the data does not have. So every transform of a proportion takes the denominator as well as the value, and shrinks toward the middle in proportion to how little was counted. That is the Haldane–Anscombe correction for a logit, the Anscombe correction for an arcsine, and a count-scaled multiplicative replacement for a CLR.

None of these are free parameters: 1/2 and 3/8 are the standard bias-minimising choices.

vdjtools.signature.transform.log10(x, floor=1.0)

log10 of a positive quantity, floored so an empty locus maps to 0 rather than -inf.

Used for counts and for Hill numbers. The floor is 1 because both are counts of things: one clonotype, one read, one effective species. Zero of them is the same as the floor for every downstream purpose, and the presence mask already records that the locus was empty.

Parameters:

floor (float)

vdjtools.signature.transform.log1p(x)

log(1+x) for a non-negative quantity that genuinely reaches zero.

Distinct from log10() in intent: this is for norms, dispersions and hit counts, where zero is a real, attainable value rather than an empty measurement.

vdjtools.signature.transform.logit(x, m)

Haldane–Anscombe logit of a proportion observed on a denominator m.

log((x·m + 1/2) / ((1−x)·m + 1/2)). Adding half an observation to each side is the standard remedy for an empty cell; its side effect is exactly the behaviour wanted here — the transform of 0 depends on how many chances there were to see something:

logit(0, m=3)   ->  -1.95      a fifth of the repertoire could hide here
logit(0, m=500) ->  -6.91      it is really absent

Parameters:

• x – Proportion(s) in [0, 1].

• m – Denominator(s) the proportion was observed on. Broadcasts against x.

Returns:

The transformed value, finite for every input including exactly 0 and exactly 1.

vdjtools.signature.transform.arcsine(x, m)

Anscombe’s variance-stabilising arcsine transform, asin(sqrt((x·m + 3/8)/(m + 3/4))).

The right transform for a sparse composition — residue and gene-usage profiles, where most cells are structurally zero at shallow depth. Unlike a CLR it is defined at zero without any replacement step, and unlike a raw proportion its variance does not collapse near the boundary. Bounded in [0, π/2], so it cannot produce the heavy tail a log-ratio would.

Parameters:

• x – Proportion(s) in [0, 1].

• m – Denominator(s) the proportion was observed on.

vdjtools.signature.transform.clr(parts, m=None, *, keys=None)

Centred log-ratio of a composition, with multiplicative zero replacement.

A CLR is the natural coordinate for a composition whose ratios carry the meaning: it is the log of each part over the geometric mean of all of them, so it is invariant to the total and a difference between two coordinates is a log-ratio of two parts.

Zeros are replaced multiplicatively, not additively: each zero part is set to delta = 0.5/m and the non-zero parts are scaled down by 1 − n_zero·delta so the composition still closes. Adding a constant to every part instead — the common shortcut — distorts the ratios among the parts that were observed, which are the only ratios the coordinate system is about.

Compute the CLR over the whole composition, then select coordinates. A CLR of a sub-composition is a different number from the corresponding coordinate of the full one, so a tier that ships four of six isotype parts must still divide by the six-part geometric mean. Doing it the other way would make the narrower tier stop being a slice of the wider one, which is the contract the layout exists to guarantee.

Parameters:

• parts – Non-negative part values as a mapping {name: value} or an array. They need not sum to 1; they are closed here.

• m – Total count the composition was observed on, setting the replacement scale. Defaults to the sum of parts when they are counts.

• keys – Part order when parts is an array. Ignored for a mapping.

Returns:

{name: clr} when parts is a mapping (or keys is given), else an array. The coordinates sum to zero by construction, which is why the layout ships all but one of them: the last is exactly determined by the rest and would make any unregularised design matrix singular.

Raises:

ValueError – If fewer than two parts are given, or any part is negative.

vdjtools.signature.transform.apply(code, x, m=None)

Apply the transform named by code to one value or array.

Parameters:

• code (str) – A transform code from vdjtools.signature.layout.TRANSFORMS.

• x – The value(s).

• m – The denominator, required by logit and arcsine.

Raises:

ValueError – If code is unknown, if a denominator is missing for a transform that needs one, or if code is "clr" (call clr() with the whole composition).