API reference#

Every public module is documented below, grouped by subpackage. See the user guide for runnable examples.

mir#

`mir`#

Package root: version, get_resource_path, re-exported TCREmp / PairedTCREmp.

mirpy — ML-oriented embeddings for immune receptor repertoires (TCR/BCR).

Version 3 is a slim, embedding-focused rewrite. Heavy machinery is delegated to sibling packages rather than duplicated:

alignment → seqtree (seqtree.gapblock),
VDJ-rearrangement / Pgen model + sampling → vdjtools.model (core dependency),
VDJdb annotation / E-values → vdjmatch (extra [annotate]),
build-time germline region annotation → arda (extra [build]).

The classical repertoire toolkit lives on the legacy-v2 branch (mirpy-lib 2.x).

mir.get_resource_path(name=None)[source]#

Return the absolute path to a bundled resource under mir/resources.

Parameters:: name (str | None) – Resource file or directory name. When None, return the sorted list of top-level resource names instead of a path.
Returns:: The absolute path as a string, or the sorted list of resource names when name is None.
Raises:: FileNotFoundError – If name does not resolve to a bundled resource.

Clonotype embedding (`mir.embedding`)#

`mir.embedding.prototypes`#

Bundled prototype loader + manifest.

Reference prototype loading for TCREMP embeddings.

Prototype files are generated once by mir/resources/prototypes/generate_prototypes.py and committed to the repository. They are never regenerated at build or test time. Each of the nine bundled files holds N_PROTOTYPES = 10 000 rows of real receptors — a uniform random sample (fixed seed=42) of germline-resolvable, unique, productive clonotypes from arda-annotated real repertoires, de-duplicated on the (v_call, j_call, junction_aa) triple. Real repertoires, not model draws: synthetic P_gen junctions have degenerate lengths and embed measurably worse (see SOURCES.md).

The default set. replicate=0 — the first n rows — is the prototype set. It is what every default, preset, bundled codec and published number uses; leave it alone unless you are deliberately measuring prototype-draw sensitivity.

Replicates. Because the file order is itself a uniform shuffle, any disjoint block of rows is an independent draw from the same pool. replicate=r returns block r (rows r*n … (r+1)*n), so replicates are disjoint — independent by construction, not a with-replacement bootstrap. You get n_replicates() = 10000 // n of them: 10 at n=1000, 5 at n=2000. Use them to answer “how much does my result depend on which prototypes I drew?”:

from mir.embedding.tcremp import TCREmp

scores = [score(TCREmp.from_defaults("human", "TRB", 1000, replicate=r).embed(df))
          for r in range(n_replicates("human", "TRB", 1000))]      # spread = draw sensitivity

Every replicate is a different coordinate system: distances within one are comparable, distances across two are not. The prototype hash covers replicate, so codecs, RepertoireSpace and DonorCohort all refuse to mix them — compare summary statistics across replicates, never raw embeddings.

Typical usage:

from mir.embedding.prototypes import load_prototypes, n_replicates

df = load_prototypes("human", "TRB")                        # all 10 000
df_small = load_prototypes("human", "TRB", n=1000)           # the default 1 000
df_alt = load_prototypes("human", "TRB", n=1000, replicate=3)  # an independent 1 000

mir.embedding.prototypes.N_PROTOTYPES: int = 10000#: Maximum number of prototypes per species/locus combination.

mir.embedding.prototypes.list_available_prototypes()[source]#

Return sorted (species, locus) pairs with a bundled prototype file.

Example

>>> ("human", "TRB") in list_available_prototypes()
True

Return type:: list[tuple[str, str]]

mir.embedding.prototypes.n_replicates(species, locus, n)[source]#

How many disjoint prototype replicates of size n the bundled pool supports.

pool_size // n — 10 at the common n=1000, 5 at n=2000. Replicate indices are range(n_replicates(...)), with 0 the default set.

Example

>>> n_replicates("human", "TRB", 1000)
10

Parameters:

species (str)
locus (str)
n (int)

Return type:

int

mir.embedding.prototypes.load_prototypes(species, locus, n=None, replicate=0)[source]#

Load prototype clonotypes for TCREMP distance-based embeddings.

Reads the pre-generated TSV in mir/resources/prototypes/ and returns n rows in the fixed generation order, so embeddings are reproducible across calls and machines.

Parameters:

species (str) – Species identifier; aliases ('hsa', 'Homo sapiens', …) are resolved via mir.aliases.normalize_species_alias().
locus (str) – Receptor locus; aliases ('beta', 'T_beta', …) are resolved via mir.aliases.normalize_locus_alias().
n (int | None) – Number of prototypes. Must be <= N_PROTOTYPES. Loads the whole pool when None.
replicate (int) – Which disjoint block of n rows to take. 0 (default) is the canonical prototype set — the first n rows, what every preset and published number uses. r > 0 returns rows r*n … (r+1)*n, an independent draw from the same pool (the file order is a uniform shuffle), for measuring prototype-draw sensitivity. Requires an explicit n; see n_replicates() for the number available.

Returns:

DataFrame with columns v_call, j_call, junction_aa in fixed generation order.

Raises:

ValueError – If n is not in [1, N_PROTOTYPES], or replicate is negative, or replicate is given without n, or the pool has no such block.
FileNotFoundError – If no prototype file exists for the species/locus.

Return type:

DataFrame

Example

>>> df = load_prototypes("human", "TRB", n=100)
>>> df.columns
['v_call', 'j_call', 'junction_aa']
>>> len(df)
100
>>> alt = load_prototypes("human", "TRB", n=100, replicate=1)
>>> set(df["junction_aa"]) & set(alt["junction_aa"])   # disjoint blocks
set()

`mir.embedding.tcremp`#

TCREmp / PairedTCREmp — the prototype (TCREMP) embedding.

TCREMP: distance-vector embeddings for TCR/BCR clonotypes.

TCREMP (T-Cell Receptor EMbedding with Prototypes) embeds each clonotype as a flat vector of distances to a fixed set of prototype clonotypes. For clonotype i and prototype k, three distances are emitted; the first two depend on the mode:

mode="vjcdr3" (default): V-germline, J-germline, junction.
mode="cdr123": CDR1, CDR2 (both germline V-gene-determined), junction.

The embedding vector (vjcdr3 shown) is:

[v_i1, j_i1, junc_i1, v_i2, j_i2, junc_i2, ..., v_iK, j_iK, junc_iK]

for K prototypes; the output matrix is (n_clonotypes, 3*K) float32.

Junction (CDR3) distances are computed with seqtree.gapblock.score_matrix() (best of contiguous gap-block placements at (3, 4, -4, -3), BLOSUM62 Gram penalty). The result d is a negative-type semimetric — a squared Hilbert distance d = ‖ψ(a)−ψ(b)‖² (symmetric, d(a,a)=0, non-negative) — not itself a metric; the induced genuine metric is ρ = √d (see mir/distances/junction.py). V/J and CDR1/CDR2 distances are looked up from baked matrices (mir.distances.germline.GermlineDistances), same Gram construction. Input is a polars DataFrame with the vdjtools.io.schema columns v_call, j_call, junction_aa.

Coordinate-system options (all default to the published v3 space): metric="sqrt" maps every block to the metric ρ=√d (benchmarked a wash, so "squared" is default); matrix= a custom seqtree.SubstitutionMatrix and alignment="sw" (paper-exact Smith-Waterman, slow) reshape the junction block — see mir.distances.junction.junction_distance_matrix().

Typical usage:

import polars as pl
from mir.embedding.tcremp import TCREmp

model = TCREmp.from_defaults("human", "TRB", n_prototypes=1000)
df = pl.DataFrame({"v_call": ["TRBV10-3*01"], "j_call": ["TRBJ2-7*01"],
                   "junction_aa": ["CASSIRSSYEQYF"]})
X = model.embed(df)          # shape: (1, 3000)

class mir.embedding.tcremp.TCREmp(species, locus, prototypes, germline, mode='vjcdr3', gap_positions=(3, 4, -4, -3), threads=0, metric='squared', matrix=None, alignment='gapblock', replicate=0)[source]#

Bases: object

Prototype distance-vector embedding for one locus.

Parameters:

species (str)
locus (str)
prototypes (pl.DataFrame)
germline (GermlineDistances)
mode (str)
gap_positions (tuple[int, ...])
threads (int)
metric (str)
alignment (str)
replicate (int)

classmethod from_defaults(species, locus, n_prototypes=None, mode='vjcdr3', replicate=0, **kwargs)[source]#

Build from the bundled prototypes and baked germline distances.

Parameters:

species (str) – Species identifier (aliases resolved).
locus (str) – Receptor locus (aliases resolved).
n_prototypes (int | None) – Prototype count; None uses the per-chain recommended preset (mir.embedding.presets.get_preset()).
mode (str) – "vjcdr3" (default) or "cdr123".
replicate (int) – Which bundled prototype draw to use. 0 (default) is the prototype set — every preset, bundled codec and published number uses it. r > 0 is an independent disjoint draw of the same size, for measuring how much a result depends on the prototype sample (see mir.embedding.prototypes.n_replicates()). Embeddings from different replicates are not mutually comparable.

Return type:

TCREmp

classmethod from_file(prototypes_path, species, locus, mode='vjcdr3', **kwargs)[source]#

Build from a custom prototype TSV (v_call/j_call/junction_aa).

Embeddings from a custom prototype set are not comparable to those from the bundled prototypes — the coordinate system differs.

Parameters:

prototypes_path (str | Path)
species (str)
locus (str)
mode (str)

Return type:

TCREmp

property n_features: int#

embed(clonotypes)[source]#

Embed a clonotype frame into (n_clonotypes, 3 * n_prototypes) float32.

Parameters:: clonotypes (DataFrame) – polars DataFrame with columns v_call, j_call, junction_aa (vdjtools.io.schema names).
Returns:: float32 array of shape (len(clonotypes), 3 * n_prototypes), interleaved per prototype as [slot0, slot1, junction].
Return type:: ndarray

class mir.embedding.tcremp.PairedTCREmp(chains)[source]#

Bases: object

Paired-chain embedding: concatenate two per-chain TCREmp embeddings.

Parameters:: chains (dict[str, TCREmp])

classmethod from_defaults(species, loci=('TRA', 'TRB'), n_prototypes=None, mode='vjcdr3', replicate=0, **kwargs)[source]#

Parameters:

species (str)
loci (tuple[str, str])
n_prototypes (int | None)
mode (str)
replicate (int)

Return type:

PairedTCREmp

property n_features: int#

embed(frames)[source]#

Embed per-chain frames (keyed by locus, row-aligned) and concatenate.

Parameters:: frames (dict[str, DataFrame]) – Mapping locus -> DataFrame; must cover every locus in the model, all with the same row count (row i is the same clonotype).
Return type:: ndarray

`mir.embedding.pca`#

PCA denoising / de-redundancy of embeddings.

PCA denoising of TCREMP embeddings.

The raw prototype-distance coordinates are highly redundant — neighbouring prototypes and correlated V/J scores span a low-rank subspace (Theory T3). A StandardScaler → PCA step smooths the embedding, removes that redundancy, and yields a compact basis that is well behaved for clustering / visualization (the paper uses n_components=50).

mir.embedding.pca.pca_denoise(X, n_components=50, scale=True, random_state=0)[source]#

Standardize then PCA-reduce a TCREMP embedding.

Parameters:

X (ndarray) – (n_clonotypes, n_features) embedding.
n_components (int) – Target dimensionality (clamped to min(n_samples, n_features)).
scale (bool) – Standardize each coordinate to zero mean / unit variance first.
random_state (int) – Seed for the (randomized) SVD solver.

Returns:

(n_clonotypes, k) float array of principal-component scores, k = min(n_components, n_samples, n_features).

Return type:

ndarray

`mir.embedding.presets`#

Per-chain recommended prototype counts + PCA dims.

Recommended per-chain presets: prototype count and PCA dimensionality.

Data-driven from the bundled prototypes (measured by the benchmarks/ scripts in the companion 2026-mirpy-analysis repo):

n_prototypes — the pairwise-distance geometry saturates (Pearson r vs a 5000-prototype reference ≥ 0.996 at 1000, ≥ 0.999 at 2000). Compact chains (IGK/IGL/TRG) saturate earlier; diverse chains (IGH/TR*) benefit from 2000. The paper’s 3000 is safe but generous.
n_components — PCA dimensions retaining ~95% variance (the paper’s clustering regime).
n_components_recon — PCA dimensions retaining ~99% variance, needed to reconstruct the longer, more diverse junctions with the inverse codec (Theory T5, arda coords: IGH exact-match 0.115 at 95% → 0.356 at 99%). Compact chains barely differ; IGH/TRD/TRA need many more.

Values are rounded up to the nearest 5; treat them as starting points, not hard limits.

class mir.embedding.presets.ChainPreset(n_prototypes, n_components, n_components_recon)[source]#

Bases: object

Recommended settings for one (species, locus).

Parameters:

n_prototypes (int)
n_components (int)
n_components_recon (int)

n_prototypes: int#

n_components: int#

n_components_recon: int#

mir.embedding.presets.get_preset(species, locus)[source]#

Return the ChainPreset for a (species, locus), aliases resolved.

Falls back to a generic diverse-chain preset for unmeasured combinations.

Example

>>> get_preset("human", "IGH").n_components_recon
300
>>> get_preset("hsa", "beta").n_prototypes
2000

Parameters:

species (str)
locus (str)

Return type:

ChainPreset

Distances (`mir.distances`)#

`mir.distances.junction`#

Junction distance via seqtree.gapblock (metric / matrix / alignment knobs).

Junction (CDR3) distance matrices.

The default backend is seqtree.gapblock — a fast contiguous-gap-block approximation to Smith-Waterman. seqtree’s blosum62() substitution matrix is already the Gram transform s_aa + s_bb − 2·s_ab, so gapblock.score_matrix returns the squared dissimilarity d directly (d(a,a)=0, symmetric, non-negative), selecting the best of contiguous gap-block placements at gap_positions.

Three knobs, all defaulting to the published v3 coordinate system:

alignment — "gapblock" (default, fast, ~5·10⁸ pairs/s) or "sw" (paper-exact Smith-Waterman via BioPython; O(n·K) pairwise, so validation / small-scale only). gap-block approximates SW in the near-neighbour regime that clustering and density use: close pairs agree exactly (e.g. CASSTTGLNTEAFF/CASSLTAMNTEAFF → 26 under both). They part company on distant pairs, because "sw" is a local alignment that simply stops extending while gap-block scores the whole junction — over all 44 850 pairs of the 300 bundled human-TRB prototypes the two rank-correlate ρ≈0.78, with equal-length outliers as far apart as 181 vs 77. The two are therefore different coordinate systems, not interchangeable scales: use "sw" to reproduce the paper’s distances, "gapblock" (default) for everything else, and never mix them in one embedding.
matrix — the substitution matrix for the gap-block backend: any seqtree.SubstitutionMatrix (blosum62() default, plus pam250(), structural(), unit(), or from_similarity for a fully custom matrix). Changing it is a coordinate-system change (the baked germline blocks stay BLOSUM62, so mix scales only deliberately).
metric — "squared" (default) returns the Gram dissimilarity d; "sqrt" returns the induced metric ρ=√d. d is a squared Hilbert distance, so it is not itself a metric (it violates the triangle inequality); ρ is. Benchmarked a wash vs d, so d is the default (the published v3 space). The gap-block placement is a monotone argmin, identical under d and √d; the sqrt is applied to the chosen value.

mir.distances.junction.junction_distance_matrix(queries, refs, gap_positions=(3, 4, -4, -3), threads=0, *, metric='squared', matrix=None, alignment='gapblock')[source]#

Return the (len(queries), len(refs)) junction distance matrix (float32).

Parameters:

queries – Query junction (CDR3) amino-acid strings.
refs – Reference/prototype junction strings.
gap_positions (tuple[int, ...]) – Candidate contiguous gap-block placements; best is chosen (gap-block only).
threads (int) – Worker threads (0 = all cores, GIL released; gap-block only).
metric (str) – "squared" (Gram dissimilarity d, default) or "sqrt" (metric ρ=√d).
matrix – A seqtree.SubstitutionMatrix for the gap-block backend; None = blosum62().
alignment (str) – "gapblock" (default, fast) or "sw" (paper-exact Smith-Waterman, slow).

Return type:

ndarray

`mir.distances.germline`#

Resource-backed V / J / CDR1 / CDR2 germline distance lookup.

Resource-backed germline distance lookup for TCREMP embeddings.

Loads the per-locus distance matrices baked by mir/resources/germline_dist/build_germline_dist.py and gathers, for a batch of query alleles against a fixed set of prototype alleles, the (n_query, n_proto) distance matrix for a component (V/J/CDR1/CDR2).

Runtime is pure numpy — no BioPython, no gene library. Unknown alleles resolve via the cascade exact allele → ``*01`` → bare gene → fallback (the fallback is the maximum observed distance for that component, so an unknown gene sits maximally far from every prototype).

class mir.distances.germline.GermlineDistances(species, locus, components)[source]#

Bases: object

Baked germline distances for one (species, locus).

Example

>>> gd = GermlineDistances.load("human", "TRB")
>>> D = gd.matrix("V", ["TRBV20-1*01"], ["TRBV20-1*01", "TRBV6-5*01"])
>>> D.shape
(1, 2)
>>> float(D[0, 0])   # identical gene -> 0
0.0

Parameters:

species (str)
locus (str)
components (dict[str, _Component])

classmethod load(species, locus)[source]#

Parameters:

species (str)
locus (str)

Return type:

GermlineDistances

has(component)[source]#

Parameters:: component (str)
Return type:: bool

matrix(component, query_alleles, proto_alleles)[source]#

Return the (len(query), len(proto)) germline distance matrix.

Parameters:

component (str) – One of V/J/CDR1/CDR2.
query_alleles – Iterable of query allele strings.
proto_alleles – Iterable of prototype allele strings.

Raises:

KeyError – If component is not available for this locus.

Return type:

ndarray

mir.distances.germline.load_germline_distances(species, locus)[source]#

Cached GermlineDistances.load(), keyed by canonical species/locus.

Parameters:

species (str)
locus (str)

Return type:

GermlineDistances

Density and background subtraction (`mir.density`)#

Graph-free continuous-density TCRNET/ALICE neighbour enrichment (Theory T6).

Continuous-density TCRNET / ALICE: enrichment in TCREMP embedding space (Theory T6).

TCRNET (Pogorelyy & Shugay, Front Immunol 2019) and ALICE (Pogorelyy et al., PLoS Biol 2019) flag antigen-driven / convergent-selected clonotypes by neighbour enrichment: count a clonotype’s near-identical (Hamming-1) neighbours and compare against a background — OLGA generation probability for ALICE, a control repertoire for TCRNET. Both are graph methods (a sequence trie).

This module reimplements the same enrichment test with neighbour counting in the TCREMP embedding space instead of on a trie — the graph-free density ratio the theory specifies (Theory T6):

E(z) = f_obs(z) / f_gen(z), f_gen = φ_# P_gen

For an observed clonotype embedded at z_i we count observed neighbours n_obs and background neighbours n_bg within a radius r (calibrated to one CDR3 substitution, so the continuous test approximates the discrete Hamming-1 one). Under the null “the observed repertoire is a background sample”, each of the other N_obs − 1 observed clonotypes lands within r of z_i with probability p_bg = n_bg / N_bg, so the expected count is (N_obs−1)·p_bg and

p-value = poisson.sf(n_obs − 1, expected) # ALICE-style (test=”poisson”)
= binom.sf(n_obs − 1, N_obs−1, p_bg) # TCRNET-style (test=”binomial”)

with Benjamini-Hochberg q-values. The fold enrichment fold = n_obs / expected is the density ratio E(z) itself. The background is either generated from the vdjtools P_gen model (generate_background(), the ALICE analog) or any supplied control repertoire (the TCRNET analog) — both just become bg_df.

Passing abundance= (clone sizes) to neighbor_enrichment() adds the clonal-depth channel (appendix §T.6 sec:dens-abund): the distinct in-ball count becomes a variance-stabilised weighted mass S=Σ g(a_j) with concave g (default log(1+a)) so a hyperexpanded clone can’t dominate, tested against a compound-Poisson Gamma tail, plus an orphan size-test P(A≥a_j) combined with breadth by Fisher. Default (abundance=None) counts each clonotype once — the shipped g≡1 behaviour.

Observed and background embeddings are comparable only in one coordinate system (same prototypes and same PCA rotation). fit_density_space() enforces that: it embeds both through one TCREmp and fits one PCA on the pooled matrix, returning a DensitySpace that projects any further frame (a control, or the radius-calibration mutants) into the same basis.

Typical usage:

from mir.density import fit_density_space, calibrate_radius, neighbor_enrichment, enriched_mask
from mir.embedding.tcremp import TCREmp
from mir.embedding.presets import get_preset

model = TCREmp.from_defaults("human", "TRB")
nc = get_preset("human", "TRB").n_components
space, obs_emb, bg_emb = fit_density_space(model, obs_df, bg_df, n_components=nc)
r = calibrate_radius(space)                       # radius ≈ one CDR3 substitution
res = neighbor_enrichment(obs_emb, bg_emb, r)     # E(z), p, q per observed clonotype
hits = obs_df.filter(enriched_mask(res))          # significantly enriched clonotypes

