Repertoire signature — the statistics half

vdjtools.signature turns one AIRR sample into a fixed, named, positional vector of repertoire statistics — vsig. Its companion is mir.signature (rsig), which covers the embedding geometry; the two are namespaced so they concatenate on sample_id without colliding, and the shared contract machinery (column layout, transform registry, frozen-reference rescaling) lives here, in vdjtools, because mirpy depends on vdjtools and not the reverse.

Quickstart — one command

No Python needed. Point vdjtools signature at your samples, pick a preset, get one row per sample:

# a metadata sheet plus a directory of samples
vdjtools signature --preset classify -m metadata.txt --base-dir samples/ -o sig.tsv

# or just pass files (AIRR, native vdjtools, Parquet, third-party — auto-detected)
vdjtools signature --preset compact sampleA.tsv sampleB.tsv.gz -o sig.tsv

# which columns am I about to get? reads no input at all
vdjtools signature --preset classify --describe

# the named feature sets, ranked, and what each one is for
vdjtools presets
vdjtools presets classify

Three presets are marked recommended: compact (the smallest vector that still describes a repertoire, usable at n = 50), classify (general-purpose, the usual random-forest / boosting input), and transfer (for a model that must work on another lab’s samples). See Feature presets — pick by intent, not by column for the full table.

Important

This command emits the vsig half only. For a classifier you usually want both halves — run mir signature --preset classify ... in mirpy, which emits vsig and rsig as one vector. A preset spanning both halves keeps only its vsig: columns here and says so on stderr.

Warning

CDR3 vs junction. The reader prefers AIRR junction_aa (conserved anchors included) and falls back to IMGT cdr3_aa (anchors excluded). A file carrying only cdr3_aa is two residues short everywhere, which shifts the length, k-mer and Pgen features. Check your headers before you trust a matrix.

The Python API

from vdjtools.signature import vsig, vsig_cohort
from vdjtools.signature import layout

v = vsig({"TRB": df})                       # {column: value}, in layout order
F = vsig_cohort(samples, tier="standard")   # one row per sample, positional
layout.columns("full")                      # the contract, computed from no data

Columns are named vsig:<block>:<locus>:<feature>. Tiers core / standard / full are exact index subsets of one frozen layout — a narrower tier is a slice of a wider one, never a differently-computed number.

Every column is transformed, and the transform is denominator-aware

A learner cannot be handed a read count, an amino-acid fraction and a Hill number in one matrix and be expected to weight them sensibly. Each block therefore declares a transform:

transform	blocks	why this one
arcsine	aa (×20)	Anscombe \(\arcsin\sqrt{(xm + 3/8)/(m + 3/4)}\) — a proportion with a known denominator
logit	qc shm div:clonality pgen:frac_atypical clon:top	Haldane–Anscombe, so 0/3 and 0/500 are different numbers
clr	iso clon (ships \(k-1\) parts)	a genuine closed composition; all \(k\) parts are linearly dependent
log10 / log1p	depth div	non-negative, right-skewed, unbounded above
none	len pchem pair pgen:mean_log10	already a moment, a log-ratio or a log-probability
none, exempt	mask	a 0/1 indicator; rescaling it against a reference would destroy it

Why arcsine and not log1p with winsorization

This comes up every time, so here are the measurements (benchmark_transform_arcsine.py in 2026-mirpy-analysis regenerates them).

Depth invariance. Hold the true proportion \(p\) fixed and vary the denominator from \(m = 91\) (the corpus median junction-residue count for a shallow RNA-seq sample) to \(m = 50{,}000\):

transform	value at m = 91	value at m = 50,000	drift
log1p(count)	0.919	6.908	7.5×
arcsine(p, m)	0.1476	0.1419	1.04×

The same biological composition sequenced twice as deep is a different feature value under log1p of counts. That alone decides it: depth is the largest nuisance axis in a public-data cohort, and a transform that encodes depth into every column hands the learner the batch label.

Variance stabilisation. Ratio of the largest to the smallest per-bin standard deviation across the 20 amino-acid columns:

denominator	log1p	arcsine	reading
m = 91	3.1–3.9×	3.1–3.9×	indistinguishable when counts are small
m = 5,000	21.6×	1.0×	only arcsine survives depth

