From c4b3c50394e3e95aec9ddb5b52b3927290ed9320 Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 19 Jul 2026 05:13:46 +0300 Subject: [PATCH] =?UTF-8?q?chore:=20src-layout,=20notebooks=E2=86=92exampl?= =?UTF-8?q?es,=20distribute=20working=20docs=20out=20of=20root?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Repo hygiene for the 3.4.0 release; no change to the installed package. - src-layout: mir/ -> src/mir/ (hatch packages=["src/mir"] maps it to `mir` at the wheel root; import name and API unchanged). Coverage/omit paths + docs sys.path updated. - notebooks/ -> examples/ (the [examples] extra + docs/examples page updated). - Removed the working result/plan docs from the repo root, distributed to the sibling repos: THEORY.md -> 2026-mirpy-ms; BENCHMARKS.md, REPERTOIRE_{EMBEDDING,LESSONS}.md, SQRT_D_MIGRATION.md, ROADMAP.md -> 2026-mirpy-analysis/benchmarks/. Root keeps README/CHANGELOG/CLAUDE/SOURCES. README/docs/CLAUDE references repointed. Verified: 131 passed / 6 skipped (fast tier) in the uv .venv; zero-warning Sphinx build; py3-none-any wheel (mir/ at root, resources bundled) + twine check PASSED. Co-Authored-By: Claude Opus 4.8 --- .gitignore | 6 +- BENCHMARKS.md | 264 -------------- CHANGELOG.md | 5 + CLAUDE.md | 14 +- README.md | 11 +- REPERTOIRE_EMBEDDING.md | 186 ---------- REPERTOIRE_LESSONS.md | 78 ----- ROADMAP.md | 91 ----- SQRT_D_MIGRATION.md | 157 --------- THEORY.md | 322 ------------------ docs/conf.py | 4 +- docs/examples.rst | 10 +- docs/index.rst | 5 +- docs/usage.rst | 6 +- .../__marimo__/session/quickstart.py.json | 0 .../__marimo__/session/theory.py.json | 0 {notebooks => examples}/density.py | 0 {notebooks => examples}/quickstart.py | 0 {notebooks => examples}/theory.py | 0 pyproject.toml | 12 +- {mir => src/mir}/__init__.py | 0 {mir => src/mir}/aliases.py | 0 {mir => src/mir}/alleles.py | 0 {mir => src/mir}/bench/__init__.py | 0 {mir => src/mir}/bench/eval.py | 0 {mir => src/mir}/bench/metrics.py | 0 {mir => src/mir}/bench/theory.py | 0 {mir => src/mir}/bench/vdjdb.py | 0 {mir => src/mir}/cli.py | 0 {mir => src/mir}/cohort.py | 0 {mir => src/mir}/density.py | 0 {mir => src/mir}/distances/__init__.py | 0 {mir => src/mir}/distances/germline.py | 0 {mir => src/mir}/distances/junction.py | 0 {mir => src/mir}/embedding/__init__.py | 0 {mir => src/mir}/embedding/pca.py | 0 {mir => src/mir}/embedding/presets.py | 0 {mir => src/mir}/embedding/prototypes.py | 0 {mir => src/mir}/embedding/tcremp.py | 0 {mir => src/mir}/explain.py | 0 {mir => src/mir}/ml/__init__.py | 0 {mir => src/mir}/ml/bundle.py | 0 {mir => src/mir}/ml/codec.py | 0 {mir => src/mir}/ml/decoder.py | 0 {mir => src/mir}/ml/encoder.py | 0 {mir => src/mir}/ml/set_encoder.py | 0 {mir => src/mir}/ml/tokenize.py | 0 {mir => src/mir}/ml/train.py | 0 {mir => src/mir}/py.typed | 0 {mir => src/mir}/repertoire.py | 0 {mir => src/mir}/resources/__init__.py | 0 .../mir}/resources/gene_library/README.md | 0 .../mir}/resources/gene_library/__init__.py | 0 .../gene_library/build_gene_library.py | 0 .../gene_library/build_region_annotations.py | 0 .../gene_library/imgt_gene_library.txt | 0 .../gene_library/olga_gene_library.txt | 0 .../gene_library/region_annotations.txt | 0 .../germline_dist/build_germline_dist.py | 0 .../resources/germline_dist/human_IGH.npz | Bin .../resources/germline_dist/human_IGK.npz | Bin .../resources/germline_dist/human_IGL.npz | Bin .../resources/germline_dist/human_TRA.npz | Bin .../resources/germline_dist/human_TRB.npz | Bin .../resources/germline_dist/human_TRD.npz | Bin .../resources/germline_dist/human_TRG.npz | Bin .../resources/germline_dist/mouse_IGH.npz | Bin .../resources/germline_dist/mouse_IGK.npz | Bin .../resources/germline_dist/mouse_IGL.npz | Bin .../resources/germline_dist/mouse_TRA.npz | Bin .../resources/germline_dist/mouse_TRB.npz | Bin .../resources/germline_dist/mouse_TRD.npz | Bin .../resources/germline_dist/mouse_TRG.npz | Bin .../mir}/resources/prototypes/.gitkeep | 0 .../prototypes/generate_prototypes.py | 0 .../mir}/resources/prototypes/human_IGH.tsv | 0 .../mir}/resources/prototypes/human_IGK.tsv | 0 .../mir}/resources/prototypes/human_IGL.tsv | 0 .../mir}/resources/prototypes/human_TRA.tsv | 0 .../mir}/resources/prototypes/human_TRB.tsv | 0 .../mir}/resources/prototypes/human_TRD.tsv | 0 .../mir}/resources/prototypes/human_TRG.tsv | 0 .../mir}/resources/prototypes/manifest.json | 0 .../mir}/resources/prototypes/mouse_TRA.tsv | 0 .../mir}/resources/prototypes/mouse_TRB.tsv | 0 85 files changed, 42 insertions(+), 1129 deletions(-) delete mode 100644 BENCHMARKS.md delete mode 100644 REPERTOIRE_EMBEDDING.md delete mode 100644 REPERTOIRE_LESSONS.md delete mode 100644 ROADMAP.md delete mode 100644 SQRT_D_MIGRATION.md delete mode 100644 THEORY.md rename {notebooks => examples}/__marimo__/session/quickstart.py.json (100%) rename {notebooks => examples}/__marimo__/session/theory.py.json (100%) rename {notebooks => examples}/density.py (100%) rename {notebooks => examples}/quickstart.py (100%) rename {notebooks => examples}/theory.py (100%) rename {mir => src/mir}/__init__.py (100%) rename {mir => src/mir}/aliases.py (100%) rename {mir => src/mir}/alleles.py (100%) rename {mir => src/mir}/bench/__init__.py (100%) rename {mir => src/mir}/bench/eval.py (100%) rename {mir => src/mir}/bench/metrics.py (100%) rename {mir => src/mir}/bench/theory.py (100%) rename {mir => src/mir}/bench/vdjdb.py (100%) rename {mir => src/mir}/cli.py (100%) rename {mir => src/mir}/cohort.py (100%) rename {mir => src/mir}/density.py (100%) rename {mir => src/mir}/distances/__init__.py (100%) rename {mir => src/mir}/distances/germline.py (100%) rename {mir => src/mir}/distances/junction.py (100%) rename {mir => src/mir}/embedding/__init__.py (100%) rename {mir => src/mir}/embedding/pca.py (100%) rename {mir => src/mir}/embedding/presets.py (100%) rename {mir => src/mir}/embedding/prototypes.py (100%) rename {mir => src/mir}/embedding/tcremp.py (100%) rename {mir => src/mir}/explain.py (100%) rename {mir => src/mir}/ml/__init__.py (100%) rename {mir => src/mir}/ml/bundle.py (100%) rename {mir => src/mir}/ml/codec.py (100%) rename {mir => src/mir}/ml/decoder.py (100%) rename {mir => src/mir}/ml/encoder.py (100%) rename {mir => src/mir}/ml/set_encoder.py (100%) rename {mir => src/mir}/ml/tokenize.py (100%) rename {mir => src/mir}/ml/train.py (100%) rename {mir => src/mir}/py.typed (100%) rename {mir => src/mir}/repertoire.py (100%) rename {mir => src/mir}/resources/__init__.py (100%) rename {mir => src/mir}/resources/gene_library/README.md (100%) rename {mir => src/mir}/resources/gene_library/__init__.py (100%) rename {mir => src/mir}/resources/gene_library/build_gene_library.py (100%) rename {mir => src/mir}/resources/gene_library/build_region_annotations.py (100%) rename {mir => src/mir}/resources/gene_library/imgt_gene_library.txt (100%) rename {mir => src/mir}/resources/gene_library/olga_gene_library.txt (100%) rename {mir => src/mir}/resources/gene_library/region_annotations.txt (100%) rename {mir => src/mir}/resources/germline_dist/build_germline_dist.py (100%) rename {mir => src/mir}/resources/germline_dist/human_IGH.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_IGK.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_IGL.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_TRA.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_TRB.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_TRD.npz (100%) rename {mir => src/mir}/resources/germline_dist/human_TRG.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_IGH.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_IGK.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_IGL.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_TRA.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_TRB.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_TRD.npz (100%) rename {mir => src/mir}/resources/germline_dist/mouse_TRG.npz (100%) rename {mir => src/mir}/resources/prototypes/.gitkeep (100%) rename {mir => src/mir}/resources/prototypes/generate_prototypes.py (100%) rename {mir => src/mir}/resources/prototypes/human_IGH.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_IGK.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_IGL.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_TRA.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_TRB.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_TRD.tsv (100%) rename {mir => src/mir}/resources/prototypes/human_TRG.tsv (100%) rename {mir => src/mir}/resources/prototypes/manifest.json (100%) rename {mir => src/mir}/resources/prototypes/mouse_TRA.tsv (100%) rename {mir => src/mir}/resources/prototypes/mouse_TRB.tsv (100%) diff --git a/.gitignore b/.gitignore index cffd24f..3aeccef 100644 --- a/.gitignore +++ b/.gitignore @@ -10,7 +10,7 @@ tests/assets/vdjdb.slim.txt.gz tests/assets/real_repertoires/ tests/assets/srx_repertoires/ *.airr.tsv -notebooks/assets/ +examples/assets/ ## VSCode @@ -184,10 +184,10 @@ cython_debug/ # and can be added to the global gitignore or merged into this file. For a more nuclear # option (not recommended) you can uncomment the following to ignore the entire idea folder. #.idea/ -mir/docs +src/mir/docs .idea public -notebooks/assets/derived/ +examples/assets/derived/ # ML run artifacts (checkpoints/metrics) — never committed experiments/runs/ diff --git a/BENCHMARKS.md b/BENCHMARKS.md deleted file mode 100644 index ae91555..0000000 --- a/BENCHMARKS.md +++ /dev/null @@ -1,264 +0,0 @@ -# Benchmark results — regression baseline - -Recorded results for the sample-level (repertoire) embedding benchmarks (`experiments/benchmark_repertoire_*.py`), -so a re-run can be checked for regression. All numbers are **repeated 50-fold CV mean±std** unless noted, on -arda-native coordinates. Theory/interpretation: `THEORY.md` T7. Data provenance: `SOURCES.md`. - -- **Env**: conda `mirpy` (Python 3.12), Apple M3 (32 GB). Seeds fixed (`seed=0`). `RUN_BENCHMARK` not required - for these (they are experiment scripts, not the guarded pytest tier). -- **Recorded**: 2026-07-14. Re-run and compare; a drop beyond ~1 std on the headline metric is a regression. -- **Reproducibility**: all load **local data first, HF fallback** — `aging`/`airr_hip` via `_hf.fetch` - (`hf_hub_download`), `airr_covid19` via `_covid.covid_path` (local `~/hf/airr_covid19` else `isalgo/airr_covid19`). - -## Sample-level embedding — established regime (airr_benchmark aging, airr_hip) - -| Benchmark | Command | Headline metric | Value | -|---|---|---|---| -| Depth-robustness (`prop:kme`) | `benchmark_repertoire_depth.py` | log–log slope `‖ΔΦ₁‖` vs `n_eff` (theory −0.5) | **−0.55** | -| Age regression (n=79) | `benchmark_repertoire_aging.py` | \|Spearman(age,pred)\|: diversity / kernel-mean / k-mer | **0.76** / 0.58 / 0.28 | -| CMV serostatus (n=240) | `benchmark_repertoire_cmvhla.py` | AUC: diversity / Φ blocks / learned | **0.83±0.05** / 0.59–0.63 / 0.49 | -| CMV age-matched control | `benchmark_repertoire_cmvhla.py` | AUC: diversity / age-only (decade-matched) | **0.83** / 0.45 | -| HLA-A\*02 (n=500, 25k) | `benchmark_repertoire_hla.py 250 25000` | AUC: second-moment / diversity | **0.623±0.048** / 0.460±0.048 | -| HLA panel (airr_hip, n=300) | `benchmark_repertoire_hla*panel` | second-moment AUC: B\*08 / B\*07 / B\*44 | **0.82** / 0.76 / 0.66 | -| Spike-in recovery | `benchmark_repertoire_spikein.py` | recall @N≤3k / FPR (clean P_gen null) | **~35–50%** / ~1.2% | - -## Aging divergence — `benchmark_repertoire_agediverge.py 500000 30` (58 donors @500k) - -Divergence rises with age at depth but is **diversity-coupled, not an independent axis**. - -| Metric | Value | Note | -|---|---|---| -| ρ(¹D, divergence) biased KME | −0.153 | biased V-stat `1/n_eff` artifact | -| ρ(¹D, divergence) **unbiased** KME | **−0.052** | `mmd_matrix(unbiased=True)`; artifact removed | -| ρ(¹D, divergence) overlap −logF | −0.680 | clonal-expansion driven | -| ρ(age, divergence overlap F) | 0.704 | (0.24 at 250k — depth-dependent) | -| **partial ρ(age, divergence \| ¹D)** | **0.067** (p=0.64) | n.s. — no independent axis (same at 40k/250k) | - -## COVID cohort (airr_covid19, 300 donors, ≤20k reads) - -### Batch cookbook — `benchmark_repertoire_covidbatch.py 300 20000` → **PASS** - -| signal | naive AUC | batch-residualized | -|---|---|---| -| batch identity (OvR) | 0.777 | **0.032** | -| COVID status (⟂̸ batch) | **0.660** | 0.409 | -| HLA-A\*02 (⟂ batch) | 0.596 | **0.609** | - -within-mixed-batch COVID (honest) **0.539**; MMD offset:biology ratio **1.05**. - -### HLA imprint class I+II — `benchmark_repertoire_covidhla.py 300 20000` (TRB) - -second-moment > diversity in direction **15/17** alleles; class-II present **8/9**. - -| allele | class | 2nd-moment AUC | Δ vs diversity | -|---|---|---|---| -| DRB1\*07:01 | II | **0.758** | +0.204 | -| B\*07:02 | I | 0.702 | +0.154 | -| DRB1\*15:01 | II | 0.702 | +0.110 | -| DQB1\*05:01 | II | 0.619 | +0.152 | -| A\*02:01 | I | 0.628 | +0.035 | - -### Paired α+β — `benchmark_repertoire_covidpaired.py 300 20000` - -**TRA carries the stronger HLA imprint**; paired concat sits between (noisier β dilutes). - -| allele | β (TRB) | α (TRA) | α+β paired | -|---|---|---|---| -| A\*02 | 0.538 | **0.636** | 0.612 | -| B\*07 | 0.646 | **0.701** | 0.664 | -| DRB1\*15 | 0.749 | **0.810** | 0.796 | - -### COVID biomarker — honest **negative** (`covidstatus.py` / `covidpaired.py`) - -| test | β | α | paired | -|---|---|---|---| -| COVID status, naive | 0.673 | 0.685 | **0.702** | -| COVID status, batch-resid | 0.488 | 0.512 | **0.522** (chance) | -| witness rediscovery of paper's clones | 0.373 | **0.450** | — | - -Convalescent (long-past) COVID leaves **no batch-robust bulk clonotype-identity signal** at RNA-seq depth; the -naive 0.66–0.72 was batch confound. Ground truth is 87% α (4393 α vs 567 β). - ---- - -## 2026-07-14 (pm) — spectral block, COVID witness, TCGA survival - -### Spectral interaction block — `benchmark_repertoire_spectral.py 300 20000` - -The opt-in top-`r` **eigenvalue** compaction of the second-moment block (new `n_eigs=` on `fit_repertoire_space`) -is **lossy for the HLA imprint**: HLA-A\*02 carriage lives in *which* public clones co-occur (directional), and -a rotation-invariant eigenvalue spectrum discards it. Default (`n_eigs=None`) keeps the full upper triangle, -unchanged. - -| representation | block dim | HLA-A\*02 AUC | -|---|---|---| -| diversity (baseline) | 4 | 0.560±0.066 | -| D₂=512 top-16 eigvals | 16 | 0.545±0.093 | -| D₂=512 top-128 eigvals | 128 | 0.536±0.092 | -| D₂=256 full upper-tri | 32896 | 0.562±0.078 | -| **D₂=512 full upper-tri** | 131328 | **0.593±0.070** | - -### COVID witness — `benchmark_repertoire_covidwitness.py 300 20000` - -Why per-clonotype Fisher-significant clones "vanish" in the bulk embedding: **cohort breadth**, not read depth. - -| level | finding | -|---|---| -| Fisher genome-wide BH (150–300 donors) | **0** clones pass at any depth (20k or 120k reads) | -| Fisher full cohort (~1137 donors, user tmp scan) | **39 β / 4 α** clones pass BH (min q≈5e-4) | -| witness β: whole-cohort → **mixed-batch** | 0.510 → **0.752** (batch control is the lever) | -| witness β: HLA-strat median / best | 0.532 / 0.712 (best = max-selection; median ≈ whole-cohort) | -| witness α: whole / mixed / strat-best | 0.463 / 0.424 / 0.473 (no recovery — α witness fails) | - -Batch control (mixed-batch contrast, where COVID+healthy share runs) recovers the paper's β clones; per-allele -HLA stratification adds more noise than signal at these donor counts. HLA+α+β is **not** the missing key — breadth is. - -### COVID motif recovery at full breadth — `benchmark_repertoire_covidmotif.py 0 60000` - -Using the tool that *works* at breadth (Emerson incidence Fisher) on the **full 1137-donor cohort** (761 COVID / -376 healthy, 60 k reads), both chains, to actually *find COVID motifs*: - -| chain | genome-wide clones q<0.05 | GT-true recovered | + HLA-restricted | -|---|---|---|---| -| **α** | 7 | **4** | +0 | -| β | 0 | 0 | +0 | - -The 4 recovered α clones are one coherent public family — **CAG·NYGGSQGNLIF** (single-residue variants), all -members of the paper's COVID-associated α **cluster 31** (27 GT clones, all `has_covid_association`). **HLA -restriction adds nothing** because these α clones are already *public* (present across HLA backgrounds) — -per-allele restriction only loses power. - -**WS1 (2026-07-15) — the α/β loop closed, and it is NOT depth (refutes the earlier "β is depth-limited" note).** -Running β at **full *native* depth** (no downsampling) over the full 1137-donor cohort -(`benchmark_repertoire_covidmotif.py 0 0`) still passes **1** feature genome-wide at any `min_incidence∈{2,3,5}` -and recovers **0** GT-β clones — full depth does not rescue β. Nor does the *restricting allele*: a per-stratum -sweep (`benchmark_repertoire_covidmotif.py native_beta`, 12 class-II/B HLA carrier groups, 117–346 carriers each -at **native depth** ~3.5 M reads/donor) recovers **0** GT-β in every stratum — refuting "β needs full depth + -the restricting allele." And it is **not absence**: **253/256 GT-β clones -are present** in the cohort (median incidence **9** donors, max **110**). The reason is that the paper's β -"COVID-associated" clones are short **public bystanders** (`CASSx…YEQYF`, TRBJ2-7 — e.g. `CASSYSYEQYF` in 110 -donors) occurring in COVID **and** healthy donors alike, so they do not discriminate status by *incidence*. -**The α/β asymmetry is which clones are status-enriched vs public bystanders — not depth or breadth:** α's -recovered `CAG·NYGGSQGNLIF` family is a genuine status-enriched public expansion; β's GT set is dominated by -non-discriminating convergent publics. *Caveat:* the incidence Fisher tests presence-breadth, not clonal- -expansion magnitude — a frequency/abundance-based test could still recover β clones enriched in *size* not -*breadth* (untested). Net: the recoverable "find-motifs" signal is the **public α compartment** (breadth- -powered, no HLA needed); the breadth-starved bulk MMD witness recovered none of it. - -### TCGA — `benchmark_repertoire_tcga.py BRCA,LUAD,KIRC 50000` - -Tumor-type separation is **depth-dependent** (deep IG light chains ≫ shallow TR); **survival adds nothing** over -clinical covariates. Base Cox C-index (age+sex+stage+log reads): BRCA 0.72, LUAD 0.66, KIRC 0.73. - -| chain | tumor-type macro-OvR AUC | survival ΔC (BRCA / LUAD / KIRC) | -|---|---|---| -| TRB (median ~24 clonotypes) | 0.523 | +0.003 / −0.003 / +0.002 | -| IGH | 0.500 | +0.002 / +0.001 / −0.004 | -| **IGK (deepest chain)** | **0.666** | −0.001 / −0.004 / −0.001 | -| concat TRB+IGH | 0.538 | −0.002 / +0.001 / +0.001 | - -IGK's 0.67 shows the method works when clonotype depth suffices (the deep IG compartment carries tissue signal); -the flat ΔC across every chain says the tumour-infiltrating repertoire adds **no prognostic value beyond clinical -covariates** at TCGA RNA-seq depth. - -### TCGA survival — biology-grounded features — `benchmark_repertoire_tcga_survival.py` - -The clonotype *embedding* adds nothing, but interpretable **biology axes** (isotype from `c_call`, infiltration / -hot-vs-cold, atypicality = gene-usage divergence from the tumour-type centroid, clonal expansion — modelled on an -internal AIRR-tissue EDA) do, cancer-specifically. Base Cox = age+sex+stage+log reads; **bold** = best gain per row. - -C-index gain over clinical (5-fold CV): - -| cancer | C base | isotype | infiltration | atypicality | clonality | all-AIRR | -|---|---|---|---|---|---|---| -| SKCM (melanoma) | 0.609 | +0.008 | +0.034 | −0.001 | +0.031 | **+0.036** | -| KIRP (renal pap.) | 0.717 | −0.008 | **+0.030** | −0.003 | +0.014 | +0.015 | -| LGG (glioma) | 0.764 | −0.006 | **+0.019** | +0.008 | +0.016 | +0.011 | -| KIRC (renal cc.) | 0.721 | +0.005 | **+0.009** | +0.002 | +0.002 | +0.006 | - -(LUAD / STAD / OV / BLCA: AIRR gain ≈ 0 or negative.) - -KM median-split log-rank p (the EDA's stratification test — surfaces threshold effects the linear C-index misses; -**bold** = p<0.05): - -| cancer | infiltration | IgA fraction | atypicality | -|---|---|---|---| -| KIRC | **0.002** | 0.089 | 0.880 | -| SKCM | **0.000** | **0.026** | 0.716 | -| KIRP | **0.020** | 0.750 | 0.555 | -| OV | **0.044** | 0.581 | 0.691 | -| LGG | **0.022** | 0.076 | **0.001** | -| BLCA | 0.676 | **0.010** | 0.243 | - -**Read:** infiltration (hot/cold) stratifies survival in **5/8** cancers (melanoma, both renal, ovarian, glioma) -and adds C-index in melanoma/KIRP; **IgA** (isotype / mucosal) stratifies **bladder** + melanoma; **atypicality** -stratifies **glioma** — each matching known immune biology. Prognosis lives in infiltration magnitude + isotype + -typicality, **not** clonotype identity, and is cancer-specific — the reason the identity embedding's ΔC was flat. - -## Repertoire embedding for TME & survival — pan-cancer (33 TCGA types) - -The reframing: the prognostic axes above are **channels of one TME-aware, multi-chain repertoire embedding** -Φ(S) (`_tcga_embedding.py`: per-chain identity ‖ diversity ‖ coverage/infiltration, + isotype + composition -+ atypicality). `benchmark_repertoire_tcga_pancancer.py ALL` fits Φ once over 9 425 OS-annotated samples -(78-dim, 5 embeddable chains + all 7 in the composition channel) and, per cancer, reports the CV **ΔC-index** -of clinical+Φ over clinical (age+sex+stage+log reads) and a **likelihood-ratio p** for the Φ block. - -Cancers where Φ is robustly prognostic (**both** LR p<0.05 **and** CV ΔC>0; bold): - -| cancer | n | events | C base | C+Φ | ΔC | LR p | top channel | -|---|---|---|---|---|---|---|---| -| **SKCM** (melanoma) | 460 | 219 | 0.609 | 0.647 | **+0.039** | **0.000** | coverage | -| **BLCA** (bladder) | 411 | 182 | 0.650 | 0.675 | **+0.025** | **0.002** | identity | -| **HNSC** (head & neck) | 500 | 218 | 0.617 | 0.639 | **+0.022** | **0.012** | coverage | -| **LGG** (glioma) | 423 | 112 | 0.764 | 0.779 | **+0.016** | **0.000** | coverage | - -Strong effect-size positives (ΔC>0, LR n.s. — power-limited): SARC +0.038 (isotype), KIRP +0.034 (atypicality), -LUAD +0.009, KIRC +0.008, BRCA +0.007. Immune-cold / small cohorts are flat-to-negative (overfit): LIHC −0.029, -STAD −0.026, PAAD −0.029, UCEC −0.072 (n=181/35ev). - -**Pan-cancer (20 evaluable):** mean ΔC ≈ 0 (median +0.003), Φ significant in 5/20 — most-informative channel -tally **coverage 5 · atypicality 5 · composition 4 · isotype 3 · identity 2 · diversity 1**. The signal is -cancer-specific and lives in the **TME channels (infiltration, atypicality, composition), not clonotype -identity**; the pan-cancer mean is ≈0 because immune-cold cohorts overfit and cancel the immune-hot wins — -so the honest claim is *Φ adds significant prognostic value in immunologically active cancers* (melanoma, -glioma, head&neck, bladder), not universally. - -### TME states (unsupervised) — `benchmark_repertoire_tcga_tme.py ALL 6` - -KMeans (k=6) on Φ's TME channels (9 425 samples) recovers coherent, interpretable microenvironment states -(z-scored channel means; **bold** = defining extreme). The UMAP (`experiments/figures/umap_tcga_tme`) shows a -clean **infiltration gradient** (hot → cold) organising the pan-cancer cohort. - -| state | n | infiltration | T-vs-B | switch | diversity | death | HR (vs ref) | enriched cancers | -|---|---|---|---|---|---|---|---|---| -| cold-humoral | 1565 | −0.53 | −0.36 B | +0.39 | −0.40 | 0.33 | **1.14** (p=0.011) | BRCA, BLCA, PRAD | -| hot-diverse | 2090 | **+1.03** | −0.70 B | +0.35 | **+0.92** | 0.31 | 0.92 (p=0.072) | BRCA, LUAD, STAD | -| T-balanced | 1277 | −0.10 | **+0.86 T** | +0.15 | +0.06 | 0.29 | 1.09 | KIRC, PRAD, BRCA | -| cold-T unsw. | 883 | **−1.58** | **+1.72 T** | **−2.20** | −1.66 | 0.30 | 1.05 | LGG, LIHC, SARC | -| cold | 871 | −1.27 | +0.74 T | −0.70 | −1.13 | 0.28 | 1.08 | THCA, LIHC, LGG | -| warm-B (ref) | 2739 | +0.48 | −0.45 B | +0.37 | +0.40 | 0.31 | — | BRCA, LUAD, HNSC | - -Stratified multivariate log-rank across states (blocking tumour type) **p=0.0038** — the states are prognostic -*beyond* cancer type. The **cold-humoral** low-infiltration state carries the **worst** outcome (HR 1.14, p=0.011); -the **hot-diverse** state trends protective (HR 0.92). Glioma (LGG) falls in the cold, strongly T-skewed states, -matching its known low-but-T-biased infiltrate. The repertoire embedding organises the TME unsupervised. - -### In-silico evolution — perturb infiltration, decode the coupled response — `benchmark_repertoire_tcga_insilico.py ALL` - -`mir.repertoire.sample_descriptor` makes every metric a **smooth derivable coordinate** (infiltration = log-mass, -diversity = log n_eff, clonality = Σw², identity = kernel mean), so the cohort's joint distribution is a -generative **manifold** and moving a sample along the infiltration axis (hot↔cold) while staying on it predicts -how the other metrics respond. Gaussian-manifold conditional slopes d(metric)/d(infiltration), pan-cancer mean: - -| response as tumour gets hotter | slope | reading | -|---|---|---| -| diversity | **+0.84** | hot ⇒ more diverse (universal, all cancers) | -| class-switch | +0.52 | hot ⇒ more class-switched B cells | -| IgG | +0.49 | hot ⇒ more IgG (mature humoral) | -| T-vs-B | −0.63 | hot ⇒ more B-skewed (TLS-like) | -| atypicality | −0.23 | hot ⇒ less atypical / more convergent | - -**In-silico "make this tumour hotter"** (move cold→hot on-manifold, read the CoxPH Δlog-HR): predicted -**protective in 12/20** cancers (SKCM −0.49, HNSC −0.63, LUAD −0.44, LAML −0.96 …) — TIL-is-good — but correctly -**adverse in glioma** (LGG **+1.19**, GBM +0.47) and renal (KIRC/KIRP +0.3), matching the known "immune-hot glioma -is worse" biology. The embedding *simulates* the hot↔cold survival axis; the couplings are learned, not imposed. diff --git a/CHANGELOG.md b/CHANGELOG.md index c624cd7..2b0dbdc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,6 +36,11 @@ Python API removed; one optional-dependency group split out. - Documentation overhaul: a use-case-driven user guide, the two CLI commands documented, an examples/notebooks page, `mir.cohort` and `mir.bench.eval` added to the API reference, a logo, and the sample-embedding schematic + real depth-robustness figure. Zero-warning Sphinx build. +- **Repo layout** (no effect on the installed package): adopted the **src-layout** (`mir/` → + `src/mir/`); renamed `notebooks/` → `examples/`; and moved the working result/plan markdown out of + the repo root — `THEORY.md` to the manuscript repo, `BENCHMARKS.md` / `REPERTOIRE_{EMBEDDING,LESSONS}.md` + / `SQRT_D_MIGRATION.md` / `ROADMAP.md` to `2026-mirpy-analysis/benchmarks/`. Root keeps + README / CHANGELOG / CLAUDE / SOURCES. ### Fixed diff --git a/CLAUDE.md b/CLAUDE.md index a7d42ee..456438f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -8,9 +8,12 @@ frozen on branch **`legacy-v2`** (`mirpy-lib` 2.x) — do not develop there. ## Repo split (2026-07-16) — three homes This repo is now **the library + its CI/test data only**. Analysis and manuscript live elsewhere: -- **`~/vcs/code/mirpy`** (here): `mir/` library, unit + fast/slow CI tests (`tests/`, on slim bundled - data), plans for new features, docs, and **slim tabular benchmark results** for README/docs. Benchmark - *scripts* and the theory *appendix* have moved out (below). +- **`~/vcs/code/mirpy`** (here): **src-layout** — the library is `src/mir/` (import `mir`); unit + + fast/slow CI tests (`tests/`, on slim bundled data), runnable `examples/` (marimo), docs. As of the + 3.4.x cleanup the result/plan markdown moved out: **`THEORY.md` → the manuscript repo**; + **`BENCHMARKS.md`, `REPERTOIRE_{EMBEDDING,LESSONS}.md`, `SQRT_D_MIGRATION.md`, `ROADMAP.md` → + `2026-mirpy-analysis/benchmarks/`**. Only `README`/`CHANGELOG`/`CLAUDE`/`SOURCES` stay at the root. + Benchmark *scripts* and the theory *appendix* had already moved out (below). - **`~/vcs/projects/2026-mirpy-analysis`**: all **benchmark scripts** (`benchmarks/` — local + aldan3), full result docs (BENCHMARKS/THEORY/…), figures, run outputs, dataset catalog. Run the mirpy library from here; refresh numbers-of-record here. @@ -53,7 +56,7 @@ negative S2) — so prototypes use real reads, not model generation. mirpy is normally read-only to the sibling repos; the `from_arda` builder + a tandem-D generation fix were added to `vdjtools` under the owner's direction (this is that owner's ecosystem). -## Layout (`mir/`) +## Layout (`src/mir/`) - `aliases.py`, `alleles.py` — species/locus + allele normalization. - `distances/junction.py` — `junction_distance_matrix` (seqtree.gapblock). Coordinate knobs (all default to the published space, threaded through `TCREmp`): `metric="squared"` (=`d`, default) | @@ -148,7 +151,8 @@ fix were added to `vdjtools` under the owner's direction (this is that owner's e GPU only in `mir.ml`: `pick_device()` = **CUDA → MPS → CPU** auto, override `device=`/`MIR_DEVICE`. ## Open loops / next steps -- **Embedding-tier roadmap** (`ROADMAP.md`, 2026-07-17) — the "vdjtools at the embedding level" audit + +- **Embedding-tier roadmap** (`2026-mirpy-analysis/benchmarks/ROADMAP.md`, moved out in the 3.4.x + cleanup) — the "vdjtools at the embedding level" audit + plan (three verbs: make / measure / generate-decode). **Phase 0** (robustness + optimization quick wins) and **Phase 1** (cohort tier: `bench/eval.py`, `repertoire.{fit_repertoire_spaces,centroid_atypicality}`, `cohort.py` digital donor) are **DONE**. Next: **Phase 2** generative loop (`generate.py` `DescriptorDensity` diff --git a/README.md b/README.md index 3d7aa59..1c4b5cd 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ **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 (see [`THEORY.md`](THEORY.md)). +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 @@ -186,8 +186,9 @@ a batch offset is first-order and cancels, while a batch-orthogonal signal (e.g. 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). See [`BENCHMARKS.md`](BENCHMARKS.md) -and [`THEORY.md`](THEORY.md) T7 (the benchmark scripts themselves live in the analysis repo). +(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. ## Reproduce the paper @@ -195,8 +196,8 @@ The self-contained theory notebooks run on bundled data: ```bash pip install "mirpy-lib[examples]" -marimo edit notebooks/theory.py # supplementary S1–S3 (distance laws, D↔d, prototype robustness) -marimo edit notebooks/quickstart.py # embed + cluster VDJdb antigens +marimo edit examples/theory.py # supplementary S1–S3 (distance laws, D↔d, prototype robustness) +marimo edit examples/quickstart.py # embed + cluster VDJdb antigens ``` The full benchmark suite (VDJdb Table S1, density, repertoire/TCGA) and result docs live in the diff --git a/REPERTOIRE_EMBEDDING.md b/REPERTOIRE_EMBEDDING.md deleted file mode 100644 index b98062d..0000000 --- a/REPERTOIRE_EMBEDDING.md +++ /dev/null @@ -1,186 +0,0 @@ -# Sample-level (repertoire) embedding — build spec - -_Self-contained guideline for a separate build session. Theory is in `appendix/tcremp_theory.tex` -§T.7 (`sec:sample`); this doc is the module + benchmark plan. Written 2026-07-14._ - -## 0. What & why - -Embed a whole **repertoire** — an order-invariant multiset of clonotypes `{(x, a_x)}` with clone -counts — into one fixed vector `Φ(S)`, so repertoires can be compared, regressed against phenotype -(age), and classified (CMV/HLA). Four requirements, each owned by one block of `Φ` (appendix §T.7): - -| # | requirement | how | appendix | -|---|---|---|---| -| a | permutation-invariant | pool a per-clonotype feature by a size-weighted sum | Prop. `prop:sampinv` | -| b | robust to sequencing depth | frequency-weighted **kernel mean embedding** (RFF); converges `O(n_eff^{-1/2})` | Prop. `prop:kme` | -| c | preserves diversity | coverage-standardized **Hill profile** `{⁰D,¹D,²D}` | Prop. `prop:antag`,`prop:simdiv` | -| d | HLA-linked interactions | **second moment** `Σ_σ w_σ ψψᵀ` (Fisher) or a learned set-network | Prop. `prop:interact`,`prop:hla` | - -**Central constraint: depth-robustness.** The real target is low-coverage bulk **RNA-seq TCR** — all -chains, but only `10²–10⁴` clonotypes/chain — of which there are 100–200k clinically-annotated samples -(default in cancer immunotherapy). `Φ` must be usable at ~100 clonotypes/chain. `~10³` receptors already -fingerprint an individual (Immprint). - -**The elegant core (answers "which K, which clustering?"): none.** `Φ₁` is the RFF **kernel mean -embedding** `Σ_σ w_σ ψ(φ(x_σ))` — the codebook-free `K→∞` soft-assignment limit of the -global-graph→cluster→histogram heuristic (VLAD/Fisher are its finite-`K` truncations; Prop. `prop:codebook`). -No `K`, no clustering rule. - -## 1. Reuse map (do NOT rebuild — confirmed by exploration) - -- **Clonotype `φ`:** `mir.embedding.tcremp.TCREmp.embed(df) -> (N, 3K) float32` (cols - `v_call,j_call,junction_aa`; counts ignored — read them separately). `PairedTCREmp` for multi-chain. - `mir.embedding.pca.pca_denoise`, `mir.embedding.presets.get_preset`. -- **Count→weight VST (lift directly):** `mir.density._WEIGHTS` = `{"distinct": 1, "log1p": log1p, - "anscombe": sqrt(a+0.375)}`; `mir.density._emp_survival` (orphan `P(A≥a)`); the compound-Poisson - dispersion `φ=E[g²]/E[g]` pattern. These already implement the concave weight `g` of §T.6.9. -- **Diversity / coverage (entire block c is reuse):** `vdjtools.stats.inext` — `inext_coverage`, - `estimate_d(base="coverage", level=C*)`, `asymptotic_diversity` (Chao1/CWJ-Shannon/MVUE-Simpson), - `sample_coverage` (Good–Turing Ĉ), `rarefaction_batch`; `vdjtools.stats.diversity.diversity_stats` - (Shannon `expH`, inverse Simpson, Chao1, Efron–Thisted, d50). -- **Depth control:** `vdjtools.preprocess.downsample(df, size, by="reads")` (hypergeometric, NOT Poisson), - `select_top`. -- **Radius / bandwidth:** `mir.density.calibrate_radius` gives `r₁` (one-substitution embedding drift) — - use as the RFF kernel length-scale so the kernel matches the sequence metric. -- **Baselines to beat:** `vdjtools.features.kmer_profile`, `physchem_profile`; - `vdjtools.overlap.{pairwise_distances, similarity_overlap}` (Morisita-Horn, sequence-weighted - Leinster–Cobbold), `cluster_samples` (MDS/hclust with metadata join). -- **Cohort loading:** `experiments/_hf.py` `fetch`/`load_repertoire` (collapses to - `junction_aa/v_call/j_call/duplicate_count`, `top=N`); `vdjtools.io.batch.read_metadata` / - `read_samples` / `iter_samples`; `io.cohort`. - -## 2. Module — `mir/repertoire.py` (torch-free core; vdjtools + optional torch lazy) - -```python -@dataclass -class RepertoireSpace: # one shared basis for a cohort (comparability invariant) - model: TCREmp | PairedTCREmp - rff: RandomFourierFeatures # omega, b; length-scale = r1 - pca: PCA | None # optional pre-RFF PCA on the pooled clonotype cloud - -def fit_repertoire_space(model, cohort_df, *, n_rff=2048, length_scale=None, - n_components=None, seed=0) -> RepertoireSpace: - """length_scale defaults to calibrate_radius(...)=r1. Fit RFF (+optional PCA) ONCE on the - pooled clonotype cloud so every sample lands in ONE coordinate system (Prop. prop:kme).""" - -def sample_embedding(space, sample_df, *, weight="log1p", blocks=("mean","diversity","second"), - coverage=None) -> np.ndarray: - """Φ(S): concat of requested blocks. Reads duplicate_count for weights. - mean: Φ1 = Σ w_σ ψ(z_σ) (block b; kernel mean embedding) - diversity: coverage-standardized {0D,1D,2D} + Ĉ (block c; via vdjtools.stats.inext) - second: upper-tri of Σ_σ w_σ ψψᵀ (or top-r eigvals) (block d; Fisher vector)""" - -def mmd_distance(phi_a, phi_b) -> float: # ||Φ1(a) − Φ1(b)|| ≈ MMD (Eq. eq:kme) -def mmd_matrix(embs) -> np.ndarray # sample×sample, feed cluster_samples / a regressor - -def hla_stratified_mmd(embs, hla) -> np.ndarray: # Prop. prop:hla: distance only within HLA-matched pairs -``` - -Design notes (all traced to a theorem — see appendix Table `tab:sample`): -- **Frequencies, not counts** (scale-freeness, Prop. `prop:kme`); concave `g` for the Zipf tail - (`weight="log1p"`/`"anscombe"`, §T.6.9). -- **One basis per cohort**: fit RFF (+PCA) on the *pooled* clonotype cloud, embed every sample through it - — same comparability invariant as `mir/ml/bundle.py` and `density.fit_density_space`. Serialize the - space (RFF `omega,b` + prototype hash) so embeddings are only ever compared within one basis. -- **Length-scale = `r₁`** (`calibrate_radius`) so the kernel resolves ~one substitution. -- **`n_eff = (Σ w²)⁻¹`** is a Hill number (Prop. `prop:antag`): report it per sample; it predicts each - sample's own depth-robustness and is the natural x-axis for the depth-robustness benchmark. -- **Multi-chain (RNA-seq):** embed each locus with `PairedTCREmp`-style per-locus `TCREmp`, concatenate - the per-locus `Φ`. Handle missing chains (some RNA-seq samples lack a locus) by zero-blocks + a - presence mask. -- **Nuisance decoupling (§T.7.8):** depth = estimation variance (frequencies + coverage-standardized - diversity handle it, `O(n_eff^{-1/2})`); **batch = a shared shift that cancels in a within-batch - contrast** (Prop. `prop:batch`, the sample-level image of the T6 differential control). So: prefer - within-batch/paired designs; else residualize `Φ` on batch indicators / stratify MMD by batch (as for - HLA); or normalize each sample by its `P_gen` pushforward. Always compare **contrasts, not raw - positions**. Record `batch`/platform per sample in the metadata and expose a `batch=`/`covariates=` arg. -- **RNA-seq normalization (§T.7.9 `sec:samp-norm`, Table `tab:norm`):** the count factorizes - `a_σ ≈ g(σ)·R·θ·c·p_σ`; **σ-independent scalars (depth `R`, infiltration `θ`, per-cell expr. `c`) cancel - under frequency weighting** (`prop:freqquotient`), the **σ-dependent kit gain `g(σ)` is the only surviving - multiplicative nuisance**. Concrete build items for the heterogeneous-cohort / tissue case: - - **Infiltration channel (`prop:infiltration`):** in tissue, `N = R·θ·c` conflates depth and T-cell content; - estimate `θ` from the **TCR read fraction** `N_TCR/N_total ∝ θ` (or a deconvolution signature, Newman 2015) - and **append it as an explicit scalar channel — never down-sample to a common depth** (that erases the - prognostic TIL signal). `Φ` on frequencies is already orthogonal to `θ`. Expose `infiltration=`/ - `read_fraction=` per-sample and concat to the embedding; flag joint-low-`R`/low-`θ` samples (coverage-limited). - - **Read-length / kit (`prop:lengthcensor`):** short reads drop long CDR3s (MNAR length-censoring). Either - restrict to the common recoverable length range across kits, or inverse-recovery-reweight `1/r_L(ℓ)`. - Chemistry bias + confounded-batch residual live in the V/J-usage+length subspace → **anchored subspace - projection** (technical replicate / reference sample = ComBat/Harmony; `P_gen`-predicted V/J usage = RUV - negative-controls), sparing the CDR3-motif directions. Expose `read_length=`/kit per sample. - -## 3. Two co-equal tracks (the learned head is `mir.ml`, torch, `[ml]` extra) - -The unsupervised backbone above and a **learned permutation-invariant set network** are co-equal, -benchmarked head-to-head (§T.7.5). - -- **Unsupervised** (default, label-free, depth-robust): blocks b+c+d as above. Ship first. -- **Learned** (`mir/ml/set_encoder.py`): a **Set Transformer** (ISAB inducing points + PMA pooling, - Lee 2019) or **DeepRC**-style attention MIL (Widrich 2020) over the clonotype `φ`-cloud, weighted by - `g(a_σ)`, supervised by age / CMV / HLA. Inducing points = learned public-cluster detectors (the block-d - signal). Bundle it like the codecs (`mir/ml/bundle.py`): serialize weights + prototype hash + RFF/PCA - basis; refuse cross-basis mixing. **Depth-robustness must be engineered in** (frequency weighting + - subsampling augmentation during training) — set-nets don't get it for free (§T.7.5). - -## 4. Benchmarks — `experiments/benchmark_repertoire_*.py` (`RUN_BENCHMARK`-guarded, subsampled first) - -Reuse the `experiments/_hf.py` fetch pattern. Fit ONE `RepertoireSpace` per cohort; embed all samples; -train a simple head (ridge / logistic / kNN) with a fixed CV split + seed; report vs baselines. - -1. **`benchmark_repertoire_aging.py` — age regression.** Use the **full-depth `aging`** cohort, NOT - `aging_lite` (see §5). Metadata `metadata_aging.txt` (`file_name, sample_id, sex, age`, 79 samples, - ages ~0–103; Britanova 2014/2016). Also `airr_hip` age. **Gate:** Spearman(age, prediction) beats the - `kmer_profile` / `physchem_profile` / diversity-only baselines; `¹D` alone is a strong baseline - (age↓diversity) so `Φ` must add over it. **Batch check (Prop. `prop:batch`, §T.7.8):** the aging cohort - has sequencing batches (sample IDs `A2/A3/A4-*`) — verify age prediction survives batch-adjustment - (leave-one-batch-out CV, or batch as a covariate); a real signal must not be a batch artifact. -2. **`benchmark_repertoire_cmvhla.py` — CMV / HLA (the key test).** `airr_hip` (Emerson 2017 HIP; 786 - subjects, CMV serostatus + HLA-A/B typing). (i) CMV⁺ vs CMV⁻ classification AUC. (ii) **HLA-stratified - MMD** (Prop. `prop:hla`): CMV⁺ samples are close *only* within an HLA-matched stratum — show - cross-CMV `mmd_distance` separates status **within** HLA-matched pairs and is uninformative across - HLA-mismatched pairs. (iii) optionally HLA-allele inference (predict presence of common HLA-A/B alleles). -3. **`benchmark_repertoire_depth.py` — depth-robustness (headline).** Take `aging` samples, `downsample` - to `N ∈ {100, 300, 1000, 3000, 10⁴}` clonotypes (and use `aging_lite` as an independent downsampled - set). Show: `‖Φ₁(sub) − Φ₁(full)‖` decays like `n_eff^{-1/2}` (Prop. `prop:kme`); age/CMV prediction - is preserved down into the RNA-seq-shallow regime (~`10²–10³`). This is the motivating use case. -4. **Multi-chain sanity** — a small all-chain (TRA+TRB) sample: per-locus concat works, missing-chain - mask behaves. - -Per progressive-scaling rules: implement + validate on ≤50 samples / subsampled repertoires, run ONE small -end-to-end per benchmark, report timing; leave full-cohort commands documented for the user. - -## 5. Datasets & SOURCES (add to `SOURCES.md`) - -- **`aging` (full depth)** — HF `isalgo/airr_benchmark`. ⚠ **Use `aging`, not `aging_lite`**: the - `vdjtools_lite/` subdir is the *downsampled* version (keep it only as a depth-robustness control). Locate - the full-depth manifest/files in the repo (age cohort, Britanova). Provenance: experimental TRBβ, ages 6–90. -- **`airr_hip`** — HF `isalgo/airr_hip` = **Emerson et al. 2017** HIP cohort (Nat Genet, PMID 28369038), - 786 subjects (666 discovery + 120 validation), TCRβ + CMV serostatus + HLA-A/B typing (+ age/sex in the - immuneACCESS metadata). Fetch the metadata TSV via `hf_hub_download(repo_id="isalgo/airr_hip", - filename=..., repo_type="dataset")`; confirm exact column names at build time (Dataset Viewer was 503 at - spec time). Provenance: experimental. -- **SRA shallow RNA-seq** — `isalgo/airr_benchmark/sra/meta.tsv` (2993 rows, `PMID Run BioProject Sample`; - all PMID 30830871, BioProject PRJNA511467). The low-coverage-RNA-seq stress set. Provenance: experimental. - -## 6. Phased build checklist - -- [ ] `mir/repertoire.py`: `RandomFourierFeatures`, `fit_repertoire_space`, `sample_embedding` (mean block - first), `mmd_distance` + `__main__` self-check → `tests/test_repertoire.py` (bundled prototypes; - injected two-cohort separation; depth-robustness of `Φ₁` under `downsample`; `n_eff` finite/positive). -- [ ] Diversity block via `vdjtools.stats.inext` (coverage-standardized) + second-moment block. -- [ ] `hla_stratified_mmd`; comparability-invariant serialization of `RepertoireSpace`. -- [ ] `benchmark_repertoire_{aging,cmvhla,depth}.py` (subsampled) + `_hf` loaders; one timed end-to-end each. -- [ ] (co-equal) `mir/ml/set_encoder.py` Set-Transformer/DeepRC + bundle; train scripts; head-to-head vs backbone. -- [ ] Docs: `THEORY.md` "Sample-level embedding" entry (already stubbed), `README` module row, `SOURCES.md` - entries; `CLAUDE.md` open-loop → done. -- [ ] `pyproject.toml`: core already has numpy/scipy/sklearn/polars; `scipy.fft`/RFF need nothing new. - Learned track under `[ml]` (torch). - -## 7. Acceptance - -- Depth-robustness curve matches `n_eff^{-1/2}` (Prop. `prop:kme`); the `n_eff = ^qD` (a Hill number) - identity holds numerically (Prop. `prop:antag`). -- HLA-stratified CMV separation (Prop. `prop:hla`): CMV signal present within HLA-matched strata, absent - across mismatched. -- `Φ` beats the `kmer_profile` / diversity-only / overlap-matrix baselines on age and CMV. -- Works at 100–1000 clonotypes/chain (the RNA-seq regime). diff --git a/REPERTOIRE_LESSONS.md b/REPERTOIRE_LESSONS.md deleted file mode 100644 index c941fa5..0000000 --- a/REPERTOIRE_LESSONS.md +++ /dev/null @@ -1,78 +0,0 @@ -# Sample-level embedding — lessons for the appendix (§T.7) - -Directions for updating `appendix/tcremp_theory.tex` §T.7 (`sec:sample`) with the 2026-07-14 empirical -findings. Each item names the proposition to add/amend, the one-line claim, and the validating script + -recorded number (see `BENCHMARKS.md`, `THEORY.md` T7). Written as directions — the appendix prose is yours. - -## 1. Unbiased MMD is mandatory when depth/diversity varies (amend the MMD definition) - -- **Claim.** The biased V-statistic `‖μ̂_a−μ̂_b‖²` carries a self-term bias `≈ 1/n_eff` (the `k(z,z)` diagonal). - For samples of unequal effective size this bias **dominates** and, worse, correlates with the phenotype - (older = lower `n_eff`), so it *manufactures* a spurious signal. Use the **diagonal-removed MMD²** - (Gretton et al. 2012): `⟨μ̂,μ̂⟩_unbiased = (‖μ̂‖²−Σw²)/(1−Σw²)` with `Σw² = 1/n_eff`. -- **Where.** Amend the `Φ` / MMD `\begin{definition}` to state both estimators and mandate the unbiased one for - cross-depth comparison. Add a one-line remark: *biased MMD is admissible only when `n_eff` is matched across - samples (e.g. a fixed downsample).* -- **Cite.** Gretton, Borgwardt, Rasch, Schölkopf, Smola, *A Kernel Two-Sample Test*, JMLR 13 (2012) 723–773. -- **Evidence.** `benchmark_repertoire_agediverge.py`: diversity↔divergence coupling biased −0.15 → unbiased −0.05. - -## 2. Aging divergence is a re-expression of the clonality decline, not an independent axis - -- **Claim (negative result worth stating).** Repertoires grow more mutually dissimilar with age at depth - (exact-overlap `−logF` vs age ρ≈0.70 at 500k, invisible at 40k — a **deep-repertoire** phenomenon in private - clonal expansions), but the component **independent of scalar diversity** is null: partial - ρ(age, divergence | ¹D) ≈ 0.07 (n.s.) across 40k/250k/500k. -- **Metric-family remark.** Sign of ρ(diversity, divergence) depends on the metric: unbiased KME ≈0 - (diverse≈naive≈shared P_gen baseline), frequency-weighted overlap **−0.68** (clonal expansion destroys - F-overlap). The naive "more diverse ⇒ less overlap" holds only for *unweighted richness* overlap; abundance - weighting flips its age-sign. Good caution against reading a single overlap number. -- **Where.** A short `\begin{remark}` after the depth-robustness proposition. Frame as: the embedding does not - add a directional aging axis beyond diversity — consistent with the "diversity for how-even" thesis. -- **Evidence.** `benchmark_repertoire_agediverge.py` (`aging` full cohort, native/500k depth). - -## 3. Batch is a first-order nuisance that cancels within-batch (new `prop:batch`, promote to a theorem+recipe) - -- **Claim.** A sequencing batch shifts a whole cohort-subset by a common offset in `Φ`; it is therefore removed - by **residualizing `Φ` on the batch indicator** (or comparing within-batch contrasts). Two consequences, - both validated on 9 real FMBA batches: a batch classifier on `Φ` is strong (OvR AUC 0.78) and collapses to - chance (0.03) after residualization; a **batch-orthogonal** biological signal (HLA) is preserved (0.60→0.61), - while a **batch-confounded** one (COVID status, some runs all-healthy) collapses to its honest within-batch - value (0.66→0.41). -- **Cookbook (add as an algorithm/COA box).** *detect* (batch one-vs-rest AUC) → *quantify* (MMD ratio - same-status-cross-batch : cross-status-within-batch; ≈1.05 here = batch as large as biology) → *correct* - (residualize `Φ` on batch, or within-batch/stratified MMD) → *verify* (batch AUC → chance; a known - batch-⟂ signal is preserved). -- **Where.** New subsection after `prop:hla`. This is the strongest new proposition — it is a clean, general - method, not COVID-specific. -- **Evidence.** `benchmark_repertoire_covidbatch.py` (airr_covid19). - -## 4. HLA imprint spans class I *and* class II, and is stronger in TRA (extend `prop:interact`/`prop:hla`) - -- **Claim.** On 4-digit-typed data the second-moment (co-occurrence) block beats diversity for HLA carriage - across loci and **both MHC classes** (15/17 alleles directional; class II 8/9; DRB1\*07:01 AUC 0.758) — a new - extension beyond the class-I β-chain literature (DeWitt 2018). On the **paired** cohort, **TRA carries the - stronger imprint** (DRB1\*15 α 0.81 vs β 0.75): TRA's lower junctional diversity makes HLA-restricted public - α clones more shared/detectable. α+β concatenation sits *between* the chains (the noisier β dilutes), so the - recommendation is **α (or α-weighted), not naive concat**, for HLA inference. -- **Where.** Extend the interaction/HLA proposition with a class-II clause and a paired-chain remark. Note the - chain asymmetry is a testable prediction of the junctional-diversity argument. -- **Evidence.** `benchmark_repertoire_covidhla.py`, `benchmark_repertoire_covidpaired.py`. - -## 5. Honest negative — a long-past exposure leaves no batch-robust bulk biomarker (a caution, state it) - -- **Claim.** Convalescent (long-resolved) SARS-CoV-2 exposure is **not** detectable from the bulk repertoire at - RNA-seq depth after batch correction: classification is chance (0.49–0.52 β/α/paired) and a from-scratch MMD - witness does not rediscover the paper's COVID clones (0.37 β / 0.45 α). The naive 0.66–0.72 was batch - confound. Detecting rare memory-phase antigen clones needs a **supervised targeted burden** of a known clone - panel (as the source paper did) or a differential control — not an unsupervised bulk contrast. -- **Where.** A `\begin{remark}` on the limits of the sample-level embedding: it reads *distributional* / - diversity structure and public-clone co-occurrence, but rare private antigen responses fall below the bulk - noise floor — deferring to the density/witness machinery (§T.6) with a proper control. -- **Evidence.** `benchmark_repertoire_covidstatus.py`, `benchmark_repertoire_covidpaired.py`. - ---- - -**One-paragraph thesis (unchanged, reinforced):** *diversity for how-even, the embedding for which-clones* — -and, added this round: **compare contrasts, not raw distances** (unbiased MMD + within-batch), the imprint is -richest in **TRA and class II**, and be willing to report **honest negatives** (aging directionality, COVID -convalescent biomarker) where the bulk signal is genuinely diversity or genuinely absent. diff --git a/ROADMAP.md b/ROADMAP.md deleted file mode 100644 index fffa630..0000000 --- a/ROADMAP.md +++ /dev/null @@ -1,91 +0,0 @@ -# mirpy ROADMAP — "vdjtools, at the embedding level" - -The goal: mirpy is the antigenomics group's **embedding / deep-learning** repertoire library — the -things vdjtools does with counts and graphs (diversity, overlap, biomarkers, generation), re-expressed -**in the language of embeddings**. This file is the durable audit + plan; per-session working notes live -in `~/.claude/plans/`. - -## The organizing idea — three verbs - -Everything the library does is one of three verbs, and the API is shaped to make that explicit: - -- **make** an embedding — `TCREmp.embed` (clonotype) · `density.fit_density_space` (cloud) · - `repertoire.sample_embedding` / `sample_descriptor` (repertoire) · `cohort.fit_donor_embeddings` - (digital donor) · `ml.SetEncoder` (learned). -- **measure across / from** embeddings — `density.neighbor_enrichment` · `repertoire.mmd_distance` / - `mmd_matrix` · `repertoire.class_witness` · `cohort.incidence_biomarkers` · `explain.channel_report`. -- **generate + decode back** — `ml.decoder.SequenceDecoder` · `repertoire.decode_metrics` · - `density.generate_background` · (planned) `generate.*`. - -**Division of labour (the house rule):** the library owns **geometry, fusion, serialization under the -prototype-hash comparability contract**; the analysis owns **which clinical/biological channels and which -statistical model**. `explain.py` ships no scorers and never sees `y`; `bench/eval.py` holds the scorers -the analysis hands it. - -## Tiers (current) - -| Tier | Module(s) | What it makes / measures | -|---|---|---| -| Clonotype | `embedding/`, `distances/`, `ml/` | prototype embedding, junction/germline distance, seq↔emb codec, Pgen regressor | -| Density (cloud) | `density.py` | TCRNET/ALICE neighbour-enrichment in embedding space, spike-in, motif denoise | -| Repertoire | `repertoire.py` | Φ(S) = kernel-mean ‖ Hill ‖ 2nd-moment; MMD; `class_witness`; descriptor | -| Explain | `explain.py` | named-channel fusion + scorer-agnostic ablation | -| **Cohort / digital donor** | `cohort.py` | multi-chain `DonorCohort`, `residualize`, `cluster_samples`, `incidence_biomarkers` | -| Scorers | `bench/eval.py` | `cv_auc` / `cv_cindex` / `km_logrank` / `kmer_matrix` (the `channel_report` closures) | -| Generative | *(planned `generate.py`)* | density over descriptors → evolve/sample → decode to clonotypes | - -## vdjtools → embedding parity - -- **DONE:** IO/schema (reuse vdjtools) · overlap→**MMD** · diversity→**Hill channel** · - TCRNET/ALICE→**density enrichment** · co-occurrence→**2nd-moment block** · Pgen→learned regressor · - biomarker-Fisher→**`cohort.incidence_biomarkers`** (delegate) · sample-clustering→**`cohort.cluster_samples`**. -- **PARTIAL:** vdjmatch annotation→nearest-epitope (benchmark only, no `annotate()`) · k-mer/physchem→named - channels (bus ready) · paired α/β (clonotype concat yes; repertoire-Φ via per-locus `DonorCohort`). -- **MISSING (planned):** clonotype-tracking→**embedding trajectory** · **preprocess** in embedding space - (downsample/batch-correct beyond `residualize`) · the **repertoire generative loop**. - -## Downstream analysis modes - -- **digital immunome / digital donor** — `cohort.fit_donor_embeddings` → `DonorCohort` (multi-chain, - hash-serialized). ✅ -- **disease / exposure classification** — `DonorCohort` + `explain.channel_report(scorer=cv_auc)`. ✅ -- **cancer prognosis & survival** — `channel_report(scorer=cv_cindex)` + `km_logrank`; the TME-aware - channels (isotype/composition/atypicality) injected by the analysis via `extra_channels`. ✅ (library); - full pan-cancer number reproduction is an analysis-repo step (needs `~/hf/airr_tcga`). -- **motif / biomarker detection** — `cohort.incidence_biomarkers` (Fisher, wins at low donor n) beside - `repertoire.class_witness` (geometry). ✅ -- **generative decode to real repertoires/clonotypes** — clonotype codec DONE; repertoire generative loop - **planned** (below). - -## Roadmap (by leverage / risk) - -- **Phase 0 — robustness + optimization quick wins. ✅ DONE** (commit `fix(phase0)`): degenerate-input - guards (unbiased-MMD singleton, empty/zero-count sample, negative-`n` prototypes, zero radius, empty - density trees); density default `backend="kdtree"`; auto-chunk + frame-sample fits; vectorized - `hla_stratified_mmd`; precomputed-witness fast-path. -- **Phase 1 — cohort tier / digital donor. ✅ DONE**: `bench/eval.py` scorers; - `repertoire.{fit_repertoire_spaces, centroid_atypicality}`; `cohort.py` - (`DonorCohort`/`fit_donor_embeddings`/`transform`/`save`/`load`, `residualize`, `cluster_samples`, - `incidence_biomarkers`). **Follow-up (analysis repo):** refactor `_tcga_embedding.build_embedding` onto - `fit_donor_embeddings` (+ `extra_channels` for isotype/composition/atypicality) and re-verify the - pan-cancer ΔC numbers (SKCM +0.039 / BLCA +0.025 / HNSC +0.022 / LGG +0.016). -- **Phase 2 — generative loop, mechanical half** (priority frontier): `CodecBundle.from_unified/from_decoder`; - `generate.py` `DescriptorDensity` (conditional Gaussian) + `evolve` + `sample`, promoting the shipped - in-silico-evolution manifold. Accept: reproduce the coupled in-silico slopes through the library. -- **Phase 3 — generative loop, research half**: RepertoireSpace-PCA ↔ codec-PCA basis bridge → - `invert_kernel_mean` (herd from real candidates + count model + germline-sampled V/J) + Pgen filter. - Accept: re-embedding the generated multiset lands near the target μ. *Fund, iterate.* -- **Phase 4 — multimodal encoders**: `modalities.isotype_fractions` / `hla_indicator` / - `fit_epitope_annotator` (nearest-epitope `annotate`) + paired-Φ recipe. (GEX / pMHC-groove deferred.) -- **Phase 5 — embedding trajectory** (research): `track.repertoire_trajectory` (Φ velocity) + - `clonotype_flux` (differential enrichment). - -## Non-goals / risks - -- **Stays in vdjtools** (delegate, never reimplement): AIRR schema/IO, `downsample`, P_gen generation, - `kmer_profile`, the `fisher_association` engine. -- **Stays analysis-local:** tissue/TME feature engineering (isotype/composition), Cox penalizer / CV scheme. -- **Comparability bites twice** in `DonorCohort`: per-locus `prototype_hash` **and** the stored identity - PCA — `load` verifies all hashes; a `residualize`d `X` is a different coordinate system. -- **Generation** can go off-manifold — herd from real candidates + Pgen filter bound (don't eliminate) the - failure rate; junction-only decoding means V/J are germline-sampled, not jointly learned. diff --git a/SQRT_D_MIGRATION.md b/SQRT_D_MIGRATION.md deleted file mode 100644 index 41fc4d4..0000000 --- a/SQRT_D_MIGRATION.md +++ /dev/null @@ -1,157 +0,0 @@ -# Migration guidelines — TCREMP coordinates from squared dissimilarity `d` to the metric `ρ = √d` - -> **DECISION (2026-07-13): NOT migrating — this plan is superseded, kept for the record.** -> The `metric="sqrt"` flag was added (commit `4e7c275`) and benchmarked: **`√d` showed no gain (a wash, -> worse on some slices)** on VDJdb retention/F1 — see the "Outcome" section below — so **v3 keeps the -> squared `d`** (default `metric="squared"`). Instead of switching, -> the theory was corrected to document `d` honestly — `(X,d)` is a *negative-type semimetric* with `ρ=√d` -> the induced metric; the coordinate map is Hölder-½ in `d` (constant `2·diam`), not 1-Lipschitz — in -> `appendix/tcremp_theory.tex` and the `junction.py`/`tcremp.py`/`bench/theory.py` docstrings. **Do not -> execute the switch below** unless a future benchmark reverses the verdict. - -_Draft plan for a separate session. Benchmark `√d` vs `d` first; switch only if it wins or ties._ -_Written 2026-07-13._ - -## Outcome (benchmarked 2026-07-13) — DO NOT SWITCH; kept `d` - -The §2 gate was prototyped and run on human TRB (VDJdb slim, arda coords): a `metric` flag added to -`junction.py` / `TCREmp` / `density` (all three blocks on one scale; `ρ=√d` is the elementwise sqrt -of the assembled squared embedding). Result — a **wash**, so `d` stays the default and the flag was -**reverted** (not shipped): - -| Metric | `d` (squared) | `√d` (sqrt) | Verdict | -|---|---|---|---| -| VDJdb mean F1 (PCA-50 preset) | 80% | 80% | tie (±1–2% jitter from seqtree threaded sums) | -| VDJdb mean F1 (best PCA dim, 100–150) | 88% | 90% | `√d` marginal, but needs ~3× the PCs | -| Retention / purity | 17% / 84% | 18% / 86% | `√d` marginal | -| **S2 `D`–`d` correlation** | **0.46** | 0.39 | **`d` clearly better** (stable across runs) | -| DBSCAN `eps` / density `r₁` | 36.6 / 9.6 | 36.6 / 9.2 | both sane | - -**Decision (per the §2 rule "if it's a wash, keep `d`"):** clustering is within noise and `√d` is -*worse* on the S2 alignment-tracking diagnostic, so the expensive switch (re-bake germline, re-fit -presets, re-train all 4 codecs, re-calibrate density, break every shipped embedding) is **not -justified**. `d` remains the published-compatible default; the `metric` flag was rolled back. It is -trivial to reinstate (sqrt of the assembled embedding, all blocks share the Gram construction — §3.2) -if the theory-cleanliness case later outweighs the wash. §5's honest `d`=Hölder‑½ / `ρ`=metric -appendix framing stands whichever way this lands. - -## 1. Why this exists - -TCREMP embeds `φ(x)_k = d(x, p_k)` where the alignment dissimilarity is the **Gram form** - -``` -d(a,b) = s(a,a) + s(b,b) − 2·s(a,b) = ‖ψ(a) − ψ(b)‖² -``` - -i.e. a **squared** Hilbert distance (Schoenberg / negative type — see appendix §T.2). A squared -Euclidean distance is **not a metric**: it violates the triangle inequality -(`ψ = 0,1,2 ⇒ d(0,2)=4 > d(0,1)+d(1,2)=2`). The genuine metric is `ρ = √d`. - -Consequences of shipping the squared `d`: - -- The appendix's clean 1‑Lipschitz / non-expansive results (`prop:lipschitz`, `prop:kuratowski` - isometry, `prop:energy` bound, `prop:samplecomplexity`, `prop:drift`, `prop:dp`) hold **verbatim for - `ρ=√d`** but only **weaken to Hölder‑½** for the squared `d` - (`D ≤ √K·d` → `D ≤ √K·2·diam·√d`; the Kuratowski isometry breaks outright). -- Euclidean distance `D` between squared-distance profiles over-weights far prototypes (the squared tail), - which **may** hurt DBSCAN / PCA. Unknown until benchmarked. - -**This is a theory-cleanliness + possible-clustering improvement, NOT a correctness fix.** The method is a -dissimilarity representation and works fine on non-metric `d`; the energy functional -`D²/K → E_p[(d(x,p)−d(y,p))²]` needs only boundedness, not the triangle inequality. - -## 2. Benchmark first — the gate (do this before touching anything else) - -Add a **coordinate-mode flag** so one pipeline computes both: - -- `mir/distances/junction.py`: `metric: Literal["squared","sqrt"] = "squared"`; when `"sqrt"`, return - `np.sqrt(np.clip(d, 0, None))`. -- Thread the flag through `TCREmp` / germline lookup so **all three blocks (V, J, junction)** are on the - same scale (see §3.2 — do not sqrt only the junction). - -Compare `d` vs `√d`, fixed seeds, per chain: - -| Metric | Where | Pass condition for switching | -|---|---|---| -| Cluster F1 / retention / purity | `mir/bench` on VDJdb (α/β/paired) | `√d` ≥ `d` (no regression) | -| `D`–`d` correlation (S2) | `mir/bench/theory.s2_dissimilarity_distance_correlation` | report both; `√d` should track alignment distance at least as well | -| DBSCAN `eps` stability | `mir/bench/metrics` (kneedle) | `√d` eps not more fragile | -| Density `r₁` calibration | `mir/density.calibrate_radius` | sane, positive; re-tune `eps_factor` if needed | - -**Decision rule:** switch only if `√d` matches or beats `d` on cluster F1/retention, or gives materially more -robust distances. If it's a wash, keep `d` (published-compatible) and just fix the appendix wording. - -## 3. Code changes (only if the benchmark says switch) - -### 3.1 Junction distance -`mir/distances/junction.py:42-51` — after `score_matrix`, clamp then sqrt: -```python -sm = np.asarray(sm, dtype=np.float32) -if metric == "sqrt": - sm = np.sqrt(np.clip(sm, 0.0, None)) # CND clamp must precede sqrt (indefinite BLOSUM ⇒ tiny negatives) -return sm -``` -- `√0 = 0` preserves the zero diagonal. The gap-block "best placement" is a **monotone** argmin, so select on - `d` (unchanged) and sqrt the chosen value — do **not** sqrt inside the placement search. -- Fix the docstring (`junction.py:3-6`): it currently calls the squared `d` a "genuine metric" — that is only - true after the sqrt. - -### 3.2 Germline distances (V / J / CDR1 / CDR2) -`mir/distances/germline.py` + the baker (`build_germline_dist.py`): the baked `.npz` blocks **must also be -`√d`** so the concatenated embedding is homogeneous across blocks. Two options: -- re-bake `√d` matrices (bump a version tag in the `.npz` / manifest), or -- sqrt at lookup time behind the same flag. -First confirm the germline blocks use the same Gram construction as the junction; if they are raw alignment -scores, reconcile the scale explicitly. - -### 3.3 Embedding + prototypes -- `mir/embedding/tcremp.py`: coordinates become `ρ`; the `embed` path is unchanged if junction + germline - return `ρ`. Assert all three interleaved slots share scale. -- `mir/resources/prototypes/manifest.json`: **bump the coordinate-system version** — `ρ` embeddings are a new, - incompatible space; consumers must check it. - -### 3.4 PCA presets -`mir/embedding/pca.py` + `presets.py`: the variance spectrum changes under `√d` → **re-fit** per-chain -`n_components` (95% / 99%) and update `CHAIN_PRESETS`. - -### 3.5 Part-2 codecs (the expensive part) -`mir/ml`: **re-train all** on `√d` embeddings — forward encoder, inverse decoder, Pgen regressor, unified -codec. Bump `CodecBundle` version + prototype hash (`bundle.py` already refuses hash mismatches — keep that -guard; old codecs become uninstallable against `ρ`, which is correct). Re-run `experiments/train_*`. - -### 3.6 Density -`mir/density.py`: `E(z)=f_obs/f_gen` is invariant to a monotone reparametrisation **in principle**, but the -balloon radius calibration `r₁` (median one-substitution drift) lives in the coordinate metric → **re-calibrate** -and re-run `experiments/benchmark_density_*`. The `ρ 0.37→−0.05` convergence and the `43%`/`1%`/`Jaccard 0.86` -numbers may shift; update THEORY.md T6. - -### 3.7 Benchmarks + docs -Re-run the bench harness and update **THEORY.md** reproduced numbers (S2 `R`, Table S1 F1/retention, T5 drift -slopes, T6 density). Update `CLAUDE.md` (coordinate-system note) and `SOURCES.md` (baked-artifact version). - -## 4. Versioning / compatibility - -- **Major coordinate-system change.** `d`-embeddings and `ρ`-embeddings are **not comparable**. Bump the - embedding schema version; the bundle's `(prototype-hash + PCA-rotation)` guard already blocks mixing — rely - on it. -- Consider keeping `metric="squared"` selectable for backward reproduction of the published-`d` results. - -## 5. Theory alignment (appendix — being done now, coordinate-agnostic) - -The appendix is being updated **now** to state the 1‑Lipschitz property correctly: it holds for the metric -`ρ=√d` (reverse triangle inequality of the genuine Hilbert metric, from negative type / Schoenberg), with the -current squared `d` documented honestly as Hölder‑½. So the manuscript is correct **whichever way the -benchmark lands**. After a switch, simply drop the squared-`d` caveats — `prop:lipschitz`, `prop:kuratowski`, -`prop:energy`, `prop:samplecomplexity`, `prop:drift`, `prop:dp` then hold verbatim. - -## 6. Checklist - -- [ ] Add `metric` flag (`junction.py`, threaded through germline + `TCREmp`) -- [ ] Benchmark `d` vs `√d`: F1/retention/purity, S2 correlation, eps stability, `r₁` — **decide** -- [ ] (if switch) sqrt junction + re-bake germline `.npz` (+ version) -- [ ] (if switch) re-fit PCA presets; bump prototype manifest version -- [ ] (if switch) re-train all `mir/ml` codecs; bump bundle version + prototype hash -- [ ] (if switch) re-calibrate density `r₁`; re-run density benchmarks -- [ ] (if switch) re-run bench harness; update THEORY.md / CLAUDE.md / SOURCES.md numbers -- [ ] Fix `junction.py` docstring "genuine metric" claim -- [ ] Drop squared-`d` caveats from the appendix once switched diff --git a/THEORY.md b/THEORY.md deleted file mode 100644 index a7d280a..0000000 --- a/THEORY.md +++ /dev/null @@ -1,322 +0,0 @@ -# TCREMP embedding — theory & experiments - -The prototype embedding rests on one claim: **Euclidean distance in embedding space -approximates pairwise alignment distance between receptors.** This note states the -propositions, points at the in-repo experiments that validate them, and records the -numbers reproduced with the v3 pipeline. The **manuscript is the source of truth for -proofs** (Kremlyakova et al., *JMB* 437 (2025) 169205, and its supplementary revision 2); -this file summarizes and links. - -Notation: `s_ij` = alignment similarity of receptors i, j; `d_ij = s_ii + s_jj − 2 s_ij` -the dissimilarity (a valid (semi)metric); `φ(i)_k = d(i, p_k)` the embedding coordinate -against prototype k; `D_ij = ‖φ(i) − φ(j)‖₂` the embedding-space Euclidean distance. - -## Propositions - -- **T1 — embedding distance ≈ alignment distance** (the core claim). With coordinates - `φ(i)_k = d(i, p_k)` over prototypes `p_k ~ P_gen`, `‖φ(i)−φ(j)‖ ≤ √K · d(i,j)` (Lipschitz - upper bound from the triangle inequality); a lower bound holds probabilistically under - prototype coverage. *Empirically:* supplementary **S2**, `Pearson(D_ij, d_ij)`. -- **T2 — optimal gap placement & substitution matrix.** The best-of-contiguous-gap-block - junction score (placements `(3,4,-4,-3)`) approximates Smith–Waterman; `d_ij ≥ 0` requires - the substitution matrix to be conditionally negative definite (seqtree's `blosum62()` is the - Gram transform `s_aa+s_bb−2s_ab`, clamped ≥ 0). Implemented via `seqtree.gapblock`. -- **T3 — PCA de-redundancy.** The V/J score blocks have rank ≤ #distinct genes ≪ K, so the - embedding lives near a low-rank subspace; `StandardScaler → PCA` recovers a compact basis - (`mir.embedding.pca`, `n_components=50`). -- **T4 — distribution laws.** `d_ij ~ Gamma` (right-skew that a Gaussian misses); nearest-in - embedding-space `D_ij ~` generalized extreme value / Fréchet (min-over-K extreme-value limit). - *Empirically:* supplementary **S1**. -- **T5 — SHM / IGH** (Part 2). Somatic hypermutation as a perturbation bounds embedding drift by - mutation load: `D_k = ‖φ(k-mutated) − φ(x)‖` is ~linear/sublinear in `k` (`mir.bench.theory. - shm_embedding_drift`; linear-R² 0.97–0.99). IGH's longer CDR3 gives the *lowest* per-mutation - drift (104 vs TRB 128) — the embedding is robust to SHM. IGH's hard *reconstruction* is instead - over-compaction: on arda coords the 95% code (95 PCs) gives exact-match 0.115, but 99% (422 PCs) - gives 0.356 (> the old 0.152; real IGH prototypes reconstruct better) — so variance retention - should be chain-adaptive, not the frame - (only 0.1% of IGH CDR3s exceed the length-40 frame). `experiments/benchmark_igh_shm.py`. -- **T6 — density space** (Part 2, `mir.density`). Enrichment as observed-density ÷ - `P_gen`-pushforward-density in embedding space, `E(z) = f_obs(z)/f_gen(z)` with `f_gen = φ_# P_gen`; - graph neighbour-enrichment (TCRNET/ALICE) is the `r→0` limit of this density ratio — the basis for - continuous, graph-free background subtraction. *Empirically* (`mir.bench.theory.tcrnet_convergence`, - real TRB repertoire vs vdjtools P_gen background, `n=1000`): the Spearman correlation between the continuous - radius-`r` neighbour count and the discrete Hamming-1 count **falls monotonically as `r` grows** past - one substitution — ρ ≈ 0.37 (½·r₁) → 0.34 (r₁) → 0.11 (2·r₁) → −0.05 (3·r₁), where r₁ is the median - one-substitution embedding drift. So the embedding test reproduces the graph test at the - one-substitution scale and generalizes it continuously beyond. The same enrichment test with a - supplied control repertoire (binomial) recovers TCRNET; with a generated `P_gen` background - (Poisson) it recovers ALICE. Benchmarks (`experiments/benchmark_density_{yfv,ankspond,tcrnet}.py`, - full repertoires from HF `isalgo/airr_{yfv19,ankspond,benchmark}`): **YFV** — day-15 vs day-0 - LLWNGPMAV(A*02) enriched hits **63 > 35** (vaccine response, ALICE regime); **AS/B27** — the public - CASSVGL[YF]STDTQYF/TRBV9/TRBJ2-3 motif appears among enriched B27+ synovial CD8 hits (**9**) and is - **absent** in B27− (0); **TCRNET vs ALICE** — a real control and a generated P_gen background agree - on the enriched clones at **Jaccard 0.86** (chance ≈ 0.30), i.e. ALICE is a special case of TCRNET. - *Lesson:* real repertoires are pervasively convergent, so P_gen enrichment flags ~40% of clones — a - **biological control** (differential) or an antigen-reference match supplies the specificity, and the - full repertoire must be used (subsampling dilutes the sparse antigen clusters). - **Full theory — appendix §T.6** (`appendix/tcremp_theory.tex`, 8 subsections T.6.1–T.6.8 with proofs): - - *Why P_gen over-flags (~40%)*: the observed law is a selection reweighting `π_obs = P_gen·Q/Z_Q` - (convergent recombination + thymic selection, "rich get richer"), so `E = q/Z_Q` and the **size-biased - bulk fold `E_obs[E] = 1 + CV²(q) > 1`** — the ~40% is null *miscalibration*, not a 40% non-null fraction. - - *Three estimators of the one ratio* `E`: **balloon** (local Poisson, exact test), **RuLSIF** (bounded - relative ratio `r_α ≤ 1/α`, valley-stable), **classifier/flow** (NCE logit `= log E + log(N/M)` — the - shared offset is *why* ALICE≡TCRNET, Jaccard 0.86); direct estimation beats plug-in at rate - `n^{-1/(2+2γ)}`. - - *Exact local test*: `μ₀ = (N−1)(n_bg+1)/(M+1)`, Poisson upper tail (ALICE) or conditional binomial - `Binom(T, N/(N+M))` (TCRNET); **ALICE is the rare-event limit of TCRNET** (`d_TV ≤ μ₀·p_bg`, Chen–Stein). - - *Water level = **Efron empirical-null median recentering*** (NOT Storey π₀): `c = max(median(n_obs)/ - median(μ₀), 1)`, robust (breakdown ½) since signal <5% barely moves the median — the DESeq median-of-ratios - size factor. Then BH under PRDS. - - *Differential control cancels convergence pointwise*: `R = f_case/f_ctrl = (1+a₁)/(1+a₀) ⊥ q` (no - water-level offset); control-absent motifs give `R→∞` (grounds B27 9-vs-0, YFV 63>35). Prefer a - biological control over P_gen whenever one exists. - - *Depth*: subsampling is Poisson thinning; detection floor `w_e^min = c/D` ⇒ process the **full** - repertoire. Ridges = density level sets → DBSCAN(ε=r₁) + per-ridge binomial epitope test. - *Parameter logic (way of action, all derived):* `λ₀∈[1,5]` default 3 (min-detectable `E* = 1 + c(q)/√λ₀`, - `c(q)=Φ⁻¹(1−q)`); `M ≥ 5N`; PCA per-chain preset (~95% var); radius `= r₁` (median one-substitution drift); - FDR `q* = 0.05`; Efron median recenter for P_gen, none for a differential control. - -- **Codec losslessness / invertibility** (Part 2, `mir.bench.theory.codec_losslessness`; appendix - §T.8 `sec:losslessness`). Three measurable levels — *geometric* (T1 distance preservation), - *informational* (`exact_ceiling = 1 − collision_rate`, decoder-independent), *reconstructive* - (decoder exact-match). On real held-out TRB the code is **injective** (collision_rate 0 ⇒ ceiling - 100% at every K/PC), so every missing exact match is decoder/data-limited, none information-limited. - Exact-match is a rate–distortion curve that saturates by ~99% var (m≈300 PCs / K≈2000; deeper is a - wash, K=10000 regresses) and is driven by *training data* — n 20k→50k→100k ⇒ 0.885→0.941→0.958, - crossing 95% on data alone, same one-shot decoder. The code is a ~10 kbit *expansion* of a ~63-bit - junction, so it is not a compressor — store the string (+ exact V/J/C) for archival recovery; the - codec earns its keep for ML/generation. Injectivity is also the linkage hazard of the privacy - section (same property, opposite sign). `experiments/benchmark_lossless_{depth,kpc,codec_losslessness}.py`. - -- **T7 — sample-level (repertoire) embedding** (v3.x, `mir.repertoire` + `mir.ml.set_encoder`; theory - **appendix §T.7** `sec:sample`). A whole repertoire is the weighted empirical measure - `ρ_S = Σ_σ w_σ δ_{φ(σ)}` on embedding space (weights `w_σ ∝ g(a_σ)`, the concave VST of T6.9). Its - fixed-vector embedding `Φ(S)` is a sketch of that measure = "the first two moments of `ψ(φ)` plus a - coverage-standardized diversity profile", three blocks each owning one requirement: - - *(a) order-invariance + (b) depth-robustness* — the **RFF kernel mean embedding** - `Φ₁ = Σ_σ w_σ ψ(φ(σ)) = μ_{ρ_S}`; converges to the population mean map at rate `n_eff^{-1/2}` with - `n_eff = (Σ w²)⁻¹`, so depth-robustness is set by the sample's spread, not its raw depth. Distance = - **MMD**. **Codebook-free** — the `K→∞` soft-assignment limit of the global-graph→cluster→histogram - recipe (VLAD/Fisher are its finite-K truncations), so there is no `K` and no clustering rule to choose. - - *(c) diversity* — coverage-standardized **Hill profile** `{⁰D,¹D,²D}` at common coverage `Ĉ*` (via - `vdjtools.stats.inext`). Depth-robustness and diversity are **mathematically antagonistic** (diversity - lives in the depth-sensitive rare tail); coverage standardization is the ecology result that reconciles - them, and `n_eff = ^qD` is itself a Hill number — one relation ties (b) to (c). `¹D` tracks the - age-related decline. Bonus: `‖μ_{ρ_S}‖²` is already an order-2 *similarity-sensitive* diversity - (Rao/Leinster–Cobbold), so a φ-aware diversity is the squared norm of the same backbone. - - *(d) HLA-linked interactions* — the compressed **second moment** `Σ_σ w_σ ψψᵀ` (codebook-free Fisher - vector), co-equal with a learned Set-Transformer/DeepRC attention head (`mir.ml`, torch). CMV clusters - are HLA-restricted, so CMV⁺ samples are close only within an **HLA-matched stratum** ⇒ use - HLA-stratified MMD. - - *(e) decoupling nuisances* — `Φ_obs = Φ_bio + ε_depth + δ_batch`: **depth** is estimation variance - `O(n_eff^{-1/2})` (handled by frequencies + coverage standardization); **batch** is a shared shift that - **cancels in a within-batch contrast** (the sample-level image of the T6 differential control) — else - residualize on batch / stratify MMD / normalize by the P_gen pushforward. Always compare *contrasts*, - not raw positions. (Variable sample length is a non-issue: the measure `ρ_S` is fixed-dimensional whatever - `|S|` is; cardinality re-enters only as richness `⁰D`.) - - *(f) RNA-seq normalization course of action* (appendix §T.7.9 `sec:samp-norm`, Table `tab:norm`) — the - raw count factorizes `a_σ ≈ g(σ)·R·θ·c·p_σ` (kit gain `g(σ)`, library size `R`, T-cell fraction/**infiltration** - `θ`, per-cell expression `c`, true frequency `p_σ`). The organizing fact: **σ-independent scalars (`R`, `θ`, `c`) - cancel under frequency normalization** (`prop:freqquotient`) — depth *and* infiltration vanish for free — while the - **σ-dependent kit gain `g(σ)` is the only surviving multiplicative nuisance**. So nuisances sort into three fates: - (1) **depth** — quotiented by frequencies, coverage-standardize diversity, never down-sample to common *depth*; - (2) **infiltration `θ`** — the tissue confound: `N = R·θ·c` conflates depth and infiltration, but `θ` is - identifiable from the **TCR read fraction** `N_TCR/N_total ∝ θ` (divides out `R`) — **carry it as an explicit - scalar channel, never normalize it away** (down-sampling erases the prognostic TIL signal); `Φ` on frequencies is - orthogonal to `θ` (`prop:infiltration`); irreducible only in the joint low-`R`, low-`θ` coverage-limited regime; - (3) **kit/read-length** — σ-dependent *shape* distortion (V/J-usage bias + long-CDR3 length-censoring, MNAR, - `prop:lengthcensor`), corrected as a **subspace projection on an anchor** (technical replicate / reference sample = - ComBat/Harmony; P_gen-predicted V/J usage = RUV negative-controls) sparing the CDR3-motif directions, or - restrict-to-common-length / inverse-recovery reweight. One line: *scalars quotiented by frequencies, shape - distortions projected out on an anchor, infiltration kept as its own channel.* **Statistical machinery** - (`rem:cmh`, refined from the BostonGene receptor-GNN whitepaper): depth acts through a per-clonotype - **detection probability** `≈1-(1-f)^{s/s̄}` (a censoring/exponent, the sample-level face of the T6 exact - point-process test) — so naive inverse-depth weighting `w∝1/s` has *no calibrated null* while the - Poisson-binomial detection model does; for covariates (age/HLA/CMV) that move marginal frequencies, **stratify - not regress** — bin by `(⌊log₁₀s⌋ × covariate band)` and use a **Cochran–Mantel–Haenszel** conditional test - (Mantel & Haenszel 1959), degrading to the exact test in small strata. Regime tiering: **blood** depth ≈ - technical (light adjustment); **tumor/tissue** depth folds in TIL fraction (10⁴-fold range → separate `θ`, - keep as channel). - *Way of action (all derived — appendix Table `tab:sample`):* frequencies not counts (scale-free); concave - `g=log1p`/Anscombe (Zipf-robust); RFF length-scale `= r₁` (one-substitution); coverage- not - depth-standardized diversity; MMD, HLA-stratified for antigen specificity. *Central use case:* low-coverage - bulk RNA-seq (all chains, `10²–10⁴` clonotypes/chain, 100–200k clinically-annotated samples). Benchmarks - (`REPERTOIRE_EMBEDDING.md`): age regression (`aging` full-depth, not `aging_lite`), CMV/HLA-stratified - (`airr_hip` = Emerson 2017), depth-robustness (`downsample` + `aging_lite`). **Build spec:** - `REPERTOIRE_EMBEDDING.md`. - *Reproduced (2026-07-14; `mir.repertoire`, `mir.ml.set_encoder`, `experiments/benchmark_repertoire_*.py`; - TRB, per-sample downsampled to the RNA-seq regime. Numbers are repeated-50-fold-CV mean±std unless noted. - **These findings were adversarially verified — two initial over-claims were caught and corrected below.**):* - - **Depth-robustness (`prop:kme`) — confirmed, but generic.** `‖Φ₁(sub)−Φ₁(full)‖` vs `n_eff` has log–log - slope **−0.55** (theory −0.5), `err·√n_eff` 1.00→0.80 across `N∈{100…10⁴}`. ⚠ *Honest reading:* this is the - **Monte-Carlo concentration rate of any weighted mean of bounded features** (`Φ₁=Σwψ`, so `err²≈Σw²·V` and - the x-axis `n_eff=(Σw²)⁻¹` is built from the same weights) — it validates the KME/MMD estimator, **not** the - TCREMB coordinate system specifically (a scrambled embedding passes too); the 20% `err·√n_eff` drift is the - only embedding-dependent content. - - **Age & CMV are clone-size (diversity) phenomena.** A coverage/Hill diversity summary dominates: age - |Spearman| **0.76** (diversity, n=79) vs 0.58 (kernel mean) vs 0.28 (k-mer); CMV AUC **0.83±0.05** - (diversity, n=240) vs 0.59–0.63 (Φ blocks) vs 0.49 (learned). For CMV this is **not an age confound** — under - decade age-matching, age-only AUC is 0.45 (chance) while diversity is 0.83: **CMV memory inflation reshapes - the clone-size distribution** (fewer, larger clones), the signal a diversity profile reads and which Φ₁'s - depth-robust concave weighting down-weights *by design*. So "the embedding adds nothing over diversity here" - is partly **by construction** (Φ₁ discards clone size) — diversity is the natural sufficient statistic for a - clone-size phenotype, not a defeat of the embedding. - - **HLA-A\*02 — clonotype identity beats diversity, established at scale (`prop:interact`).** A pure identity - signal (public A\*02-restricted clones; HLA leaves diversity unchanged): diversity is at **chance - (0.46±0.05, as predicted)**, while the **second-moment co-occurrence block separates** at **0.623±0.048** - (n=500 donors, 25k reads, `n_rff_second=256`, repeated 50-fold CV — the intervals clear each other). ⚠ - *Scale + depth + resolution were required*: an early single 70/30 split reported an inflated 0.64 (noise, - n_test≈30), the n=240 CV gave a borderline 0.535±0.08 (overlapping), and only at n=500/25k does the - second-moment CI separate cleanly from chance. The kernel-mean Φ₁ (0.58) and k-mer (0.54) also rise above - diversity, but the second moment is strongest — the HLA signal lives in clonotype **co-occurrence**, exactly - where the clone-size distribution is blind. Grounded in DeWitt et al. 2018 (*eLife* 7:e38358: TCR occurrence - patterns encode HLA on this same Emerson HIP cohort). - - **Finding motifs (`prop:witness`).** `class_witness` (`w=μ_A−μ_B`, score `s(σ)=⟨w,ψ(φ(σ))⟩`) surfaces - coherent A\*02-associated `CASS…EQYF` clones (TRBV4/6/7); the injected-motif unit test recovers a planted - public clone. On real YF data (`benchmark_repertoire_yfv.py`, day-15 vs day-0) it ranks LLWNGPMAV/A\*02 clones - at mean AUC **0.57** (vs naive fold-change 0.53) — *marginal* (n=1 sample/group per donor, ≤4 donors, no CI); - the witness is essentially a kernel-smoothed fold-change, so the **density-ratio** recovers convergent - clusters far better. - - **Spike-in recovery from RNA-seq depth (corrected)** (`benchmark_repertoire_spikein.py`; VDJdb ground truth): - plant a real VDJdb epitope's **CDR3-Hamming-selected motif family** (NLVPMVATV/A\*02 CMV, clonally expanded) - into naive P_gen backgrounds of depth `N` and recover it with `mir.density`. ⚠ *The selection metric (Hamming) - is independent of the BLOSUM-gapblock detection embedding — a cross-metric test; an earlier version that - selected the core in the same embedding it detected in was circular (72%→50% once de-circularised).* Honest - result: **recall ~35–50% at RNA-seq depth (N≤3k), FPR ~1.2%**; breadth dilutes at bulk depth but the - abundance/clonal-depth channel holds it to ~25% at N=10k. Caveat: FPR is vs a *clean* P_gen null; a real - repertoire has its own convergent clusters, so a **biological differential control** (T6) is the honest - false-positive test. Two robust lessons: antigen specificity ≠ sequence convergence (spike a real motif - family, not a diffuse epitope sample → 0% recall), and shallow depth is *favorable* for a fixed response. - - **Aging = clonality, not an independent divergence axis (unbiased MMD).** Repertoires *do* grow more - dissimilar with age at depth (exact-overlap `−logF` divergence vs age ρ **0.70** at 500k, vs 0.24 at 250k — - the signal lives in deep private expansions, invisible at the old 40k), but it is a **re-expression of the - clonality/diversity decline**, not a directional axis the embedding adds. ⚠ *The raw pairwise KME-MMD used - the biased V-statistic, whose `1/n_eff` self-term inflates distances for low-diversity (old) samples — an - artifact masquerading as signal. `mmd_matrix(unbiased=True)` removes it (diagonal-removed MMD², Gretton - 2012); the diversity↔divergence coupling drops from biased −0.15 to unbiased −0.05.* Metric-family sign - check: unbiased KME ≈0 (diverse≈naive≈shared P_gen baseline), frequency-weighted overlap **−0.68** (old - private expansions destroy F-overlap — the "more diverse ⇒ less overlap" intuition holds only for exact - overlap, and even then abundance-weighting flips its age-sign via clonal expansion). The age-divergence - *independent of scalar diversity* is n.s. across 40k/250k/500k (partial ρ(age, divergence|¹D) ≈ **0.07**, - p≈0.6). Net: deeper depth *reveals* the divergence, its content *is* the diversity decline. - (`benchmark_repertoire_agediverge.py`, `aging` full cohort at native/500k depth.) - - **Batch effects & correction (`prop:batch`) — clean validation on 9 real batches.** On `airr_covid19` - (Vlasova 2026; 9 FMBA sequencing runs), Φ strongly encodes batch (one-vs-rest AUC **0.78**) and - **residualizing Φ on the batch indicator cancels it** (→ **0.03**, chance). Natural experiment: **HLA ⟂ - batch** (donor genetics) so its signal **survives** (A\*02 0.60→0.61), while **COVID status ⟂̸ batch** (some - runs are ~all-healthy) so its naive AUC **collapses to the honest within-batch value** (0.66→0.41 - residualized, 0.54 within-mixed-batch) — the naive number rode the batch confound. MMD decomposition: - same-status-cross-batch (offset) ≈ cross-status-within-batch (biology), ratio **1.05**. **Cookbook:** - *detect* (batch OvR AUC) → *quantify* (MMD offset:biology ratio) → *correct* (residualize / within-batch - stratify) → *verify* (batch→chance, batch-⟂ signal preserved). (`benchmark_repertoire_covidbatch.py`.) - - **HLA imprint across loci and both MHC classes; TRA > TRB (`prop:interact`, DeWitt 2018).** On the - 4-digit-typed `airr_covid19`, second-moment beats diversity in direction for **15/17** class-I+II alleles, - **class II present in 8/9** (a new extension beyond class-I-only airr_hip) — **DRB1\*07:01 0.758** (Δ+0.20) - strongest. Using the **paired** cohort (α+β, 1258/1258), **TRA carries the stronger HLA imprint** (α vs β: - A\*02 0.64/0.54, B\*07 0.70/0.65, **DRB1\*15 0.81/0.75**) — TRA's lower junctional diversity makes - HLA-restricted public α clones more shared; α+β concat sits *between* the chains (noisier β dilutes), so - **use α, not paired, for HLA**. (`benchmark_repertoire_covidhla.py`, `benchmark_repertoire_covidpaired.py`.) - - **COVID-exposure biomarker — honest negative at RNA-seq depth.** A convalescent (long-past) SARS-CoV-2 - exposure leaves **no batch-robust bulk clonotype-identity signal**: COVID⁺-vs-healthy classification is naive - 0.67–0.70 but **chance after batch correction** (0.49–0.52 across β/α/paired), and the from-scratch MMD - witness does **not** rediscover the paper's COVID-associated clones (AUC **0.45 α / 0.37 β**, below chance) — - the ground truth is **87% α** (4393 vs 567), so a β-only test is doubly wrong. Paired α+β lifts the *naive* - AUC (0.67→0.70) but not the corrected one. *Reading:* the paper's biomarkers need the supervised clone panel - as a targeted burden (or deeper/larger cohorts), not an unsupervised bulk contrast — memory-phase antigen - clones are too rare in 20k-read bulk to surface without a differential control. - (`benchmark_repertoire_covidstatus.py`, `benchmark_repertoire_covidpaired.py`.) - *Lesson (verified):* two complementary regimes. **Clone-size phenotypes (age, CMV)** are diversity's turf — - the embedding discards clone size by design, so a Hill/coverage summary wins (and CMV's clonality is real - biology, not an age confound). **Clonotype-identity phenotypes (HLA)** are the embedding's turf — diversity is - at chance while the **second-moment co-occurrence block separates at scale (A\*02: 0.62 vs 0.46, n=500)**, and - the supervised witness / density recover the underlying public motifs (strongest in **TRA** and **class II** — - DRB1\*07:01 0.76, DRB1\*15 α 0.81 — a novel extension beyond class-I β work). The first-moment kernel mean's own - CI-backed value is depth-robustness (a generic KME property) + being a fixed fusion modality. Net: *diversity - for how-even, the embedding for which-clones.* Two cross-cutting cautions, both from `airr_covid19`: **(1) - batch is a first-order nuisance** — it can be as large as the biology (MMD ratio 1.05) and *inflates any - batch-confounded contrast* (naive COVID 0.66→0.41 corrected), so always detect it (OvR AUC) and compare - *within-batch / residualized* contrasts (`prop:batch`); **(2) honest negatives** — a long-past (convalescent) - antigen exposure leaves no batch-robust bulk biomarker at RNA-seq depth (COVID status = chance after - correction; the paper's clones need a supervised targeted burden, not an unsupervised bulk witness). Also: - pairwise repertoire MMD must use the **unbiased** estimator (`mmd_matrix(unbiased=True)`) when samples differ - in depth/diversity — the biased V-statistic's `1/n_eff` self-term otherwise fakes a divergence signal. - -- **2026-07-14 (pm) additions** (`BENCHMARKS.md` "2026-07-14 (pm)"; three more generalizable lessons): - **(3) the interaction block's information is directional, not spectral.** The second-moment block - `Σ w ψ₂ψ₂ᵀ` carries the HLA imprint in *which* public clones co-occur (its eigen-*vectors*), so an opt-in - top-`r` **eigenvalue** compaction (`fit_repertoire_space(n_eigs=r)`) — rotation-invariant — is lossy: HLA-A\*02 - AUC ≤0.55 vs 0.593 for the full upper triangle. Keep the upper triangle; the spectral path is a compact - fingerprint only. **(4) clone-level significance is cohort-breadth-limited.** Per-clonotype Fisher (Emerson - incidence) needs the full cohort: 0 COVID clones clear BH FDR at 150–300 donors *at any read depth*, but ~1137 - donors yield 39β/4α — so a bulk sample-level embedding cannot recover a signal that is not yet significant at - its own donor count; the tractable sample-level lever is **batch control** (mixed-batch witness β 0.51→0.75), - not per-allele HLA stratification. **(5) tissue signal is depth-gated.** On TCGA (bulk RNA-seq, IG-dominant), - the *deepest* chain IGK separates 33 tumour types at AUC 0.67 while shallow TR (~tens of clonotypes/sample) - sits at chance — but the embedding adds **no** survival value over clinical covariates (ΔC-index≈0, every - chain), a clean prognostic negative at RNA-seq depth (`benchmark_repertoire_tcga.py`). - -### Repertoire embeddings for the tumour microenvironment & survival (pan-cancer) - -The clonotype-*identity* embedding being flat for tumour survival is not the end of the paradigm — it is a -statement about *which channel* carries the signal. A repertoire embedding is a **multi-channel object**, and -different questions engage different channels: - -| channel | what it measures | question it answers | -|---|---|---| -| kernel-mean + second-moment (identity) | *which* clones / public structure | antigen specificity, HLA imprint, infection/vaccine response | -| diversity (Hill) | clonality / richness | immune ageing, exhaustion | -| coverage / infiltration | receptor read load (magnitude Φ₁ normalises away) | **tumour infiltration — hot vs cold** | -| isotype / C-gene (IGH class-switch) | plasma / mucosal humoral state | **B-cell / TLS-driven TME** | -| composition (T-vs-B, 7-chain balance) | TME cell-type mix | **immune contexture** | -| atypicality (Φ-distance to tumour-type centroid) | selection / divergence from the typical repertoire | antigen-driven remodelling | - -Making these channels explicit (`_tcga_embedding.py`) turns the "biology features" into one **TME-aware, -multi-chain repertoire embedding** Φ(S). Pan-cancer (33 TCGA types, 9 425 OS-annotated samples, -`benchmark_repertoire_tcga_pancancer.py`): over a clinical Cox (age+sex+stage+log reads), Φ is **robustly -prognostic** (LR p<0.05 *and* CV ΔC>0) in **SKCM +0.039, BLCA +0.025, HNSC +0.022, LGG +0.016** — the -immunologically active cancers — with strong power-limited positives in SARC (isotype +0.038) and KIRP -(atypicality +0.034). The most-informative channel is **coverage (infiltration) / atypicality / composition**, -essentially never clonotype identity; immune-cold small cohorts overfit (LIHC/PAAD/UCEC negative), so the -pan-cancer mean ΔC≈0 masks a clear immune-hot-cancer effect. **Lesson (6):** the same Φ(S) that fingerprints -infection and HLA in blood also stratifies the tumour microenvironment and survival in tissue — one just reads -the *non-identity* channels (infiltration magnitude + humoral class + contexture), which the depth-normalised, -single-chain kernel mean discards by construction. Unsupervised clustering of Φ recovers interpretable TME -states (hot-diverse / cold-humoral / T-skewed-unswitched), characterised by these channels, enriched in the -expected cancers, and prognostic **beyond** tumour type (stratified multivariate log-rank p=0.004; the -low-infiltration cold-humoral state carries the worst outcome, HR 1.14) — `benchmark_repertoire_tcga_tme.py`. - -**Lesson (7) — derivable metrics ⇒ a simulatable embedding.** Because the mass-preserving descriptor -(`mir.repertoire.sample_descriptor`) makes infiltration (log-mass), diversity (log n_eff) and clonality (Σw²) -**smooth analytic coordinates** — the "smoother form" of the Hill block, keeping the total mass Φ normalises -away — the cohort's joint coordinate distribution is a generative **manifold**. Perturbing the infiltration -coordinate and conditioning the rest (Gaussian manifold) is *in-silico evolution*: pan-cancer, a hotter tumour -is predicted to become more diverse (+0.84), more class-switched (+0.52) and more B-skewed (−0.63 T-vs-B), and -the CoxPH survival response is protective in 12/20 cancers but adverse in glioma (LGG Δlog-HR +1.19) — learned -couplings that recover known immunobiology (`benchmark_repertoire_tcga_insilico.py`). The embedding is thus not -only a descriptor but a *model* one can sample and evolve. - -## Reproduced numbers (v3 pipeline) - -Run `python experiments/reproduce_supplementary.py` (S1–S3) and -`python experiments/benchmark_vdjdb.py` (Table S1). `n = 3000` CDR3β from the bundled -`human_TRB` prototypes; SW = the paper's Smith-Waterman/BLOSUM62 metric, gapblock = the v3 -pipeline metric. - -Prototypes are **arda-native** (2026-07): arda-annotated real repertoires (`isalgo/airr_model_read` -functional reads → `arda rnaseq map`) over arda-baked germline distances — one IMGT allele frame -shared with query data. - -| Claim | Paper | mirpy SW | mirpy gapblock (v3) | -|---|---|---|---| -| **S2 / T1** `Pearson(D_ij, d_ij)` | 0.56 | **0.575** | 0.404 | -| **S1 / T4** `d_ij` law | Gamma > Normal | Gamma ≈ Normal (AIC tie) | Gamma ≈ Normal (AIC tie) | -| **S1 / T4** `D_ij` law | GEV/Fréchet ≫ Normal, ξ=+0.11 | **GEV wins** (KS .021 vs .078), ξ≈−0.03 | **GEV wins** (KS .017 vs .082), ξ≈0.00 | -| **S3** real vs model prototypes | 0.96 | — | **0.940** | -| **Table S1** VDJdb TRB antigen clustering | mean F1 91%, retention 18% | — | mean F1 81%, retention 17% | - -The core claims reproduce on the arda-native coordinate system: the embedding distance tracks -alignment distance (T1: SW 0.575 ≈ paper 0.56), `D_ij` is extreme-value distributed (T4), and the -prototype *source* barely matters (S3, 0.940 ≈ 0.96). `d_ij` is now a Gamma/Normal near-tie (the -real-repertoire junction-length distribution is tighter than the previous set). The v3 embedding is -a deliberately new, versioned coordinate system; switching prototypes to arda-annotated real -repertoires (from the earlier real/OLGA mix) shifts the exact numbers but preserves every law. diff --git a/docs/conf.py b/docs/conf.py index 84fec2f..542b635 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -4,8 +4,8 @@ import sys from importlib.metadata import PackageNotFoundError, version as _pkg_version -# mir is installed (pip install -e .), but add the repo root so autodoc works from a bare checkout. -sys.path.insert(0, os.path.abspath("..")) +# mir is installed (pip install -e .), but add the src dir so autodoc works from a bare checkout. +sys.path.insert(0, os.path.abspath("../src")) project = "mirpy" author = "ISALGO laboratory" diff --git a/docs/examples.rst b/docs/examples.rst index 6d95430..48fd486 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -1,14 +1,14 @@ Examples ======== -mirpy ships three runnable `marimo `_ notebooks under ``notebooks/``. They are +mirpy ships three runnable `marimo `_ notebooks under ``examples/``. They are self-contained (they run on the bundled prototypes / test assets — no downloads) and double as living documentation for the three tiers. Install the extra and open one: .. code-block:: bash pip install "mirpy-lib[examples]" # marimo, matplotlib, umap-learn - marimo edit notebooks/quickstart.py # interactive; or `marimo run …` for read-only + marimo edit examples/quickstart.py # interactive; or `marimo run …` for read-only .. list-table:: :header-rows: 1 @@ -16,13 +16,13 @@ living documentation for the three tiers. Install the extra and open one: * - Notebook - What it shows - * - ``notebooks/quickstart.py`` + * - ``examples/quickstart.py`` - Clonotype embedding end-to-end: ``TCREmp.embed`` a VDJdb-style set, PCA-denoise, cluster antigen-specific TCRs, and a UMAP coloured by epitope. - * - ``notebooks/density.py`` + * - ``examples/density.py`` - Background subtraction (Theory T6): fit a density space, run balloon ``neighbor_enrichment`` against a P_gen / control background, and pull out the enriched convergent family. - * - ``notebooks/theory.py`` + * - ``examples/theory.py`` - Reproduces the supplementary results S1–S3 on bundled data — the distance laws (Gamma / extreme-value), the D↔d correlation, and prototype-source robustness. diff --git a/docs/index.rst b/docs/index.rst index 32068bd..28fb46a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -29,8 +29,9 @@ sequences. The classical v1.x/v2 toolkit is frozen on the ``legacy-v2`` branch ( New here? Start with the :doc:`user guide ` (install → a real result, one task per section), skim the runnable :doc:`examples `, and reach for the :doc:`API reference ` for -signatures. The mathematical theory (T1–T7) lives in ``THEORY.md`` and recorded benchmark numbers -in ``BENCHMARKS.md``. +signatures. The mathematical theory (T1–T7) and recorded benchmark numbers live in the companion +`2026-mirpy-analysis `_ repo (``benchmarks/{THEORY,BENCHMARKS}.md``), +with the LaTeX theory appendix in ``2026-mirpy-ms``. 30-second quickstart -------------------- diff --git a/docs/usage.rst b/docs/usage.rst index 1e0ae81..48d3264 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -87,7 +87,7 @@ scores per-antigen F1 / retention against known labels (e.g. a VDJdb dump): metrics = cluster_metrics(labels, epitopes) # {epitope: AntigenMetric} DBSCAN (default) is tightest/purest; HDBSCAN trades precision for ~3× coverage on variable-density -data. See ``notebooks/quickstart.py`` for the end-to-end VDJdb example with a UMAP. +data. See ``examples/quickstart.py`` for the end-to-end VDJdb example with a UMAP. Density and background subtraction ---------------------------------- @@ -109,7 +109,7 @@ adaptive-bandwidth balloon estimator with a Poisson/binomial test and BH q-value Prefer a biological control (e.g. pre/post-vaccination) over the P_gen background — differential enrichment cancels generic public convergence and isolates the antigen-specific response. At whole-repertoire scale, pass ``backend="kdtree"`` (exact, multicore) or ``backend="ann"`` -(approximate, ``[ann]`` extra). See ``notebooks/density.py``. +(approximate, ``[ann]`` extra). See ``examples/density.py``. Repertoire embedding Φ(S) + MMD ------------------------------- @@ -261,7 +261,7 @@ Benchmark harness & reproducing the paper ``mir.bench`` provides the VDJdb clustering benchmark (F1 / retention / purity) and the reproduced theory experiments (S1–S3, T5–T6, codec losslessness). The self-contained theory runs on bundled -data via ``notebooks/theory.py``; the full benchmark suite lives in the ``2026-mirpy-analysis`` repo. +data via ``examples/theory.py``; the full benchmark suite lives in the ``2026-mirpy-analysis`` repo. .. code-block:: python diff --git a/notebooks/__marimo__/session/quickstart.py.json b/examples/__marimo__/session/quickstart.py.json similarity index 100% rename from notebooks/__marimo__/session/quickstart.py.json rename to examples/__marimo__/session/quickstart.py.json diff --git a/notebooks/__marimo__/session/theory.py.json b/examples/__marimo__/session/theory.py.json similarity index 100% rename from notebooks/__marimo__/session/theory.py.json rename to examples/__marimo__/session/theory.py.json diff --git a/notebooks/density.py b/examples/density.py similarity index 100% rename from notebooks/density.py rename to examples/density.py diff --git a/notebooks/quickstart.py b/examples/quickstart.py similarity index 100% rename from notebooks/quickstart.py rename to examples/quickstart.py diff --git a/notebooks/theory.py b/examples/theory.py similarity index 100% rename from notebooks/theory.py rename to examples/theory.py diff --git a/pyproject.toml b/pyproject.toml index c532642..c2f5430 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -73,7 +73,7 @@ Documentation = "https://antigenomics.github.io/mirpy" Issues = "https://github.com/antigenomics/mirpy/issues" [tool.hatch.build.targets.wheel] -packages = ["mir"] +packages = ["src/mir"] [tool.hatch.build.targets.sdist] exclude = [ @@ -81,7 +81,7 @@ exclude = [ ".idea/**", ".vscode/**", "docs/_build/**", - "notebooks/**", + "examples/**", "benchmarks/**", "tests/**", "tmp/**", @@ -99,12 +99,12 @@ markers = [ # Coverage is opt-in (not in addopts) so the default fast suite stays fast. # Measure with: pytest -m "not benchmark and not integration" --cov=mir --cov-report=term-missing [tool.coverage.run] -source = ["mir"] +source = ["src/mir"] branch = true omit = [ - "mir/resources/gene_library/build_gene_library.py", - "mir/resources/gene_library/build_region_annotations.py", - "mir/resources/prototypes/generate_prototypes.py", + "src/mir/resources/gene_library/build_gene_library.py", + "src/mir/resources/gene_library/build_region_annotations.py", + "src/mir/resources/prototypes/generate_prototypes.py", "*/__init__.py", ] diff --git a/mir/__init__.py b/src/mir/__init__.py similarity index 100% rename from mir/__init__.py rename to src/mir/__init__.py diff --git a/mir/aliases.py b/src/mir/aliases.py similarity index 100% rename from mir/aliases.py rename to src/mir/aliases.py diff --git a/mir/alleles.py b/src/mir/alleles.py similarity index 100% rename from mir/alleles.py rename to src/mir/alleles.py diff --git a/mir/bench/__init__.py b/src/mir/bench/__init__.py similarity index 100% rename from mir/bench/__init__.py rename to src/mir/bench/__init__.py diff --git a/mir/bench/eval.py b/src/mir/bench/eval.py similarity index 100% rename from mir/bench/eval.py rename to src/mir/bench/eval.py diff --git a/mir/bench/metrics.py b/src/mir/bench/metrics.py similarity index 100% rename from mir/bench/metrics.py rename to src/mir/bench/metrics.py diff --git a/mir/bench/theory.py b/src/mir/bench/theory.py similarity index 100% rename from mir/bench/theory.py rename to src/mir/bench/theory.py diff --git a/mir/bench/vdjdb.py b/src/mir/bench/vdjdb.py similarity index 100% rename from mir/bench/vdjdb.py rename to src/mir/bench/vdjdb.py diff --git a/mir/cli.py b/src/mir/cli.py similarity index 100% rename from mir/cli.py rename to src/mir/cli.py diff --git a/mir/cohort.py b/src/mir/cohort.py similarity index 100% rename from mir/cohort.py rename to src/mir/cohort.py diff --git a/mir/density.py b/src/mir/density.py similarity index 100% rename from mir/density.py rename to src/mir/density.py diff --git a/mir/distances/__init__.py b/src/mir/distances/__init__.py similarity index 100% rename from mir/distances/__init__.py rename to src/mir/distances/__init__.py diff --git a/mir/distances/germline.py b/src/mir/distances/germline.py similarity index 100% rename from mir/distances/germline.py rename to src/mir/distances/germline.py diff --git a/mir/distances/junction.py b/src/mir/distances/junction.py similarity index 100% rename from mir/distances/junction.py rename to src/mir/distances/junction.py diff --git a/mir/embedding/__init__.py b/src/mir/embedding/__init__.py similarity index 100% rename from mir/embedding/__init__.py rename to src/mir/embedding/__init__.py diff --git a/mir/embedding/pca.py b/src/mir/embedding/pca.py similarity index 100% rename from mir/embedding/pca.py rename to src/mir/embedding/pca.py diff --git a/mir/embedding/presets.py b/src/mir/embedding/presets.py similarity index 100% rename from mir/embedding/presets.py rename to src/mir/embedding/presets.py diff --git a/mir/embedding/prototypes.py b/src/mir/embedding/prototypes.py similarity index 100% rename from mir/embedding/prototypes.py rename to src/mir/embedding/prototypes.py diff --git a/mir/embedding/tcremp.py b/src/mir/embedding/tcremp.py similarity index 100% rename from mir/embedding/tcremp.py rename to src/mir/embedding/tcremp.py diff --git a/mir/explain.py b/src/mir/explain.py similarity index 100% rename from mir/explain.py rename to src/mir/explain.py diff --git a/mir/ml/__init__.py b/src/mir/ml/__init__.py similarity index 100% rename from mir/ml/__init__.py rename to src/mir/ml/__init__.py diff --git a/mir/ml/bundle.py b/src/mir/ml/bundle.py similarity index 100% rename from mir/ml/bundle.py rename to src/mir/ml/bundle.py diff --git a/mir/ml/codec.py b/src/mir/ml/codec.py similarity index 100% rename from mir/ml/codec.py rename to src/mir/ml/codec.py diff --git a/mir/ml/decoder.py b/src/mir/ml/decoder.py similarity index 100% rename from mir/ml/decoder.py rename to src/mir/ml/decoder.py diff --git a/mir/ml/encoder.py b/src/mir/ml/encoder.py similarity index 100% rename from mir/ml/encoder.py rename to src/mir/ml/encoder.py diff --git a/mir/ml/set_encoder.py b/src/mir/ml/set_encoder.py similarity index 100% rename from mir/ml/set_encoder.py rename to src/mir/ml/set_encoder.py diff --git a/mir/ml/tokenize.py b/src/mir/ml/tokenize.py similarity index 100% rename from mir/ml/tokenize.py rename to src/mir/ml/tokenize.py diff --git a/mir/ml/train.py b/src/mir/ml/train.py similarity index 100% rename from mir/ml/train.py rename to src/mir/ml/train.py diff --git a/mir/py.typed b/src/mir/py.typed similarity index 100% rename from mir/py.typed rename to src/mir/py.typed diff --git a/mir/repertoire.py b/src/mir/repertoire.py similarity index 100% rename from mir/repertoire.py rename to src/mir/repertoire.py diff --git a/mir/resources/__init__.py b/src/mir/resources/__init__.py similarity index 100% rename from mir/resources/__init__.py rename to src/mir/resources/__init__.py diff --git a/mir/resources/gene_library/README.md b/src/mir/resources/gene_library/README.md similarity index 100% rename from mir/resources/gene_library/README.md rename to src/mir/resources/gene_library/README.md diff --git a/mir/resources/gene_library/__init__.py b/src/mir/resources/gene_library/__init__.py similarity index 100% rename from mir/resources/gene_library/__init__.py rename to src/mir/resources/gene_library/__init__.py diff --git a/mir/resources/gene_library/build_gene_library.py b/src/mir/resources/gene_library/build_gene_library.py similarity index 100% rename from mir/resources/gene_library/build_gene_library.py rename to src/mir/resources/gene_library/build_gene_library.py diff --git a/mir/resources/gene_library/build_region_annotations.py b/src/mir/resources/gene_library/build_region_annotations.py similarity index 100% rename from mir/resources/gene_library/build_region_annotations.py rename to src/mir/resources/gene_library/build_region_annotations.py diff --git a/mir/resources/gene_library/imgt_gene_library.txt b/src/mir/resources/gene_library/imgt_gene_library.txt similarity index 100% rename from mir/resources/gene_library/imgt_gene_library.txt rename to src/mir/resources/gene_library/imgt_gene_library.txt diff --git a/mir/resources/gene_library/olga_gene_library.txt b/src/mir/resources/gene_library/olga_gene_library.txt similarity index 100% rename from mir/resources/gene_library/olga_gene_library.txt rename to src/mir/resources/gene_library/olga_gene_library.txt diff --git a/mir/resources/gene_library/region_annotations.txt b/src/mir/resources/gene_library/region_annotations.txt similarity index 100% rename from mir/resources/gene_library/region_annotations.txt rename to src/mir/resources/gene_library/region_annotations.txt diff --git a/mir/resources/germline_dist/build_germline_dist.py b/src/mir/resources/germline_dist/build_germline_dist.py similarity index 100% rename from mir/resources/germline_dist/build_germline_dist.py rename to src/mir/resources/germline_dist/build_germline_dist.py diff --git a/mir/resources/germline_dist/human_IGH.npz b/src/mir/resources/germline_dist/human_IGH.npz similarity index 100% rename from mir/resources/germline_dist/human_IGH.npz rename to src/mir/resources/germline_dist/human_IGH.npz diff --git a/mir/resources/germline_dist/human_IGK.npz b/src/mir/resources/germline_dist/human_IGK.npz similarity index 100% rename from mir/resources/germline_dist/human_IGK.npz rename to src/mir/resources/germline_dist/human_IGK.npz diff --git a/mir/resources/germline_dist/human_IGL.npz b/src/mir/resources/germline_dist/human_IGL.npz similarity index 100% rename from mir/resources/germline_dist/human_IGL.npz rename to src/mir/resources/germline_dist/human_IGL.npz diff --git a/mir/resources/germline_dist/human_TRA.npz b/src/mir/resources/germline_dist/human_TRA.npz similarity index 100% rename from mir/resources/germline_dist/human_TRA.npz rename to src/mir/resources/germline_dist/human_TRA.npz diff --git a/mir/resources/germline_dist/human_TRB.npz b/src/mir/resources/germline_dist/human_TRB.npz similarity index 100% rename from mir/resources/germline_dist/human_TRB.npz rename to src/mir/resources/germline_dist/human_TRB.npz diff --git a/mir/resources/germline_dist/human_TRD.npz b/src/mir/resources/germline_dist/human_TRD.npz similarity index 100% rename from mir/resources/germline_dist/human_TRD.npz rename to src/mir/resources/germline_dist/human_TRD.npz diff --git a/mir/resources/germline_dist/human_TRG.npz b/src/mir/resources/germline_dist/human_TRG.npz similarity index 100% rename from mir/resources/germline_dist/human_TRG.npz rename to src/mir/resources/germline_dist/human_TRG.npz diff --git a/mir/resources/germline_dist/mouse_IGH.npz b/src/mir/resources/germline_dist/mouse_IGH.npz similarity index 100% rename from mir/resources/germline_dist/mouse_IGH.npz rename to src/mir/resources/germline_dist/mouse_IGH.npz diff --git a/mir/resources/germline_dist/mouse_IGK.npz b/src/mir/resources/germline_dist/mouse_IGK.npz similarity index 100% rename from mir/resources/germline_dist/mouse_IGK.npz rename to src/mir/resources/germline_dist/mouse_IGK.npz diff --git a/mir/resources/germline_dist/mouse_IGL.npz b/src/mir/resources/germline_dist/mouse_IGL.npz similarity index 100% rename from mir/resources/germline_dist/mouse_IGL.npz rename to src/mir/resources/germline_dist/mouse_IGL.npz diff --git a/mir/resources/germline_dist/mouse_TRA.npz b/src/mir/resources/germline_dist/mouse_TRA.npz similarity index 100% rename from mir/resources/germline_dist/mouse_TRA.npz rename to src/mir/resources/germline_dist/mouse_TRA.npz diff --git a/mir/resources/germline_dist/mouse_TRB.npz b/src/mir/resources/germline_dist/mouse_TRB.npz similarity index 100% rename from mir/resources/germline_dist/mouse_TRB.npz rename to src/mir/resources/germline_dist/mouse_TRB.npz diff --git a/mir/resources/germline_dist/mouse_TRD.npz b/src/mir/resources/germline_dist/mouse_TRD.npz similarity index 100% rename from mir/resources/germline_dist/mouse_TRD.npz rename to src/mir/resources/germline_dist/mouse_TRD.npz diff --git a/mir/resources/germline_dist/mouse_TRG.npz b/src/mir/resources/germline_dist/mouse_TRG.npz similarity index 100% rename from mir/resources/germline_dist/mouse_TRG.npz rename to src/mir/resources/germline_dist/mouse_TRG.npz diff --git a/mir/resources/prototypes/.gitkeep b/src/mir/resources/prototypes/.gitkeep similarity index 100% rename from mir/resources/prototypes/.gitkeep rename to src/mir/resources/prototypes/.gitkeep diff --git a/mir/resources/prototypes/generate_prototypes.py b/src/mir/resources/prototypes/generate_prototypes.py similarity index 100% rename from mir/resources/prototypes/generate_prototypes.py rename to src/mir/resources/prototypes/generate_prototypes.py diff --git a/mir/resources/prototypes/human_IGH.tsv b/src/mir/resources/prototypes/human_IGH.tsv similarity index 100% rename from mir/resources/prototypes/human_IGH.tsv rename to src/mir/resources/prototypes/human_IGH.tsv diff --git a/mir/resources/prototypes/human_IGK.tsv b/src/mir/resources/prototypes/human_IGK.tsv similarity index 100% rename from mir/resources/prototypes/human_IGK.tsv rename to src/mir/resources/prototypes/human_IGK.tsv diff --git a/mir/resources/prototypes/human_IGL.tsv b/src/mir/resources/prototypes/human_IGL.tsv similarity index 100% rename from mir/resources/prototypes/human_IGL.tsv rename to src/mir/resources/prototypes/human_IGL.tsv diff --git a/mir/resources/prototypes/human_TRA.tsv b/src/mir/resources/prototypes/human_TRA.tsv similarity index 100% rename from mir/resources/prototypes/human_TRA.tsv rename to src/mir/resources/prototypes/human_TRA.tsv diff --git a/mir/resources/prototypes/human_TRB.tsv b/src/mir/resources/prototypes/human_TRB.tsv similarity index 100% rename from mir/resources/prototypes/human_TRB.tsv rename to src/mir/resources/prototypes/human_TRB.tsv diff --git a/mir/resources/prototypes/human_TRD.tsv b/src/mir/resources/prototypes/human_TRD.tsv similarity index 100% rename from mir/resources/prototypes/human_TRD.tsv rename to src/mir/resources/prototypes/human_TRD.tsv diff --git a/mir/resources/prototypes/human_TRG.tsv b/src/mir/resources/prototypes/human_TRG.tsv similarity index 100% rename from mir/resources/prototypes/human_TRG.tsv rename to src/mir/resources/prototypes/human_TRG.tsv diff --git a/mir/resources/prototypes/manifest.json b/src/mir/resources/prototypes/manifest.json similarity index 100% rename from mir/resources/prototypes/manifest.json rename to src/mir/resources/prototypes/manifest.json diff --git a/mir/resources/prototypes/mouse_TRA.tsv b/src/mir/resources/prototypes/mouse_TRA.tsv similarity index 100% rename from mir/resources/prototypes/mouse_TRA.tsv rename to src/mir/resources/prototypes/mouse_TRA.tsv diff --git a/mir/resources/prototypes/mouse_TRB.tsv b/src/mir/resources/prototypes/mouse_TRB.tsv similarity index 100% rename from mir/resources/prototypes/mouse_TRB.tsv rename to src/mir/resources/prototypes/mouse_TRB.tsv