Metadata-Version: 2.5
Name: mirpy-lib
Version: 4.1.0
Summary: ML-oriented embeddings for immune receptor repertoires (TCR/BCR)
Project-URL: Homepage, https://github.com/antigenomics/mirpy
Project-URL: Repository, https://github.com/antigenomics/mirpy
Project-URL: Documentation, https://antigenomics.github.io/mirpy
Project-URL: Issues, https://github.com/antigenomics/mirpy/issues
Author-email: ISALGO lab <mikhail.shugay@gmail.com>
License: GPL-3.0-or-later
License-File: LICENSE
Keywords: AIRR,BCR,TCR,VDJdb,embedding,immunosequencing,prototypes
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.11
Requires-Dist: numpy
Requires-Dist: polars
Requires-Dist: scikit-learn
Requires-Dist: scipy
Requires-Dist: seqtree>=1.0.0
Requires-Dist: vdjtools>=4.1.0
Provides-Extra: ann
Requires-Dist: numba>=0.60; extra == 'ann'
Requires-Dist: pynndescent; extra == 'ann'
Provides-Extra: annotate
Requires-Dist: vdjmatch; extra == 'annotate'
Provides-Extra: bench
Requires-Dist: huggingface-hub; extra == 'bench'
Requires-Dist: kneed; extra == 'bench'
Requires-Dist: lifelines; extra == 'bench'
Requires-Dist: matplotlib; extra == 'bench'
Requires-Dist: seaborn; extra == 'bench'
Provides-Extra: build
Requires-Dist: biopython; extra == 'build'
Requires-Dist: vdjtools[model]>=4.1.0; extra == 'build'
Provides-Extra: dev
Requires-Dist: psutil>=5; extra == 'dev'
Requires-Dist: pytest-cov; extra == 'dev'
Requires-Dist: pytest>=8; extra == 'dev'
Provides-Extra: docs
Requires-Dist: pydata-sphinx-theme<1,>=0.15; extra == 'docs'
Requires-Dist: sphinx<9,>=7.2; extra == 'docs'
Provides-Extra: examples
Requires-Dist: marimo; extra == 'examples'
Requires-Dist: matplotlib; extra == 'examples'
Requires-Dist: umap-learn; extra == 'examples'
Provides-Extra: ml
Requires-Dist: torch; extra == 'ml'
Description-Content-Type: text/markdown

<p align="center">
  <picture>
    <source media="(prefers-color-scheme: dark)" srcset="assets/mirpy_dark.svg">
    <source media="(prefers-color-scheme: light)" srcset="assets/mirpy_light.svg">
    <!-- Absolute PNG fallback: PyPI strips <picture>/<source> and cannot render a relative or
         raw-served SVG, so the logo must be an absolute-URL raster here. GitHub uses the SVG sources. -->
    <img alt="mirpy" src="https://raw.githubusercontent.com/antigenomics/mirpy/master/assets/mirpy_light.png" width="360">
  </picture>
</p>

<h1 align="center">mirpy — ML embeddings for immune repertoires</h1>