Zeros carry information. arcsine(0, m) depends on \(m\) — a zero from 91 residues and a zero from 50,000 residues are different evidence, and the transform says so. log1p(0) = 0 always, which asserts that never observing a residue in a shallow sample and never observing it in a deep one are the same observation.

Winsorization is a fitted parameter. A 5% clip has to be estimated from a corpus, so it inherits that corpus’s depth distribution — and at the top of the distribution it clips the dominant V genes and the most abundant residues, which is where the biology usually is. Arcsine has no free parameter: the \(3/8\) is Anscombe’s, not ours.

Can you compose them? arcsine(log1p(x)) is undefined for \(x \ge 2\) — log1p(2) = 1.10 and \(\arcsin\) needs \([0,1]\). The useful composition is the one the pipeline already does: log the counts into a clone weight, normalise to a proportion, then arcsine that proportion against the true count denominator. Under heavy clonal expansion the log2p1 clone weight cuts the column standard deviation about 4× (0.0319 → 0.0080) and restores arcsine’s 1.14× stabilisation, which raw frequency weighting loses.

Holes are nan, never 0

A locus that was not sequenced, or a statistic the sample is too shallow to estimate, yields nan plus a vsig:mask: column. A model that reads “absent” as “zero” reads an unsequenced chain as a biological finding.

Ambiguous V calls resolve to the first gene

Real data does not hand you one V gene. Adaptive/immunoSEQ material realigned with MiXCR against the IMGT reference — the better reference, but junction plus short flanks rather than full primered reads — routinely calls a comma-separated tie. In one store slice that produced 1,296 distinct “genes”, of which 1,235 were comma-strings rather than genes at all, which silently shatters any V-keyed feature space into singleton columns.

vdjtools.io.schema.resolve_gene() takes the first gene in the list and strips the allele — MiXCR orders by alignment score, so the first is the best-supported call:

from vdjtools.io.schema import resolve_gene, strip_allele

df.with_columns(resolve_gene(pl.col("v_call")))   # "TRBV5-1*01,TRBV5-5*01" -> "TRBV5-1"
df.with_columns(strip_allele(pl.col("v_call")))   # keeps the tie, sorted, for reporting

The two are deliberately different functions. strip_allele sorts the tie so that reporting is order-insensitive; composing it with “take the first” would therefore give you the alphabetically first gene (TRBV19 where the aligner said TRBV5), which is a plausible-looking wrong answer. Use resolve_gene for anything keyed on V.

The V + k-mer space

vdjtools.features.kmer_space is the heavy feature: a junction k-mer profile keyed jointly on the V gene, TF-IDF scaled, and projected onto a frozen truncated-SVD basis. The k-mer counting is C++ (src/kmer.cpp); a Python implementation is not viable at corpus scale.

from vdjtools.features.kmer_space import fit_kmer_space, save_kmer_spaces
from vdjtools.signature.kmer import register_kmer

sp = fit_kmer_space(frames, pattern="xxxx", n_groups=8, flank=4,
                    min_df=0.02, max_df=0.80, n_components=32)
register_kmer({"TRB": sp})            # adds vsig:kmer:TRB:PC01.. to the layout
sp.transform(df, weight="freq")       # project one repertoire

Patterns. "xxxx" is a contiguous 4-mer; "xx.x" is gapped (the . is a skipped position). Gapped patterns exist because a conserved motif with one variable position is otherwise split across every substitution at that position.

Alphabets. n_groups=20 is the plain amino-acid alphabet. Smaller values cluster residues by BLOSUM62: seqtree.SubstitutionMatrix.blosum62().penalty() is \(s(a,a) + s(b,b) - 2s(a,b)\), which is a squared distance in the Gram sense, so classical MDS recovers the exact embedding and Ward linkage on it gives the groups. No hand-picked chemistry classes.

The sparsity trade, measured. On TRB with V-keying and a 4-mer:

alphabet	code space	surviving columns	effect
A = 20	3,360,000	1,022	most codes never occur; the survivors are common motifs
A = 8	86,016	25,402	dense, but distinct motifs collapse into one column