class mir.density.DensitySpace(model, space, scaler, pca)[source]#

Bases: object

A single fitted TCREMP → PCA coordinate system for density comparison.

Holds the embedding model plus the standardization + PCA fitted on the pooled observed+background matrix. transform() projects any frame into this same basis so that observed, background, control and calibration points are all directly comparable — the invariant that makes E(z) meaningful.

Parameters:

model (object)
space (str)
scaler (StandardScaler)
pca (PCA)

model: object#

space: str#

scaler: StandardScaler#

pca: PCA#

transform(df)[source]#

Embed df and project it into this fitted coordinate system.

Parameters:: df (DataFrame)
Return type:: ndarray

class mir.density.EnrichmentResult(n_obs, n_bg, expected, fold, pvalue, qvalue, radius, score=None, pvalue_breadth=None, pvalue_size=None)[source]#

Bases: object

Per-observed-clonotype neighbour enrichment (input order preserved).

Parameters:

n_obs (ndarray)
n_bg (ndarray)
expected (ndarray)
fold (ndarray)
pvalue (ndarray)
qvalue (ndarray)
radius (float)
score (ndarray | None)
pvalue_breadth (ndarray | None)
pvalue_size (ndarray | None)

n_obs: ndarray#

n_bg: ndarray#

expected: ndarray#

fold: ndarray#

pvalue: ndarray#

qvalue: ndarray#

radius: float#

score: ndarray | None = None#

pvalue_breadth: ndarray | None = None#

pvalue_size: ndarray | None = None#

mir.density.fit_density_space(model, obs_df, bg_df, *, n_components, space='full', seed=0, pca_fit_cap=None, chunk_size=None)[source]#

Embed observed + background frames into one shared PCA coordinate system.

Both frames are embedded through the same model and a single PCA (fit on the pooled matrix) is applied to both, guaranteeing the two point clouds live in one comparable basis — the invariant that makes E(z) meaningful.

Parameters:

model – A fitted single-chain TCREmp.
obs_df (DataFrame) – Observed clonotypes (v_call, j_call, junction_aa).
bg_df (DataFrame) – Background clonotypes (generated P_gen sample or a control repertoire).
n_components (int) – PCA dimensionality (clamped to min(n_samples, n_features)); use get_preset(species, locus).n_components.
space (str) – "full" (whole V+J+junction vector — V/J restriction is soft/continuous) or "junction" (CDR3 sub-block only, mimicking pure-CDR3 TCRNET/ALICE and a third of the memory at whole-repertoire scale).
seed (int) – RNG / PCA solver seed.
pca_fit_cap (int | None) – Fit the StandardScaler + PCA on at most this many randomly-sampled pooled rows (then transform all rows). None fits on everything. NB on its own this caps the fit, not the memory: without chunk_size the full raw matrix for both frames is materialized before the PCA is fitted at all.
chunk_size (int | None) – Embed and project in batches of this many rows, so the full raw matrix is never resident (peak ≈ chunk_size × n_features instead of len(df) × n_features). Required for whole-cohort clouds: 4.2M clonotypes × 1000 prototypes is ~51 GB raw, versus ~2.4 GB at chunk_size=200_000. When set and pca_fit_cap is None, the fit is capped at 200k pooled rows — fitting on everything would defeat the purpose. None (default) keeps the original single-shot path exactly.

Returns:

(density_space, obs_emb, bg_emb) — the fitted DensitySpace and the two reduced float arrays, row-aligned to obs_df / bg_df.

Return type:

tuple[DensitySpace, ndarray, ndarray]

mir.density.calibrate_radius(space, *, sample_df=None, sample=2000, seed=0, quantile=0.5)[source]#

Radius corresponding to one CDR3 substitution, in space’s coordinate system.

Mutates one interior residue of each sampled junction and measures the embedding drift ‖φ(seq) − φ(mutated)‖ through the same fitted transform used for enrichment (reusing the Theory-T5 SHM-drift idea). The requested quantile of that drift is the neighbourhood radius, so the continuous test approximates the discrete Hamming-1 one. V/J are held fixed, isolating the junction contribution.

Parameters:

space (DensitySpace) – A DensitySpace from fit_density_space().
sample_df (DataFrame | None) – Frame to draw calibration sequences from; defaults to the model’s own prototype set (the coordinate anchor).
sample (int) – Number of prototypes to use when sample_df is None.
seed (int) – RNG seed for the substitutions.
quantile (float) – Drift quantile to return (0.5 = median).

Returns:

The calibrated radius (a positive float).

Return type:

float

mir.density.neighbor_enrichment(obs_emb, bg_emb, radius=None, *, lambda0=3.0, test='poisson', pseudocount=1.0, calibrate='median', abundance=None, weight='log1p', orphan=True, backend='kdtree', k_max=96, seed=0)[source]#

Continuous neighbour-enrichment test in one embedding coordinate system.

Two neighbourhood-scale modes (appendix §T.6):

balloon (default, radius=None) — adaptive bandwidth. Each observed point’s radius is the distance to its k-th background neighbour, with k chosen so the expected background occupancy equals lambda0. The adaptive radius encodes the P_gen density shape (small where generation is dense, large where sparse); it is robust to the distance concentration that makes a single global radius fragile in high dimensions, and keeps the test well-powered everywhere (minimum detectable E ≳ 1 + c/√lambda0).
fixed (radius given) — one global radius (e.g. from calibrate_radius()).

A real repertoire’s background is systematically denser than a generative P_gen model (gene usage, sampling structure, thymic selection) — the “water level” π of §T.6. calibrate="median" rescales the null so the bulk of clones sits at fold ≈ 1 and only clones exceeding the repertoire’s own typical local density are called. Signal (< 5 % in naive repertoires) barely moves the median, so this is robust; pass calibrate=None to test against the raw P_gen null instead.

Clonal abundance (appendix §T.6 sec:dens-abund). By default each clonotype counts once (abundance=None), scoring only convergence breadth. Supplying per-clonotype clone sizes abundance adds the depth channel: the in-ball count is replaced by a variance-stabilised weighted mass S(z) = Σ g(a_j) over neighbours, with g non-decreasing and concave (weight="log1p" g=log(1+a) or "anscombe" g=√(a+3/8)) so a hyperexpanded clone contributes a bounded increment, not a thousand-fold one. Under H0 the mass is compound-Poisson with dispersion index φ=E[g²]/E[g] (Prop. abund), so S is tested against a moment-matched Gamma(μ_S/φ, φ) upper tail (which collapses to the Poisson count test when g≡1). With orphan=True each clone additionally gets a size p-value P(A≥a_j) and the breadth/depth channels are combined by Fisher (χ²₄), recovering a hyperexpanded orphan (depth, no breadth) that the count test alone misses.

Parameters:

obs_emb (ndarray) – (N_obs, d) observed embedding.
bg_emb (ndarray) – (N_bg, d) background embedding, same basis as obs_emb (produce both with fit_density_space()). Use M ≥ 5N for a stable ratio.
radius (float | None) – Fixed neighbourhood radius; None selects the balloon estimator.
lambda0 (float) – Target expected background occupancy for the balloon estimator ([1, 5]).
test (str) – "poisson" (ALICE analog) or "binomial" (TCRNET control analog). Applies to the distinct-count path; the abundance-weighted path always uses the Gamma tail.
pseudocount (float) – Added to the background count to stabilize p_bg.
calibrate (str | None) – "median" empirical-null water-level calibration (default), or None.
abundance (ndarray | None) – Optional (N_obs,) clone sizes (e.g. duplicate_count), row-aligned to obs_emb; enables the weighted/orphan channels.
weight (str) – Concave size transform g — "log1p" (default), "anscombe", or "distinct" (g≡1, ignore sizes even if abundance is given).
orphan (bool) – When abundance-aware, combine the size p-value with the breadth p-value by Fisher.
backend (str) – neighbour engine. "kdtree" (default, scipy cKDTree — exact, multithreaded, 5–9× faster than BallTree), "exact" (BallTree — exact, single-threaded; the historical default, kept for bit-reproducing the BallTree baseline), or "ann" (approximate pynndescent kNN for whole-repertoire scale — ~30× faster at ≥40k but recall < 1 undercounts, biasing enrichment conservatively; see _ann_neighbors()). "kdtree" is bit-identical to "exact" at a fixed radius; in balloon mode counts differ by at most ±1 at the ball boundary (float-epsilon in the computed k-th-neighbour distance), which is negligible — pass backend="exact" only to reproduce the BallTree baseline exactly.
k_max (int) – backend="ann" only — kNN-graph degree. A clonotype whose ball holds all k_max neighbours is undercounted (and warned about); raise this when that warning fires.
seed (int) – backend="ann" only — pynndescent random_state.

Returns:

An EnrichmentResult; fold is the (calibrated) density ratio E(z). In balloon mode radius is the median adaptive radius used. When abundance-aware, score holds the weighted mass S and pvalue is the combined breadth×depth test.

Return type:

EnrichmentResult

mir.density.enriched_mask(res, *, alpha=0.05, min_fold=1.0, min_neighbors=2)[source]#

Boolean hit mask: q < alpha & fold > min_fold & ≥ min_neighbors (self + neighbours).

min_neighbors=2 reproduces the legacy “self plus at least one neighbour” criterion (n_obs excludes self, so the threshold is n_obs ≥ min_neighbors − 1).

For an abundance-aware result (res.score set) the min_fold/min_neighbors breadth gates are dropped: the combined breadth×depth q already governs, so a hyperexpanded orphan (significant depth, no neighbours) is kept rather than filtered out.

Parameters:

res (EnrichmentResult)
alpha (float)
min_fold (float)
min_neighbors (int)

Return type:

ndarray

mir.density.exposure_score(res, *, abundance=None, alpha=0.05, min_fold=1.0, min_neighbors=2)[source]#

Aggregate one repertoire’s per-clonotype EnrichmentResult into 3 exposure scalars.

neighbor_enrichment() answers “which clones are enriched” (a per-clonotype hit list, the TCRNET/ALICE readout); this answers “how exposed does the whole repertoire look” — a sample-level summary suitable as an explainable channel (mir.explain) alongside identity/diversity/ coverage, rather than a clonotype list. Exposure detection moves from clone-level (density.py’s core job) to repertoire-level (a scalar a cohort model can score) with no new statistics: it is a read-off of the enrichment already computed.

Parameters:

res (EnrichmentResult) – One donor’s enrichment result (obs_emb = that donor’s repertoire, bg_emb shared across the cohort so results are comparable).
abundance (ndarray | None) – Optional clone sizes (duplicate_count), row-aligned to res — enables the mass-weighted column. None falls back to unweighted (distinct-count) breadth.
alpha (float) – Forwarded to enriched_mask().
min_fold (float) – Forwarded to enriched_mask().
min_neighbors (int) – Forwarded to enriched_mask().

Returns:

breadth — fraction of distinct clonotypes flagged enriched.
mass_fraction — abundance-weighted fraction of repertoire mass flagged (equals breadth when abundance is None).
mean_log2_fold — mean log2(fold) among flagged clones (0.0 if none flagged).

Return type:

[breadth, mass_fraction, mean_log2_fold]

mir.density.exposure_channel(results, *, abundances=None, **kwargs)[source]#

Stack exposure_score() across a cohort into one ready-to-fuse channel matrix.

Parameters:

results (list[EnrichmentResult]) – One EnrichmentResult per donor (same shared background embedding).
abundances (list[ndarray | None] | None) – Optional per-donor clone-size arrays, row-aligned to results (None per donor falls back to unweighted breadth for that donor).
**kwargs – Forwarded to exposure_score() (alpha, min_fold, min_neighbors).

Returns:

(n_donors, 3) array — columns [breadth, mass_fraction, mean_log2_fold] — pass directly to mir.explain.ChannelBuilder.add() (e.g. builder.add("exposure", exposure_channel(results))) or return it from a mir.cohort.fit_donor_embeddings() extra_channels closure.

Return type:

ndarray

mir.density.denoise_and_cluster(obs_emb, res, *, alpha=0.05, min_fold=1.0, min_neighbors=2, **cluster_kw)[source]#

Noise-filter then cluster: DBSCAN the enriched subset, non-enriched → label -1.

Background subtraction (enriched_mask()) removes the naive-repertoire noise, then mir.bench.metrics.cluster() groups the surviving hits into convergent motifs.

Returns:

(labels, mask) — labels is length N_obs (-1 for non-enriched or DBSCAN-noise points); mask is the enriched-hit boolean.

Parameters:

obs_emb (ndarray)
res (EnrichmentResult)
alpha (float)
min_fold (float)
min_neighbors (int)

Return type:

tuple[ndarray, ndarray]

mir.density.generate_background(locus, n, *, species='human', source='learned', seed=0, productive_only=True)[source]#

Draw n synthetic background clonotypes from the bundled vdjtools P_gen model.

Lazy vdjtools.model import (vdjtools is a core dependency; imported here only to defer the cost). Returns a frame with junction_aa, v_call, j_call — ready for TCREmp.embed — sampled from the vdjtools VDJ-rearrangement model, i.e. a Monte-Carlo estimate of the f_gen = φ_# P_gen pushforward (the ALICE analog background).

A biological control repertoire is the better background whenever you have one (see the module docstring): a P_gen sample flags the pervasive convergence of any real repertoire.

Parameters:

locus (str) – Bundled model locus (e.g. "TRB").
n (int) – Number of sequences to generate.
species (str) – "human" (default) or "mouse". Mouse requires source="arda" — it is the only bundled set with a non-human organism.
source (str) – Bundled model source — "learned" (vdjtools EM-inferred, the default), "olga" (legacy OLGA-parameter bootstrap, retained only for comparison), or "arda". Prefer "arda" when the observed data is arda-annotated: it is the only set in the arda IMGT allele namespace, the same frame as mirpy’s prototypes and baked germline distances, so generated V/J calls resolve exactly instead of taking the allele-cascade fallback.
seed (int) – Sampling seed.
productive_only (bool) – Reject out-of-frame / stop-codon rearrangements.

Raises:

ValueError – If the installed vdjtools predates the arda/organism support that source="arda" and non-human species need.

Return type:

DataFrame

Sample-level (repertoire) embedding (`mir.repertoire`)#

One vector per repertoire: RFF kernel mean ‖ Hill diversity ‖ second moment; MMD (Theory T7). Optionally a sub-probability measure: missing_mass (Good-Turing / bias-corrected Chao1), naive_reference (the germline location of the unseen) and contrast_embedding (Ψ = mass·(Φ − naive), signed and magnitude-carrying). Also the measure algebra around it: rao_q (Rao’s quadratic entropy = 1 − ‖Φ₁‖²), depth_threshold (the κ below which Φ is mostly sampling noise), sample_statistics / cohort_statistics, the compartment decomposition band_frames / band_embeddings / mixture_weights, and rarefy_embedding. Derivations: Mathematical foundations.

Sample-level (repertoire) embedding: one fixed vector per repertoire (Theory §T.7).

A repertoire is an order-invariant multiset of clonotypes with clone counts, S = {(σ, a_σ)}. We embed it as the empirical measure ρ_S = Σ_σ w_σ δ_{z_σ} on the clonotype embedding space (z_σ = φ(σ) from TCREmp), with clone-size weights w_σ = g(a_σ) / Σ_τ g(a_τ) (g = log2p1 — log2(1+a) — by default, so one hyperexpanded clone can’t dominate; "duplicate_count" weights linearly by clone size and "distinct" ignores size entirely, g≡1). Φ(S) is a sketch of ρ_S in three blocks, each owning one requirement (appendix §T.7):

