{ "cells": [ { "cell_type": "markdown", "id": "646944b2", "metadata": { "language": "markdown" }, "source": [ "# Precursor frequency for an epitope\n", "\n", "How much of a repertoire can see a given pMHC? The estimand is\n", "\n", "$$F(e) = \\sum_{\\tau \\in C_e} \\pi(\\tau)$$\n", "\n", "— the probability that a random naive-repertoire junction recognises epitope $e$. This is the\n", "continuous quantity behind the word \"immunogenic\".\n", "\n", "This notebook walks the three things that make it harder than summing `Pgen` over the TCRs a\n", "database happens to hold:\n", "\n", "1. **The values span orders of magnitude.** A mean is meaningless; the sum is set by a handful of\n", " public clonotypes.\n", "2. **The neighbourhoods overlap.** Cognate junctions are near-duplicates by construction, so adding\n", " their per-sequence ball masses double-counts the shared region. `union_mass` is exact and does\n", " not enumerate.\n", "3. **Most of the cognate set was never catalogued.** Two answers, because they are different\n", " questions: a finite census of the neighbourhood, and a Horvitz–Thompson extrapolation whose\n", " *mass* converges but whose *count* does not.\n", "\n", "Needs the optional extra:\n", "\n", "```\n", "pip install 'vdjmatch[precursor]'\n", "```" ] }, { "cell_type": "code", "execution_count": 1, "id": "e95ba646", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:26:50.246520Z", "iopub.status.busy": "2026-08-17T15:26:50.244714Z", "iopub.status.idle": "2026-08-17T15:26:50.389048Z", "shell.execute_reply": "2026-08-17T15:26:50.371947Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "vdjmatch 0.1.2\n", "vdjtools 3.9.3\n", "seqtree 0.7.0\n", "polars 1.41.2\n" ] } ], "source": [ "# Environment: record the versions every number below depends on.\n", "from importlib.metadata import version\n", "\n", "import polars as pl\n", "\n", "import vdjmatch\n", "from vdjmatch import db, precursor as P\n", "\n", "SEED = 0 # nothing here samples, but the VDJdb release pin is the reproducibility knob\n", "VDJDB_TAG = \"2026-06-11-ZENODO\"\n", "for pkg in (\"vdjmatch\", \"vdjtools\", \"seqtree\", \"polars\"):\n", " print(f\"{pkg:10s} {version(pkg)}\")" ] }, { "cell_type": "markdown", "id": "5ac7d515", "metadata": { "language": "markdown" }, "source": [ "## 1. Load a cognate set\n", "\n", "**CDR3 is not junction.** VDJdb's column is *named* `cdr3` but holds junctions — Cys104 and\n", "Phe118/Trp included — which is what the recombination model wants. An anchor-stripped IMGT CDR3\n", "scores exactly `0.0` with no error, so `check_junctions` exists to make that failure loud." ] }, { "cell_type": "code", "execution_count": 2, "id": "ef43eff7", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:26:50.458475Z", "iopub.status.busy": "2026-08-17T15:26:50.456647Z", "iopub.status.idle": "2026-08-17T15:26:53.849764Z", "shell.execute_reply": "2026-08-17T15:26:53.816078Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GILGFVFTL: 6,629 records, 6,627 distinct junctions, 2 failed the anchor check\n" ] } ], "source": [ "# One well-sampled HLA-A*02:01 epitope: influenza A M1 58-66.\n", "EPITOPE = \"GILGFVFTL\"\n", "\n", "vdj = db.load(asset=\"slim\", species=\"HomoSapiens\", pin=VDJDB_TAG)\n", "beta = vdj.filter((pl.col(\"gene\") == \"TRB\") & (pl.col(\"epitope\") == EPITOPE))\n", "\n", "junctions, suspect = P.check_junctions(beta[\"cdr3\"].unique().to_list())\n", "junctions = list(dict.fromkeys(junctions))\n", "print(f\"{EPITOPE}: {beta.height:,} records, {len(junctions):,} distinct junctions, \"\n", " f\"{len(suspect)} failed the anchor check\")" ] }, { "cell_type": "markdown", "id": "3b415557", "metadata": { "language": "markdown" }, "source": [ "## 2. The spread is the point\n", "\n", "`observed_mass` is the sum of `Pgen` over the recorded junctions — a **strict lower bound** on\n", "$F(e)$, and a biased one: a TCR enters a specificity database roughly in proportion to its\n", "repertoire frequency, so the recorded members are systematically the high-`Pgen` ones." ] }, { "cell_type": "code", "execution_count": 3, "id": "3377919c", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:26:53.875069Z", "iopub.status.busy": "2026-08-17T15:26:53.874681Z", "iopub.status.idle": "2026-08-17T15:26:56.700058Z", "shell.execute_reply": "2026-08-17T15:26:56.662042Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "observed_mass 3.305e-04\n", "median Pgen 1.724e-09\n", "spread 30.0 orders of magnitude\n", "top-1 share of sum 0.010\n", "top-10 share of sum 0.080\n", "Pgen == 0 under model 3\n" ] } ], "source": [ "# Per-junction Pgen, then how far apart the values are and how concentrated their sum is.\n", "import numpy as np\n", "\n", "model = P.load_model(\"TRB\")\n", "pg = np.array(P.pgen(model, junctions))\n", "nz = np.sort(pg[pg > 0])[::-1]\n", "\n", "print(f\"observed_mass {P.observed_mass(model, junctions):.3e}\")\n", "print(f\"median Pgen {np.median(nz):.3e}\")\n", "print(f\"spread {np.log10(nz[0] / nz[-1]):.1f} orders of magnitude\")\n", "print(f\"top-1 share of sum {nz[0] / nz.sum():.3f}\")\n", "print(f\"top-10 share of sum {nz[:10].sum() / nz.sum():.3f}\")\n", "print(f\"Pgen == 0 under model {(pg == 0).sum()}\")" ] }, { "cell_type": "markdown", "id": "ab5ba53d", "metadata": { "language": "markdown" }, "source": [ "## 3. The union, not the sum\n", "\n", "One substitution is part of the estimator, not a tuned parameter: at radius 0 a repertoire is too\n", "sparse for either route to be estimable. But two cognate junctions one substitution apart have balls\n", "sharing 20 sequences, and adding their masses counts that region twice.\n", "\n", "`union_mass` is exact without enumerating the union. Every member is counted $\\mathrm{cov}(x)$ times\n", "by the naive sum, so\n", "\n", "$$m\\Big(\\bigcup_a B_r(a)\\Big) = \\sum_a m\\big(B_r(a)\\big) - \\sum_{x\\,:\\,\\mathrm{cov}(x)\\ge 2} \\big(\\mathrm{cov}(x)-1\\big)\\,P_{\\mathrm{gen}}(x)$$\n", "\n", "with no inclusion–exclusion and hence no truncation error. Only centres inside one connected\n", "component of the $2r$ graph can contribute, so singleton components cost nothing." ] }, { "cell_type": "code", "execution_count": 4, "id": "2d4bea67", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:26:56.720060Z", "iopub.status.busy": "2026-08-17T15:26:56.719798Z", "iopub.status.idle": "2026-08-17T15:27:36.374820Z", "shell.execute_reply": "2026-08-17T15:27:36.371019Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "naive per-sequence sum 1.257e-02\n", "exact union 7.432e-03\n", "double-counting avoided 40.9%\n", "components 2,968 (3,897 of 6,627 junctions have a close neighbour)\n", "Pgen calls for the correction: 42,928 (the union itself holds 1,709,201) [39.5s]\n" ] } ], "source": [ "# The exact union, and the double-counting a naive sum would have invented.\n", "import time\n", "\n", "t0 = time.perf_counter()\n", "u = P.union_mass(model, junctions, r=1)\n", "dt = time.perf_counter() - t0\n", "\n", "print(f\"naive per-sequence sum {u['naive_sum']:.3e}\")\n", "print(f\"exact union {u['union']:.3e}\")\n", "print(f\"double-counting avoided {u['overlap']:.1%}\")\n", "print(f\"components {u['n_components']:,} \"\n", " f\"({u['n_clustered']:,} of {u['n_seqs']:,} junctions have a close neighbour)\")\n", "print(f\"Pgen calls for the correction: {u['n_multiply_covered']:,} \"\n", " f\"(the union itself holds {u['n_union']:,}) [{dt:.1f}s]\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "db6ea4c3", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:27:36.383236Z", "iopub.status.busy": "2026-08-17T15:27:36.382925Z", "iopub.status.idle": "2026-08-17T15:27:44.360154Z", "shell.execute_reply": "2026-08-17T15:27:44.359583Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "union_mass 2.278723889e-04\n", "ball_mass 2.278723889e-04\n", "relative 1.19e-16\n" ] } ], "source": [ "# Regression against the oracle: enumerate the union and score every member. Slow, and exact.\n", "small = junctions[:150]\n", "fast, oracle = P.union_mass(model, small, r=1), P.ball_mass(model, small, r=1)\n", "rel = abs(fast[\"union\"] - oracle[\"union\"]) / oracle[\"union\"]\n", "print(f\"union_mass {fast['union']:.9e}\\nball_mass {oracle['union']:.9e}\\nrelative {rel:.2e}\")" ] }, { "cell_type": "markdown", "id": "c7c94b03", "metadata": { "language": "markdown" }, "source": [ "## 4. The ball is a smoother, not a coverage correction\n", "\n", "It is tempting to read the neighbourhood mass as \"the observed mass plus what sampling missed\". It\n", "is not. The ball is a statement about **cognacy**: a junction one substitution from a cognate TCR is\n", "itself cognate with some probability, measured by Mayer & Callan (*PNAS* 2023;120:e2213264120) to\n", "fall about ten-fold per unit of edit distance. `shell_profile` applies that per shell.\n", "\n", "Shells are obtained by *differencing unions*, which is exact because min-distance shells partition\n", "the ball — so the shell masses inherit `union_mass`'s exactness and nothing is enumerated to get\n", "them. What still has to be enumerated is the multiply-covered set, and only within a connected\n", "component. **That is where the real ceiling is**, and on a well-sampled epitope at radius 2 it is\n", "reached: this cognate set has one component of a couple of thousand junctions whose radius-2 union\n", "runs to tens of millions of sequences. The estimator refuses loudly rather than thrashing, and says\n", "what to do about it." ] }, { "cell_type": "code", "execution_count": 6, "id": "462b6765", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:27:44.362745Z", "iopub.status.busy": "2026-08-17T15:27:44.362381Z", "iopub.status.idle": "2026-08-17T15:30:51.542376Z", "shell.execute_reply": "2026-08-17T15:30:51.538535Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " shell r=0 weight 1.000 n 6,627 mass 3.305e-04\n", " shell r=1 weight 0.100 n 1,702,574 mass 7.101e-03\n", "\n", "raw union (alpha=1) 7.432e-03\n", "retained (alpha=0.1) 1.041e-03\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "observed bound 3.305e-04\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "r=2 on all 6,627 junctions ->\n", " one connected component of 1199 junctions has a radius-2 union of 43,726,919 sequences (~8.3 GB as Python strings), above max_members=20,000,000. Raise max_members if you have the memory, split this group, or lower r; components are independent, so their union masses add exactly and splitting by component loses nothing.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "radius-2 profile on 200 of them:\n", " shell r=0 weight 1.000 n 200 mass 8.251e-06\n", " shell r=1 weight 0.100 n 53,380 mass 3.327e-04\n", " shell r=2 weight 0.010 n 6,707,831 mass 4.134e-03\n", " retained 8.287e-05 against a raw union of 4.475e-03\n" ] } ], "source": [ "# Shell-resolved mass with cognacy retention alpha = 0.1 per edit, at the recommended radius.\n", "prof = P.shell_profile(model, junctions, r=1, alpha=0.1)\n", "for s in prof[\"shells\"]:\n", " print(f\" shell r={s['r']} weight {s['alpha']:.3f} n {s['n']:>9,} mass {s['mass']:.3e}\")\n", "print(f\"\\nraw union (alpha=1) {prof['union']:.3e}\")\n", "print(f\"retained (alpha=0.1) {prof['retained']:.3e}\")\n", "print(f\"observed bound {P.observed_mass(model, junctions):.3e}\")\n", "\n", "# Radius 2 on the whole set hits the component ceiling. That is the documented behaviour, not a\n", "# crash: components are independent, so their union masses add exactly and splitting loses nothing.\n", "try:\n", " P.shell_profile(model, junctions, r=2, alpha=0.1)\n", "except MemoryError as e:\n", " print(f\"\\nr=2 on all {len(junctions):,} junctions ->\\n {e}\")\n", "\n", "# On a subset small enough to fit, the radius-2 profile runs and shell 2 is the dominant volume.\n", "prof2 = P.shell_profile(model, junctions[:200], r=2, alpha=0.1)\n", "print(\"\\nradius-2 profile on 200 of them:\")\n", "for s in prof2[\"shells\"]:\n", " print(f\" shell r={s['r']} weight {s['alpha']:.3f} n {s['n']:>9,} mass {s['mass']:.3e}\")\n", "print(f\" retained {prof2['retained']:.3e} against a raw union of {prof2['union']:.3e}\")" ] }, { "cell_type": "markdown", "id": "6291224c", "metadata": { "language": "markdown" }, "source": [ "## 5. What was never catalogued\n", "\n", "Two answers, and they are not interchangeable.\n", "\n", "**The ball census** is finite and exact given the ball: `n_union - n_observed` sequences sit in the\n", "neighbourhood of a known cognate TCR and appear in no database.\n", "\n", "**Horvitz–Thompson** extrapolates instead, using the fact that `Pgen` *is* the sampling probability,\n", "so the inclusion probability $\\pi(j) = 1 - e^{-N p_j}$ is known rather than fitted. Good–Turing is\n", "the wrong tool here — Laydon et al. (*PLoS Comput Biol* 2014;10:e1003646) measure 61.7% median\n", "error for it on real TCR abundance data, because the capture-probability distribution is far too\n", "heterogeneous for the uniform-multinomial assumption behind it.\n", "\n", "The **mass** converges, because the weight $p/\\pi \\to 1/N$ as $p \\to 0$. The **count** does not,\n", "because $1/\\pi$ diverges — so `richness_reliable` says when the extrapolated count has stopped\n", "meaning anything." ] }, { "cell_type": "code", "execution_count": 7, "id": "35fb550f", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:30:51.548081Z", "iopub.status.busy": "2026-08-17T15:30:51.547912Z", "iopub.status.idle": "2026-08-17T15:30:52.930164Z", "shell.execute_reply": "2026-08-17T15:30:52.929276Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ball census: 1,702,574 uncatalogued candidates carrying mass 7.101e-03\n", "Horvitz-Thompson: declined -- every junction is a singleton: the capture curve is unidentified (theta -> 0, the Horvitz-Thompson sum diverges)\n" ] } ], "source": [ "# Capture units: how many distinct studies re-reported each junction for this epitope.\n", "mult = (beta.group_by(\"cdr3\").agg(pl.col(\"reference_id\").n_unique().alias(\"m\"))\n", " .filter(pl.col(\"cdr3\").is_in(junctions)))\n", "\n", "un = P.unseen_junctions(model, mult[\"cdr3\"].to_list(), mult[\"m\"].to_list(),\n", " n_units=beta[\"reference_id\"].n_unique())\n", "print(f\"ball census: {u['n_union'] - u['n_seqs']:,} uncatalogued candidates \"\n", " f\"carrying mass {prof['union'] - P.observed_mass(model, junctions):.3e}\")\n", "if un[\"degenerate\"]:\n", " print(f\"Horvitz-Thompson: declined -- {un['reason']}\")\n", "else:\n", " print(f\"Horvitz-Thompson: {un['n_unseen']:,.0f} unseen carrying {un['unseen_mass']:.3e}; \"\n", " f\"an average unseen junction is {un['rarity_ratio']:,.0f}x rarer than an observed one\")\n", " print(f\" richness_reliable = {un['richness_reliable']} \"\n", " f\"(min inclusion {un['min_inclusion']:.3g})\")" ] }, { "cell_type": "markdown", "id": "8e220d19", "metadata": { "language": "markdown" }, "source": [ "## 6. From a mass to cells in a person\n", "\n", "`q` is the selection constant carrying a generation probability to a post-selection repertoire\n", "frequency. It defaults to `1`, i.e. **uncalibrated** — the returned `F` is then a raw model mass\n", "whose ranking, not scale, is meaningful. ALICE's published TRB value is `9.41`; the matched\n", "comparison against the model-free event ratio puts the same factor at a median 14.8×, so the two\n", "are the same order and not the same number. Whatever you pass is echoed back, so a reported number\n", "always carries its calibration.\n", "\n", "`F(e)` answers \"is there a precursor\". A detectable response needs more than one, which is what\n", "`p_ge_k` is for — and once $N_{\\mathrm{eff}} F$ is of order one the two stop being a monotone\n", "reparametrisation of each other." ] }, { "cell_type": "code", "execution_count": 8, "id": "3d26611d", "metadata": { "execution": { "iopub.execute_input": "2026-08-17T15:30:52.932102Z", "iopub.status.busy": "2026-08-17T15:30:52.931939Z", "iopub.status.idle": "2026-08-17T15:31:09.648220Z", "shell.execute_reply": "2026-08-17T15:31:09.647347Z" }, "language": "python" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "F(e) 9.792e-03 (q = 9.41, alpha = 0.1)\n", "expected cells 293,772,385 of 1e11 x 0.3 CD8\n", "lambda 9.79e+05\n", "P(>= 1 precursor) 1.0000\n", "P(>= 10) 1.0000\n" ] } ], "source": [ "# CD8 compartment of a 1e11-cell pool, and precursor probabilities at a realistic clonotype count.\n", "out = P.precursor_frequency(model, junctions, r=1, q=P.ALICE_Q,\n", " n_cells=1e11, compartment=0.3, n_eff=1e8)\n", "print(f\"F(e) {out['F']:.3e} (q = {out['q']}, alpha = {out['alpha']})\")\n", "print(f\"expected cells {out['cells']:,.0f} of 1e11 x 0.3 CD8\")\n", "print(f\"lambda {out['lambda']:.3g}\")\n", "print(f\"P(>= 1 precursor) {out['p_ge_1']:.4f}\")\n", "print(f\"P(>= 10) {out['p_ge_10']:.4f}\")" ] }, { "cell_type": "markdown", "id": "63688019", "metadata": { "language": "markdown" }, "source": [ "## 7. The same thing from the command line\n", "\n", "```console\n", "$ vdjmatch precursor --vdjdb --mhc-class MHCI --min-junctions 10 --q 9.41 -o precursor.txt\n", "```\n", "\n", "writes one row per (epitope, chain) with every column above. `--group-by` and `--chain-col` do the\n", "same for your own table.\n", "\n", "## What this does not answer\n", "\n", "`F(e)` is the probability a precursor **exists**, not that a response is **mounted**: that also\n", "needs priming, help and an absence of tolerance, none of which is modelled here. And every set\n", "total — the observed mass, the union, the event ratio alike — grows with how much of the true\n", "cognate set the database holds, because that enters the numerator alone. A claim that ranks\n", "epitopes by a set total needs that control stated alongside it." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.7" } }, "nbformat": 4, "nbformat_minor": 5 }