Neither is universally right. max_df matters as much as the alphabet: at A=8, max_df=0.80 a known disease motif was cut from the vocabulary entirely (0/4 windows present), while at A=20, max_df=0.99 all 4 survived. If you are looking for a rare discriminative motif, keep the document-frequency ceiling high and do not project — see below.

Do not select components by variance

A truncated SVD keeps the directions of greatest variance in the fitting corpus — which, for a repertoire matrix, are sequencing depth, V-gene usage and batch. A disease motif carried by a handful of clonotypes in a handful of donors is not one of those, so projecting onto the leading components discards it by construction.

Measured on ankylosing spondylitis vs healthy, both HLA-B27+, with the space fitted on a disjoint TRB cohort (benchmark_kmer_ankspond.py):

read-out	AUC	perm. p	verdict
sum of the 17 columns the published motif occupies (A = 20)	0.769	0.031	the representation carries it
same, A = 8	0.731	0.061	the reduced alphabet blurs it
best of 64 SVD components	0.841	0.20	chance — null median 0.297, 95th 0.385
best single vocabulary column	0.813	—	a label-selected maximum, not an estimate

The two “best-of-N” rows are the point. They are the larger numbers and they are the ones that mean nothing: a maximum taken over 64 components reaches \(|\mathrm{AUC} - 0.5| = 0.34\) under a label permutation null. Any best-of-N read-out must be nulled or not quoted.

How many components, then? Measured on the emitted signature matrix (14,553 samples × 1,369 columns, 182 studies), mean AUC over the four largest tasks with study-disjoint folds and the rotation refit inside every fold:

components	8	16	64	256	512	all 1,369
mean AUC	0.575	0.581	0.589	0.591	0.584	0.555

Flat from 16 to 256, then a cliff — the full matrix scores worse than 16 components. Use 16–64: 16 buys 98% of the achievable AUC at a quarter of the width, 64 is the plateau. Keep all columns only for a regularised learner (L1, boosting) or a rare sparse signal.

The practical consequence for feature selection:

• rare, discriminative signal — keep the un-projected TF-IDF columns and use an L1 model. The SVD is the wrong tool; it optimises for the wrong thing.

• broad compositional shift — project, and keep only components that survive a study-disjoint split-half refit (per-component score correlation ≥ 0.95), not components that reach a cumulative-variance threshold.

• either way — report the permutation null alongside the number.

Regenerating these numbers

Every table on this page is produced by a script in the companion 2026-mirpy-analysis repo, so it can be re-measured rather than trusted:

script	what it re-measures
benchmark_transform_arcsine.py	the arcsine vs log1p depth-drift and variance-stabilisation tables
benchmark_kmer_ankspond.py	the k-mer transfer AUCs and their permutation nulls
benchmark_signature_dimension.py	how many components reproduce across a study-disjoint refit
fit_kmer_spaces.py	the fitted per-locus spaces and their vocabulary sizes

Measurements on this page were last taken 2026-08-13.

Feature presets — pick by intent, not by column

The signature is over 1,400 columns. Almost nobody wants all of them, and which subset is right depends on the question — a model that must run on another lab’s samples wants different columns from one scoring samples inside a single study. vdjtools.signature.presets names those choices, documents each, and ranks it:

recommended

Use this unless you have a reason not to.

specific

Correct for a stated purpose and wrong outside it.

avoid

A control, a baseline, or a measured dead end. Named so that choosing it is deliberate.

preset	rank	columns	what it is
compact	recommended	152	The smallest vector that still describes a repertoire. Start here.
transfer	recommended	550	For models that must work on another lab’s samples. Drops the columns whose level moves most between studies.
classify	recommended	615	The general-purpose set. Best measured task performance when train and test come from comparable cohorts.
statistics	specific	101	Classical repertoire statistics only. Needs no embedding, so vdjtools alone suffices.
bcell	specific	286	B-cell receptor work: the immunoglobulin loci with somatic hypermutation and isotype.
geometry	specific	514	Embedding coordinates only — no count statistics at all.
full	specific	1403	Every contract column. For feature selection, not for fitting.
nuisance	avoid	73	Sequencing protocol only. A control, not a feature set.