[![PyPI](https://img.shields.io/pypi/v/mirpy-lib.svg)](https://pypi.org/project/mirpy-lib/)
[![Python](https://img.shields.io/pypi/pyversions/mirpy-lib.svg)](https://pypi.org/project/mirpy-lib/)
[![License](https://img.shields.io/badge/license-GPLv3-green)](LICENSE)
[![Docs](https://img.shields.io/badge/docs-docs.isalgo.dev-blue)](https://docs.isalgo.dev/mirpy/)

**mirpy v3** turns T-/B-cell receptor sequences into fixed-length numeric vectors you can
cluster, visualize, and feed to ML models. It implements **TCREMP** — embedding each receptor
by its alignment distances to a fixed set of *prototype* sequences — so that Euclidean distance
in embedding space approximates pairwise alignment distance (Theory T1).

> v3 is a slim, embedding-focused rewrite. The classical repertoire toolkit (parsing, overlap,
> diversity, TCRnet, GLIPH, …) lives on the **`legacy-v2`** branch (`mirpy-lib` 2.x) and in the
> sibling tools [`vdjtools`](https://github.com/antigenomics/vdjtools) /
> [`vdjmatch`](https://github.com/antigenomics/vdjmatch).

## Install

```bash
pip install mirpy-lib            # core: numpy, polars, scipy, scikit-learn, seqtree, vdjtools
pip install "mirpy-lib[bench]"   # + benchmark / theory experiments
```

Pure-Python wheel; the heavy lifting (alignment, Pgen, sampling) is reused from
[`seqtree`](https://github.com/antigenomics/seqtree) and `vdjtools`.

## Where to start

| You want to | Go to |
|---|---|
| Turn receptors into vectors you can cluster or classify | [Quick start](#quick-start) · [User guide](https://docs.isalgo.dev/mirpy/usage.html) |
| Pick prototype counts and PCA dimensions | [Recommended presets](#recommended-presets) |
| Find enriched / antigen-driven neighbourhoods | [`mir.density`](#background-subtraction-and-clustering-mirdensity) |
| Compare whole repertoires, not clonotypes | [`mir.repertoire`](#sample-level-repertoire-embedding-mirrepertoire) |
| One fixed feature vector per sample, for a classifier | [Repertoire signatures](#repertoire-signatures-extended) · [Signature](https://docs.isalgo.dev/mirpy/signature.html) |
| Name which part of a vector carries a signal | [`mir.explain`](https://docs.isalgo.dev/mirpy/channels.html) |
| Worked examples as notebooks | [Notebook gallery](https://docs.isalgo.dev/mirpy/notebooks.html) |

## Quick start

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

model = TCREmp.from_defaults("human", "TRB", n_prototypes=3000)   # mode="vjcdr3" | "cdr123"
df = pl.DataFrame({
    "v_call":      ["TRBV10-3*01", "TRBV20-1*01"],
    "j_call":      ["TRBJ2-7*01",  "TRBJ1-2*01"],
    "junction_aa": ["CASSIRSSYEQYF", "CSARVSGYYGYTF"],
})
X = model.embed(df)          # (2, 9000) float32 — 3 distances × 3000 prototypes
```

Downstream (cluster antigen-specific TCRs — needs `[bench]` for the kneedle `eps`, and a real
set of clonotypes rather than the two above):

```python
from mir.embedding.pca import pca_denoise
from mir.bench.metrics import cluster, cluster_metrics
labels = cluster(pca_denoise(model.embed(vdjdb_df), n_components=50))
```

Paired chains concatenate per-chain embeddings via `PairedTCREmp`. Input/output are AIRR polars
frames keyed by `vdjtools.io.schema` column names.

## Command line

`pip install mirpy-lib` also installs a `mir` command for the two embedding scales — no Python
needed. Inputs are any format `vdjtools.io` reads (AIRR TSV, vdjtools, MiXCR, immunoSEQ, parquet).

```bash
# one repertoire  ->  per-clonotype embedding table (e0…), the input to clustering / ML
mir embed clonotypes sample.tsv --pca 50 -o clonotypes.parquet

# a dataset of repertoires  ->  one fingerprint Φ(S) per sample, per chain (phi0…), on one
# shared basis so the rows are mutually comparable; --mmd also writes the pairwise MMD matrix
mir embed repertoires cohort/*.tsv.gz -o phi.tsv --mmd mmd.tsv

# the portable signature  ->  one fixed, named, standardised feature vector per sample
# One tool per half: mirpy emits the geometry, vdjtools the statistics. Join on sample_id.
mir corpus --corpus synthetic-blood -o rsig_synthetic-blood.npz  # fit one; no cohort needed
mir signature --corpus synthetic-blood cohort/*.tsv.gz -o rsig.parquet  # geometry half
mir signature --corpus synthetic-blood --components 32 --describe       # exactly what you get
```

`mir embed clonotypes -h` / `mir embed repertoires -h` list every flag (species, locus,
prototype count, weight, Φ blocks, …). Sample id defaults to the filename stem; the locus is
inferred per file (or restrict with `--locus`, which takes aliases — `beta`, `T-alpha` — and errors
on anything it can't resolve). An MMD matrix is per chain, so across several loci `--mmd mmd.tsv`
writes `mmd.TRB.tsv`, `mmd.TRA.tsv`, …; with one locus the name is used as given.

Both `embed` commands drop non-coding clonotypes (stop codon / legacy out-of-frame markers in
`junction_aa`) before embedding, and there is **no flag to turn it off**: a stop codon is in
seqtree's alphabet, so an unfiltered frame does not crash — it embeds to a finite, meaningless
distance and contaminates the geometry silently. `--no-filter-functional` is refused with a
pointer to `vdjtools filter --nonproductive`, which is what to use when the non-productive
fraction is the thing you want.

## Recommended presets

`TCREmp.from_defaults(species, locus)` uses the per-chain preset when `n_prototypes` is
omitted. Values are data-driven from the bundled prototypes (prototype geometry saturates by
these counts; PC columns are the PCA dims retaining ~95% / ~99% variance):

| chain | n_prototypes | PCs (95%, clustering) | PCs (99%, reconstruction) |
|---|--:|--:|--:|
| human TRA | 2000 | 65 | 220 |
| human TRB | 2000 | 65 | 260 |
| human TRG | 1000 | 25 | 100 |
| human TRD | 2000 | 65 | 280 |
| human IGH | 2000 | 65 | 300 |
| human IGK | 1000 | 20 | 65 |
| human IGL | 1000 | 20 | 65 |
| mouse TRA | 2000 | 50 | 150 |
| mouse TRB | 2000 | 55 | 225 |

Use **95%** PCs for clustering/visualization (the paper's regime); use **99%** PCs when
*reconstructing* sequences with the neural inverse codec (diverse chains like IGH/TRD/TRA lose
too much sequence detail at 95%). Programmatically: `from mir.embedding import get_preset`.

```python
from mir.embedding import get_preset
from mir.embedding.pca import pca_denoise
p = get_preset("human", "IGH")
Xc = pca_denoise(X, n_components=p.n_components)          # clustering
Xr = pca_denoise(X, n_components=p.n_components_recon)    # codec reconstruction
```

## Prototypes — which receptors, and how much do they matter?

Every embedding is *distances to prototypes*, so the prototype set **is** the coordinate system.
mirpy ships one per chain, and you get them without downloading anything:

| | |
|---|---|
| What | **10 000 real receptors** per chain — a uniform random sample (fixed `seed=42`) of unique, productive, germline-resolvable clonotypes from arda-annotated real repertoires |
| Why real | Model-generated junctions have degenerate lengths and embed measurably worse (negative self-prototype distance correlation); real repertoires give a tight, well-behaved manifold |
| The default | `replicate=0` — the first `n` rows. This is *the* set: every preset, bundled codec, and published number uses it. Don't change it unless you're deliberately testing sensitivity |
| Chains | human TRA/TRB/TRG/TRD/IGH/IGK/IGL, mouse TRA/TRB (`list_available_prototypes()`) |

**Is my result an artefact of which prototypes I drew?** Take a replicate. The file order is itself a
uniform shuffle, so each disjoint block of `n` rows is an independent draw from the same pool —
`n_replicates()` of them, **10 at `n=1000`**, 5 at `n=2000`:

```python
from mir.embedding.prototypes import n_replicates
from mir.embedding.tcremp import TCREmp

scores = [my_metric(TCREmp.from_defaults("human", "TRB", 1000, replicate=r).embed(df))
          for r in range(n_replicates("human", "TRB", 1000))]   # 10 draws; spread = sensitivity
```

Same from the shell: `mir embed clonotypes sample.tsv --n-prototypes 1000 --replicate 3`.

**How much *does* it matter?** Usually very little, and you can check for yourself —
`bench.theory.prototype_source_correlation(queries, protos_a, protos_b)` correlates the pairwise
junction-distance geometry under two prototype sets. Two independent draws, 400 held-out human-TRB
queries:

| prototypes `n` | 100 | 250 | 500 | 1000 (default) | 2000 |
|---|--:|--:|--:|--:|--:|
| R between two draws | 0.922 | 0.971 | 0.990 | 0.993 | **0.997** |

So the geometry is essentially draw-independent from `n≈500` up: at the default counts, *which*
prototypes you drew is not what your result rests on. Below `n≈250` it starts to be.

> Each replicate is a **different coordinate system**. Distances *within* one are comparable;
> distances *across* two are not. The prototype hash covers the replicate index, so codecs,
> `RepertoireSpace` and `DonorCohort` all refuse to mix them — compare summary statistics across
> replicates (AUC, F1, cluster counts), never raw embeddings. Sweeping `n_prototypes` instead is a
> *nested* comparison (draw `r=0` at `n=500` is a prefix of `n=1000`), which answers "how many do I
> need", not "does it matter which".

The regenerate command is `src/mir/resources/prototypes/generate_prototypes.py` (needs `[build]`
and `ARDA_HOME`); the shipped TSVs are the versioned reference and need no rebuild. Per-artifact
provenance travels in `manifest.json` beside each one.

## What's inside

Grouped by what you are doing. Full API in
[`skills/mirpy/SKILL.md`](skills/mirpy/SKILL.md) and the
[API reference](https://docs.isalgo.dev/mirpy/api.html).

**Receptors to vectors**

| Module | What it does |
|---|---|
| `mir.embedding` | `TCREmp` / `PairedTCREmp`, PCA denoising, the per-chain presets |
| `mir.distances` | junction distance via `seqtree.gapblock`, plus baked germline distances |

**Repertoires to vectors**

| Module | What it does |
|---|---|
| `mir.repertoire` | one vector per sample: RFF kernel mean, Hill diversity, second moment; MMD, motif witness, sub-probability measures |
| `mir.signature` | the **portable signature** — fixed, named, already-standardised columns you can hand to a collaborator |
| `mir.explain` | which named channel carries a signal, and which clonotypes drive it |
| `mir.density` | continuous TCRNET / ALICE: neighbourhood enrichment and noise filtering |

**Cohorts, time and generation**

| Module | What it does |
|---|---|
| `mir.cohort` | the **digital donor** — multi-chain donor embeddings, residualisation, incidence biomarkers |
| `mir.track` | **exposure trajectory** — a latent progression axis disentangled from a known covariate |
| `mir.generate` | the **generative loop**: sample new synthetic donor states, or evolve one along a coordinate |
| `mir.twin` | the **digital twin**: perturb or resample one donor's state through a generator |

**Supporting**

| Module | What it does |
|---|---|
| `mir.cli` | the `mir` console script — `embed clonotypes`, `embed repertoires`, `signature`, `presets` |
| `mir.bench` | VDJdb loader, clustering metrics, theory experiments, cohort scorers |
| `mir.ml` | neural codecs, learned set encoders, diffusion generator — experimental, `[ml]` extra |
| `mir.aliases`, `mir.alleles` | species / locus aliases and allele-name normalisation |

## Background subtraction and clustering (`mir.density`)

TCRNET/ALICE find antigen-driven convergent clusters by *neighbour enrichment*. `mir.density`
does the same test with neighbour-counting in the **embedding space** instead of on a sequence
graph (Theory T6): the enrichment `E(z) = f_obs(z)/f_gen(z)` is estimated by an adaptive-bandwidth
**balloon** estimator with a per-clonotype Poisson/binomial significance test and BH q-values —
no graph, and it scales to whole repertoires.

```python
from mir.density import fit_density_space, neighbor_enrichment, enriched_mask, denoise_and_cluster
from mir.embedding.tcremp import TCREmp

model = TCREmp.from_defaults("human", "TRB", n_prototypes=1000)
# background = a control repertoire (TCRNET) or generate_background(...) (ALICE, P_gen)
space, obs_emb, bg_emb = fit_density_space(model, obs_df, control_df, n_components=20, space="full")
res  = neighbor_enrichment(obs_emb, bg_emb, test="binomial")   # balloon + water-level calibration
hits = obs_df.filter(enriched_mask(res, alpha=0.05))            # background-subtracted clones
labels, mask = denoise_and_cluster(obs_emb, res)               # noise-filter + DBSCAN the hits
```

Use a **biological control** as the background when you have one (e.g. pre- vs post-vaccination,
patient vs healthy) — differential enrichment cancels generic public convergence and isolates the
antigen-specific response. No control of your own? Pooled healthy-donor repertoires are one fetch
away (HF [`isalgo/airr_control`](https://huggingface.co/datasets/isalgo/airr_control), read with
`vdjtools.io.read`; the dataset card carries the caveats).
Failing that, `generate_background(locus, n)` samples the vdjtools P_gen model (the ALICE regime);
the "water level" of a naive repertoire is handled by the empirical-null calibration. Pass
`source="arda"` there when your data is arda-annotated (same allele namespace as the prototypes),
and `species="mouse"` for mouse — both need a vdjtools shipping the
bundled `arda` model set. The density benchmarks (YFV, ankylosing-spondylitis B27, TCRNET)
live in the companion [`2026-mirpy-analysis`](https://github.com/antigenomics) repo.

The default backend is `"kdtree"` (exact scipy cKDTree, all cores). At whole-repertoire scale pass
`backend="ann"` (pynndescent, ~30× faster past ~10⁵ clones; `pip install "mirpy-lib[ann]"`). Only
the **observed** side is approximate there — recall < 1 undercounts the observed ball, which biases
enrichment *down* (conservative). The **background** occupancy is always exact, because
undercounting it would shrink the expected count and inflate fold and significance, which is the one
direction an enrichment test must never err in.

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

One fixed vector `Φ(S)` per **repertoire** — an order-invariant multiset of clonotypes with clone
sizes — depth-robust into the low-coverage bulk-RNA-seq regime (Theory §T.7). `Φ(S)` sketches the
empirical measure `ρ_S = Σ_σ w_σ δ_{φ(σ)}` in three blocks: an RFF **kernel mean** (depth-robust,
codebook-free — no `K`, no clustering), a coverage-standardized **Hill diversity** profile, and a
**second-moment** Fisher vector carrying clonotype co-occurrence (HLA-linked public structure).
Repertoire distance is the **MMD** `‖Φ₁(S) − Φ₁(S')‖`.

The per-clonotype weights `w_σ = g(a_σ)/Σ_τ g(a_τ)` come from a clone-size transform `g` (`weight=`
on `sample_embedding`/`fit_repertoire_space`/`mir embed repertoires --weight`): `"log2p1"` —
`g=log2(1+a)` — is the **default**, concave so one hyperexpanded clone can't dominate;
`"duplicate_count"` weights linearly by clone size (`g=a`); `"distinct"` ignores size entirely
(`g≡1`, presence only). `"log1p"` (natural log) and `"anscombe"` remain available.

```python
from mir.repertoire import fit_repertoire_space, sample_embedding, mmd_matrix, class_witness
from mir.embedding.tcremp import TCREmp
import polars as pl

model  = TCREmp.from_defaults("human", "TRB", n_prototypes=1000)
space  = fit_repertoire_space(model, pl.concat(samples))   # ONE basis for the whole cohort
embs   = [sample_embedding(space, s) for s in samples]     # Φ(S): mean ‖ diversity ‖ second moment
D      = mmd_matrix(embs, unbiased=True)                    # pairwise repertoire distance (unbiased MMD²)
motifs = class_witness(space, pos_samples, neg_samples, candidates)   # public clones separating two groups
```

**Comparability invariant** (as with the codecs / density): 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 and `RepertoireSpace` refuses a prototype-hash mismatch.

Use the **unbiased** MMD (`unbiased=True`) whenever samples differ in depth/diversity — the biased
V-statistic's `1/n_eff` self-term otherwise inflates low-diversity samples and fakes a signal. When a
nuisance batch is present, compare *within-batch* contrasts (residualize `Φ` on the batch indicator):
a batch offset is first-order and cancels, while a batch-orthogonal signal (e.g. HLA) survives. The
empirical rule of thumb — **diversity for how-even, the embedding for which-clones**: clone-size
phenotypes (age, CMV) are a diversity summary's turf, while clonotype identity (HLA — strongest in
TRA and class II) lives in the second moment / witness. A learned co-equal set encoder
(Set-Transformer / DeepRC) is in `mir.ml.set_encoder` (`[ml]` extra). Recorded results and theory
(T7) live in the companion [`2026-mirpy-analysis`](https://github.com/antigenomics) repo
(`benchmarks/{BENCHMARKS,THEORY}.md`) alongside the benchmark scripts.

### Sub-probability embeddings: the deficient measure

`Φ(S)` above is the kernel mean embedding of a **probability** measure — the weights sum to 1, so
every sample asserts one full unit of confidence. At RNA-seq depth that premise fails. Measured: the
**median tissue TRB sample holds 21 unique clonotypes** (blood TRB
254, 1st percentile 1), so `w_σ = a_σ/Σa` is `1/n` for a *technical draw size*, not a clonal
frequency — the true frequencies live at 1e-5…1e-8, and one singleton's weight spans **21,454×**
across blood TRB purely from sample size. Worse, normalising to 1 *forces* a 5-clonotype tumour to
assert full confidence, so it lands somewhere arbitrary on the unit sphere instead of where it
belongs, and callers respond with a minimum-clonotype floor which in tumour deletes the **immune
desert** — the phenotype of interest. (A floor once cut 7,179 labelled donors to 2,129.)

Let the measure be **sub-probability** instead:

```python
from mir.repertoire import contrast_embedding, missing_mass, naive_reference, sample_embedding

emb = sample_embedding(space, sample, missing_mass="chao")   # or "turing"; "none" = old behaviour
emb.mass                                    # retained mass 1 − M₀ ∈ [0, 1]
ref = naive_reference(space)                 # kernel mean of 20k naive V(D)J recombinations (~8 s)
psi = contrast_embedding(emb, ref)           # Ψ = mass·(Φ − naive): signed, magnitude-carrying
```

* **`missing_mass`** estimates the mass `M₀` of the clonotypes that were never drawn — Good–Turing
  (`f₁/N`) or **bias-corrected Chao1** (`S_u/(N+S_u)`, `S_u = f₁(f₁−1)/(2(f₂+1))`; never the
  classical `f₁²/2f₂`, undefined when no clone was seen exactly twice, which is common here).
  `missing_mass=` only sets `.mass`; the blocks are untouched, and the `"none"` default is
  bit-identical to before.
* **Not** a negative measure. `Φ`'s value is that `‖Φ_P − Φ_Q‖` *is* the MMD and that a convex
  combination of two `Φ`'s is the `Φ` of a real pooled repertoire (what makes `mir.twin` and
  trajectory interpolation mean anything); a measure allowed to go negative on a set is not a
  probability measure and loses both. A sub-probability measure costs neither.
* **`naive_reference`** gives the unseen block a principled location — the germline recombination
  model (`vdjtools.model.generate`), not the corpus centroid. This is the load-bearing choice:
  shrinking toward the centroid is James–Stein toward the mean and it measurably **hurt** (it piles
  shallow samples into a dense ball that is itself depth-correlated), while the germline draw 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.
* **`contrast_embedding`** is where legitimate negativity lives: a signed *difference of two
  probability measures*, negative wherever the sample is depleted relative to unselected
  recombination, still an ordinary RKHS element with `‖Ψ_S‖ = MMD(S, naive)`. Magnitude then reads
  as **confidence × deviation-from-naive**: an immune desert has `M₀ → 1` and lands at the
  **origin**, the right place for "no infiltrate detected", and a shallow blood sample says so by
  its norm instead of being dropped. Give the caller `mass` to weight with; don't add a floor.

> WARNING: **Scale a magnitude-carrying block with one global scalar, never per column.** Per-column
> standardisation forces every coordinate to unit variance across samples, so a matrix where half the
> rows sit at the origin comes out looking exactly like one where none do — it deletes the deficiency
> it was built to preserve. Use `ChannelBuilder.add(..., preserve_magnitude=True)`, which applies one
> pooled RMS per channel and fills holes with `0`. (`stack_embeddings` warns if it is handed
> deficient-mass embeddings, since `Φ.vector` does not carry the mass.)

### Functional diversity, compartments, and depth

Four more consequences of the same measure algebra, all of them cheap once `Φ` is a kernel mean:

```python
from mir.repertoire import (band_embeddings, depth_threshold, mixture_weights,
                            rao_q, rarefy_embedding, sample_statistics)

rao_q(emb)                          # 1 − ‖Φ₁‖² — Rao's quadratic entropy, exactly
depth_threshold(embs).kappa         # the size below which Φ is mostly sampling noise
bands = band_embeddings(space, sample)              # singleton / expanded / top (or IGH isotypes)
mixture_weights(emb, bands)["weights"]              # π per compartment, by NNLS
rarefy_embedding(space, sample, depth=20_000).v_rep # matched-depth Φ + its replicate noise
```

* **`rao_q`** — Rao's quadratic entropy is `1 − ‖Φ₁‖²` *exactly* (verified against an explicit
  Gram to ~1e-16), so the **norm** of the kernel mean is a diversity statistic and no Gram matrix
  is needed.
  It 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, while Rao's Q weights each pair by how different the receptors are. Measured: this one
  scalar recovers R² 0.74–0.85 of classical diversity, and embedding derivatives reach R² 0.974–0.994
  for Shannon. Valid only on the **uncentred** Φ₁ — centring preserves differences (MMD) but not norms.
* **`depth_threshold`** — the damage depth does to a kernel mean is not bias but **variance ∝ 1/n**.
  Regressing `‖Φ_S − Φ̄‖²` on `1/n` splits the spread into between-sample signal `τ²` and sampling
  noise `σ²`, so **κ = σ²/τ²** is the size at which they are equal. Measured κ ≈ 40–70 clonotypes
  across four independent views, with 23–69% of samples below it. Report κ for *your* cohort instead
  of importing a cutoff — and note the library still applies no floor.
* **`band_frames` / `band_embeddings` / `mixture_weights`** — Φ₁ is an *average*, which is the right
  operation for a population mean and the wrong one for a minority signal: writing the repertoire as
  `(1−π)ρ_naive + π ρ_expanded` shows a compartment-confined effect reaching Φ₁ attenuated to `πΔ`
  while the naive compartment supplies the noise. Bands (`singleton` / `expanded` / `top 1%`, or IGH
  isotypes from `c_call`) are embedded through the **same frozen space** — never refit, or band-to-band
  distances stop meaning anything — and bands under 5 clonotypes are recorded absent (`None`), not
  upsampled. Because mixture linearity `Φ(S) = Σ π_c Φ(c)` is exact, a non-negative least
  squares recovers each compartment's share (exact in exact arithmetic; float32 embeddings put the
  realised residual near 1e-5). Measured on IGH isotypes: class-switched **IgG carries
  π = 0.070** of Φ₁(IGH) (IgM 0.230, IgA 0.176, 0.520 uncalled) — the dilution factor measured rather
  than assumed, and a power calculation before you spend compute: a subset with π ≈ 0.001 is not
  detectable by any aggregate distance, so use the per-clonotype witness. Straight about the negative:
  on survival endpoints the bands did **not** beat a diversity reference (0/22 pre-registered cells),
  though `singleton` reproduced a clinical-only score and `expanded` reproduced the whole-repertoire
  score exactly as the mixture argument predicts, and banding did win on kNN entropy in tissue IGH
  (0.3875 vs 0.4686).
* **`rarefy_embedding`** — averaging Φ over multinomial subsamples is itself a kernel mean (of the
  mixture over subsamples), so it is the **only** depth correction that keeps MMD, Rao and mixture
  linearity exactly (~1e-15); an orthogonal projection breaks the norm identity and a per-coordinate
  rescale breaks both. It also hands you a free per-sample noise estimate, because
  `Rao(Φ̄) = mean_r Rao(Φ_r) + v_rep` holds **exactly** — the excess diversity of the average *is* the
  replicate variance. Not a default: rarefying a cohort to its shallowest useful depth throws away the
  deep samples' advantage.

Also `sample_statistics` / `cohort_statistics` (the sampling fingerprint: f₁/f₂/f₃₊, singleton
fraction, missing mass, top-clone fraction, library size — the `stats=` input below, and candidate
biology in their own right) and `mir.cohort.depth_report`, which regresses the leading PCs on that
fingerprint. Read it beside `explained_variance`: with the deficient measure, 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".

> WARNING: **`mir.cohort.residualize(..., shrink=True)`** for batch offsets fitted from few samples in many
> dimensions. Plain per-group centring can make the batch *easier* to read, measured out-of-sample:
> batch-identity AUC 0.863 raw → **0.985** after centring, 0.978 after ComBat. The mechanism is
> estimation error — `‖µ̂ − µ‖ ≈ √(σ²d/n)` was ≈16 against a true offset of 7–24, so subtracting
> it injects a batch-constant vector as large as the one it removes, invisible in-sample by
> construction.
> Positive-part James–Stein shrinkage recovered most of it: 0.985 → **0.889**.

The matching evaluation criterion is **recoverability, not competition**:
`mir.bench.recovery_report(X, stats, groups)` runs a grouped-CV ridge from the embedding's PCs back
to each basic repertoire statistic (richness, Shannon, top-clone fraction, singleton fraction, Chao
unseen fraction, library size) and reports R² — high means the statistic is *carried inside* the
embedding and nothing has to be bolted on beside it. Renormalising to mass 1 deletes the magnitude,
so coverage and richness are unrecoverable from `Φ` **by construction**; the deficient measure should
win that question as a design consequence. It sits beside `mir.cohort.missingness_report` as the
other "is this object honest" check.

## Repertoire signatures (extended)

An **optional** layer on top of everything above, for when the deliverable is a *table a
collaborator can join* rather than an embedding. `Φ(S)` is a fingerprint but not a portable one —
its basis is fitted on *your* cohort, so two labs get incomparable vectors. The **signature** fixes
the basis and the scale: a fixed-width, name-addressed, already-standardised vector that anyone who
`pip install mirpy-lib` can compute from their own AIRR files and drop into PCA, logistic
regression, boosting or an MLP with no scaler of their own.

Command line and library both, and they emit the same columns:

```bash
mir signature --corpus synthetic-blood cohort/*.tsv.gz -o rsig.parquet
mir corpus    --corpus synthetic-tissue -o rsig_synthetic-tissue.npz   # or build your own
```

```python
from mir.signature import Corpus, rsig_cohort

corpus = Corpus.load("rsig_synthetic-blood.npz")   # REQUIRED -- there is no default
F = rsig_cohort(samples, corpus, n_jobs=0)         # one row per sample
```

**Four corpora ship, all synthetic.** `synthetic-blood` and `synthetic-tissue` draw each repertoire
as a naive/memory mixture across three quantile ladders measured per locus on that compartment --
clonotype richness, reads per expanded clone, and the singleton fraction that stands in for the naive
share -- so the corpus spans the depth and clone-size range real samples have (blood TRB richness 74
to 3,162 clonotypes, n = 34,365 reference samples). `naive` and `memory` remain as the pure-regime
references. Every receptor comes from vdjtools' bundled recombination models, so no cohort is needed
to build or use one.

Two halves, joined on `sample_id` and namespaced so they never collide: `vsig` (statistics of the
clone-size vector, from [vdjtools](https://github.com/antigenomics/vdjtools)) and `rsig` (geometry —
functionals of `Φ`). A column is `<sig>:<block>:<locus>:<feature>`. What comes out per locus is
`<sig>:pc:<locus>:PCnn` plus **channels**, which are carried in their own units and never rotated:
`rsig:div:*:rao` (sequence-aware diversity) and `rsig:qc:-:winsor_frac` (how much of the row the
corpus's bounds clamped). See [**Channels**](https://docs.isalgo.dev/mirpy/channels.html).

**A corpus is required, and it fixes everything that is fitted.** Bounds, centre, scale, rotation
and per-PC scaling all come out of one pass over one matrix of repertoires. There is no default: a
signature is comparable to another one only if both were rotated through the same corpus, and
nothing about the numbers would say otherwise.

Until 4.0 the rotation was *fit-free* — the PCA of the bundled prototype panel, 10,000 individual
**clonotypes** — while every one of the 399 PC columns it produced was a **repertoire** statistic,
and its centre and scale were fitted on zero rows of that artifact and arrived from a separate
corpus. Two independent fits, stitched; one shipped reference paired a centre of exactly `0.0` with
a scale plainly fitted from data and put a corpus-typical sample **81 robust deviations** out.

All four shipped corpora need no cohort at all and are reproducible by anyone who installs the
library — `synthetic-blood` and `synthetic-tissue` (the naive/memory mixture drawn across a
compartment's measured per-locus ladders), plus the pure regimes `naive` (every clone size 1, what
the recombination model emits) and `memory` (Zipf rank-abundance clone sizes). `deep-tcr`, `blood`
and `tissue` will be fitted on real cohorts and are not shipped yet.

**Holes are never zeros.** An unsequenced locus, a compartment below its clonotype floor, or a
statistic the sample is too shallow to estimate is `nan` plus a `mask:` column, because a model that
reads "absent" as "zero" reads an unsequenced chain as biology.

Full documentation — the corpora, the supports that decide which tail is trimmed, what is fitted and
what is not: [**Signature**](https://docs.isalgo.dev/mirpy/signature.html).

## Exposure trajectory, generative loop, digital twin

A repertoire cohort often has a **known covariate** (HLA, batch, vaccine arm) but an **unknown or
noisy progression axis** (days since exposure, response severity). `mir.track.fit_exposure_trajectory`
recovers that latent trajectory `tau` while disentangling it from the covariate — a PhenoPath-style
model (Campbell & Yau 2018, *Nat. Commun.*
[10.1038/s41467-018-04696-6](https://doi.org/10.1038/s41467-018-04696-6), adapted from genes×cells to
repertoire-channels×samples) over *any* per-sample channel matrix (a stacked `Φ`, a `ChannelBuilder`
build, or a raw embedding block):

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

X, spec = stack_embeddings(embs)                 # channels: mean ‖ diversity ‖ second
fit = fit_exposure_trajectory(X, hla_indicator, channel_names=spec.names)
fit.tau                    # inferred exposure/progression pseudotime, one per sample
fit.top_interactions(5)    # which channels respond to progression differently by covariate
```

Complementing that, `mir.generate.DescriptorDensity` fits a (optionally class-conditional) density
over `RepertoireDescriptor` vectors — `sample` draws brand-new synthetic donor states, `evolve`
perturbs one donor along a coordinate ("what if hotter") and propagates the coupled shift through
every other coordinate via the fitted covariance's conditional mean. `mir.ml.diffusion` (`[ml]`
extra) is the non-linear alternative — a compact conditional DDPM/DDIM generator with
classifier-free guidance, sharing the same `sample(n, condition=…)` call shape so it drops in
unchanged. `mir.twin.DonorTwin` glues a donor's descriptor + trajectory position + covariate into
one object you perturb or resample through, instead of threading the three APIs together by hand:

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

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)
synthetic = twins[0].simulate(density, n=20)     # 20 new synthetic peers of donor 0
```

See `examples/trajectory_and_twin.py` for a runnable end-to-end demo.

## Reproduce the paper

The self-contained theory notebooks run on bundled data:

```bash
pip install "mirpy-lib[examples]"
marimo edit examples/theory.py          # supplementary S1–S3 (distance laws, D↔d, prototype robustness)
marimo edit examples/quickstart.py      # embed + cluster (epitope colours need a local VDJdb dump)
```

The full benchmark suite (VDJdb Table S1, density, repertoire/TCGA) and result docs live in the
companion analysis repo [`2026-mirpy-analysis`](https://github.com/antigenomics) — this repo is the
library + CI tests only.

Method: Kremlyakova *et al.*, *TCREMP: a bioinformatic pipeline for efficient embedding of
T-cell receptor sequences*, **J Mol Biol** 437 (2025) 169205.

## Performance & parallelism

mirpy is CPU-parallel by default and uses the GPU for the neural codecs. Knobs, by hot path:

| Stage | Knob | Default | Notes |
|---|---|---|---|
| Embedding (junction distance) | `TCREmp(..., threads=N)` | `0` = **all cores** | The C++ `seqtree.gapblock` scorer; releases the GIL, ~530 M pairs/s @16 cores. `threads=1` for a serial run. |
| Density kNN / balloon | `neighbor_enrichment(..., backend=…)` | `"kdtree"` (scipy cKDTree, **all cores**) | Exact and multithreaded (`workers=-1`), 5–9× faster than the BallTree baseline. `backend="ann"` = pynndescent, auto all-core, ~30× at ≥1e5 (observed side approximate/conservative, background exact); `backend="exact"` = the 1-core BallTree baseline, for reproducing older runs. |
| Clustering | `cluster(..., n_jobs=-1)` | sklearn default (1) | forwarded to DBSCAN/OPTICS/HDBSCAN via `**kwargs`; parallelizes the neighbour search. |
| BLAS (PCA, RFF, matmul) | `OMP_NUM_THREADS` / `OPENBLAS_NUM_THREADS` env | all cores | numpy/sklearn use the platform BLAS; cap via env if oversubscribed. |
| Neural codecs (`mir.ml`) | `pick_device()` / `device=` / `MIR_DEVICE` env | **CUDA → MPS → CPU**, auto | every `train_*` / codec / bundle takes `device=`; e.g. `MIR_DEVICE=cuda:1` pins the second GPU. Torch-free paths (`density`, `repertoire`) never touch the GPU. |

Rule of thumb: leave `threads=0` (all cores) for embedding; leave density on the default
`backend="kdtree"` (exact, multicore) and switch to `"ann"` only at whole-repertoire scale; the GPU
is used only by `mir.ml`.

## Development

Repo-local `.venv` via [uv](https://docs.astral.sh/uv/) (bash/zsh): `bash setup.sh` — add
`--dev-parents` to editable-install the sibling `seqtree` / `vdjtools` / `vdjmatch` checkouts and
`--tests` to run the fast suite. Tests: `python -m pytest tests/ -q`. See [`CLAUDE.md`](CLAUDE.md)
for the architecture and reuse map.