mean Φ₁ = Σ_σ w_σ ψ(z_σ) — the random-Fourier-feature kernel mean embedding (order-invariant + depth-robust: converges to the population mean map at rate n_eff^{-1/2}, n_eff = (Σ w²)⁻¹). Distance ‖Φ₁(S)−Φ₁(S')‖ ≈ MMD. Codebook-free — no K, no clustering (Prop. prop:kme/prop:codebook).
diversity Φ₂ = {⁰D, ¹D, ²D} (Hill numbers), optionally coverage-standardized to a common Good–Turing coverage Ĉ* via vdjtools.stats.inext; n_eff is itself a Hill number, n_eff ∈ [²D, ⁰D] (Prop. prop:antag).
second Σ_S = Σ_σ w_σ ψ₂(z_σ) ψ₂(z_σ)ᵀ — the codebook-free Fisher vector carrying clonotype co-occurrence / HLA-linked public structure (Prop. prop:interact); computed on a small RFF so its upper triangle stays cheap.

The comparability invariant (as in mir.ml.bundle.CodecBundle and mir.density.fit_density_space()): every sample in a cohort must be embedded through one prototype set and one PCA + RFF basis, or the measures are incomparable. fit_repertoire_space() fits that basis once on the pooled clonotype cloud; RepertoireSpace serializes it and refuses a prototype-hash mismatch on load.

Sub-probability embeddings. Φ₁ above is the kernel mean of a probability measure — the weights are normalised to sum to 1, so every sample asserts one full unit of confidence. At RNA-seq depth that premise fails: a tissue sample holding 21 unique clonotypes has w_σ = 1/n for a technical draw size, not a clonal frequency, and normalising forces a 5-clonotype tumour to land somewhere arbitrary on the unit sphere. sample_embedding(missing_mass="turing"|"chao") instead estimates the probability mass of the clonotypes that were never drawn (missing_mass()) and records the retained mass on SampleEmbedding.mass ≤ 1 — a sub-probability measure, which keeps everything that makes Φ useful (‖Φ_P − Φ_Q‖ = MMD, convex combinations are real pooled repertoires) while a signed measure would not. The unseen block gets a principled location from naive_reference() (the germline recombination model, not the corpus centroid — measured worse), and contrast_embedding() composes the two into Ψ_S = mass·(Φ_S − naive), where magnitude = confidence × deviation-from-naive and an immune desert lands at the origin instead of being filtered out. See mir.bench.recovery_report() for the “is this object honest” check.

Torch-free (numpy / sklearn; vdjtools + the learned set-network in mir.ml are lazy).

Typical usage:

from mir.repertoire import fit_repertoire_space, sample_embedding, mmd_distance
from mir.embedding.tcremp import TCREmp
import polars as pl

model = TCREmp.from_defaults("human", "TRB", n_prototypes=1000)
pooled = pl.concat(samples)                    # all clonotypes across the cohort
space = fit_repertoire_space(model, pooled)    # ONE basis for the cohort
embs = [sample_embedding(space, s) for s in samples]
d = mmd_distance(embs[0], embs[1])             # ≈ MMD between the two repertoires

mir.repertoire.missing_mass(counts, method='chao')[source]#

Estimated probability mass M₀ of the clonotypes that were never drawn.

Both estimators feed the same identity with a different M₀; they rank samples nearly identically (measured r = 0.93 blood TRB, 0.76 tissue IGH) but differ in level — Chao runs above Turing when singletons dominate doubletons (f₁ ≫ 2f₂) and below it when doubletons are abundant. Measured means: blood TRB 0.552 (Turing) vs 0.649 (Chao), tissue IGH 0.189 vs 0.458:

none    M0 = 0                        the sample is treated as complete
turing  M0 = f1 / N                   Good-Turing, read-weighted
chao    M0 = S_u / (N + S_u)          S_u = f1(f1-1) / (2(f2+1))

f₁/f₂ are the singleton / doubleton clonotype counts and N the total reads. The chao form rests on one modelling assumption: an unseen clonotype is rare, so had it been drawn it would carry at most one read, and unseen clonotypes do not overlap. That pins each unseen clone’s frequency at the detection boundary 1/(N + S_u) and is what makes the units legal — N counts reads and S_u counts clonotypes, and adding them is only defensible because each unseen clone contributes exactly one read. The observed term then collapses to exactly a_σ/(N + S_u).

S_u is the bias-corrected Chao1 f1(f1-1)/(2(f2+1)), never the classical f1²/(2f2): at these depths a sample with no clonotype seen exactly twice is common, and the classical form is undefined there.

Parameters:

counts – Per-clonotype read counts (duplicate_count); rounded to integers to count singletons and doubletons.
method (str) – "none" | "turing" | "chao".

Returns:

M₀ ∈ [0, 1] — the missing (unobserved) mass. The retained mass is 1 − M₀.

Raises:

ValueError – If method is not recognised.

Return type:

float

Example

>>> missing_mass([5000, 3000, 2000], "chao")     # deep, no singletons
0.0

mir.repertoire.sample_statistics(sample_df)[source]#

The sampling fingerprint of one repertoire — the basic statistics, from counts alone.

Two uses. (a) The targets mir.bench.recovery_report() asks the embedding to carry: if a statistic is recoverable from Φ, nothing has to be bolted on beside it. (b) Features in their own right — the abundance classes f1/f2/f3plus and the top-clone fraction are clonality and expansion, i.e. candidate biology, not only depth nuisance.

Parameters:: sample_df (DataFrame) – One sample’s clonotypes with duplicate_count.
Returns:: {name: float} with n_reads, log_reads, richness, log_richness, f1, f2, f3plus, singleton_fraction, top_clone_fraction, shannon, missing_mass_turing, missing_mass_chao.
Raises:: ValueError – On an empty or all-zero-count repertoire (via _sample_weights()).
Return type:: dict

mir.repertoire.cohort_statistics(frames)[source]#

Per-statistic arrays over a cohort — sample_statistics() transposed.

Parameters:: frames – One clonotype frame per sample, in the same order as the embedding rows.
Returns:: {name: (n_samples,) array}, ready as mir.bench.recovery_report()’s stats.
Raises:: ValueError – If frames is empty.
Return type:: dict

class mir.repertoire.RandomFourierFeatures(omega, b, length_scale)[source]#

Bases: object

Random Fourier features for a Gaussian kernel k(z,z')=exp(−‖z−z'‖²/2ℓ²).

ψ(z) = √(2/D)·cos(Ωzᵀ + b) with Ω ~ N(0, 1/ℓ²), b ~ U[0, 2π], so that E[ψ(z)·ψ(z')] = k(z,z') (Rahimi–Recht 2007). Bandwidth ℓ is set to the one-substitution embedding scale r₁ so the kernel resolves ~one CDR3 mutation.

Parameters:

omega (ndarray)
b (ndarray)
length_scale (float)

omega: ndarray#

b: ndarray#

length_scale: float#

property dim: int#

transform(Z)[source]#

Map (n, p) embedding rows to (n, D) random features.

Parameters:: Z (ndarray)
Return type:: ndarray

class mir.repertoire.RepertoireSpace(clono, rff, rff2, meta, _naive=<factory>)[source]#

Bases: object

One fitted clonotype→PCA→RFF basis shared by a whole cohort.

clono projects a clonotype frame into the shared PCA coordinates; rff lifts those into the kernel-mean feature space (mean block); rff2 is a smaller RFF for the second-moment block. meta records the embedding-space identity (prototype hash + knobs) so load() can refuse an incomparable basis.

Parameters:

clono (DensitySpace)
rff (RandomFourierFeatures)
rff2 (RandomFourierFeatures | None)
meta (dict)
_naive (dict)

clono: DensitySpace#

rff: RandomFourierFeatures#

rff2: RandomFourierFeatures | None#

meta: dict#

transform_clonotypes(df)[source]#

Project a clonotype frame into the shared PCA coordinate system.

Parameters:: df (DataFrame)
Return type:: ndarray

sample_cloud(df, *, weight='log2p1')[source]#

A sample as (Z, w): PCA-coord clonotypes Z and normalized weights w=g(a)/Σg.

The raw material for both the fixed kernel mean (sample_embedding()) and the learned set network (mir.ml.set_encoder) — same basis, so their embeddings are comparable.

Parameters:

df (DataFrame)
weight (str)

save(path)[source]#

Pickle the basis (scaler + PCA + RFF + meta); the model is reconstructed on load.

Raises:: ValueError – If the model uses a non-default matrix / alignment. Neither is recorded in meta, so load() would silently rebuild a different coordinate system while the prototype-hash check still passed.
Return type:: None

classmethod load(path, *, verify=True)[source]#

Load a basis, rebuilding the model and verifying prototype comparability.

Parameters:: verify (bool)
Return type:: RepertoireSpace

mir.repertoire.fit_repertoire_space(model, cohort_df, *, n_rff=2048, n_rff_second=128, n_eigs=None, length_scale=None, n_components=None, space='full', pca_fit_cap=200000, seed=0)[source]#

Fit ONE PCA + RFF basis on the pooled clonotype cloud of a cohort.

Parameters:

model – A fitted single-chain TCREmp.
cohort_df (DataFrame) – Pooled clonotypes across all samples (v_call/j_call/junction_aa); counts are ignored here — only the geometry of the cloud sets the basis.
n_rff (int) – Mean-block RFF dimension D (~1–4k; §T.7 tab:sample).
n_rff_second (int) – Second-moment-block RFF dimension D₂ (kept small — the block stores D₂(D₂+1)/2 upper-triangle entries, or its top-n_eigs eigenvalues).
n_eigs (int | None) – If set, the second-moment block keeps the top-``n_eigs`` eigenvalues of the weighted covariance Σ_σ w_σ ψ₂ψ₂ᵀ (a compact, rotation-invariant spectral signature) instead of its full D₂(D₂+1)/2 upper triangle. None (default) keeps the upper triangle — unchanged behaviour. Must satisfy 0 < n_eigs ≤ n_rff_second.
length_scale (float | None) – Gaussian-kernel bandwidth ℓ. None → mir.density.calibrate_radius() (the one-substitution scale r₁), so the kernel resolves ~one CDR3 mutation.
n_components (int | None) – PCA dimensionality; None → get_preset(species, locus).n_components.
space (str) – "full" (V+J+junction) or "junction" (CDR3 sub-block only).
pca_fit_cap (int | None) – Fit the scaler + PCA on at most this many randomly-sampled pooled rows.
seed (int) – RNG seed (PCA solver + RFF draws).

Returns:

A RepertoireSpace ready for sample_embedding().

Return type:

RepertoireSpace

mir.repertoire.fit_repertoire_spaces(models, cohort_frames, *, min_clonotypes=50, **kwargs)[source]#

Fit one RepertoireSpace per locus — the multi-chain (digital-donor) basis.

Each locus gets an independent basis (its own prototype set + PCA + RFF), so per-chain kernel means are only ever compared within their own chain, honouring the comparability contract. Loci whose pooled cloud is too small to fit a stable PCA are skipped (returned dict omits them).

Parameters:

models (dict) – {locus: TCREmp} — one fitted single-chain model per locus.
cohort_frames (dict) – {locus: pooled clonotype frame} (v_call/j_call/junction_aa), the pooled cloud for that locus across the whole cohort.
min_clonotypes (int) – Skip a locus whose pooled frame has fewer rows (PCA would be degenerate).
**kwargs – Forwarded to fit_repertoire_space() (n_rff, n_components, seed …).

Returns:

{locus: RepertoireSpace} for the loci that could be fit.

Return type:

dict

class mir.repertoire.SampleEmbedding(mean, diversity, second, n_eff, mass=1.0)[source]#

Bases: object

One repertoire’s fixed-width embedding, kept block-wise so MMD uses only the mean block.

mass is the retained probability mass 1 − M₀: 1.0 (the default) means the sample is treated as a complete probability measure, < 1 that missing_mass() estimated part of the repertoire was never drawn. The blocks themselves are unchanged by mass — the deficiency is applied where it is meaningful, in contrast_embedding().

Parameters:

mean (ndarray)
diversity (ndarray | None)
second (ndarray | None)
n_eff (float)
mass (float)

mean: ndarray#

diversity: ndarray | None#

second: ndarray | None#

n_eff: float#

mass: float = 1.0#

property vector: ndarray#

The concatenated feature vector (the multimodal-fusion tensor).

Only the blocks that were computed are concatenated — blocks=("diversity",) is a valid (mean-less) Φ, even though MMD then is not (see mmd_distance()).

mir.repertoire.sample_embedding(space, sample_df, *, weight='log2p1', blocks=('mean', 'diversity', 'second'), coverage=None, missing_mass='none')[source]#

Embed one repertoire into Φ(S) (mean ‖ diversity ‖ second).

Parameters:

space (RepertoireSpace) – A RepertoireSpace from fit_repertoire_space().
sample_df (DataFrame) – One sample’s clonotypes with duplicate_count (counts drive the weights).
weight (str) – Clone-size weight g — "log2p1" (default, g=log2(1+a)) / "duplicate_count" (g=a, linear) / "distinct" (g≡1, presence) / "log1p" (g=ln(1+a)) / "anscombe" (g=√(a+3/8)). Frequencies are w = g(a)/Σg (scale-free).
blocks (tuple[str, ...]) – Which blocks to compute/return.
coverage (float | None) – Common Good–Turing coverage Ĉ* for the diversity block; None uses the sample’s observed Hill numbers.
missing_mass (str) – Estimator for the mass of the never-drawn clonotypes (missing_mass()): "none" (default — a complete probability measure, mass=1, bit-identical to pre-3.8 output), "turing" or "chao". Only .mass changes; the blocks do not.

Returns:

A SampleEmbedding; .vector is the concatenated fixed-width tensor.

Return type:

SampleEmbedding

mir.repertoire.rao_q(emb)[source]#

Rao’s quadratic entropy of the repertoire, read straight off Φ₁ as 1 − ‖Φ₁‖².

Rao’s Q = Σ_{σ,τ} w_σ w_τ d(σ,τ) with the kernel dissimilarity d = 1 − k collapses, for weights summing to 1 and a normalised kernel (k(z,z)=1), to exactly 1 − ‖Φ₁‖² — so the norm of the kernel mean is a diversity statistic and no Gram matrix is needed. Verified against an explicit Gram to machine precision (max relative error ~1e-16).

This is the diversity the Hill block cannot express. Every Hill number is a functional of the clone-size distribution alone, hence invariant to permuting which receptor carries which abundance — it describes the shape of the abundance distribution and is blind to composition. Rao’s Q weights each pair by how different the receptors are, so it is a functional diversity: sequence-aware, and already carried inside Φ. Measured: this one scalar recovers R² 0.74–0.85 of classical diversity, while embedding derivatives reach R² 0.974–0.994 for Shannon and 0.985–0.9999 for richness.

Two caveats. The kernel is a random-Fourier approximation, so k(z,z)=1 holds only to O(D^{-1/2}) and a single-clone sample reads Q of order 1e-2 rather than exactly 0 — widen n_rff if the absolute level matters, or read differences, which are unaffected. And it is valid only for the true, uncentred Φ₁: Centring or PCA-projecting Φ₁ preserves differences (so MMD survives) but not norms, so Q cannot be read off a centred/reduced identity block — check which object you hold before quoting any norm-based quantity. With a deficient mass this is the Q of the observed (renormalised) part, which is the honest reading: it is the diversity you measured.

Parameters:: emb (SampleEmbedding) – A SampleEmbedding with its mean block (uncentred, as returned).
Returns:: Q ∈ [0, 1] — 0 for a single clone, rising as the receptors present grow more dissimilar.
Raises:: ValueError – If emb has no mean block.
Return type:: float

class mir.repertoire.DepthThreshold(kappa, tau2, sigma2, fraction_below, sizes)[source]#

Bases: object

Where sampling noise stops being smaller than between-sample signal (see depth_threshold()).

Parameters:

kappa (float)
tau2 (float)
sigma2 (float)
fraction_below (float)
sizes (ndarray)

kappa: float#

tau2: float#

sigma2: float#

fraction_below: float#

sizes: ndarray#

mir.repertoire.depth_threshold(embs, sizes=None)[source]#

Estimate κ, the sample size below which a repertoire’s Φ is mostly sampling noise.

The damage depth does to a kernel mean is not bias — it is variance ∝ 1/n, over an n spanning four orders of magnitude, and in a neighbour graph that heteroscedasticity is a depth axis. So regress each sample’s squared distance from the cohort centroid on 1/n:

E‖Φ_S − Φ̄‖² ≈ τ² + σ²/n

The intercept τ² is between-sample (biological) spread, the slope σ² is within-sample sampling noise, and ``κ = σ²/τ²`` is the size at which the two are equal — below it a sample’s Φ carries more noise than signal. Measured across four independent views, κ came out at 40–70 clonotypes every time, with 23–69% of samples below it.

This is the estimable replacement for a hand-picked clonotype floor: report κ and fraction_below for the cohort in front of you rather than importing someone else’s cutoff. Note the library still does not apply a floor — in blood a low clonotype count is shallow sequencing, in a tumour it is low infiltration (the phenotype of interest). Pair κ with missing_mass() / contrast_embedding() and weight instead of filtering.

Parameters:

embs – Per-sample SampleEmbedding (mean block present), one cohort, one basis.
sizes – Per-sample size to use as n. None → each embedding’s n_eff. Pass observed clonotype richness to read κ in clonotypes rather than effective clonotypes.

Returns:

A DepthThreshold. kappa is inf when the fit finds no between-sample spread (τ² ≤ 0) and 0.0 when it finds no depth dependence (σ² ≤ 0).

Raises:

ValueError – If fewer than three samples are given, or sizes has the wrong length.

Return type:

DepthThreshold

mir.repertoire.naive_reference(space, *, n=20000, seed=0, sequences=None)[source]#

Kernel mean of n naive V(D)J recombinations — where the unseen clonotypes live.

Once a sample’s mass is deficient (missing_mass()), something has to occupy the unseen block, and the choice decides whether the object is a shrinkage estimator toward a meaningful point. Three priors were measured:

the sample’s own singletons — begs the question, adds no information;
the corpus mean — James–Stein toward the centroid, and it measurably hurt: it piles shallow samples into a dense ball that is itself depth-correlated;
the germline draw — the one that works. Filling the unseen block with it dropped R²(PC1, depth) from 0.259 to 0.001 (blood TRB) and 0.067 to 0.006 (tissue IGH), with kNN label entropy unchanged or better and PC1’s explained variance unchanged (so PC1 is a different direction, not a collapsed one).

The draw comes from vdjtools.model.generate() over the bundled germline model for this space’s locus (~8 s for 20,000 sequences), and is cached per (n, seed) on the space. Pass sequences= to supply the reference repertoire yourself (an uncached path) if you do not want the vdjtools generative model in the loop.

Parameters:

space (RepertoireSpace) – The shared RepertoireSpace — supplies the locus and the coordinate system, so the reference is automatically comparable to every Φ fit through it.
n (int) – Number of naive recombinations to draw.
seed (int) – Generation seed; the result is deterministic in (n, seed).
sequences (DataFrame | None) – Optional pre-built clonotype frame to use instead of generating (v_call / j_call / junction_aa). Not cached.

Returns:

(n_rff,) unweighted kernel mean of the naive draw, in space’s RFF coordinates.

Raises:

ValueError – If the space is non-human (only human germline models are bundled).

Return type:

ndarray

Example

>>> ref = naive_reference(space)
>>> psi = contrast_embedding(emb, ref)

mir.repertoire.contrast_embedding(emb, reference)[source]#

Ψ_S = mass · (Φ₁(S) − reference) — the deficient measure as a signed contrast.

This is the answer to “can the missing mass be a negative probability”. It cannot: Φ’s value is that ‖Φ_P − Φ_Q‖ = MMD(P,Q) and that a convex combination of two Φ’s is the Φ of a real pooled repertoire (what makes mir.twin and trajectory interpolation mean anything), and a measure allowed to go negative on a set is not a probability measure. But a signed difference of two probability measures already gives signed coordinates — negative wherever the sample is depleted relative to unselected recombination — and stays an ordinary RKHS element with ‖Ψ_S‖ = MMD(S, naive) intact.

Note Ψ_S is exactly the reference-centred sub-probability embedding: with Φ_true = mass·Φ₁ + M₀·reference, Φ_true − reference = mass·(Φ₁ − reference).

Combined with a deficient mass, magnitude = confidence × deviation-from-naive. An immune desert has M₀ → 1 and lands at the origin, which is the correct place for “no infiltrate detected”; a vague shallow blood sample lands there too and says so by its norm rather than by being dropped from the cohort. Scale such a block with one global scalar — see mir.explain.ChannelBuilder.add()’s preserve_magnitude.

Parameters:

emb (SampleEmbedding) – A SampleEmbedding with a mean block (embed with missing_mass= for mass < 1; at the "none" default this is just the centred kernel mean).
reference (ndarray) – (n_rff,) unseen-block location, normally naive_reference().

Returns:

(n_rff,) signed contrast vector.

Raises:

ValueError – If emb has no mean block.

Return type:

ndarray

mir.repertoire.band_frames(sample_df, *, bands=('singleton', 'expanded', 'top'), top_fraction=0.01, top_clip=(10, 500), min_clonotypes=5)[source]#

Partition one sample’s clonotypes into compartments — abundance bands or IGH isotypes.

Φ₁ is a clone-size-weighted average, which is the right operation for estimating a population mean and the wrong one for detecting a minority signal. Writing the repertoire as a mixture ρ_S = (1−π)ρ_naive + π ρ_expanded gives Φ₁(S) = (1−π)Φ₁(naive) + π Φ₁(expanded), so a difference confined to the expanded compartment reaches Φ₁ attenuated to πΔ while the noise is supplied by the naive compartment that owns most of the clonotypes. Embedding each compartment separately occupies exactly the blind spot of the Hill block, which is invariant to which receptor carries which abundance.

Bands (on duplicate_count, never on a frequency column):

band	definition	compartment
`singleton`	`count == 1`	naive-dominated background; also where sequencing error concentrates
`expanded`	`count >= 2`	clones that divided at least once
`top`	top `top_fraction` by count, clipped to `top_clip`	the dominant clonal compartment (chronic exposure history)
`igm`	`c_call` in IGHM/IGHD	unswitched (IGH only)
`igg`	`c_call` in IGHG1–4	class-switched (IGH only)
`iga`	`c_call` in IGHA1–2	mucosal (IGH only)

top and expanded overlap deliberately — “what does the dominant compartment look like” and “what does everything that expanded look like” are different questions. The isotype cut partitions on an irreversible molecular event rather than an abundance threshold, so IgM/IgG/IgA are genuinely different cell populations sharing one locus. Rows with a null c_call are excluded from every isotype band, never defaulted to IgM: an uncalled constant region is missing data, and roughly 43% of IGH reads have none at RNA-seq coverage.

Parameters:

sample_df (DataFrame) – One sample’s clonotypes with duplicate_count (and c_call for isotypes).
bands (tuple[str, ...]) – Which bands to build. Any of singleton/expanded/top/igm/igg/iga.
top_fraction (float) – Fraction of clonotypes in the top band before clipping.
top_clip (tuple[int, int]) – (min, max) clonotype count for the top band.
min_clonotypes (int) – A band with fewer clonotypes is returned as None — recorded absent, never embedded. A 3-clone “repertoire” is not an estimate; None is the same hole convention ChannelBuilder and mir.cohort.align_loci() already understand.

Returns:

{band: frame or None} in the requested order.

Raises:

ValueError – On an unknown band name, or an isotype band without a c_call column.

Return type:

dict

mir.repertoire.band_embeddings(space, sample_df, *, bands=('singleton', 'expanded', 'top'), min_clonotypes=5, **kwargs)[source]#

Embed each compartment of one sample through the same basis (see band_frames()).

A band is never a reason to refit: refitting would put each compartment in its own geometry and make band-to-band and cohort-to-cohort distances meaningless. One space, every band.

Parameters:

space (RepertoireSpace) – The shared RepertoireSpace.
sample_df (DataFrame) – One sample’s clonotypes.
bands (tuple[str, ...]) – Which bands to embed (band_frames()).
min_clonotypes (int) – Bands smaller than this are None (absent, not embedded).
**kwargs – Forwarded to sample_embedding() (weight, blocks, missing_mass…).

Returns:

{band: SampleEmbedding or None}.

Return type:

dict

Note

To compare two compartments’ directions, centre the cohort first (X - X.mean(0), or centroid_atypicality() for the per-group version) and cosine the residuals. Do not read a raw cosine between two TCREmp-derived vectors, and do not think subtracting each vector’s own scalar mean fixes it: on mean distance profiles, measured raw cos 0.9993, self-centred 0.9985, cohort-centred −0.14. Only cohort centring removes the shared profile. Measured compartment overlaps of order 0.01 (RFF kernel means, whose null is already ≈0) are real signal above that baseline.

mir.repertoire.mixture_weights(whole, parts)[source]#

Non-negative mixture weights π_c of compartments inside a whole-repertoire Φ₁.

Φ is linear in the clone-weight measure, so Φ(S) = Σ_c π_c Φ(c) exactly for any partition, which makes a non-negative least squares of the whole on its compartments well-posed rather than a heuristic fit. (Exact in exact arithmetic: the realised agreement is limited by the float32 clonotype embedding, whose per-row values depend slightly on the thread batching, so expect a relative residual around 1e-5 rather than machine epsilon.) The weights are the answer to “how much of what I am measuring does this compartment actually own” — the dilution factor π of band_frames(), measured instead of assumed.

Two things it is for. Interpretation: measured on IGH, the class-switched IgG compartment — affinity-matured, T-dependent, the fraction most likely to carry antigen-specific signal — carries a median π of only 0.070 of Φ₁(IGH), with unswitched IgM at 0.230, IgA at 0.176 and 0.520 unexplained (the uncalled-isotype share, agreeing in ballpark with the ~0.43 counted from reads — different denominators, embedding weight vs reads). Power: π is a prior estimate of how far a subset can move Φ, so a subset carrying π ≈ 0.001 will not be detectable by any aggregate distance on Φ, and the per-clonotype route (class_witness()) is the sensitive one. Check π before spending compute on the aggregate.

Parameters:

whole – The whole repertoire’s Φ₁ — a SampleEmbedding or a (D,) array.
parts (dict) – {name: SampleEmbedding | array | None}; None entries are skipped (an absent band is not a zero-weight band).

Returns:

{"weights": {name: π}, "residual": 1 − Σπ, "r2": reconstruction R²}. residual is the share of the whole no compartment accounts for (an unmeasured or excluded compartment), and can be negative if the parts overlap (top ⊂ expanded).

Raises:

ValueError – If no usable parts were given, or a part’s width differs from whole’s.

Return type:

dict

class mir.repertoire.RarefyResult(embedding, v_rep, depth, n_replicates)[source]#

Bases: object

A depth-standardised embedding plus the replicate dispersion (see rarefy_embedding()).

Parameters:

embedding (SampleEmbedding)
v_rep (float)
depth (int)
n_replicates (int)

embedding: SampleEmbedding#

v_rep: float#

depth: int#

n_replicates: int#

rao_gap()[source]#

v_rep again, under the name of the identity it satisfies.

Rao(Φ̄) = mean_r Rao(Φ_r) + v_rep holds exactly: Φ̄ embeds the mixture over subsamples, which is genuinely more diverse than any single subsample, so the excess diversity of the average is the replicate variance. Averaging commutes with Φ but not with nonlinear functionals of it — mean(Rao) ≠ Rao(mean) — and this is the exact gap.

Return type:: float

mir.repertoire.rarefy_embedding(space, sample_df, depth, *, n_replicates=10, weight='log2p1', seed=0)[source]#

Embed a repertoire rarefied to a common read depth, averaged over replicate draws.

The only depth correction that preserves the kernel-mean semantics exactly. Φ is linear in the clone-weight measure, so the mean of Φ over independent multinomial subsamples is itself a kernel mean embedding — of the mixture distribution over subsamples — and MMD, Rao’s Q and mixture linearity all survive (verified to ~1e-15). An orthogonal projection breaks Rao; a per-coordinate location-scale rescale breaks MMD and Rao both.

Not a default, and not a substitute for carrying depth explicitly: rarefying a whole cohort to the depth of its shallowest useful samples throws away the deep samples’ entire advantage, which is why depth_threshold() + missing_mass() are the first-line tools. Reach for this when two groups must be compared at matched depth and you need the comparison to stay an MMD.

Only the mean block is averaged — see RarefyResult.rao_gap().

Parameters:

space (RepertoireSpace) – The shared RepertoireSpace.
sample_df (DataFrame) – One sample’s clonotypes with duplicate_count.
depth (int) – Reads to draw per replicate. Must not exceed the sample’s total.
n_replicates (int) – Independent multinomial draws to average (v_rep needs ≥ 2).
weight (str) – Clone-size weight applied to the rarefied counts.
seed (int) – RNG seed.

Returns:

A RarefyResult.

Raises:

ValueError – If depth exceeds the sample’s read total, depth < 1, or n_replicates < 1.

Return type:

RarefyResult

class mir.repertoire.RepertoireDescriptor(log_mass, log_neff, simpson, mean)[source]#

Bases: object

Smooth, mass-preserving repertoire descriptor — every summary metric is a derivable coordinate.

sample_embedding()’s Φ frequency-normalises, which discards the total mass (= infiltration). The descriptor instead keeps the mass as a coordinate alongside diversity and clonality, so the whole object is:

decodable — metrics() reads infiltration / diversity / clonality off analytically;
smooth — log mass, log n_eff and Simpson λ are continuous (no integer richness ⁰D), the “smoother form” of the Hill block;
simulatable — vector is a fixed-width continuous vector you can perturb and generate: fit a density over a cohort, move along the infiltration coordinate → in-silico “hotter / colder” (decode_metrics() turns any perturbed vector back into named metrics).

mean is the normalised kernel mean μ/G (clonotype identity / composition); the scalar coordinates are the count-distribution summaries. Distances/generation live in the concatenated vector.

Parameters:

log_mass (float)
log_neff (float)
simpson (float)
mean (ndarray)

log_mass: float#

log_neff: float#

simpson: float#

mean: ndarray#

property scalar: ndarray#: The derivable-metric coordinates [infiltration, log n_eff, clonality].

property vector: ndarray#: The full descriptor [scalar ‖ mean] — the continuous object to perturb / generate.

metrics()[source]#

Named metrics read off the coordinates (all analytic, no recomputation from counts).

Return type:: dict

mir.repertoire.sample_descriptor(space, sample_df, *, weight='log2p1')[source]#

Mass-preserving smooth descriptor of one repertoire (see RepertoireDescriptor).

The scale (log_mass = infiltration) is retained rather than normalised away, so infiltration, diversity and clonality are all smooth coordinates of the same object — the representation the in-silico-evolution / embedding-simulation workflow perturbs.

Parameters:

space (RepertoireSpace)
sample_df (DataFrame)
weight (str)

Return type:

RepertoireDescriptor

mir.repertoire.decode_metrics(vector)[source]#

Read named metrics off a (possibly perturbed or generated) descriptor vector — the inverse used for in-silico evolution: perturb RepertoireDescriptor.vector, decode the new metrics.

Parameters:: vector (ndarray)
Return type:: dict

mir.repertoire.mmd_distance(a, b, *, unbiased=False)[source]#

MMD between two repertoires, ‖Φ₁(a) − Φ₁(b)‖ (Eq. eq:kme).

The default (unbiased=False) is the biased V-statistic ‖μ̂_a − μ̂_b‖ — simple, but its self-terms carry a positive bias ≈ 1/n_eff (from the k(z,z) diagonal), so a low-diversity (small n_eff) sample gets its distances inflated by construction. When diversity is itself the variable of interest (e.g. divergence-vs-age), that bias masquerades as signal with the wrong sign (distance tracks low diversity, not high). unbiased=True removes the diagonal analytically using the stored n_eff (Gretton et al. 2012, unbiased MMD²) — the estimator to trust when comparing samples of unequal depth/diversity.

Parameters:

a (SampleEmbedding)
b (SampleEmbedding)
unbiased (bool)

Return type:

float

mir.repertoire.mmd_matrix(embs, *, unbiased=False)[source]#

Symmetric sample×sample MMD matrix (feeds a regressor / cluster_samples).

unbiased=True uses the diagonal-removed MMD² (see mmd_distance()) — necessary whenever samples differ in depth/n_eff, else the 1/n_eff self-bias confounds the comparison.

Parameters:

embs (list[SampleEmbedding])
unbiased (bool)

Return type:

ndarray

mir.repertoire.class_witness(space, pos, neg, candidates, *, weight='log2p1', top=30, witness=None)[source]#

Rank clonotypes by how much they drive the pos vs neg group difference (Prop. prop:witness).

The MMD witness is the mean-embedding difference w = μ_pos − μ_neg in RFF feature space; a clonotype σ scores s(σ) = ⟨w, ψ(φ(σ))⟩, so the top-scoring candidates are the discriminative public clones / motifs separating the two groups — the supervised way to find motifs that the unsupervised bulk kernel mean cannot surface (it is swamped by the naive background).

Parameters:

space (RepertoireSpace) – The shared RepertoireSpace.
pos (list[DataFrame]) – Per-group lists of sample clonotype frames (with duplicate_count).
neg (list[DataFrame]) – Per-group lists of sample clonotype frames (with duplicate_count).
candidates (DataFrame) – Clonotype frame to score (e.g. all clonotypes seen in pos).
weight (str) – Clone-size weight for the per-sample kernel means.
top (int) – Number of top motifs to return.
witness (ndarray | None) – Optional precomputed witness direction μ_pos − μ_neg ((D,)). When given, pos/neg are not re-embedded — reuse it to score several candidate sets cheaply.

Returns:

candidates with a witness_score column, sorted descending, truncated to top.

Return type:

DataFrame

mir.repertoire.hla_stratified_mmd(embs, hla)[source]#

MMD restricted to HLA-matched sample pairs; nan where two samples share no allele.

Response ridges are presented only by a specific HLA type, so cross-sample overlap is informative only within an HLA-matched stratum (Prop. prop:hla). Pairs that share no HLA allele are masked to nan rather than compared.

Parameters:

embs (list[SampleEmbedding])
hla (list[set])

Return type:

ndarray

mir.repertoire.centroid_atypicality(X, groups)[source]#

Per-sample cosine distance of its identity vector to its group’s centroid — a Φ-geometry op.

A sample far from its group’s mean identity has an atypical clonotype composition (selection / divergence). Grouping is caller-supplied (tumour type, cohort, …); the geometry — centroid then 1 − cos — is the library concern. Feeds an atypicality channel of a digital-donor embedding.

Parameters:

X (ndarray) – (n, d) identity block (a per-sample kernel-mean or its PCA reduction).
groups (ndarray) – (n,) group label per row; centroids are computed within each group.

Returns:

(n,) atypicality in [0, 2] (0 = on the group centroid direction).

Return type:

ndarray

mir.repertoire.correct_batch(X, batch, *, covariates=None, n_clusters=8, theta=1.0, sigma=0.1, max_iter=10, ridge=1.0, seed=0)[source]#

Harmony-like cluster-aware batch correction on a stacked sample×feature Φ matrix.

Removes a batch offset per soft cluster rather than globally, so a batch confounded with a biological cluster is corrected without erasing that biology — the failure mode of plain per-group mean subtraction (mir.cohort.residualize()). Follows Harmony (Korsunsky et al. 2019, Nat. Methods): soft-cluster the samples with a batch-diversity penalty theta (clusters pushed toward batch-balanced membership), then subtract each cluster’s membership-weighted batch offset (covariates retained), and iterate. Reduces exactly to residualize at n_clusters=1 or theta=0.

Parameters:

X (ndarray) – Stacked (n_samples, n_features) matrix (e.g. mir.explain.stack_embeddings()).
batch (ndarray) – Length-n_samples batch label per row.
covariates (ndarray | None) – Optional (n_samples, k) biological covariates to retain (never removed).
n_clusters (int) – Number of soft clusters (Harmony’s K); <=1 disables clustering.
theta (float) – Batch-diversity penalty strength; 0 disables clustering (→ residualize).
sigma (float) – Soft-assignment temperature.
max_iter (int) – Correction iterations (converges in a few — the batch fit vanishes once removed).
ridge (float) – L2 ridge on the per-cluster batch regression (stabilises small clusters).
seed (int) – KMeans seed for the soft-cluster initialisation.

Returns:

The corrected (n_samples, n_features) matrix — a new coordinate system (as with residualize, never compare a corrected X to an uncorrected one).

Return type:

ndarray

Explainable readouts (`mir.explain`)#

Which named channel of Φ carries the signal, and which clonotypes drive it (Theory T7). ChannelBuilder.add(..., preserve_magnitude=True) scales a magnitude-carrying (sub-probability) block by one global scalar instead of per-column z-scoring.

Explainable readouts over a repertoire feature matrix: which channel carries the signal (§T.7).

vdjtools-era repertoire analysis answers “which summary statistic separates my groups?” by running a fixed menu of named statistics side by side — diversity, clonality, CDR3 length, V-usage, overlap — one test per statistic, and reading off which one moved. The menu is the explanation: every number has a name, so a result is a sentence (“the groups differ in clonality”).

The embedding replaces the menu with one vector Φ(S) (mir.repertoire), which scores better and explains worse: Φ.vector is an anonymous concatenation, so “the classifier found something” has no noun. This module restores the noun without giving up the vector. A channel is a named group of Φ columns (the mean/identity block, the Hill block, the second-moment block, and any extra blocks a study bolts on — per-chain blocks of the same kind merge under one name). A ChannelSpec is the name→column-index map that Φ.vector alone does not carry; ChannelBuilder assembles one alongside the matrix, so the layout is declared where the blocks are stacked, not rediscovered later.

Given the map, “which part of Φ carries this signal” is an ablation against a user-supplied scorer f: X_block → float (higher is better):

delta_in(g)  = f(X[:, cols(g)]) − base        # does channel g carry signal on its own?
delta_out(g) = f(X) − f(X[:, ¬cols(g)])       # is channel g's signal anywhere else?

The scorer is a closure over the study’s labels, so the same call works for a Cox C-index, a cross-validated AUC, or anything else — this module never sees y and ships no scorers of its own (model choices belong to the analysis, not the library). delta_in is the default and is marginal: it is inflated by correlation between channels, so two redundant channels both look important. delta_out is conditional and deflated by the same correlation. Reported together they separate the two claims — high in / high out = irreplaceable; high in / ~zero out = redundant, its signal is duplicated elsewhere. Significance, when asked for, is a row permutation of the block: the scorer holds y in row order, so shuffling the block’s rows breaks the association exactly as shuffling y would, and it is the only scorer-agnostic null available.

Leave-one-in is the default for a concrete reason: the scorers in practice reduce their input (an in-fold PCA to n_pc components), so dropping one narrow channel out of a wide block leaves the reduction to re-mix the rest and reconstruct nearly the same components — leave-one-out is structurally near-blind to exactly the 1-column channels that often win. Ask for mode="both" when the question is whether the winner is necessary rather than sufficient.

The last hop, channel → clonotypes, exists only for channels with a clonotype pre-image: a kernel-mean block, where mir.repertoire.class_witness() reads the MMD witness s(σ)=⟨μ_pos−μ_neg, ψ(φ(σ))⟩. A Hill number or a read-count fraction has no such pre-image — asking which clones drive the diversity channel is a category error, and channel_drivers() raises rather than answer it. Attributability is therefore declared per channel at build time, never guessed from its name.

Torch-free (numpy / polars); scorer-agnostic (no sklearn / lifelines import here).

Typical usage:

from mir.explain import ChannelBuilder, channel_drivers, channel_report, stack_embeddings

embs = [sample_embedding(space, s) for s in samples]
X, spec = stack_embeddings(embs)                     # channels: mean ‖ diversity ‖ second
rep = channel_report(X, spec, lambda B: cv_auc(B, y)[0], base=0.5)
print(rep.frame())                                   # one row per channel: score, delta, rank
rep.best                                             # -> "second"   (the HLA imprint)

# multi-source assembly (per-chain blocks merge by name) + the channel -> clonotypes hop
b = ChannelBuilder()
for c in chains:
    b.add("identity", ident[c], attributable=True).add("diversity", hill[c])
X, spec = b.add("coverage", log_reads).build()       # median-impute + z-score
rep = channel_report(X, spec, lambda B: cv_cindex(dur, evt, base=clin, block=B, n_pc=8),
                     base=cv_cindex(dur, evt, base=clin), mode="both")
channel_drivers(rep, space=space, pos=pos, neg=neg, candidates=cands)   # if rep.best is a KME

class mir.explain.ChannelSpec(columns_by_name, attributable=frozenset({}))[source]#

Bases: object

Name → column-index map for a feature matrix — the label Φ.vector does not carry.

spec[name] is an alias for columns(), so a spec drops into code written against the plain dict[str, list[int]] it replaces (X[:, spec["coverage"]], or spec["a"] + spec["b"] for a union of columns).

Variables:

columns_by_name (dict[str, list[int]]) – Channel name → its column indices in X, in assembly order.
attributable (frozenset[str]) – Channels with a clonotype pre-image (kernel-mean blocks) — the only ones channel_drivers() will attribute. Declared at build time, never inferred.

Parameters:

columns_by_name (dict[str, list[int]])
attributable (frozenset[str])

columns_by_name: dict[str, list[int]]#

attributable: frozenset[str] = frozenset({})#

property names: list[str]#: Channel names in assembly order.

property width: int#: Total column count — must equal X.shape[1].

columns(name)[source]#

Column indices of name.

Parameters:: name (str) – Channel name.
Returns:: The channel’s column indices in X.
Raises:: ValueError – If name is not a channel.
Return type:: list[int]

class mir.explain.ChannelBuilder(_blocks=<factory>, _cols=<factory>, _attr=<factory>, _mag=<factory>, _col=0)[source]#

Bases: object

Accumulate named blocks into one matrix + its ChannelSpec.

Blocks added under the same name merge into one channel (per-chain identity blocks become a single identity channel), which is the granularity ablation asks about. For per-chain resolution, add under distinct names (f"identity:{chain}").

Parameters:

_blocks (list[ndarray])
_cols (dict[str, list[int]])
_attr (set[str])
_mag (set[str])
_col (int)

add(name, mat, *, attributable=False, preserve_magnitude=False)[source]#

Append a block. Returns self for chaining.

Parameters:

name (str) – Channel name; repeat a name to merge blocks into one channel.
mat (ndarray) – (n_samples,) or (n_samples, k). 1-D is treated as a single column.
attributable (bool) – Mark the channel clonotype-attributable (a kernel-mean block). Sticky: one attributable block makes the whole merged channel attributable.
preserve_magnitude (bool) – Scale this channel by one global scalar (its pooled RMS over observed entries, no centring) instead of per-column z-scoring, and fill holes with 0. Required for any block whose magnitude carries information — a sub-probability / contrast block from mir.repertoire.contrast_embedding(), where a near-zero row means “no confident deviation from naive”. Per-column standardisation forces every coordinate to unit variance across samples, so a matrix where half the rows are near-zero comes out looking exactly like one where none are: it deletes the deficiency it was meant to preserve. (This already invalidated one experiment, whose deficient arm scored 0.6481 against a 0.6179 baseline purely from the rescaling.) Sticky, as with attributable.

Returns:

self.

Raises:

ValueError – On a row-count mismatch with the blocks already added.

Return type:

ChannelBuilder

build(*, standardize=True, impute=True)[source]#

Assemble (X, spec).

Parameters:

impute (bool) – Replace non-finite entries with the column median (0.0 if a column is all non-finite) — chains a sample lacks leave holes.
standardize (bool) – Z-score each column. Default True because the builder exists for heterogeneous assembly, where a log-read-count column and a read-fraction column differ by orders of magnitude and any distance/PCA/penalised fit is otherwise dominated by the loudest unit. Pass False to keep raw values.

Return type:

tuple[ndarray, ChannelSpec]

Both statistics come from the observed entries of each column, before imputation. The order matters: imputing first and then z-scoring computes sd over a column whose holes have all been set to one constant, which deflates it in proportion to how sparse the column is — so a locus present in 30% of samples ends up with its real values scaled up ~1.8x against a fully-covered one, giving the least-observed chain the most weight in every downstream distance, PCA and penalised fit. Standardising on observed entries makes an imputed hole land at exactly 0 (the column mean, i.e. no information) in a column whose scale the hole did not set.

Channels added with preserve_magnitude=True are exempt from both: they get one global scalar and no centring, so their row norms survive, and their holes are filled with 0.

Returns:

(X, spec).

Raises:

ValueError – If no blocks were added.

Parameters:

standardize (bool)
impute (bool)

Return type:

tuple[ndarray, ChannelSpec]

mir.explain.stack_embeddings(embs)[source]#

Stack Φ(S) embeddings into (X, spec) with channels named after their blocks.

X is row-wise identical to np.stack([e.vector for e in embs]) — this only attaches the names, it does not transform. No z-scoring: the blocks are one homogeneous Φ and rescaling them would silently change every MMD/PCA downstream.

Parameters:: embs (Sequence[SampleEmbedding]) – Per-sample embeddings from mir.repertoire.sample_embedding(), all fitted from one RepertoireSpace.
Returns:: (X, spec) with X shape (len(embs), Φ_dim) and channels mean / diversity / second (whichever are present). mean is marked attributable.
Raises:: ValueError – If embs is empty or the embeddings disagree on present blocks / widths.
Warns:: UserWarning – If any embedding carries a deficient mass (mass < 1). Φ.vector does not encode it — the blocks are renormalised measures either way — so a deficient cohort stacked this way silently discards the deficiency. Build the mean block with mir.repertoire.contrast_embedding() and add it under ChannelBuilder.add(..., preserve_magnitude=True) instead.
Return type:: tuple[ndarray, ChannelSpec]

class mir.explain.ChannelReport(channels, score, delta, rank, base, full, n_samples, spec, mode, delta_out=None, pvalue=None)[source]#

Bases: object

Per-channel ablation readout — the explainable answer to “which part of Φ carries this”.

Variables:

channels (list[str]) – Channel names in report order.
score (numpy.ndarray) – scorer(X[:, cols(g)]) — the leave-one-in score of each channel alone.
delta (numpy.ndarray) – score − base. nan if no base was given.
rank (numpy.ndarray) – 1 = best, by the primary statistic of mode.
base (float) – Reference score with no channel features (Cox clinical-only, or 0.5 for chance); nan if not supplied.
full (float) – scorer over all reported channels — the “C+Φ” headline.
n_samples (int) – Rows scored.
spec (mir.explain.ChannelSpec) – The ChannelSpec scored (provenance; also carries attributable).
mode (str) – "in" | "out" | "both".
delta_out (numpy.ndarray | None) – full − scorer(X[:, ¬cols(g)]) — the drop caused by removing g, so larger = more uniquely necessary (same sign convention as delta). None unless requested.
pvalue (numpy.ndarray | None) – One-sided row-permutation p per channel, add-one smoothed (≥ 1/(1+n_permutations), never 0). None unless requested.

Parameters:

channels (list[str])
score (ndarray)
delta (ndarray)
rank (ndarray)
base (float)
full (float)
n_samples (int)
spec (ChannelSpec)
mode (str)
delta_out (ndarray | None)
pvalue (ndarray | None)

channels: list[str]#

score: ndarray#

delta: ndarray#

rank: ndarray#

base: float#

full: float#

n_samples: int#

spec: ChannelSpec#

mode: str#

delta_out: ndarray | None = None#

pvalue: ndarray | None = None#

property best: str#: The top-ranked channel name.

property gain: float#: full − base — the whole embedding’s gain over the reference (ΔC / ΔAUC).

frame()[source]#

Long-form readout, one row per channel, sorted by rank.

Returns:: Columns channel, n_columns, score, delta, rank, attributable, plus delta_out / pvalue when present. Pivot for the wide one-value-per-cell table; this returns the tidy form.
Return type:: DataFrame

mir.explain.channel_report(X, spec, scorer, *, base=None, mode='in', channels=None, n_permutations=0, seed=0)[source]#

Score every channel of X under scorer and rank them (§T.7).

Parameters:

X (ndarray) – (n_samples, spec.width) feature matrix.
spec (ChannelSpec) – Channel map for X’s columns.
scorer (Callable[[ndarray], float]) – X_block → float, higher is better, closing over the study’s labels (lambda B: cv_auc(B, y)[0], lambda B: cv_cindex(dur, evt, base=clin, block=B, n_pc=8)). Any block-width-dependent choice — in-fold PCA, penaliser — belongs inside the closure; this function only slices columns.
base (float | None) – Reference score with no channel features. Cox: cv_cindex(dur, evt, base=clin). Classification: 0.5. None → delta is nan and ranking falls back to score (an identical ordering — delta is a constant shift of score).
mode (str) – "in" (default; each channel alone vs base — the marginal question), "out" (each channel removed vs full — the uniqueness question), or "both".
channels (Sequence[str] | None) – Restrict / order the report. None → spec order. For "out"/"both", full and the complements are taken over the union of the reported channels.
n_permutations (int) – If > 0, a one-sided permutation p per channel. Costs n_permutations × len(channels) extra scorer calls; off by default.
seed (int) – Permutation RNG seed.

Returns:

A ChannelReport; .best is the winning channel, .frame() the tidy table.

Raises:

ValueError – If X.shape[1] != spec.width, channels names an unknown channel, the channel list is empty, or mode is not one of "in"/"out"/"both".

Return type:

ChannelReport

mir.explain.channel_drivers(report, *, space, pos, neg, candidates, channel=None, weight='log2p1', top=30)[source]#

Clonotypes driving a channel — the channel → clones hop (Prop. prop:witness).

Delegates to mir.repertoire.class_witness() (the MMD witness s(σ)=⟨μ_pos−μ_neg, ψ(φ(σ))⟩) and guards it: only a channel declared attributable — a kernel-mean block — has a clonotype pre-image. A Hill number, a read-count fraction or a centroid distance is a summary; its “drivers” are the count distribution, not any clone, and asking for them is a category error this function refuses rather than answers.

NB the witness is recomputed from pos/neg through space’s raw kernel mean, so it attributes the channel’s source, not the exact (typically reduced) columns that were scored. That is the intended reading — the reduction is a fitting convenience, the mean is the object — but a driver list explains μ_pos−μ_neg, not the specific components.

Parameters:

report (ChannelReport) – The readout naming the channel.
space (RepertoireSpace) – The shared RepertoireSpace the channel was built from.
pos (list[DataFrame]) – Per-sample clonotype frames for the positive group (with duplicate_count).
neg (list[DataFrame]) – Per-sample clonotype frames for the negative group.
candidates (DataFrame) – Clonotype frame to score (e.g. all clonotypes seen in pos).
channel (str | None) – Which channel to attribute. None → report.best.
weight (str) – Clone-size weight for the per-sample kernel means.
top (int) – Number of driving clonotypes to return.

Returns:

candidates with witness_score (descending, truncated to top) plus a channel column, so a stacked driver frame stays self-describing.

Raises:

ValueError – If channel is unknown, or is not clonotype-attributable.

Return type:

DataFrame

Cohort / digital donor (`mir.cohort`)#

Fuse per-chain repertoire embeddings into one hash-verified, serialisable donor matrix; batch residualisation (with optional James–Stein shrinkage), depth and missingness diagnostics, sample clustering, and incidence biomarkers (Theory T7).

The digital donor: one fixed, explainable descriptor per subject, fused across chains (§T.7).

Where mir.repertoire embeds one repertoire (one chain of one sample) and mir.explain names the pieces of a feature matrix, this module assembles the cohort object the clinical pipelines actually consume — a digital donor whose vector concatenates, per locus, the identity (kernel mean), diversity (Hill) and coverage (log receptor load) channels into one row, ready for survival / classification via mir.explain.channel_report().

The library owns the geometry, fusion and serialization; the study owns which extra channels (isotype, composition, HLA, …) and which scorer. So fit_donor_embeddings() builds the per-chain geometry and fuses it through mir.explain.ChannelBuilder, and takes an extra_channels hook for the analysis to inject its own tissue/clinical blocks into the same matrix — the channels then flow through channel_report indistinguishably.

Comparability bites twice here (two bases, not one): each locus carries its own prototype hash and a cross-sample identity-PCA rotation. DonorCohort stores both; DonorCohort.load() verifies every locus’s prototype hash, and DonorCohort.transform() is the only way to project held-out donors into the fitted basis. A batch-residualize``d ``X is likewise incomparable to the raw one — record the correction in provenance.

Torch-free (numpy / sklearn / polars); vdjtools is lazy (only incidence_biomarkers()).

class mir.cohort.LocusAlignment(ids, blocks, mask, loci)[source]#

Bases: object

Per-locus matrices aligned onto one sample axis, with holes where a locus is missing.

Variables:

ids (list) – The sample ids, sorted — the row order of every block and of mask.
blocks (dict) – {locus: (n_samples, k) array}; rows for samples lacking that locus are nan.
mask (numpy.ndarray) – (n_samples, n_loci) boolean presence, column order loci.
loci (list) – Locus names, in the order given.

Parameters:

ids (list)
blocks (dict)
mask (ndarray)
loci (list)

ids: list#

blocks: dict#

mask: ndarray#

loci: list#

property n_loci: ndarray#: (n_samples,) count of loci present per sample.

property pattern: ndarray#: (n_samples,) presence pattern as a string, e.g. "1101100" — the missingness stratum.

build(*, standardize=True, attributable=True)[source]#

Fuse to (X, spec) through ChannelBuilder, one channel per locus.

The holes are imputed and standardized by the builder, which takes both statistics from observed entries — so a locus present in a third of samples is not upweighted by its own sparsity, and a hole lands at exactly 0.

Parameters:

standardize (bool) – Z-score each column on its observed entries.
attributable (bool) – Mark the per-locus blocks clonotype-attributable (they are kernel means).

Returns:

(X, spec).

missingness_report(labels)[source]#

Whether a grouping tracks which loci a sample has rather than what is in them.

See missingness_report().

Return type:: dict

mir.cohort.align_loci(blocks, *, how='union', require=None, min_loci=1)[source]#

Align per-locus embedding matrices, keyed by sample id, onto one sample axis.

The step between “one matrix per locus, each over its own samples” and “one matrix over one sample set”. Loci are sequenced to different depths, so their sample sets differ — and an inner join across seven loci is bound by the thinnest two. Measured in a 62,293-sample tissue cohort: intersecting all seven left 19,346 samples (31%), while the union keeps every one. That is complete-case deletion wearing a join’s clothes, so how="union" is the default.

A missing locus becomes nan rows, which is the hole convention ChannelBuilder already understands: it imputes them to the observed centre and standardizes on observed entries, so an absent locus contributes no information and does not set any channel’s scale.

Do not zero-fill before calling this. A literal 0 is a value — after a naive global z-score it becomes -mean/std, a large shared constant that every sample missing that locus carries, and the cohort then clusters by coverage. (Measured: cluster-vs-locus-count AMI 0.741 and read-depth eta-squared 0.42 cohort-wide against 0.01 inside the fully-observed subset.) nan says “not observed”; 0 says “observed to be zero”.

Parameters:

blocks (dict) – {locus: (ids, matrix)} — ids a sequence of sample ids, matrix an (len(ids), k) array. Per-locus k may differ.
how (str) – "union" keeps every sample seen in any locus (holes where absent); "inner" keeps only samples present in all of them. "inner" exists to be compared against, not as a default.
require (list | None) – Loci a sample must have to be kept at all — for when one chain genuinely carries the question and the others are context.
min_loci (int) – Drop samples carrying fewer than this many of the loci.

Returns:

A LocusAlignment.

Raises:

ValueError – If blocks is empty, a matrix’s row count disagrees with its ids, how is not recognised, require names an absent locus, or nothing survives the filters.

Return type:

LocusAlignment

Example

>>> al = align_loci({"TRB": (trb_ids, TRB), "IGH": (igh_ids, IGH)})
>>> X, spec = al.build()
>>> al.missingness_report(clusters)["ami_vs_pattern"] < 0.1
True

mir.cohort.missingness_report(labels, mask)[source]#

How much of a grouping is explained by which blocks each sample has.

A cohort assembled from partly-observed blocks can cluster by coverage rather than by content, and nothing in X says which happened. This is the check to run before believing a clustering, an enrichment or a trajectory built on holed data.

Parameters:

labels – (n_samples,) cluster or group assignment to test.
mask – (n_samples, n_blocks) boolean presence, e.g. LocusAlignment.mask.

Returns:

ami_vs_pattern — adjusted mutual information between labels and the presence pattern; eta2_n_present — correlation ratio of the per-sample present-block count; n_patterns; mean_present. Near zero on the first two is what you want. High means the grouping is a missingness stratification: re-read it inside a fully-observed subset before reporting it as biology.

Return type:

dict

Example

>>> missingness_report(clusters, al.mask)["ami_vs_pattern"]
0.041

mir.cohort.depth_report(X, stats, *, n_pc=5, subset=None)[source]#

How much of an embedding’s leading geometry is explained by sampling depth.

The companion to missingness_report(): that one asks whether a grouping tracks which blocks a sample has, this one whether the embedding’s dominant axes track how deeply it was sequenced. Each of the top n_pc principal components is regressed on the sampling fingerprint (mir.repertoire.cohort_statistics()) and reported as R².

Report the leading block, not only PC1 — a correction can move depth off PC1 and leave it in PC2. Measured with the deficient measure (mir.repertoire.contrast_embedding()), R²(PC1, depth) fell 0.259 → 0.001 and best-of-PC1–5 0.253 → 0.047, while PC1’s explained variance was unchanged, which is what distinguishes “a different direction” from “a collapsed one” — so read explained_variance alongside the R².

Parameters:

X – (n_samples, d) embedding matrix (uncentred is fine; PCA centres internally).
stats (dict) – {name: (n_samples,) values} sampling fingerprint. Non-finite entries are mean-filled so one missing library size cannot drop a sample.
n_pc (int) – Number of leading components to test.
subset – Optional boolean mask of a trusted subset (e.g. samples above mir.repertoire.depth_threshold()’s κ). When given, the same statistics are recomputed inside it and returned under *_subset: high overall and low within the subset means the structure was between coverage strata.

Returns:

r2_pc1, r2_best, pc_best (1-indexed), r2_by_pc, explained_variance (per-PC ratio), n_pc, stats (names, in fit order), n_samples — plus r2_pc1_subset / r2_best_subset / n_samples_subset when subset is given.

Raises:

ValueError – If stats is empty or a statistic’s length disagrees with X.

Return type:

dict

class mir.cohort.DonorCohort(X, spec, spaces, identity_pca, rows, meta=<factory>)[source]#

Bases: object

A fitted cohort of digital donors: the fused matrix plus everything needed to reproduce it.

Variables:

X (numpy.ndarray) – (n_donors, width) fused, imputed, z-scored feature matrix.
spec (mir.explain.ChannelSpec) – ChannelSpec — the channel → column map for X.
spaces (dict) – {locus: RepertoireSpace} — the per-locus clonotype bases (hash-verified on load).
identity_pca (dict) – {locus: fitted StandardScaler→PCA} — the cross-sample identity reducer, stored (not thrown away) so held-out donors project comparably.
rows (list | None) – Per-donor metadata, row-aligned to X (None if not supplied).
meta (dict) – Provenance — per-locus prototype hashes, id_pca, block order, standardization stats.

Parameters:

X (ndarray)
spec (ChannelSpec)
spaces (dict)
identity_pca (dict)
rows (list | None)
meta (dict)

X: ndarray#

spec: ChannelSpec#

spaces: dict#

identity_pca: dict#

rows: list | None#

meta: dict#

missingness_report(labels)[source]#

Whether a grouping tracks which chains each donor has, rather than what is in them.

See missingness_report(); the presence matrix is the one recorded at fit, one column per assembled block.

Parameters:: labels – (n_donors,) cluster or group assignment to test.
Returns:: {ami_vs_pattern, eta2_n_present, n_patterns, mean_present}.
Raises:: ValueError – If the cohort predates block-presence recording.
Return type:: dict

Example

>>> cohort.missingness_report(clusters)["ami_vs_pattern"] < 0.1
True

transform(donor_frames, *, extra=None)[source]#

Project held-out donors into this fitted basis (the only comparable path for new donors).

Rebuilds the per-locus identity/diversity/coverage blocks through the stored spaces and identity PCAs, re-adds any extra channels the caller supplies (same names/widths as at fit — the study owns those), then applies the fit-cohort impute medians and z-scores.

Parameters:

donor_frames (list[dict]) – Per held-out donor, {locus: chain clonotype frame}.
extra (dict | None) – {name: (m, k) array} for the non-core channels used at fit (isotype/…). Required iff the cohort was fit with extra_channels; the widths must match.

Returns:

(len(donor_frames), width) matrix in the fitted basis.

Return type:

ndarray

save(path)[source]#

Pickle the cohort: X/spec/rows/meta/identity PCAs + each locus’s basis.

Raises:: ValueError – If any locus basis uses a non-default matrix / alignment — load could not rebuild that coordinate system (see mir.repertoire._check_rebuildable()).
Return type:: None

classmethod load(path, *, verify=True)[source]#

Load a cohort, rebuilding every locus basis and verifying all prototype hashes.

Parameters:: verify (bool)
Return type:: DonorCohort

mir.cohort.fit_donor_embeddings(spaces, donor_frames, *, rows=None, id_pca=8, min_clones=5, coverage=True, extra_channels=None, standardize=True, impute=True, seed=0)[source]#

Fuse per-chain repertoire embeddings into one digital-donor matrix (§T.7).

For each donor and each fitted locus, computes the kernel-mean identity (cross-sample PCA-reduced), Hill diversity and log receptor-load coverage; per-chain blocks of the same kind merge under one channel name. Missing chains leave holes that ChannelBuilder imputes.

Parameters:

spaces (dict) – {locus: RepertoireSpace} from mir.repertoire.fit_repertoire_spaces().
donor_frames (list[dict]) – Per donor, {locus: chain clonotype frame} (already canonical-filtered / grouped by the study; the library only embeds). Row-aligned to rows.
rows (list | None) – Optional per-donor metadata, carried onto the cohort.
id_pca (int) – Cross-sample identity-PCA dimensionality per locus.
min_clones (int) – Skip a chain for a donor with fewer clonotypes (its blocks become holes).
coverage (bool) – Include the per-chain log-receptor-load channel.
extra_channels – Optional (rows, identity_concat) -> {name: (n, k) array} — the study’s own blocks (isotype / composition / atypicality …), fused into the same matrix. All are treated as non-attributable. identity_concat is the (n, Σ id_pca) per-chain identity, so e.g. mir.repertoire.centroid_atypicality() can be computed against it.
impute (bool) – Passed through to the final assembly (z-score / median-impute).
seed (int) – RNG seed for the identity PCA.
standardize (bool)
impute

Returns:

A DonorCohort; .X / .spec feed mir.explain.channel_report().

Return type:

DonorCohort

mir.cohort.residualize(X, group, *, shrink=False)[source]#

Subtract each group’s mean vector — remove the first-order batch offset (Prop. prop:batch).

The detect→correct→verify batch cookbook’s correction step, applied to a stacked donor matrix. NB the corrected matrix is a different coordinate system: never compare a residualized X to an uncorrected one.

⚠ The plain correction can make the batch easier to read, not harder, and this is measured rather than hypothetical: evaluated out-of-sample (leave-one-batch-out for any shared basis, plus a donor-level split inside each batch), batch-identity AUC went 0.863 raw → 0.985 after per-group centring, and 0.978 after ComBat. The mechanism is estimation error, not a bug — with a mean fitted from as few as 8–15 samples in ~1,280 dimensions, ‖µ̂ − µ‖ ≈ √(σ²d/n) was ≈16 against a true offset of norm 7–24, so subtracting it injects a batch-constant vector as large as the one it removes. In-sample this vanishes by construction, which is why only an out-of-sample evaluation sees it. shrink=True (positive-part James–Stein) recovered most of the damage in the same experiment: 0.985 → 0.889.

Parameters:

X (ndarray) – (n_samples, n_features) stacked matrix.
group (ndarray) – Length-n_samples batch label per row.
shrink (bool) – Scale each group’s offset by the positive-part James–Stein factor max(0, 1 − (d−2)·(σ̂²/n_g)/‖µ̂_g‖²) — where σ̂² is the pooled within-group per-coordinate variance — so a group whose apparent offset is no bigger than its own estimation error is left alone. Recommended whenever a group has few samples relative to the feature count; False keeps the exact previous behaviour.

Returns:

The corrected (n_samples, n_features) matrix.

Return type:

ndarray

mir.cohort.cluster_samples(embs, *, unbiased=True, method='dbscan', eps=None, min_samples=3, k=4, **kwargs)[source]#

Cluster repertoires by MMD — the repertoire-level analog of mir.bench.metrics.cluster() (TME states etc.).

Builds the MMD distance matrix (unbiased by default, since cohorts differ in depth) and clusters it with a precomputed-metric density estimator (mir.bench.metrics.cluster()).

Parameters:

embs – Per-sample SampleEmbedding (mean block present).
unbiased (bool) – Use the diagonal-removed MMD² (recommended for unequal-depth cohorts).
method (str) – "dbscan" / "hdbscan" / "optics".
eps (float | None) – DBSCAN radius; None → median of each sample’s k-th nearest MMD (self excluded).
k (int) – density / eps-estimation neighbours.
**kwargs – forwarded to the estimator.
min_samples (int)
k

Returns:

Cluster labels (-1 = noise), length len(embs).

mir.cohort.incidence_biomarkers(cohort, phenotype, *, pheno_col, match='1mm', min_incidence=3, **kwargs)[source]#

Cohort biomarker detection: per-clonotype subject-incidence Fisher test (Emerson 2017).

A thin delegate to vdjtools.biomarker.fisher.fisher_association() — the presence/absence biomarker that complements the geometry witness (mir.repertoire.class_witness()) and, at realistic donor counts, recovers public motifs the witness misses. match="1mm" groups single-mismatch metaclonotypes (the paper’s method; exact-match typically returns ~0 hits).

Parameters:

cohort – Pooled clonotype frame with a sample-id column (vdjtools cohort form).
phenotype – Per-sample frame carrying the binary pheno_col label.
pheno_col (str) – Phenotype column name in phenotype.
match (str) – "1mm" (metaclonotypes) or "exact".
min_incidence (int) – Minimum subject incidence to test a feature.
**kwargs – forwarded to fisher_association (key, scope, alternative, …).

Returns:

incidence, odds-ratio, p_value, BH q_value, direction.

Return type:

One row per feature

Exposure trajectory (`mir.track`)#

A PhenoPath-style (Campbell & Yau 2018) covariate-disentangled latent trajectory over a per-sample channel matrix — repertoire-level exposure detection, complementing mir.density’s clone-level TCRNET/ALICE enrichment.

Covariate-disentangled repertoire trajectory: a PhenoPath-style exposure pseudotime.

Exposure detection so far in mirpy is clone-level: mir.density.neighbor_enrichment() answers “which clonotypes look antigen-driven” against a background null. This module answers a different question at the repertoire level: given a cohort of samples with a known covariate (HLA type, vaccine arm, batch — anything already recorded) but an unknown or noisy progression axis (days since exposure, response stage), can we recover a shared latent trajectory while explicitly separating out which channels of the embedding respond to that trajectory differently depending on the covariate?

This is the model of PhenoPath (Campbell & Yau 2018, Nat. Commun. 9:2442, doi:10.1038/s41467-018-04696-6), adapted from genes × cells to repertoire channels × samples. For sample n and channel g:

Y[n, g] = c_g + alpha_g . x_n + (kappa_g + gamma_g . x_n) * tau_n + eps_{n,g}

x_n is the known covariate vector, tau_n the latent trajectory (pseudo-exposure/progression), kappa_g the channel’s baseline sensitivity to the trajectory, and gamma_g the covariate × trajectory interaction — the object of interest: a channel with large ‖gamma_g‖ evolves along the trajectory differently depending on the covariate (e.g. an HLA-restricted response channel that contracts over the response only in carriers). alpha_g is the covariate’s ordinary main effect, independent of progression.

Inference here is a simplified, closed-form alternative to PhenoPath’s coordinate-ascent variational inference: rather than a full mean-field posterior over every parameter, this alternates (a) per-channel ridge regression for (c, alpha, kappa, gamma) given tau, (b) a generalized-least-squares update of tau given the channel loadings, and (c) an iteratively- reweighted (empirical-Bayes) precision update on gamma that shrinks small interaction coefficients toward zero — the fast/approximate limit of PhenoPath’s ARD prior, not a literal reimplementation of its CAVI engine. Both blocks are ordinary linear-Gaussian conditionals, so every step is closed-form; this keeps the module torch-free (numpy only), consistent with mir.repertoire / mir.density.

Typical usage:

from mir.explain import stack_embeddings
from mir.track import fit_exposure_trajectory

embs = [sample_embedding(space, s) for s in samples]     # one shared RepertoireSpace
X, spec = stack_embeddings(embs)                         # channels: mean | diversity | second
covariates = hla_indicator                                # (n_samples, n_alleles) known covariate
fit = fit_exposure_trajectory(X, covariates, channel_names=spec.names)
fit.tau                                                    # inferred exposure/progression pseudotime
fit.top_interactions(top=5)                                # which channels are covariate-modulated

class mir.track.TrajectoryFit(tau, intercept, alpha, kappa, gamma, sigma2, ard_precision, channel_names, covariate_names, n_iter, converged, loglik_trace)[source]#

Bases: object

A fitted covariate-disentangled repertoire trajectory (PhenoPath-style, see module docstring).

Variables:

tau (numpy.ndarray) – (n_samples,) latent trajectory, standardized (mean 0, unit variance). Sign is fixed so tau correlates positively with the first principal component of the (standardized) input channels, so it is not arbitrary between refits.
intercept (numpy.ndarray) – (n_channels,) per-channel intercept c_g.
alpha (numpy.ndarray) – (n_channels, n_covariates) direct covariate effect (independent of tau).
kappa (numpy.ndarray) – (n_channels,) baseline trajectory loading — the channel’s progression sensitivity shared across covariate levels.
gamma (numpy.ndarray) – (n_channels, n_covariates) covariate x trajectory interaction — the channel’s progression sensitivity changes by this much per unit covariate. The quantity top_interactions() ranks channels by.
sigma2 (numpy.ndarray) – (n_channels,) residual variance per channel.
ard_precision (numpy.ndarray) – (n_channels, n_covariates) shrinkage precision applied to gamma (large -> that interaction is forced toward 0). Only meaningful when fit with ard=True.
covariate_names (channel_names /) – Labels carried through for reporting, if supplied.
n_iter (int) – Iterations actually run.
converged (bool) – Whether the relative change in (tau, gamma) fell below tol.
loglik_trace (list[float]) – Per-iteration Gaussian log-likelihood (up to the additive normalizing constant) — should increase monotonically-ish; a non-monotone trace signals a fitting problem.

Parameters:

tau (ndarray)
intercept (ndarray)
alpha (ndarray)
kappa (ndarray)
gamma (ndarray)
sigma2 (ndarray)
ard_precision (ndarray)
channel_names (list[str] | None)
covariate_names (list[str] | None)
n_iter (int)
converged (bool)
loglik_trace (list[float])

tau: ndarray#

intercept: ndarray#

alpha: ndarray#

kappa: ndarray#

gamma: ndarray#

sigma2: ndarray#

ard_precision: ndarray#

channel_names: list[str] | None#

covariate_names: list[str] | None#

n_iter: int#

converged: bool#

loglik_trace: list[float]#

property interaction_score: ndarray#: (n_channels,) covariate-trajectory interaction magnitude, ‖gamma_g‖ over covariates.

top_interactions(top=10)[source]#

Rank channels by covariate x trajectory interaction strength (descending).

Parameters:: top (int) – Number of channels to return.
Returns:: Columns channel, interaction_score (‖gamma_g‖), kappa (the channel’s baseline trajectory sensitivity, for context), sorted by interaction_score descending.
Return type:: DataFrame

mir.track.fit_exposure_trajectory(X, covariates, *, channel_names=None, covariate_names=None, ard=True, max_iter=200, tol=0.0001, seed=0)[source]#

Fit a covariate-disentangled latent trajectory over a per-sample channel matrix.

Parameters:

X (ndarray) – (n_samples, n_channels) per-sample feature matrix — any stacked repertoire channel matrix (mir.explain.stack_embeddings(), a ChannelBuilder build, or a raw embedding block). Standardized internally; the input scale doesn’t matter.
covariates (ndarray) – (n_samples,) or (n_samples, n_covariates) known covariate(s) — e.g. an HLA indicator, batch, vaccine arm. Continuous or one-hot; standardized internally.
channel_names (list[str] | None) – Optional length-n_channels labels, carried onto TrajectoryFit for top_interactions().
covariate_names (list[str] | None) – Optional length-n_covariates labels.
ard (bool) – Iteratively-reweighted shrinkage on gamma (see module docstring) — the empirical-Bayes approximation to PhenoPath’s ARD prior. False fits an unregularized interaction term (a fixed light ridge only, no shrinkage), useful for comparison / small cohorts where ARD’s extra sparsity isn’t wanted.
max_iter (int) – Maximum alternating-update iterations.
tol (float) – Convergence tolerance on the relative change of (tau, gamma) between iterations.
seed (int) – RNG seed for a tiny initialization jitter that breaks exact PC1 ties (e.g. a rank-deficient X); the fit is otherwise deterministic.

Returns:

A fitted TrajectoryFit.

Raises:

ValueError – On a row-count mismatch between X and covariates, or fewer than 3 samples.

Return type:

TrajectoryFit

mir.track.fit_exposure_trajectory_from_samples(space, samples, covariates, *, weight='log2p1', blocks=('mean', 'diversity'), **kwargs)[source]#

Convenience: build the channel matrix from repertoire samples, then fit the trajectory.

Parameters:

space – A shared RepertoireSpace (mir.repertoire.fit_repertoire_space()).
samples (list) – Per-sample clonotype frames (duplicate_count present), one SampleEmbedding per sample via mir.repertoire.sample_embedding().
covariates (ndarray) – (len(samples),) or (len(samples), n_covariates) known covariate(s).
weight (str) – Forwarded to mir.repertoire.sample_embedding().
blocks (tuple[str, ...]) – Forwarded to mir.repertoire.sample_embedding().
**kwargs – Forwarded to fit_exposure_trajectory() (ard, max_iter, tol, seed, covariate_names).

Returns:

A fitted TrajectoryFit, with channel_names set from the embeddings’ blocks.

Return type:

TrajectoryFit

Generative loop (`mir.generate`)#

A fitted (optionally class-conditional) density over RepertoireDescriptor vectors: sample new synthetic donor states, or evolve one along a coordinate via the fitted covariance’s conditional mean — the mechanical half of the generative loop (ROADMAP Phase 2).

The generative loop: a fitted density over repertoire descriptors, sample from it, evolve one.

Promotes the ad-hoc “in-silico evolution” pattern (a bare np.cov + a hardcoded conditional-mean slope over a hand-picked metric list, prototyped analysis-side over TCGA) into a real, reusable library object built on mir.repertoire.RepertoireDescriptor — the mass-preserving, decodable, smooth coordinate system mir.repertoire already ships ([infiltration, log n_eff, clonality] ‖ identity mean). DescriptorDensity is the (optionally class-conditional) Gaussian fit over that coordinate; sample() draws brand-new synthetic donor states from it (a digital twin generator — see mir.track for the complementary trajectory axis), and evolve() perturbs one coordinate of a real descriptor and propagates the coupled shift through every other coordinate via the fitted covariance’s multivariate- normal conditional mean — the same move as “hotter ⇒ diversity ↑ / switch ↑ / T-vs-B ↓” but as a library primitive instead of a one-off script.

This is the mechanical half of the generative loop (ROADMAP Phase 2): a linear-Gaussian density is what the shipped in-silico-evolution numbers were fit with, so it stays the default here. The research half — a genuinely non-linear, learned generative model over the full embedding (not just the compact descriptor) — is mir.ml.diffusion, which this module’s sample/evolve API deliberately mirrors so the two are drop-in alternatives.

What stays analysis-local (ROADMAP non-goals): fitting a survival/classification model on the perturbed metrics (e.g. “does hotter predict better survival”) is the study’s scorer, not this module’s job — mir.explain’s channel_report/the analysis’s own CoxPH fit consume the output of evolve(), this module never sees an outcome y.

Torch-free (numpy / sklearn).

Typical usage:

from mir.repertoire import sample_descriptor
from mir.generate import fit_descriptor_density, evolve

descriptors = [sample_descriptor(space, s) for s in samples]
density = fit_descriptor_density(descriptors, labels=tumor_type)   # one Gaussian per label

synthetic = density.sample(50, condition="SKCM")                   # 50 new synthetic donor states
hotter = evolve(density, descriptors[0], coordinate="infiltration", delta=2.0, condition="SKCM")
hotter.metrics()                                                   # named metrics after the move

class mir.generate.DescriptorDensity(mean=<factory>, cov=<factory>, dim=0)[source]#

Bases: object

A fitted (optionally class-conditional) Gaussian density over descriptor vectors.

Variables:

mean (dict) – {label: (dim,) mean} — a single key None when unconditional.
cov (dict) – {label: (dim, dim) shrinkage covariance}.
dim (int) – Vector width (len(RepertoireDescriptor.vector)).

Parameters:

mean (dict)
cov (dict)
dim (int)

mean: dict#

cov: dict#

dim: int = 0#

classmethod fit(X, labels=None)[source]#

Fit a Gaussian (or, with labels, one Gaussian per label) over descriptor vectors.

Parameters:

X (ndarray) – (n_samples, dim) stacked descriptor vectors (e.g. np.stack([d.vector for d in descriptors])).
labels (list | None) – Optional length-n_samples group labels (tumor type, batch, condition, …) — fits one mean/covariance per distinct label instead of a single pooled Gaussian.

Returns:

A fitted DescriptorDensity.

Raises:

ValueError – If X has fewer than 2 rows.

Return type:

DescriptorDensity

sample(n=1, *, condition=None, seed=0)[source]#

Draw n synthetic descriptor vectors from the fitted (conditional) Gaussian.

Parameters:

n (int) – Number of synthetic vectors to draw.
condition – Label to condition on (must be one of the labels fit saw); None for an unconditional density, or the sole label when only one was fit.
seed (int) – RNG seed.

Returns:

(n, dim) array — decode with mir.repertoire.decode_metrics() per row.

Return type:

ndarray

evolve(vector, *, coordinate, delta, condition=None)[source]#

Shift vector by delta along one coordinate, propagating the coupled response.

The multivariate-normal conditional mean: fixing coordinate d to move by delta shifts every other coordinate j by Sigma[j,d]/Sigma[d,d] * delta in expectation — this is the “hotter -> diversity up / switch up / T-vs-B down” move, generalized from a hand-picked metric list to the full descriptor (and, via RepertoireDescriptor.mean, the identity block too).

Parameters:

vector (ndarray) – A (dim,) descriptor vector to perturb (e.g. one real donor’s).
coordinate (int | str) – Vector index, or one of "infiltration" / "log_neff" / "clonality" for the named scalar coordinates.
delta (float) – Amount to shift coordinate by (in the descriptor’s own units — e.g. log-mass).
condition – Which fitted Gaussian’s covariance defines the coupling (see sample()).

Returns:

The perturbed (dim,) vector — decode with mir.repertoire.decode_metrics().

Return type:

ndarray

mir.generate.fit_descriptor_density(descriptors, labels=None)[source]#

Fit a DescriptorDensity directly from RepertoireDescriptor objects.

Parameters:

descriptors (list[RepertoireDescriptor]) – One RepertoireDescriptor per sample (mir.repertoire.sample_descriptor()).
labels (list | None) – Optional per-sample group labels — see DescriptorDensity.fit().

Returns:

A fitted DescriptorDensity.

Return type:

DescriptorDensity

mir.generate.evolve(density, descriptor, *, coordinate, delta, condition=None)[source]#

Perturb one RepertoireDescriptor and rebuild the result as one.

Thin wrapper around DescriptorDensity.evolve() that unpacks/repacks RepertoireDescriptor instead of a bare vector, so the result’s .metrics() reads off directly.

Parameters:

density (DescriptorDensity) – A fitted DescriptorDensity (the coupling comes from its covariance).
descriptor (RepertoireDescriptor) – The descriptor to perturb.
coordinate (int | str) – See DescriptorDensity.evolve().
delta (float) – See DescriptorDensity.evolve().
condition – See DescriptorDensity.evolve().

Returns:

The perturbed RepertoireDescriptor.

Return type:

RepertoireDescriptor

Digital twin (`mir.twin`)#

One donor’s perturbable, simulatable state — glues RepertoireDescriptor + an optional mir.track trajectory position + a mir.generate/mir.ml.diffusion generator into one object.

The digital twin: one donor’s compact, perturbable, simulatable state.

Where DonorCohort fuses a cohort’s measured state into one comparable matrix for scoring, this module wraps one donor’s state as an object that can be perturbed (mir.generate.evolve(), a “what if” move) or used to seed brand-new synthetic realizations from a fitted generator — mir.generate.DescriptorDensity (linear-Gaussian) or mir.ml.diffusion.DiffusionModel (non-linear, needs [ml]) — both of which share the same sample(n, *, condition=None, seed=0) call shape, so either drops in unchanged. A twin optionally carries this donor’s position on a fitted TrajectoryFit (where along an inferred exposure/progression axis it sits) and its known covariate, so simulate() can resample condition-matched to the donor’s own group by default.

This closes the loop the rest of the library already opened: mir.cohort measures a donor, mir.track locates it on a trajectory, mir.generate / mir.ml.diffusion generate new states — DonorTwin is the one object a caller perturbs or resamples through, instead of threading three APIs together by hand for every donor.

Torch-free itself; the diffusion generator (only if actually passed to simulate()) is imported lazily by the caller, not by this module.

Typical usage:

from mir.twin import make_twins
from mir.generate import fit_descriptor_density

density = fit_descriptor_density(descriptors, labels=tumor_type)
twins = make_twins(descriptors, conditions=tumor_type, donor_ids=sample_ids)

hotter = twins[0].perturb(density, coordinate="infiltration", delta=2.0)  # this donor, "what if hotter"
synthetic = twins[0].simulate(density, n=20)                              # 20 new synthetic peers

class mir.twin.DonorTwin(descriptor, tau=None, condition=None, donor_id=None)[source]#

Bases: object

One donor’s digital twin: measured state + optional trajectory position + covariate.

Variables:

descriptor (mir.repertoire.RepertoireDescriptor) – The donor’s measured RepertoireDescriptor.
tau (float | None) – This donor’s position on a fitted TrajectoryFit, if known (typically fit.tau[donor_index]).
condition (object | None) – The donor’s known covariate/group label (tumor type, HLA, batch, …) — the default conditioning value for perturb() / simulate().
donor_id (str | None) – Optional identifier, carried through unchanged by every method here.

Parameters:

descriptor (RepertoireDescriptor)
tau (float | None)
condition (object | None)
donor_id (str | None)

descriptor: RepertoireDescriptor#

tau: float | None = None#

condition: object | None = None#

donor_id: str | None = None#

metrics()[source]#

Named metrics of the current descriptor (metrics()).

Return type:: dict

perturb(density, *, coordinate, delta, condition=None)[source]#

Move this twin along one descriptor coordinate (mir.generate.evolve()).

Parameters:

density – A fitted mir.generate.DescriptorDensity — the coupling between coordinates comes from its covariance.
coordinate – Vector index, or "infiltration" / "log_neff" / "clonality".
delta (float) – Amount to shift coordinate by (in the descriptor’s own units).
condition – Which of density’s fitted groups defines the coupling; defaults to this twin’s own condition.

Returns:

A new DonorTwin with the perturbed descriptor — tau, condition and donor_id carry over unchanged (a perturbation is a “what if”, not a re-measurement of a different donor).

Return type:

DonorTwin

simulate(generator, n=1, *, condition=None, seed=0, **kwargs)[source]#

Draw n new synthetic donor states from a fitted generator, decoded to named metrics.

Parameters:

generator – A fitted mir.generate.DescriptorDensity or mir.ml.diffusion.DiffusionModel — both share the sample(n, *, condition, seed) call shape; either drops in unchanged.
n (int) – Number of synthetic states to draw.
condition – Which fitted group to sample from; defaults to this twin’s own condition.
seed (int) – RNG seed.
**kwargs – Forwarded to generator.sample (e.g. a DiffusionModel’s steps=/guidance_scale=; unused by DescriptorDensity).

Returns:

n decoded metric dicts (mir.repertoire.decode_metrics()) — synthetic “what a donor like this typically looks like” states, not perturbations of this specific twin (see perturb() for that).

Return type:

list[dict]

mir.twin.make_twins(descriptors, *, tau=None, conditions=None, donor_ids=None)[source]#

Zip a cohort’s descriptors (+ optional trajectory/condition/id) into one DonorTwin each.

Parameters:

descriptors (list[RepertoireDescriptor]) – One RepertoireDescriptor per donor.
tau (ndarray | None) – Optional (n_donors,) trajectory position per donor (a fitted TrajectoryFit’s .tau), row-aligned to descriptors.
conditions (list | None) – Optional per-donor covariate/group labels, row-aligned.
donor_ids (list[str] | None) – Optional per-donor identifiers, row-aligned.

Returns:

One DonorTwin per donor, in the same order as descriptors.

Raises:

ValueError – If a supplied tau / conditions / donor_ids length disagrees with descriptors.

Return type:

list[DonorTwin]

Benchmark harness (`mir.bench`)#

`mir.bench.vdjdb`#

VDJdb loader + antigen subsets.

Load VDJdb dumps into the AIRR polars frame TCREmp expects.

VDJdb “slim” dumps use dotted column names (cdr3, v.segm, j.segm, antigen.epitope, gene, complex.id, mhc.class). This maps them to the vdjtools.io.schema names plus an epitope label column.

mir.bench.vdjdb.load_vdjdb(path)[source]#

Return an AIRR frame with v_call, j_call, junction_aa, locus, epitope.

Parameters:: path (str) – VDJdb slim TSV (optionally gzipped).
Return type:: DataFrame

mir.bench.vdjdb.antigen_subset(df, chain, min_records)[source]#

Rows for chain restricted to epitopes with >= min_records records.

Parameters:

df (DataFrame)
chain (str)
min_records (int)

Return type:

DataFrame

`mir.bench.metrics`#

Clustering (DBSCAN/HDBSCAN/OPTICS) + F1 / retention / purity.

Antigen-cluster metrics for the TCREMP benchmark (paper Table 1 / S1).

DBSCAN clustering with a kneedle-selected eps on embedded TCRs, then per-antigen F1 (on clustered TCRs, via cluster→majority-label assignment) and retention (fraction of an antigen’s TCRs that land in a cluster), following the decoupling used in the paper: F1 measures clustering quality, retention measures coverage.

mir.bench.metrics.estimate_dbscan_eps(X, k=4)[source]#

Kneedle eps: the knee of the sorted k-th nearest-neighbour distance curve.

Parameters:

X (ndarray)
k (int)

Return type:

float

mir.bench.metrics.cluster(X, eps=None, min_samples=3, k=4, eps_factor=0.4, method='dbscan', **kwargs)[source]#

Density-cluster an embedding; noise = -1 for every method.

Parameters:

method (str) – "dbscan" (default), "hdbscan", or "optics". All three are density estimators that emit -1 noise, so cluster_metrics() and mir.density.denoise_and_cluster() work unchanged across them. "dbscan" is the default for reproducibility of the published Table-S1 numbers.
k (int) – DBSCAN only. When eps is None it is kneedle_knee × eps_factor — the paper’s eps := (distance at knee) × (dataset-specific factor) (Fig 1); the raw knee over-merges in the PCA embedding, so ~0.3–0.4 recovers tight antigen clusters.
min_samples (int) – core-point / density threshold, shared by all three methods.
**kwargs – passed to the underlying estimator (e.g. min_cluster_size for HDBSCAN, cluster_method / xi for OPTICS).
X (ndarray)
eps (float | None)
k
eps_factor (float)

Return type:

ndarray

HDBSCAN is a variable-density estimator — the persistent cluster tree (Hartigan level sets) rather than DBSCAN’s single global eps slice — so it is the natural choice when local density spans orders of magnitude across antigen ridges (appendix sec:dens-depth).

class mir.bench.metrics.AntigenMetric(epitope, n, f1, precision, recall, retention)[source]#

Bases: object

Per-antigen clustering quality: F1 of the recovered cluster plus retention.

Variables:

epitope (str) – The antigen epitope this metric scores.
n (int) – Number of clonotypes labelled with this epitope.
f1 (float) – F1 of the best-matching predicted cluster against the epitope’s clonotypes.
precision (float) – Precision of that best-matching cluster.
recall (float) – Recall of that best-matching cluster.
retention (float) – Fraction of the epitope’s clonotypes that were clustered (not noise).

Parameters:

epitope (str)
n (int)
f1 (float)
precision (float)
recall (float)
retention (float)

epitope: str#

n: int#

f1: float#

precision: float#

recall: float#

retention: float#

mir.bench.metrics.cluster_metrics(labels, antigens)[source]#

Per-antigen F1 (on clustered TCRs) and retention.

Parameters:

labels (ndarray) – DBSCAN cluster labels (-1 = noise/unclustered).
antigens – True antigen (epitope) label per TCR, same length as labels.

Returns:

epitope -> AntigenMetric.

Return type:

dict[str, AntigenMetric]

`mir.bench.theory`#

Reproduced supplementary theory (S1–S3, T5–T6, codec losslessness).

Theory-validation experiments (reproduce TCREMP supplementary S1–S3).

These numerically check the properties the embedding is claimed to have, using the actual v3 pipeline distances (seqtree.gapblock junction dissimilarity):

S2 / Theory T1 — Euclidean distance D_ij in embedding space tracks the pairwise junction dissimilarity d_ij (s2_dissimilarity_distance_correlation(); paper Pearson R ≈ 0.56). Using the sequences as their own prototypes, the embedding of sequence i is row i of the dissimilarity matrix and D_ij = ‖d_i − d_j‖₂.
S1 / Theory T4 — d_ij follows a Gamma law (Gaussian fails); D_ij follows a Fréchet (GEV, shape ξ>0) law (fit_distributions()).
S3 — distances from real vs model prototypes agree (prototype_source_correlation(); paper Pearson R ≈ 0.96), i.e. the prototype source barely matters.

mir.bench.theory.junction_dissimilarity(cdr3s, threads=0)[source]#

Symmetric (n, n) junction dissimilarity d_ij (v3 gapblock; a negative-type semimetric / squared Hilbert distance, not a metric — ρ=√d is).

Parameters:: threads (int)
Return type:: ndarray

mir.bench.theory.junction_dissimilarity_sw(cdr3s)[source]#

Paper-exact d_ij via Smith-Waterman BLOSUM62 (supplementary S1/S2).

d_ij = s_ii + s_jj − 2·s_ij from a local BLOSUM62 alignment (linear gap). Requires BioPython (the [build] extra). O(n²) alignments — use for the theory validation on a few thousand sequences, not at scale.

Return type:: ndarray

class mir.bench.theory.CorrelationResult(n: 'int', n_pairs: 'int', pearson: 'float', d: 'np.ndarray', D: 'np.ndarray')[source]#

Bases: object

Parameters:

n (int)
n_pairs (int)
pearson (float)
d (ndarray)
D (ndarray)

n: int#

n_pairs: int#

pearson: float#

d: ndarray#

D: ndarray#

mir.bench.theory.s2_dissimilarity_distance_correlation(cdr3s, threads=0, dissimilarity='gapblock')[source]#

Correlate junction dissimilarity d_ij with embedding distance D_ij.

Uses the self-prototype construction from supplementary Fig. S2: the embedding of each sequence is its row of the dissimilarity matrix, and D_ij is the Euclidean distance between rows. Returns the flat i>j vectors and Pearson R.

Parameters:

dissimilarity (str) – "gapblock" (the v3 pipeline metric) or "sw" (paper-exact Smith-Waterman BLOSUM62).
threads (int)

Return type:

CorrelationResult

mir.bench.theory.fit_distributions(d_flat, D_flat, sample=200000, seed=0)[source]#

Fit S1 distributions: Gamma vs Normal for d, GEV/Fréchet vs Normal for D.

Returns a nested dict of {ks, aic, params} per fit plus the GEV shape xi (xi > 0 ⇒ Fréchet). Large inputs are subsampled to sample points.

Parameters:

d_flat (ndarray)
D_flat (ndarray)
sample (int)
seed (int)

Return type:

dict

mir.bench.theory.shm_embedding_drift(cdr3s, prototypes, max_mut=8, n_rep=3, seed=0, threads=0)[source]#

Embedding drift vs mutation load (Theory T5, the IGH/SHM case).

Applies k random interior substitutions (a somatic-hypermutation proxy) to each CDR3 and measures the junction-embedding Euclidean distance to the original, D_k = ‖φ(mutated) − φ(original)‖. Returns k -> (mean D_k, std). The claim: the drift is bounded by the mutation load — D_k grows ~linearly in k — so heavily hypermutated IGH sequences sit controllably far from germline in embedding space.

Parameters:

max_mut (int)
n_rep (int)
seed (int)
threads (int)

Return type:

dict[int, tuple[float, float]]

mir.bench.theory.prototype_source_correlation(query_cdr3s, real_prototypes, model_prototypes, threads=0)[source]#

Correlate embedding distances built from real vs model prototypes (S3).

Embeds query_cdr3s against each prototype set (junction only), takes pairwise Euclidean distances under each embedding, and returns their Pearson R (paper ≈ 0.96).

Parameters:: threads (int)
Return type:: dict

mir.bench.theory.tcrnet_convergence(obs_cdr3s, bg_cdr3s, prototypes, *, n_components=25, scales=(0.5, 1.0, 1.5, 2.0, 3.0), seed=0, threads=0)[source]#

T6: continuous embedding enrichment converges to discrete Hamming-1 enrichment.

Embeds observed and background CDR3s against the prototype set (junction only), reduces with one shared StandardScaler → PCA, and at radius scale × r₁ (r₁ = the median one-substitution embedding drift) counts each observed clonotype’s neighbours in embedding space. Returns the Spearman correlation between those continuous counts and the discrete Hamming-1 neighbour counts. The correlation is high at small radii and fades as the radius grows past one substitution — numerically confirming that graph neighbour-enrichment (TCRNET/ALICE) is the r→0 limit of the density ratio E(z).

Parameters:

obs_cdr3s – Observed junction sequences.
bg_cdr3s – Background junction sequences (P_gen sample / control).
prototypes – Prototype junctions defining the embedding.
n_components (int) – PCA dimensionality of the shared coordinate system.
scales – Radius multiples of the one-substitution scale r₁ to evaluate.
seed (int) – RNG seed for the calibration mutations and PCA solver.
threads (int)

Returns:

{radius_1sub, hamming1_mean, spearman_by_scale, spearman_at_1sub}.

Return type:

dict

mir.bench.theory.codec_losslessness(codes, seqs, recon=None, *, eps=1e-06)[source]#

Measurable losslessness of a codec, split into decoder-independent vs -dependent parts.

Informational ceiling (decoder-independent, from codes alone). The compact code is injective iff distinct sequences map to distinct codes; two distinct sequences whose codes coincide can never both be recovered, so exact_ceiling = 1 − collision_rate upper-bounds any decoder’s exact-match — the information the code retains, separate from how well a decoder reads it. Encoding is deterministic (identical sequences ⇒ identical codes), so duplicate sequences are collapsed first and a collision is a near-coincident code pair between distinct sequences (nearest-code distance < eps). nn_dist_{median,min} report the code-resolution scale so eps can be read in context.

Reconstructive loss (needs recon, the decoded strings row-aligned to seqs): exact_match, mean edit distance (a graded distortion, unlike binary exact-match), and per-position token accuracy over the length-40 frame — localizing where loss falls (conserved C…[FW] anchors vs the specificity-bearing variable middle). Per-position accuracy re-encodes recon into the canonical frame, so it measures positional token agreement.

Parameters:

codes – (n, m) compact codes, one per sequence (e.g. whitened-PCA junction code).
seqs – the n true junction strings.
recon – optional n decoded strings (e.g. InverseDecoder.decode(codes)). When None, only the decoder-independent ceiling is returned.
eps (float) – code-distance below which two distinct sequences count as colliding.

Returns:

Always {n, n_unique, collision_rate, exact_ceiling, nn_dist_median, nn_dist_min}; and when recon is given, additionally {exact_match, mean_edit, token_acc, anchor_acc, middle_acc, pos_acc} (pos_acc a length-40 array).

Return type:

dict

`mir.bench.eval`#

Scorers for the explainable readout: cross-validated AUC, Cox C-index, log-rank — plus recovery_report, the grouped-CV ridge that asks whether each basic repertoire statistic is carried inside the embedding.

Scorers for the explainable cohort readout (mir.explain) — the closures you hand to mir.explain.channel_report().

mir.explain is deliberately scorer-agnostic: it slices channels and never sees y, so the model choice lives here, in the caller’s closure. This module collects the recurring choices — cross-validated classification AUC and Cox survival C-index — as small, reusable functions so every cohort benchmark stops re-gluing them. A scorer takes a feature block and returns a float where higher is better:

from mir.bench.eval import cv_auc, cv_cindex
from mir.explain import stack_embeddings, channel_report

X, spec = stack_embeddings(embs)
rep = channel_report(X, spec, lambda B: cv_auc(B, y)[0], base=0.5)              # classification
rep = channel_report(X, spec, lambda B: cv_cindex(dur, evt, base=C, block=B, n_pc=8),
                     base=cv_cindex(dur, evt, base=C, block=None), mode="both")  # survival

Needs the [bench] extra (scikit-learn always; lifelines for the survival scorers; vdjtools for the k-mer baseline). Everything is lazily imported so importing this module stays cheap.

mir.bench.eval.held_out_auc(Xtr, ytr, Xte, yte, *, pca_cols=0)[source]#

Held-out ROC-AUC of a logistic head; the first pca_cols columns get in-fold PCA(0.9 var).

Wide blocks (kernel mean, second moment, k-mer) overfit unless reduced inside the fold; the few diversity/coverage features pass through raw so their signal isn’t buried. pca_cols=0 = no PCA.

Parameters:: pca_cols (int)
Return type:: float

mir.bench.eval.cv_auc(X, y, *, pca_cols=0, n_splits=5, n_repeats=10, seed=0)[source]#

Repeated stratified k-fold AUC as (mean, std) — a CI, not a single-split point estimate.

A single 70/30 split at small n has AUC SD ≈ 0.5/√n_test (~0.1 for n_test≈30), so point estimates are near-meaningless; repeated CV exposes whether two methods’ intervals separate. Pass cv_auc(B, y)[0] as the channel_report scorer.

Parameters:

pca_cols (int)
n_splits (int)
n_repeats (int)
seed (int)

mir.bench.eval.cv_cindex(durations, events, *, base=None, block=None, n_pc=0, n_splits=5, seed=0, penalizer=0.1)[source]#

5-fold CV Cox C-index of base covariates + an optional feature block.

The survival analog of cv_auc(), and the scorer for survival channel reports. base is the clinical design (e.g. age+sex+stage+log-reads — whatever the study built; this function is schema-agnostic) and block is the feature block being scored; when block is wider than n_pc it is PCA-reduced to n_pc components inside each fold (train-only fit) so a wide kernel-mean block doesn’t overfit. Returns the mean concordance across folds (nan-robust).

Use as lambda B: cv_cindex(dur, evt, base=C, block=B, n_pc=8) with base score cv_cindex(dur, evt, base=C, block=None).

Parameters:

n_pc (int)
n_splits (int)
seed (int)
penalizer (float)

Return type:

float

mir.bench.eval.km_logrank(durations, events, groups)[source]#

Multivariate log-rank p-value across groups — do the KM survival curves differ?

The test behind a Kaplan–Meier stratification (e.g. TME states from cluster_samples). Returns the p-value.

Return type:: float

mir.bench.eval.recovery_report(X, stats, groups=None, *, n_pc=20, n_splits=5, alpha=1.0, seed=0)[source]#

Grouped-CV ridge R² from an embedding’s PCs back to each basic repertoire statistic.

The honest scoring rule for a repertoire embedding is not “does Φ beat the one-number marker” — that is a competition between an object and one of its own coordinates. It is recoverability: is the statistic carried inside the embedding, so nothing has to be bolted on beside it? A high R² says yes; a low one says the embedding threw that information away and a downstream model will need the statistic as a separate feature.

This is where a sub-probability embedding earns its keep by construction. Renormalising every sample to mass 1 (mir.repertoire.sample_embedding()’s default) deletes the magnitude, so coverage- and richness-like statistics are unrecoverable from Φ no matter how good the model — whereas the deficient measure (mir.repertoire.contrast_embedding()) keeps them. Sits beside mir.cohort.missingness_report() as the other “is this object honest” check.

Report it together with the increment: the endpoint score for the embedding, for the statistic alone, and for both together (cv_auc() / cv_cindex()).

Parameters:

X (ndarray) – (n_samples, d) embedding matrix.
stats (dict) – {name: (n_samples,) values} — the basic statistics to recover (clonotype richness, Shannon, top-clone fraction, singleton fraction, Chao unseen fraction, library size, AIRR-per-million …). Non-finite entries are dropped per statistic.
groups – Optional (n_samples,) grouping (subject / batch / cohort) — folds are then grouped so a subject never appears in both train and test. None → plain shuffled k-fold.
n_pc (int) – PCA components fit inside each training fold; 0 uses the raw columns.
n_splits (int) – Number of folds.
alpha (float) – Ridge penalty.
seed (int) – Fold RNG seed.

Returns:

{name: R²} from out-of-fold predictions (can be negative — worse than the mean).

Raises:

ValueError – If a statistic’s length disagrees with X’s row count.

Return type:

dict

Example

>>> recovery_report(X, {"richness": rich, "top_clone": top}, groups=subject)
{'richness': 0.81, 'top_clone': 0.44}

mir.bench.eval.kmer_matrix(frames, k=3)[source]#

Sample × k-mer frequency matrix (vdjtools kmer_profile long-form → wide) — the classic baseline.

Parameters:

frames – One clonotype polars.DataFrame per sample (with junction_aa + counts).
k (int) – k-mer length.

Returns:

(len(frames), n_kmers) frequency matrix over the pooled vocabulary.

Return type:

ndarray

Neural codecs and learned encoders (`mir.ml`)#

`mir.ml.tokenize`#

CDR3 tokenisation for the neural codecs.

Fixed-length-40 one-hot tokenization of CDR3/junction sequences.

Mirrors the irrm-codec representation: every sequence is placed into a length-40 frame by anchoring the conserved ends and inserting a contiguous gap block in the variable middle — the N-terminal 4 residues left-align at the start, the C-terminal 3 residues right-align at the end, the middle left-aligns after position 4, and gaps fill the remainder (mirrors the TCREMP gap placement). Long IGH middles that overflow are truncated (the hard case). The alphabet is the 20 amino acids plus a gap symbol (21 tokens); unknown characters map to gap.

mir.ml.tokenize.encode_indices(cdr3s)[source]#

Return (N, 40) int64 token indices (gap-padded, ends anchored).

Return type:: ndarray

mir.ml.tokenize.encode_onehot(cdr3s)[source]#

Return (N, 40, 21) float32 one-hot encoding.

Return type:: ndarray

`mir.ml.encoder`#

Forward encoder: sequence → compact embedding code.

Forward codec: a CNN that maps a CDR3 sequence to its TCREMP embedding.

The GPU-friendly approximation of the prototype mapping (irrm-codec forward model). Trained with free supervision — targets are TCREMP embeddings computed by mir.embedding.tcremp.TCREmp, so labelled data is never the bottleneck.

class mir.ml.encoder.SequenceEncoder(*args, **kwargs)[source]#

Bases: Module

1-D CNN over the length-40 one-hot CDR3 → a fixed embedding vector.

Parameters:

embed_dim (int) – Output dimensionality (the target embedding width).
channels (tuple[int, ...]) – Conv channel sizes.
kernel (int) – Conv kernel size.
hidden (int) – Width of the dense head.
dropout (float) – Dropout on the dense head.

forward(x)[source]#

Parameters:: x (torch.Tensor)
Return type:: torch.Tensor

`mir.ml.decoder`#

Inverse decoder: code → sequence.

Inverse codec: reconstruct a CDR3 sequence from its (compact) embedding.

A parallel decoder (irrm-codec style): the compact PCA code is expanded and refined by 1-D convolutions into length-40 per-position token logits in one shot (no autoregression). Trained with cross-entropy against the fixed-length-40 tokenization; the reconstructed string drops gap tokens.

class mir.ml.decoder.SequenceDecoder(*args, **kwargs)[source]#

Bases: Module

Compact embedding code → (B, 40, 21) per-position token logits.

Parameters:

code_dim (int) – Input code width (e.g. #PCs of the compact embedding).
hidden (int) – Width of the expansion MLP.
channels (tuple[int, ...]) – Conv refinement channels (first is the reshaped depth).
kernel (int) – Conv kernel size.

forward(z)[source]#

Parameters:: z (torch.Tensor)
Return type:: torch.Tensor

mir.ml.decoder.tokens_to_seq(idx_row)[source]#

Map a length-40 token-index row back to an amino-acid string (gaps dropped).

Return type:: str

`mir.ml.train`#

Training loops + device selection (CUDA → MPS → CPU).

Train the forward codec (CDR3 → TCREMP embedding) with free supervision.

Targets are TCREMP embeddings computed by mir.embedding.tcremp.TCREmp, so training data is unlimited — sample sequences from vdjtools.model.generate (or a real repertoire) and embed them. Reports test-set mean cosine similarity (irrm-codec’s forward metric).

mir.ml.train.pick_device(device=None)[source]#

Best available accelerator: explicit device override, else CUDA, else Apple mps, else cpu.

Set MIR_DEVICE (e.g. cuda:1) to override without threading device= through every call.

Parameters:: device (str | None)
Return type:: torch.device

mir.ml.train.seed_everything(seed)[source]#

Parameters:: seed (int)
Return type:: None

class mir.ml.train.ForwardEncoder(model, transform, device, is_pca)[source]#

Bases: object

A trained forward codec: CDR3 strings → embedding vectors.

transform maps the model output back to the original embedding space — either an inverse StandardScaler (raw target) or an inverse whitened PCA (compact 95%-variance target). code returns the compact code.

Parameters:

model (nn.Module)
device (torch.device)
is_pca (bool)

code(cdr3s, batch=1024)[source]#

Return the compact code the model predicts (PCA coords, or normalized target).

Parameters:: batch (int)
Return type:: ndarray

encode(cdr3s, batch=1024)[source]#

Predict embeddings in the original space (inverse PCA / scaler).

Parameters:: batch (int)
Return type:: ndarray

mir.ml.train.train_forward_encoder(cdr3s, targets, *, target_pca=0.95, epochs=40, batch=256, lr=0.001, val_frac=0.1, test_frac=0.1, seed=0, device=None, verbose=True)[source]#

Train a SequenceEncoder to predict targets from cdr3s.

Parameters:

target_pca (float | None) – If a float in (0, 1], compact the (redundant) target embedding with a whitened PCA keeping that fraction of variance — fit on the train split only. None trains on the raw standardized target.
targets (ndarray)
epochs (int)
batch (int)
lr (float)
val_frac (float)
test_frac (float)
seed (int)
device (str | None)
verbose (bool)

Return type:

tuple[ForwardEncoder, dict]

Returns the fitted ForwardEncoder and a metrics dict (test_cosine — reconstruction cosine in the original embedding space — val_mse, n, n_components). No leakage: PCA / scaler are fit on train only.

class mir.ml.train.InverseDecoder(model, code_scaler, device)[source]#

Bases: object

A trained inverse codec: embedding codes → CDR3 strings.

Parameters:

model (nn.Module)
code_scaler (StandardScaler)
device (torch.device)

decode(codes, batch=1024)#

Parameters:: batch (int)
Return type:: list[str]

mir.ml.train.train_inverse_decoder(codes, cdr3s, *, epochs=50, batch=256, lr=0.001, val_frac=0.1, test_frac=0.1, seed=0, device=None, verbose=True)[source]#

Train a SequenceDecoder to reconstruct cdr3s from codes.

Returns the fitted InverseDecoder and metrics (exact_match — full sequence reconstructed correctly — token_acc, n). Codes are standardized on the train split only.

Parameters:

codes (ndarray)
epochs (int)
batch (int)
lr (float)
val_frac (float)
test_frac (float)
seed (int)
device (str | None)
verbose (bool)

Return type:

tuple[InverseDecoder, dict]

class mir.ml.train.PgenRegressor(model, scaler, device)[source]#

Bases: object

A trained regressor: CDR3 strings → log10 Pgen.

Parameters:

model (nn.Module)
scaler (StandardScaler)
device (torch.device)

predict(cdr3s, batch=1024)#

Parameters:: batch (int)
Return type:: ndarray

mir.ml.train.train_pgen_regressor(cdr3s, log_pgen, *, epochs=40, batch=256, lr=0.001, val_frac=0.1, test_frac=0.1, seed=0, device=None, verbose=True)[source]#

Train the shared sequence encoder to predict log10 Pgen from CDR3.

Returns the fitted PgenRegressor and metrics (pearson, rmse in log10 units, n). Targets standardized on the train split only.

Parameters:

epochs (int)
batch (int)
lr (float)
val_frac (float)
test_frac (float)
seed (int)
device (str | None)
verbose (bool)

Return type:

tuple[PgenRegressor, dict]

`mir.ml.codec`#

Unified encoder+decoder codec with a geometry-anchor term.

Unified codec: jointly train encoder + decoder, anchored to the embedding geometry.

The encoder (seq → code) and decoder (code → seq) are trained together, but a geometry-preservation term keeps the encoder’s code close to the true (compact PCA) embedding — so the codec round-trips sequences without the code drifting away from the embedding space that makes distances meaningful (Theory T1). lambda_embed controls the anchor strength: large → geometry preserved, decoder adapts to it.

class mir.ml.codec.UnifiedCodec(encoder, decoder, device)[source]#

Bases: object

Trained encoder+decoder: encode (seq→code), decode (code→seq), roundtrip.

Parameters:

encoder (nn.Module)
decoder (nn.Module)
device (torch.device)

encode(cdr3s, batch=1024)#

Parameters:: batch (int)
Return type:: ndarray

decode(codes, batch=1024)#

Parameters:: batch (int)
Return type:: list[str]

roundtrip(cdr3s, batch=1024)[source]#

Parameters:: batch (int)
Return type:: list[str]

mir.ml.codec.train_unified_codec(cdr3s, codes, *, lambda_embed=1.0, epochs=60, batch=256, lr=0.001, val_frac=0.1, test_frac=0.1, seed=0, device=None, verbose=True)[source]#

Jointly train encoder+decoder to reconstruct cdr3s while matching codes.

Parameters:

codes (ndarray) – The true compact embedding codes (whitened PCA), one per sequence — the geometry anchor and the encoder target.
lambda_embed (float) – Weight of the encoder geometry term vs the reconstruction term.
epochs (int)
batch (int)
lr (float)
val_frac (float)
test_frac (float)
seed (int)
device (str | None)
verbose (bool)

Return type:

tuple[UnifiedCodec, dict]

Returns the UnifiedCodec and metrics: encode_cosine (encoder code vs true code — geometry preservation), roundtrip_exact (seq→code→seq), and decode_true_exact (decode the true code — reference upper bound).

`mir.ml.bundle`#

CodecBundle — prototype-hash-verified shipping of a trained codec.

Ship a trained codec together with the things that define its embedding space.

Two embeddings are only comparable if they share the same prototype set and the same PCA rotation — different prototypes or a different PCA give an incomparable coordinate system. A CodecBundle therefore serializes, as one artifact:

the PCA transform (the rotation/whitening — the exact coordinate system),
a prototype hash (identity of the prototype set the embedding was built on),
the model weights and the metadata needed to reproduce/verify it.

Reconstructing a usable encoder (CodecBundle.forward_encoder()) verifies the prototype hash against the currently-bundled prototypes, so a codec can never be silently used to produce embeddings in a different, incomparable space. CodecBundle.load() only deserializes — call matches_current_prototypes() if you touch transform / model directly.

mir.ml.bundle.prototype_hash(species, locus, n, replicate=0)[source]#

Stable 16-hex-char hash of the ordered prototype junction set (its identity).

replicate is part of the identity: two draws of the same size are different coordinate systems, so they must hash differently or the comparability guards would let them mix.

Parameters:

species (str)
locus (str)
n (int)
replicate (int)

Return type:

str

class mir.ml.bundle.CodecBundle(meta, transform, model)[source]#

Bases: object

A trained codec plus its embedding-space identity (PCA + prototype hash).

Parameters:

meta (dict)
transform (object)
model (object)

meta: dict#

transform: object#

model: object#

classmethod from_forward(encoder, species, locus, n_prototypes, gap_positions=(3, 4, -4, -3), kind='forward', replicate=0)[source]#

Parameters:

species (str)
locus (str)
n_prototypes (int)
gap_positions (tuple[int, ...])
kind (str)
replicate (int)

Return type:

CodecBundle

save(path)[source]#

Parameters:: path (str | Path)
Return type:: None

classmethod load(path)[source]#

Parameters:: path (str | Path)
Return type:: CodecBundle

matches_current_prototypes()[source]#

True iff the current bundled prototypes reproduce this codec’s embedding space.

Return type:: bool

forward_encoder(device=None, verify=True)[source]#

Reconstruct the ForwardEncoder, verifying prototype comparability.

Parameters:

device (str | None)
verify (bool)

`mir.ml.set_encoder`#

Learned repertoire track (Set-Transformer / DeepRC attention pooling).

Learned permutation-invariant repertoire encoder — the co-equal learned track (Theory §T.7/§T.8).

The fixed backbone (mir.repertoire) pools clonotypes by a size-weighted sum of fixed RFF features. This module replaces that with a trained set-pooling network (Deep Sets form, Prop. prop:sampinv): a per-clonotype MLP ψ_θ followed by pooling-by-multihead-attention (PMA, Lee et al. 2019 Set Transformer) whose learned seed vectors act as inducing points — i.e. trained metaclonotype / public-cluster detectors (the interaction signal of Prop. prop:interact, the alternative to the fixed second-moment Σ_S). It reads the same PCA-coordinate clonotype cloud as the fixed backbone (so the two are comparable) and is fit end-to-end to a label (age / CMV).

Depth-robustness is engineered in, not free (§T.7.5): clonotypes enter the attention weighted by their frequency g(a_σ) (a log-weight bias on the attention logits), and each epoch every sample is randomly subsampled so the network learns to read shallow repertoires — the RNA-seq regime.

Torch ([ml] extra). Shipped as a SetEncoderBundle (weights + PCA transform + prototype hash), refusing a prototype mismatch on load, exactly like mir.ml.bundle.CodecBundle.

class mir.ml.set_encoder.SetEncoder(*args, **kwargs)[source]#

Bases: Module

Per-clonotype MLP + frequency-weighted PMA pooling → fixed vector / prediction.

forward(Z, w) maps one sample’s clonotype cloud Z (n×p PCA coords) and weights w (n, summing to 1) to an out_dim output, permutation-invariantly. The n_seeds learned queries are the inducing points (metaclonotype detectors); frequency enters as an additive log w bias on the attention logits so abundant clones weigh more.

Parameters:

p (int)
d (int)
n_seeds (int)
heads (int)
out_dim (int)

forward(Z, w)[source]#

Parameters:

Z (torch.Tensor)
w (torch.Tensor)

Return type:

torch.Tensor

class mir.ml.set_encoder.SetEncoderModel(model, device, task)[source]#

Bases: object

A trained SetEncoder plus its device — turns a sample cloud into a prediction.

Parameters:

model (torch.nn.Module)
device (torch.device)
task (str)

model: torch.nn.Module#

device: torch.device#

task: str#

predict(clouds)#

Predict for a list of (Z, w) clonotype clouds (probabilities if classification).

Return type:: ndarray

mir.ml.set_encoder.train_set_encoder(clouds, y, *, task='regression', epochs=60, lr=0.001, n_seeds=16, d=128, val_frac=0.2, augment=True, seed=0, device=None, verbose=True)[source]#

Train the set encoder on (clouds, y) (clouds = list of (Z, w) per sample).

Parameters:

clouds – per-sample (Z: (n_i, p) float, w: (n_i,) float) from RepertoireSpace.sample_cloud.
y (ndarray) – per-sample target (float for regression, {0,1} for classification).
task (str) – "regression" (MSE) or "classification" (BCE-with-logits).
augment (bool) – subsample each sample every epoch to engineer depth-robustness (§T.7.5).
epochs (int)
lr (float)
n_seeds (int)
d (int)
val_frac (float)
seed (int)
device (str | None)
verbose (bool)

Returns:

(SetEncoderModel, metrics) — metrics carry the held-out score (Spearman / AUC).

Return type:

tuple[SetEncoderModel, dict]

class mir.ml.set_encoder.SetEncoderBundle(meta, transform, model)[source]#

Bases: object

A trained set encoder plus its embedding-space identity (PCA transform + prototype hash).

Parameters:

meta (dict)
transform (object)
model (object)

meta: dict#

transform: object#

model: object#

classmethod from_model(sem, space, task)[source]#

Parameters:

sem (SetEncoderModel)
task (str)

Return type:

SetEncoderBundle

save(path)[source]#

Parameters:: path (str | Path)
Return type:: None

classmethod load(path, *, verify=True)[source]#

Parameters:

path (str | Path)
verify (bool)

Return type:

SetEncoderBundle

`mir.ml.diffusion`#

Conditional diffusion generator (DDPM/DDIM + classifier-free guidance) over a compact descriptor/code space — the research half of the generative loop, complementing mir.generate’s linear Gaussian.

Diffusion generator over a compact embedding/descriptor space (Part 2, needs [ml]).

The research half of the generative loop (ROADMAP Phase 2/3), complementing mir.generate’s linear-Gaussian DescriptorDensity: a small conditional denoising diffusion model (DDPM, Ho et al. 2020) trained directly on a compact coordinate system — a RepertoireDescriptor vector, a codec’s PCA-compacted code (mir.ml.bundle), or any other (N, dim) matrix — that can express a genuinely non-linear generative manifold where the Gaussian baseline can only fit an ellipsoid. Sampling uses DDIM (Song et al. 2021) rather than the full ancestral chain: deterministic given a seed and needs far fewer steps (tens, not the training horizon T), the standard modern choice over vanilla DDPM sampling. Classifier-free guidance (Ho & Salimans 2022) — random label dropout during training, so one model samples both unconditionally and conditionally, with guidance_scale trading diversity for how strongly a sample matches its condition — lets one model serve as the digital-twin generator mir.generate’s sample/evolve API mirrors: draw a synthetic donor state unconditionally, or conditioned on a covariate (tumor type, HLA, batch), with the same call shape.

The coordinate space here is already compact (tens to a few hundred dimensions after PCA/descriptor reduction), so the epsilon-predictor is a small MLP (DiffusionMLP) — matching mir.ml’s existing compact-network style (SequenceDecoder), not a U-Net or transformer built for pixel/token grids that would be over-parameterized for this data.

DiffusionModel ships the same way CodecBundle does — meta + weights via save()/load() — so a trained generator carries its own provenance; pass a prototype_hash (mir.ml.bundle.prototype_hash()) in meta when the coordinate space came from a specific prototype set, and check it with matches_prototype_hash() before sampling.

Typical usage:

from mir.ml.diffusion import train_diffusion

model, metrics = train_diffusion(X, tumor_type, epochs=300)     # X: (n, dim) descriptor/code matrix
synth = model.sample(50, condition="SKCM", steps=50, guidance_scale=2.0)  # 50 synthetic SKCM-like
model.save("skcm_diffusion.pt")

mir.ml.diffusion.cosine_beta_schedule(T, s=0.008)[source]#

Cosine noise schedule (Nichol & Dhariwal 2021) — smoother SNR decay than linear DDPM betas, the modern default over the original linear schedule.

Parameters:

T (int)
s (float)

Return type:

torch.Tensor

class mir.ml.diffusion.DiffusionMLP(*args, **kwargs)[source]#

Bases: Module

Epsilon-predictor: noised vector + timestep (+ optional class) -> predicted noise.

Class index 0 is reserved as the “null” conditioning token (classifier-free guidance): real classes are indices 1..n_classes, so y=None (or the null token during training-time label dropout) samples the unconditional model.

Parameters:

dim (int) – Data dimensionality (the descriptor/code width).
time_dim (int) – Sinusoidal timestep embedding width.
hidden (int) – Hidden layer width.
n_classes (int) – Number of real classes; 0 disables conditioning entirely.
class_dim (int) – Class embedding width (unused if n_classes == 0).

forward(x_t, t, y=None)[source]#

Parameters:

x_t (torch.Tensor)
t (torch.Tensor)
y (torch.Tensor | None)

Return type:

torch.Tensor

class mir.ml.diffusion.DiffusionModel(model, betas, meta=<factory>)[source]#

Bases: object

A trained diffusion generator: noised-vector-space -> new samples via DDIM.

Variables:

model (torch.nn.Module) – The trained DiffusionMLP (eval mode).
betas (torch.Tensor) – (T,) noise schedule the model was trained with.
meta (dict) – Provenance dict — always carries dim, n_classes, T, and the per-dimension mu/sd the training data was standardized with (sample() un-standardizes with these before returning); carries classes (the original label values, index i -> class token i+1) when class-conditional; carries a prototype_hash iff the caller supplied one at train time.

Parameters:

model (torch.nn.Module)
betas (torch.Tensor)
meta (dict)

model: torch.nn.Module#

betas: torch.Tensor#

meta: dict#

matches_prototype_hash(current)[source]#

True iff meta has no prototype_hash recorded, or it equals current.

Parameters:: current (str)
Return type:: bool

save(path)[source]#

Serialize model + schedule + meta (mirrors mir.ml.bundle.CodecBundle.save()).

Parameters:: path (str | Path)
Return type:: None

classmethod load(path, device=None)[source]#

Rebuild a saved DiffusionModel (weights on CPU unless device given).

Parameters:

path (str | Path)
device (str | None)

Return type:

DiffusionModel

sample(n=1, *, condition=None, steps=50, guidance_scale=1.0, seed=0, device=None)#

Draw n synthetic vectors via deterministic DDIM (Song et al. 2021, eta=0).

Same call shape as mir.generate.DescriptorDensity.sample() (n, condition, seed) so the two generators are drop-in alternatives.

Parameters:

n (int) – Number of samples to draw.
condition – Class label to condition on (must be one of the labels train_diffusion saw); None samples unconditionally.
steps (int) – DDIM steps — far fewer than the training horizon T (state-of-art vs. the full ancestral chain); clipped to T if larger.
guidance_scale (float) – Classifier-free guidance weight. 1.0 = plain conditional/ unconditional prediction; >1.0 sharpens toward the condition (requires condition set and the model to be class-conditional).
seed (int) – RNG seed.
device (str | None) – Override device (default: wherever the model currently lives).

Returns:

(n, dim) array.

Raises:

ValueError – If condition is unknown, guidance is requested without condition, or the model isn’t class-conditional but condition was given.

Return type:

ndarray

mir.ml.diffusion.train_diffusion(X, condition=None, *, T=200, epochs=200, batch=64, lr=0.001, hidden=256, time_dim=64, class_dim=32, class_drop_prob=0.1, val_frac=0.1, seed=0, device=None, verbose=True, meta=None)[source]#

Train a conditional (or unconditional) diffusion model over X.

Parameters:

X (ndarray) – (n_samples, dim) data — a descriptor/code matrix.
condition – Optional length-n_samples class labels (any hashable value); None trains unconditionally.
T (int) – Diffusion horizon (training timesteps). Sampling can use far fewer via DDIM.
epochs (int) – Standard training hyperparameters.
batch (int) – Standard training hyperparameters.
lr (float) – Standard training hyperparameters.
hidden (int) – DiffusionMLP architecture knobs.
time_dim (int) – DiffusionMLP architecture knobs.
class_dim (int) – DiffusionMLP architecture knobs.
class_drop_prob (float) – Per-example probability of replacing the true label with the null token during training (classifier-free guidance) — ignored if condition is None.
val_frac (float) – Held-out fraction for the reported validation loss (best-epoch checkpointing).
seed (int) – RNG seed (data split, weight init, training stochasticity).
device (str | None) – Override device (default: mir.ml.train.pick_device()’s auto-pick).
verbose (bool) – Print progress.
meta (dict | None) – Extra provenance to carry on the returned model (e.g. {"prototype_hash": ...}).

Returns:

(model, metrics) — the fitted DiffusionModel (best-validation-loss weights) and a metrics dict (val_loss, n, n_classes).

Raises:

ValueError – If X has fewer than 4 rows (too few to hold out a validation split).

Return type:

tuple[DiffusionModel, dict]

Utilities#

`mir.aliases`#

Species / locus aliases.

Species and locus alias normalization for user-facing entry points.

Users, notebooks, and other tools spell species and loci many ways (beta, T_beta, hsa, Homo sapiens). Prototype resource files and germline lookups are keyed by the canonical IMGT locus (TRB) and canonical species (human / mouse), so normalize once at the boundary.

mir.aliases.normalize_species_alias(species)[source]#

Return the canonical species ('human' / 'mouse') for an alias.

Parameters:: species (str)
Return type:: str

mir.aliases.normalize_locus_alias(locus)[source]#

Return the canonical IMGT locus (e.g. 'TRB') for an alias.

Parameters:: locus (str)
Return type:: str

`mir.alleles`#

Allele normalisation with default-allele cascade.

Scalar V/D/J/C allele-notation helpers used by the germline distance lookup.

For frame-level allele operations use vdjtools.io.schema.strip_allele (a polars expression). These string helpers drive the per-clonotype resolution cascade (exact allele → *01 → bare gene → fallback) when gathering baked germline distances.

mir.alleles.strip_allele(gene_name)[source]#

Return the gene base without allele suffix (TRBV6-5*02 → TRBV6-5).

Parameters:: gene_name (str | None)
Return type:: str

mir.alleles.allele_with_default(gene_name, default_allele='01')[source]#

Append *01 when no allele is given; preserve an explicit allele.

TRBV6-5 → TRBV6-5*01; TRBV6-5*02 → TRBV6-5*02.

Parameters:

gene_name (str | None)
default_allele (str)

Return type:

str

mir.alleles.allele_to_major(gene_name)[source]#

Normalize to major-allele form *01 (TRBV6-5*02 → TRBV6-5*01).

Parameters:: gene_name (str | None)
Return type:: str

Command-line interface (`mir`)#

pip install mirpy-lib installs a mir console script (also python -m mir.cli) with two commands, one per embedding scale:

Both commands drop non-coding clonotypes (stop codon / legacy out-of-frame junction_aa) before embedding by default via vdjtools.preprocess.filter_functional — pass --no-filter-functional to disable.

mir embed clonotypes INPUT: One repertoire’s clonotype table → a per-clonotype TCREMP embedding table (columns e0…). Flags: --species, --locus (inferred when the file has one locus), --n-prototypes, --mode {vjcdr3,cdr123}, --replicate R (prototype draw), --pca K (compact the table), --filter-functional/--no-filter-functional, --threads, -o (.tsv / .parquet; default stdout TSV).
mir embed repertoires INPUT...: A dataset of clonotype tables → one repertoire vector Φ(S) per sample per chain on one shared basis (columns phi0…; sample id = filename stem). Flags: --species, --locus (restrict), --n-prototypes, --replicate R, --weight {log2p1,duplicate_count,distinct,log1p,anscombe} (default log2p1), --blocks mean,diversity[,second], --n-rff, --n-rff-second, --n-components, --mmd OUT (also write the per-chain pairwise unbiased-MMD matrix), --filter-functional/--no-filter-functional, --threads, --seed, -o.

mir command-line interface — turn receptor tables into embeddings.

Two commands cover the two scales mirpy embeds at:

mir embed clonotypes SAMPLE — one repertoire’s clonotype table → a per-clonotype TCREMP embedding table (e0…), the input to clustering / ML.
mir embed repertoires SAMPLE… — a dataset of clonotype tables → one repertoire vector Φ(S) (phi0…) per sample, per chain, on one shared basis (so the rows are mutually comparable / MMD-able).

Inputs are any format vdjtools.io.read sniffs (AIRR TSV, vdjtools, MiXCR, immunoSEQ, parquet, …). Output is TSV (default / .tsv) or Parquet (.parquet — recommended for the wide raw embedding); -o - (or no -o) writes TSV to stdout.

Run mir embed clonotypes -h / mir embed repertoires -h for the full flag list.

mir.cli.build_parser()[source]#

Return the mir argument parser (embed clonotypes / embed repertoires).

Return type:: ArgumentParser

mir.cli.main(argv=None)[source]#

Parse argv (default sys.argv[1:]) and run the requested mir subcommand.

Parameters:: argv (list[str] | None)
Return type:: None

API reference#

mir#

mir#

Clonotype embedding (mir.embedding)#

mir.embedding.prototypes#

mir.embedding.tcremp#

mir.embedding.pca#

mir.embedding.presets#

Distances (mir.distances)#

mir.distances.junction#

mir.distances.germline#

Density and background subtraction (mir.density)#

Sample-level (repertoire) embedding (mir.repertoire)#

Explainable readouts (mir.explain)#

Cohort / digital donor (mir.cohort)#

Exposure trajectory (mir.track)#

Generative loop (mir.generate)#

Digital twin (mir.twin)#

Benchmark harness (mir.bench)#

mir.bench.vdjdb#

mir.bench.metrics#

mir.bench.theory#

mir.bench.eval#

Neural codecs and learned encoders (mir.ml)#

mir.ml.tokenize#

mir.ml.encoder#

mir.ml.decoder#

mir.ml.train#

mir.ml.codec#

mir.ml.bundle#

mir.ml.set_encoder#

mir.ml.diffusion#

Utilities#

mir.aliases#

mir.alleles#

Command-line interface (mir)#

`mir`#

Clonotype embedding (`mir.embedding`)#

`mir.embedding.prototypes`#

`mir.embedding.tcremp`#

`mir.embedding.pca`#

`mir.embedding.presets`#

Distances (`mir.distances`)#

`mir.distances.junction`#

`mir.distances.germline`#

Density and background subtraction (`mir.density`)#

Sample-level (repertoire) embedding (`mir.repertoire`)#

Explainable readouts (`mir.explain`)#

Cohort / digital donor (`mir.cohort`)#

Exposure trajectory (`mir.track`)#

Generative loop (`mir.generate`)#

Digital twin (`mir.twin`)#

Benchmark harness (`mir.bench`)#

`mir.bench.vdjdb`#

`mir.bench.metrics`#

`mir.bench.theory`#

`mir.bench.eval`#

Neural codecs and learned encoders (`mir.ml`)#

`mir.ml.tokenize`#

`mir.ml.encoder`#

`mir.ml.decoder`#

`mir.ml.train`#

`mir.ml.codec`#

`mir.ml.bundle`#

`mir.ml.set_encoder`#

`mir.ml.diffusion`#

`mir.aliases`#

`mir.alleles`#

Command-line interface (`mir`)#