Every preset resolves to a column list from the frozen layout alone — block names, loci, tier. No corpus, no fitted artifact and no private data is involved, so two people selecting the same preset get the same columns in the same order.

vdjtools presets                      # the table above
vdjtools presets transfer             # one preset in full: features, how, use cases, caveats
vdjtools signature *.tsv --preset transfer --describe    # the exact columns it selects

from vdjtools.signature import presets

presets.get("transfer").rank        # 'recommended'
cols = presets.columns("compact")   # a concrete, ordered column list
presets.table()                     # the whole registry as a DataFrame

Where the rankings come from

A benchmark over a public multi-study AIRR corpus — 182 study groups over 198 accessions, 14,553 samples — scored with study-disjoint folds: fit on some studies, predict on studies the fit never saw. Under that split a column that merely encodes sequencing protocol scores at chance, which is the point. Three findings shaped the presets:

• A nuisance floor of depth + presence masks + call quality is a surprisingly strong predictor on many contrasts. Any feature set worth using has to beat its own floor, which is why nuisance ships as a named control rather than being hidden.

• Projection did not help. Plain robust or asinh scaling beat PCA at every rank tested, so no preset projects by default and full is documented as a feature-selection tool rather than a model input.

• The two halves have opposite nuisance profiles. The embedding geometry carries several times less study-to-study variance than the count statistics and the most donor-to-donor variance, and is nearly unaffected by whether a sample is blood or tissue — but wins fewer supervised tasks outright. Hence transfer and geometry for robustness, classify and statistics for raw accuracy.

Anyone with a comparable SRA/AIRR corpus can reproduce this; none of it depends on a private dataset.

API

vdjtools.signature — the repertoire-signature contract, and the statistics half of it.

See vdjtools.signature.layout for what a signature is and how columns are named.

class vdjtools.signature.Block(sig, name, features, loci=None, attributable=False, exempt=False, magnitude=False)

One named feature family.

Parameters:

• sig (str) – Owning signature — "vsig" (statistics) or "rsig" (geometry).

• name (str) – Block name; may be declared by several Block entries as long as their features are disjoint (e.g. the all-loci and IGH-only halves of mask).

• features (dict[str, tuple[str, str]]) – {feature_name: (minimum_tier, transform)}, most easily built with feats(). A feature appears in every tier at or above its own; insertion order is the emitted column order.

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

• attributable (bool) – Whether a column has a clonotype pre-image, i.e. whether asking “which clones drive this” is a well-posed question. Declared here at build time, never inferred from the name — a Hill number and a read fraction are summaries and have no pre-image, so asking is a category error rather than an unanswered question.

• exempt (bool) – Skip the frozen reference rescaling entirely (masks, which are already 0/1).

• magnitude (bool) – Rescale by one frozen scalar for the whole block, with no centring, because the block’s magnitude is its signal. Per-column standardisation would force every coordinate to unit variance and make a near-zero sample indistinguishable from a typical one — it deletes exactly the deficiency the block exists to carry.

sig: str

name: str

features: dict[str, tuple[str, str]]

loci: tuple[str, ...] | None = None

attributable: bool = False

exempt: bool = False

magnitude: bool = False

transform(feature)

The transform declared for one of this block’s features.

Parameters:

feature (str)

Return type:

str

property emitted_loci: tuple[str, ...]

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

columns(tier='full')

Column names this block contributes at tier, in emitted order.

Parameters:

tier (str)

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.registry(sig=None)

Registered blocks, optionally filtered to one signature.

Parameters:

sig (str | None)

Return type:

list[Block]

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.columns(tier='standard', sig=None, blocks_=None)

The column names of a signature, in emitted order.

Parameters:

• tier (str) – "core", "standard" or "full". Tiers are nested and each is an exact index subset of "full".

• sig (str | None) – "vsig", "rsig", or None for both concatenated.

• blocks – Restrict to these block names (both signatures unless sig is given).

• blocks_ (tuple[str, ...] | None)

Returns:

Column names as <sig>:<block>:<locus>:<feature>.

Raises:

ValueError – If tier is unknown, or blocks_ names a block that is not registered.

Return type:

list[str]

vdjtools.signature.describe(tier='standard', sig=None)

The column dictionary — one row per column, so a collaborator can read the contract.

Returns:

A pl.DataFrame with column, sig, block, locus, feature, tier, transform, attributable, exempt, magnitude.

Parameters:

• tier (str)

• sig (str | None)

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

{name: (tier, transform)} for a run of features that share both.

Blocks are frequently heterogeneous — a clonality block carries CLR-transformed count fractions beside a logit-transformed top-clone share — so the transform is a property of the feature, not of the block. This keeps the homogeneous runs terse anyway.

Parameters:

• tier (str)

• transform (str)

• names (str)

Return type:

dict[str, tuple[str, str]]

vdjtools.signature.index(tier='core', sig=None)

Positions of tier’s columns within the full-width layout of the same sig.

This is what makes the tiers cheap: emit full once, then slice.

Parameters:

• tier (str)

• sig (str | None)

Return type:

list[int]

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.magnitude_scale(block, rms)

Rescale a whole block by one frozen scalar, with no centring.

For a block whose magnitude is its meaning — the signed contrast to an unselected reference, where a sample with no immune deviation should land at the origin. Standardising such a block per column would force every coordinate to unit variance across samples, which makes a near-zero sample look exactly like a typical one and deletes the deficiency the block exists to carry.

vdjtools.signature.parse(column)

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

Raises:

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

Parameters:

column (str)

Return type:

tuple[str, str, str, str]

vdjtools.signature.reference_z(x, loc, scale, clip=8.0)

Rescale against the frozen reference: (x − loc) / scale, clipped.

loc and scale are the reference median and 1.4826·MAD — robust, so a handful of pathological samples in the reference corpus cannot set the scale for everyone. They are frozen: a collaborator does not fit them, which is what makes their vector comparable to ours rather than merely internally consistent.

A zero scale means the reference never saw this column vary; the value is passed through centred but unscaled rather than divided by zero.

Parameters:

clip (float)

vdjtools.signature.register(block)

Add a block to the registry (used by mir.signature for late-bound vocabularies).

Raises:

ValueError – If a block with the same sig and name already declares any of the same features — a duplicate column would silently shift every index after it.

Parameters:

block (Block)

Return type:

None

vdjtools.signature.robust_loc_scale(x, axis=0)

Reference (median, 1.4826·MAD) from observed values only.

Non-finite entries are ignored rather than imputed. Computing the statistics before any imputation matters: filling first and then measuring deflates the scale in proportion to how sparse the column is, so the least-observed locus ends up with the largest apparent values and dominates every distance and every principal component.

vdjtools.signature.vsig(sample, *, tier='standard', cstar=0.2, weight='log2p1', pgen_q05=None, kmer_spaces=None, threads=0)

The vsig half of one sample’s signature, as {column_name: value}.

Parameters:

• sample – {locus: clonotype frame}, or a single frame with a locus column.

• tier (str) – "core", "standard" or "full".

• cstar (float | dict[str, float]) – Coverage level for the standardised Hill numbers; a scalar or {locus: level}.

• weight (str) – Clone-size weight g (see work_frame()).

• pgen_q05 (dict[str, float] | None) – Per-locus frozen 5th-percentile log10 Pgen for pgen:*:frac_atypical; that column stays nan without it, since “atypical” is meaningless without a reference to be atypical against.

• kmer_spaces (dict | None) – Per-locus frozen KmerSpace. The kmer block is emitted only for loci present here AND only if the block has been registered (see vdjtools.signature.kmer.register_kmer()) – the columns do not exist in the layout otherwise, so passing spaces without registering them is a no-op rather than a silent width change.

• threads (int) – Worker threads for the Pgen batch; 0 = auto. Off by default so tier="full" still runs; on when a caller needs the guarantee that every declared column was actually computed.

Returns:

Every column vdjtools.signature.layout.columns() lists for tier and "vsig", in that order, with nan where the sample could not support one.

Raises:

ValueError – If tier is unknown.

Return type:

dict[str, float]

vdjtools.signature.vsig_cohort(samples, *, tier='standard', **kw)

Assemble a whole cohort into one frame: sample_id plus the vsig columns.

Parameters:

• samples – {sample_id: sample} or an iterable of (sample_id, sample).

• tier (str) – Passed to vsig().

• **kw – Passed to vsig().

Returns:

A pl.DataFrame, one row per sample, columns in layout order.

A frozen (V gene x k-mer) feature space: reduced alphabets, TF-IDF, and a sparse basis.

vdjtools.features.kmer produces tidy per-sample profiles. This module produces the thing a model wants: a fixed-width vector per sample, on a vocabulary and a rotation that were decided once and then never move, so two collaborators’ matrices are comparable.

Three problems stand between a k-mer count and that vector, and each has a specific answer here.

Sparsity. A V-pinned 4-mer over 20 residues is a code space of ~8e6, of which one repertoire touches ~4e5 and any two repertoires share far fewer. Two independent reductions apply: the alphabet (reduced_alphabet() collapses the 20 residues into BLOSUM62 groups, so chemically equivalent substitutions stop being different features) and the vocabulary (a document-frequency window, below).

Scale. Raw k-mer weight tracks sequencing depth and CDR3 length before it tracks anything immunological. TF-IDF plus an L2 row norm removes both: IDF down-weights the germline-adjacent k-mers that appear in every repertoire, and the row norm makes a deep and a shallow sample comparable.

Noise. The surviving matrix is still tens of thousands of columns wide and mostly zero, so it is projected onto a small basis by truncated SVD. scipy.sparse.linalg.svds operates on the sparse matrix directly and never densifies it – and, deliberately, scikit-learn is not used: scipy is a base dependency of vdjtools while sklearn is an extra, and a portable feature block that imports an extra fails late, after every other block of a sample has already been computed.

The document-frequency window is what makes the vocabulary robust rather than merely small. A k-mer present in nearly every sample carries no contrast (it is germline or near it), and one present in a handful cannot be estimated; both ends are cut, which is the same reasoning that bounds the public-clonotype panel’s incidence window.

vdjtools.features.kmer_space.reduced_alphabet(n_groups=8)

Partition the 20 amino acids into n_groups by BLOSUM62 similarity.

seqtree’s SubstitutionMatrix.penalty(a, b) = s(a,a) + s(b,b) - 2 s(a,b) is the Gram transform of the log-odds matrix – that is, a squared Euclidean distance. So the spectral embedding needs no Laplacian and no kernel choice: classical MDS (double-centre, eigendecompose) recovers coordinates whose pairwise distances are exactly that penalty, and the clustering happens there. Building an affinity matrix and a normalised Laplacian would be the same computation with an arbitrary bandwidth bolted on.

Ward linkage on those coordinates, via scipy.cluster.hierarchy – a base dependency, and for 20 points the choice of clusterer is not where the answer comes from.

Parameters:

n_groups (int) – Number of groups. 20 returns the identity partition (each residue its own group), which is the un-reduced alphabet.

Returns:

{residue: group index} over the 20 standard amino acids, group indices contiguous from 0. Deterministic: no random initialisation anywhere in the path.

Return type:

dict[str, int]

vdjtools.features.kmer_space.alphabet_table(groups)

A 256-entry char -> group table for the native kernel, and the group count.

Everything not in groups maps to -1, which voids any window containing it. That includes X, * and lowercase: a k-mer spanning an unknown residue is not a k-mer, and mapping it to a wildcard group would invent a count that was never observed.

Parameters:

groups (dict[str, int])

Return type:

tuple[ndarray, int]

vdjtools.features.kmer_space.pattern_of(spec)

Parse a k-mer shape: "xxx" is an ungapped 3-mer, "xx.x" a gapped 3-mer over 4.

x keeps a position, . (or _) skips it. A pattern must start and end with x – a leading or trailing gap is the same feature as the shorter pattern, shifted, and admitting both would put two names on one column.

Parameters:

spec (str)

Return type:

list[int]

class vdjtools.features.kmer_space.KmerSpace(pattern, groups, v_genes, flank, codes, idf, components=None, meta=<factory>, _lookup=None)

A frozen (V, k-mer) vocabulary with its IDF and optional SVD basis.

Everything needed to turn a repertoire into the same vector next year, in another lab. Fitted once by fit_kmer_space(); applied by transform().

Parameters:

• pattern (list[int])

• groups (dict[str, int])

• v_genes (list[str])

• flank (int)

• codes (ndarray)

• idf (ndarray)

• components (ndarray | None)

• meta (dict)

• _lookup (ndarray | None)

pattern: list[int]

groups: dict[str, int]

v_genes: list[str]

flank: int

codes: ndarray

idf: ndarray

components: ndarray | None = None

meta: dict

property n_columns: int

property n_components: int

lookup()

code -> column index (-1 outside the vocabulary), built once and reused.

Held as one int32 buffer and passed to the kernel by pointer. Materialised lazily because it is the size of the whole code space – ~8e6 int32 = 33 MB for a V-pinned 4-mer over 20 residues – which is cheap once per corpus and ruinous once per sample.

Return type:

ndarray

counts(df, *, weight='freq')

Raw clone-weighted counts on this vocabulary, one entry per column.

Parameters:

• df (DataFrame)

• weight (str)

Return type:

ndarray

transform(df, *, weight='freq', sublinear=True, residual=False)

TF-IDF vector for one repertoire, L2-normalised, projected if a basis was fitted.

Parameters:

• df (DataFrame) – One repertoire’s clonotype frame, already filtered. Pass the work frame if the caller uses one – weight_expr reads whatever frequency holds.

• weight (str) – Clone weight ladder, as elsewhere in vdjtools.

• sublinear (bool) – log1p the term frequency before IDF. On by default: a clone expanded 1000-fold is not 1000 times more informative about which k-mers a repertoire contains, and without it one clonal expansion dominates the whole vector.

• residual (bool) – Append ||x - V^T V x||, the norm of what the basis discarded. Only meaningful when a basis was fitted.

Returns:

(n_components [+1],) if a basis was fitted, else (n_columns,).

Return type:

ndarray

vdjtools.features.kmer_space.fit_kmer_space(frames, *, pattern='xxxx', n_groups=8, flank=4, v_genes=None, weight='freq', min_df=0.02, max_df=0.8, max_columns=50000, n_components=32, sublinear=True, threads=0)

Fit the vocabulary, IDF and SVD basis on a reference corpus.

This is the only corpus-fitted part of the k-mer block, and it is fitted once, on studies disjoint from anything it will later be scored on. Nothing here is re-estimated per cohort; that is the whole point.

Parameters:

• frames – Iterable of per-sample clonotype frames (one locus).

• pattern (str) – Shape spec – "xxx"/"xxxx" ungapped, "xx.x"/"x.x.x" gapped.

• n_groups (int) – BLOSUM62 alphabet groups; 20 disables the reduction.

• flank (int) – Residues trimmed from each end. 4 is the measured germline core (an N-terminal 4-mer is shared by 31% of clonotypes, a central one by 0.08%).

• v_genes – V vocabulary; inferred from the corpus when None.

• min_df (float) – Drop k-mers seen in fewer than this fraction of samples – too rare to estimate.

• max_df (float) – Drop k-mers seen in more than this fraction – present everywhere is no contrast, which for a CDR3 core means germline or near it.

• max_columns (int) – Keep at most this many columns, highest document frequency first among the survivors. A cap, and reported when it binds rather than silently truncating.

• n_components (int) – SVD rank; 0 ships the TF-IDF columns unprojected.

• sublinear (bool) – log1p the term frequency (see KmerSpace.transform()).

• threads (int) – Native worker threads across samples; 0 = auto.

• weight (str)

Returns:

A fitted KmerSpace.

Return type:

KmerSpace

vdjtools.features.kmer_space.save_kmer_spaces(spaces, path)

Freeze a per-locus set of spaces to one .npz.

A space that cannot be written down is not frozen, whatever the docstring says: the whole claim is that a collaborator recomputes the same columns, which requires the vocabulary, the IDF and the basis to travel. Stored per locus with the locus in the key.

The alphabet is stored as the group id of each of the 20 residues in AMINO_ACIDS order, not as the n_groups that produced it – a future change to the clustering must not silently re-partition an already-frozen space.

Parameters:

spaces (dict[str, KmerSpace])

Return type:

None

vdjtools.features.kmer_space.load_kmer_spaces(path)

Read back what save_kmer_spaces() wrote.

Return type:

dict[str, KmerSpace]