Designing a cassette: what goes in, and what it is worth ========================================================= A vaccine cassette is a **set**, and the quantity that decides whether it works is not how good its units are on average but whether *several* of them elicit a response. Sorting a candidate list and keeping the top *m* answers that question correctly only if the units respond independently. They do not, and this page is what the difference costs and what to do about it. Two commands: .. code-block:: bash mhcmatch cassette select --candidates pool.tsv -k 20 --tol 3 --out cassette.tsv mhcmatch cassette score --cassettes cassette.tsv --pool pool.tsv :doc:`safety` is the step before — which units to withdraw before capacity is spent on them. :doc:`portfolio` is the geometry underneath: the response model, the Pareto/reachability results, and the measured over-dispersion. This page is the operational middle. .. contents:: On this page :local: :depth: 2 The objective, and where it comes from -------------------------------------- It is **derived from the design goal**, not fitted to an outcome cohort. Write :math:`R_i` for unit *i*'s response indicator and :math:`p_i = E[R_i]` for its calibrated probability. The breadth of a cassette :math:`S` is :math:`B(S) = \sum_{i \in S} R_i`, with .. math:: E[B] &= \sum_i p_i \\ Var[B] &= \sum_i s_i^2 + 2 \sum_{i` records it. .. note:: **The sequence channel counted exact shared 3-mers until v2, and that was measured to be a duplicate detector rather than a similarity.** Over the eight TESLA and eleven HiTIDE donors only **4,053 of 150,994 within-donor pairs (2.68 %)** share any 3-mer at all, and **17.6 % of those that do are the same peptide window**. It is also blind to chemistry: ``GILGFVFTL`` against ``GILGFVFTV`` and against ``GILGFVFTW`` share the same six 3-mers, though one substitution is conservative and the other is not. :func:`~mhcmatch.cassette.sequence_overlap` scores those two pairs 6 and 19 on the BLOSUM distance. The k-mer form is kept, as ``overlap(..., features=None)`` with :data:`~mhcmatch.cassette.KMER`, so results recorded under it reproduce. .. _cassette-v2: Selecting on the degeneracy (``rule="v2"``) ------------------------------------------- ``p_i`` is a probability, so the number of units that respond is a **random variable**, and many size-*k* sets are indistinguishable in it. A sort already maximises the expected count; the sets it cannot tell apart are not a nuisance, they are the design freedom. ``mhcmatch cassette select --rule v2`` returns, among all cassettes that are — with stated probability — no worse than the ranked list, the one whose units share the fewest ways of failing. :func:`~mhcmatch.cassette.not_worse` computes that probability, and it is cheap for an exact reason: units in both sets are the **same random variable**, not merely identically distributed, so they cancel and only the symmetric difference carries variance. ``--not-worse 1.0`` returns the reference exactly; lower values buy diversity and say how often you are willing to be wrong. .. warning:: **``--not-worse`` is a per-donor guarantee.** It bounds ``P(this donor's cassette catches at least as much as this donor's sort)`` and says nothing about a sum over donors. At ``0.5`` every donor independently accepts a coin-flip, so a cohort-level count is worse than the sort most of the time. A pooled comparison needs a tighter floor than intuition suggests. The reference is the top-*k* sort unless ``reference=`` names another set. That matters: v2 only ever trades capture *away* from its reference, so the reference is a floor the rule cannot fall below by more than the stated probability, and never a rival it can beat. ``cassette select`` ------------------- .. code-block:: bash mhcmatch rank fasta windows.fa --alleles "$HLA" --cls mhc1 --out ranked.tsv mhcmatch cassette select --candidates ranked.tsv -k 20 --tol 3 -v Four steps, in order: 1. **The offset** is fitted once over the pool by :func:`~mhcmatch.cassette.prob_offset` at ``--prevalence``, and held. Fitting it over the chosen set instead would pin every donor's cassette to the same mean and destroy the comparison the score exists to make (below). 2. **``rho``** is the measured background, or yours. 3. :func:`~mhcmatch.cassette.goal_energy` turns ``(p, overlap, rho, gamma)`` into ``(h, J)``. 4. :func:`~mhcmatch.cassette.greedy` takes ``k + tol`` units in :math:`O(kN)` — about 4,000 operations for twenty of two hundred — :func:`~mhcmatch.cassette.refine` swaps until no single exchange raises ``H``, and the reported size is the one with the largest ``H`` in ``[k - tol, k + tol]`` --- with the lower end raised to the coverage floor, the number of ``universe`` allotypes the pool can supply, since no smaller cassette can hold them. Greedy plus the swap pass reaches the **brute-force optimum** on every pool small enough to enumerate; that is the only warrant the :math:`O(kN)` rule has and it is a test rather than a claim. .. important:: **Give it the whole candidate pool, not a shortlist.** ``binder`` and ``expr_lvl`` are the two largest positive coefficients in the shipped model and ``expr_norm`` is positive too — run ``mhcmatch rank --coefficients`` for the sizes, which move at every refit — so a pool that has already been cut on binding and expression has no range left along them. This is measurable rather than arguable: on the 46-patient half of the NCI gastrointestinal screen held out of the EPIC fit, an **exhaustive** exome screen responding at 0.0144 per mutation, selection lifts captured responses to **3.31× the base rate** at *k* = 5 (11 of 58 positives against 3.3 expected). On TESLA's *nominated* list — the same disease question, but candidates a consortium's pipelines had already put forward, responding at 0.0612, **4.25×** the NCI rate — every rule sits at the base rate, because the selection had already been done. Full table in ``bench/results/cassette_select.md``. .. note:: ``bench/results/...`` paths on this page resolve in the benchmark repository, ``2026-mhcmatch-code`` (private; released with the manuscript), not in the library repo. **``--tol`` is spent on the objective, not on the largest size that fits.** A mean–variance objective has an internal optimum size, and where it falls moves with the prevalence and with ``rho``, so ``-k 20 --tol 5`` returns whichever size in 15–25 carries the largest *H* and says so on stderr. With the per-unit ``gamma`` this is a per-donor answer: on the eight TESLA pools ``-k 20 --tol 5`` returns sizes 19 to 25 and on the eleven HiTIDE pools 20 to 23. With ``gamma`` passed undivided it returned 15 — the floor of the window — for every donor of both, which is what a :math:`k^\star` below the requested size looks like from the outside. With ``--tol 0`` the size is exactly *k*, which is what a fixed manufacturing budget wants. **When the budget is a confidence rather than a count, ask for the size.** :func:`~mhcmatch.cassette.size_for` returns the smallest cassette reaching :math:`\Pr(\ge m \text{ responses}) \ge C` for one donor's own pool, and ``--confidence C`` --- with ``-k`` read as the manufacturing ceiling --- asks for it from the command line. ``--block-live`` reaches the probe as well as the selection, so a cassette that can lose a whole allotype at once is sized for that rather than against it. .. code-block:: bash mhcmatch cassette select --candidates pool.tsv -k 40 --confidence 0.90 \ --prevalence 0.026 --out cassette.tsv A donor whose head of list is genuinely strong reaches 0.90 in ten units; a donor whose is not needs thirty, and one who cannot reach it inside the ceiling is reported at the ceiling with ``reached = False`` rather than rounded down into a cassette that claims the target. **``--prevalence`` is what lets it see that, and the default is one pool's number.** The map from score to probability is :math:`\sigma(s + b)` with a single additive offset: the *slope* is measured and is right — :math:`\alpha = 1.0004 \pm 0.0364` over 339,598 labelled rows, likelihood ratio 0.0 against 1 — and the *level* is stated, because EPIC carries one unpenalised intercept per screen and no global one, so it is not identifiable from the fit. :data:`~mhcmatch.rank.POOL_PREVALENCE` is 0.0602; TESLA's own candidates respond at 0.0462 and HiTIDE's at 0.0263, so that one default over-states predicted yield by 1.3× on the first pool and 2.3× on the second. Pass the pool's own expected base rate. ``cassette score`` ------------------ .. code-block:: bash mhcmatch cassette score --cassettes manufactured.tsv --pool candidates.tsv Group rows by a ``donor`` column and one file may hold many cassettes of different sizes. Returned per cassette: =================== ===================================================================== ``yield`` ``sum p`` — the expected number of responding units. **A level, not a probability** ``p_at_least`` ``P(X >= target)`` under the block model, exactly, via a convolution of per-block Poisson binomials — no Monte Carlo and no :math:`2^B` enumeration ``n_effective`` how many *independent* shots the cassette is worth ``lam`` nats above a uniform random subset of that donor's own pool ``rho_hla`` / ``rho_seq`` / ``rho_dom`` the three pairwise statistics, each by its exact closed form ``yield_loh`` / ``lost_allotype`` expected responding units left after the **worst single** allotype is lost, and which one that is. ``yield_loh / yield`` is the share of expected response that does not depend on any one allotype ``coverage`` allotype counts, Gini, and share of maximum entropy --- against ``--universe`` when given, which is what makes an allotype holding zero units visible =================== ===================================================================== ``lam`` is the one that crosses donors **and** sizes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. math:: \lambda(S) = H(S) - \log \sum_{|S'| = k} e^{H(S')} + \log \binom{N}{k} The middle term is the exact log partition function over every size-*k* subset of the donor's pool, computed by the elementary-symmetric recurrence in log space (:func:`~mhcmatch.cassette.log_ek`) — so it never enumerates a subset, and :math:`\binom{5000}{20}` is not a number anybody was going to sum over. Adding :math:`\log \binom{N}{k}` back makes the comparison against the *average* subset rather than their sum, so **zero is a cassette exactly as good as a uniformly random one from the same pool**, positive is better, and the units are nats. Dividing by the donor's own pool is what removes both pool depth and *k*. Measured on 3,064 TCGA donors: a cassette built by sorting the candidate list on the ranker scores a median :math:`\lambda = -0.408` nats — *below* a uniform random subset of the same pool — against **+3.164** for the greedy argmax of ``H``, a gain of **+3.490** nats. .. note:: ``score`` does **not** report ``H``. :func:`~mhcmatch.cassette.goal_energy` renormalises the overlap to the set it is handed, and the dominance channel is scaled by that set's range — so an ``H`` computed on a cassette alone is not the ``H`` ``select`` maximised over the pool, and a rule that spent expected count on non-overlapping units would score identically to one that did not. To compare two rules on the objective, build ``(h, J)`` once over the pool and evaluate both index sets with :func:`~mhcmatch.cassette.energy`. That is five lines and it is exact. Allotype coverage, and why it is the HLA-loss question ------------------------------------------------------ ``coverage`` looks like a tidiness metric and is not. It is the readout of the one failure mode that takes a whole group of units at once. **What it measures.** Given the units' allotype labels, :func:`~mhcmatch.portfolio.coverage` returns the per-allotype counts, a Gini index (0 = every allotype equally covered, → 1 = every unit on one), the share of maximum entropy, and ``n_covered`` of ``n_allotypes``. **Pass ``--universe`` --- the donor's *distinct* allotypes --- or the index answers a different question.** Computed over the labels the cassette happens to carry, an allotype holding **zero** units is invisible, and a zero is exactly the inequality the index exists to report. The same argument runs the other way for a homozygous donor: a patient homozygous at *B* has five distinct class-I allotypes, not six, so an even cassette over five is perfectly even and scoring it against a denominator of six would report a genotype as a design flaw. .. code-block:: bash mhcmatch cassette select --candidates pool.tsv -k 20 \ --universe "$HLA" --block-live 0.8 --max-share 0.5 --out cassette.tsv Why it matters: an allotype is a group of units that fail together ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Every unit credited to one class-I molecule shares that molecule's presentation, its precursor niche, and its fate. If the tumour loses the allele --- or downregulates it, or the typing was wrong --- **all of them go at once**. A cassette of twenty units on two allotypes is two shots, not twenty, and no per-unit score can see that, because it is a property of the set. :func:`~mhcmatch.portfolio.survival` has modelled this since it was written: a unit responds only if its block is live *and* its own term fires, :math:`R_i = B_b \varepsilon_i` with :math:`B_b \sim \mathrm{Bern}(q_b)`. ``--block-live`` is that :math:`q`. What it buys the objective is not a heuristic but a covariance --- for two units on one allotype, .. math:: \mathrm{Cov}(R_i, R_j) = q_b r_i r_j - q_b^2 r_i r_j = (1 - q_b)\, p_i p_j / q_b and zero across allotypes. So losing an allele contributes exactly :math:`\gamma (1 - q_b) p_i p_j / q_b` to :math:`J_{ij}` and nothing anywhere else. No ``rho``, no overlap channel, and no parameter that is not the loss rate the designer stated. At :math:`q = 1` the term is identically zero and every cassette built before it existed is reproduced bit for bit. That is worth entering as itself. :func:`~mhcmatch.cassette.overlap` returns the *mean* of its two or three channels, so with all three populated a same-allotype pair reaches :math:`J` at one third weight, diluted by whether the two peptides happen to share 3-mers. What it is worth, measured ~~~~~~~~~~~~~~~~~~~~~~~~~~ On the six TESLA donors (605 nominated candidates, 37 validated-immunogenic, pools 73–144) at *k* = 20, scored genotype-free through the identical path ``cassette_select.md`` uses so the only difference between arms is the selection rule: ============== ========== ================ =========== arm captured ``captured_loh`` ``rho_hla`` ============== ========== ================ =========== ``sort`` 7 **1** 0.457 ``select`` 8 2 0.305 ``select+loh`` **10** **4** 0.290 ============== ========== ================ =========== ``captured`` is validated units in the cassette, pooled over the six donors; ``captured_loh`` is how many are left after the **worst single allotype** is lost. The worst case rather than an average over losses, because LOH takes a specific allele and a designer asking to be protected is asking about the bad draw. Ranking the list and taking the head keeps **1 of its 7** captured units through that draw; pricing the loss at :math:`q = 0.8` keeps **4 of 10**. Full table, per donor and at *k* = 5/10/20, in ``bench/results/cassette_tesla_donors.md``. .. note:: ``select`` already spread without being told to --- ``rho_hla`` 0.305 against the sort's 0.457 --- because the allotype channel of the overlap was always one of the three. Naming the loss rate is what turns that from a side effect into a stated design parameter with a number on it. The floor is a constraint, not an objective term ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``--max-share`` caps any one allotype's share of the cassette; ``--universe`` --- or ``--floor``, which takes the floor from the allotypes the donor's own pool carries rather than from a stated genotype --- gives every allotype the pool can supply a unit before the free slots are filled. Both are **manufacturing constraints** and are deliberately outside :math:`H`: the loss coupling already prefers spread, and stacking a second diversity term inside the objective double-counts unless you mean it --- the argument :func:`~mhcmatch.portfolio.compose` already makes for ``weight_evenness``. An infeasible pair (a share cap too tight to fill *k* across the allotypes the floor demands) raises with the arithmetic rather than quietly returning a cassette that breaks one of the two. ``--weight-coverage`` is the third option and it is neither of those two: a **stated** exchange rate, in expected responding units, paid to put a unit on an allotype the cassette does not yet reach. Coverage is a property of the set rather than of a unit, so it cannot be a field term and enters the greedy marginal gain instead; the bonus is monotone and submodular, so the :math:`1-1/e` bound above is untouched. It ships at ``0.0``, where every path is bit-identical and coverage remains what it has always been here — **reported, never optimised**. .. note:: **A constraint that does not bind changes nothing, which is why the weight exists.** At *k* = 20 over four to six allotypes the floor is already satisfied by the unconstrained argmax and a 0.34 share cap is seven slots where no allotype holds seven — measured on 19 donors of two pools, ``--universe`` plus ``--max-share 0.34`` returned numbers *identical* to plain ``select`` on every one of them. Use the constraints when a floor must hold; use the weight to trade spread against expected response continuously. Quote it with its counterfactual at zero, the contract ``--weight-escape`` already follows — the CLI prints both. An allotype the pool cannot supply is skipped rather than raising. That is a fact about the donor's candidates, and it shows up as ``n_covered`` below ``n_allotypes`` where a caller can act on it. Tumour selectivity: a stated preference, not a refit ----------------------------------------------------- "High in the tumour, low in healthy tissue" is a design goal, and the shipped ranker does not share it. EPIC fits **both** expression terms positive --- v12 puts ``expr_lvl`` at **+0.5000** and ``expr_norm`` at **+0.2222** log-odds per standard deviation, the first being the largest coefficient after presentation itself --- so *as fitted, high normal-tissue expression is rewarded*. (``mhcmatch rank --coefficients`` prints the set an install actually scores with.) That is not a defect: the model was fitted on **will this respond**, and a gene transcribed everywhere responds more often. Selectivity is a different question, and it is a **safety** question. So it enters as a declared exchange rate, the way ``gamma`` does: .. math:: h_i = p_i - \tfrac{\gamma}{2} s_i^2 + w \,(\mathrm{expr\_lvl}_i - \mathrm{expr\_norm}_i) ``w`` is in expected responding units per **log2-fold** of tumour-over-normal abundance --- both terms are :math:`\log_2(1 + \mathrm{TPM}/c)` on one floor, so their difference is a log2 ratio. .. code-block:: bash mhcmatch rank fasta windows.fa --alleles "$HLA" --out ranked.tsv # emits both terms mhcmatch cassette select --candidates ranked.tsv -k 20 --selectivity 0.05 -v Three properties, and each is the reason for a design choice: * **Charged to the objective, never to** ``p``. ``p`` is a calibrated marginal that :func:`~mhcmatch.portfolio.survival` reads literally, so discounting it would silently restate the response model as well as the preference. Same rule ``compose``'s ``weight_cost`` follows. * **Nothing is asserted about the fit.** Both coefficients stay as measured and both terms stay reported. Imposing the tumour/normal *ratio* on the model --- equal and opposite coefficients --- would assert an answer the data rejects. * **The run reports its own trade**: what the same pool would have built at ``w = 0``, the expected units given up, and the mean log2-fold bought. A stated weight that does not report its cost is a knob, not a preference. A candidate missing either term takes a delta of **0**, not ``nan`` --- ``nan`` would reach the argmax and delete the candidate, where 0 leaves it ranked on everything else. Escape: what the tumour has to give up to lose a unit ----------------------------------------------------- The objective has priced one escape route since it was written and did not price the other. :func:`~mhcmatch.cassette.goal_energy` charges same-allotype pairs for the loss-of-heterozygosity event that takes both --- derived from the block model, not fitted, and identically zero at ``block_live=1.0``. But a tumour that cannot lose an allotype can still lose the *antigen*: delete the mutation, silence the locus, or let the subclone carrying it be outgrown by one that never had it. Two arguments, and they are separate. ``escape`` with ``weight_escape`` is the **field**: a per-unit cost in ``[0, 1]`` --- the tumour's own price for deleting that candidate --- and a stated exchange rate on it, charged as ``h_i += weight_escape * escape_i``. A clonal driver in a gene the tumour cannot silence is expensive to lose; a subclonal passenger is free. ``genes`` is the **coupling**: two units from one source gene fall together under one deletion, the antigen-loss twin of the allotype channel. Both are needed, because a set can be made of expensive units that all sit in one locus. :func:`~mhcmatch.cassette.escape_cost` builds the field for you, so the cost is one definition rather than one per caller:: eps = CA.escape_cost(is_clonal, detected=above_floor, gene_driver=is_cancer_gene, residue_driver=is_hotspot, expr=tpm) **Clonality enters as two states, not as a fraction, and that is deliberate.** A cancer-cell fraction near 0.05 is consistent with a great many small subclones, and what the units carried by any one of them do to each other is not something bulk sequencing observes --- so treating 0.05 and 0.25 as a five-fold difference in escape cost asserts a resolution the assay does not have. What bulk data does support is two claims: that a variant is **above the detection floor**, and that it sits at the **heterozygous-clonal or homozygous frequency** with the coverage to say so, which is the case where a known initiating driver reads at or near half the reads. Hence ``1.0`` for clonal, :data:`~mhcmatch.cassette.C_SUBCLONAL` for detected-but-not-clonal, ``0`` for neither. A caller holding a CCF rather than a call can threshold it at :data:`~mhcmatch.cassette.CLONAL_CCF`; a caller holding the depositor's own binary call should pass that and ignore the constant. The two driver flags are **graded, not conjoined**: :data:`~mhcmatch.cassette.D_PASSENGER` where neither fires, :data:`~mhcmatch.cassette.D_ONE_SIDED` where one does, ``1.0`` where both do. Measured on 465,343 TCGA units the strict conjunction fires on 3,269, and 4,709 of 7,261 donors carry none at all, so a binary term would leave the weight selecting on clonality and expression alone for two donors in three while appearing to select on drivers. ``expr`` becomes a **within-pool percentile**, because a locus already near the floor is one the tumour can silence for nothing. Nothing in it is fitted to an outcome: durability is a preference over a horizon no one-timepoint response screen observes, so both constants are stated once and held. A missing annotation becomes a number rather than a ``nan`` --- a ``nan`` reaches the argmax and silently deletes the candidate. .. code-block:: console $ mhcmatch cassette select --candidates pool.tsv -k 20 \ --escape-column eps --weight-escape 0.25 --gene-channel -v # -: escape w=0.25 traded yield 1.671 -> 1.250 unit(s) for mean escape cost # 0.390 -> 0.900; 6 of 20 slot(s) changed The counterfactual is printed, not optional --- it is the ``--selectivity`` contract, for the same reason. A stated weight is auditable only if what it gave up is on the record beside what it bought, and a cassette quoting one without the other has not said what it cost. **Stated, not fitted**, like ``gamma`` and ``selectivity``. Durability is a preference over a horizon no response screen observes: a screen reading out at one timepoint cannot supply an exchange rate between catching a response now and keeping it later. A non-finite cost contributes **0**, the same contract :func:`~mhcmatch.cassette.selectivity_delta` follows --- for a driver annotation covering under half its mutations, that is the difference between ranking a candidate on everything else and deleting it silently. :func:`~mhcmatch.cassette.size_for` takes the same three arguments and must be given them. It walks its own greedy order, so without them a ``--confidence`` size and a weighted selection are answers to two different questions. Joining the metadata an escape cost needs ------------------------------------------ The clone assignment, cancer-cell fraction, driver calls and expression that ``escape_cost`` reads live in the variant caller's output and the RNA table, never in a candidate list. ``cassette select`` joins them, and builds the cost from the joined columns, so nothing has to be precomputed:: $ mhcmatch cassette select --candidates pool.tsv -k 20 \ --metadata clones.tsv --metadata-on donor,gene \ --ccf-column ccf --driver-genes cgc.txt --escape-expr-column tpm \ --weight-escape 0.5 --gene-channel ``--metadata`` is a **left join and never drops a candidate**: an unmatched row keeps its own columns and takes no escape cost, because the pool is what defines the background the choice is made against and a silently narrowed pool is a different question. A key column missing from either side raises, and a join matching zero rows raises --- a silent no-op returns the table unchanged and looks like it worked. ``--metadata-on`` takes one or more comma-separated keys; ``donor,gene`` and ``donor,peptide`` are the usual ones. The escape flags mirror :func:`~mhcmatch.cassette.escape_cost` argument for argument: ``--clonal-column`` a boolean call, or ``--ccf-column`` with ``--clonal-ccf`` for a caller holding only a fraction; ``--detected-column``; ``--driver-genes`` a file of gene symbols, one per line; ``--driver-column`` for residue-level evidence; ``--escape-expr-column`` for the abundance whose within-pool percentile becomes the floor. They are mutually exclusive with ``--escape-column``, which takes a cost the caller already computed, and passing both raises rather than silently preferring one. One design per group inside a donor ------------------------------------ ``--group-column`` designs a separate cassette for each value of a column *within* each donor --- one per tumour clone, one per variant class --- and emits the group beside the chosen units:: $ mhcmatch cassette select --candidates pool.tsv -k 20 --group-column clone **The calibration offset is the donor's, not the group's**, fitted once on that donor's whole pool and passed to every group. That is the whole reason the flag exists rather than being a shell loop: :func:`~mhcmatch.cassette.prob_offset` anchors the mean of the batch it is handed, so calibrating each group on itself pins all of them to the same declared prevalence and deletes exactly the difference between them. The next section is that trap in full. In the library the same thing is ``CA.select(..., offset=base.offset)``. A page for the person who has to defend the design -------------------------------------------------- Every counterfactual a selection run computes goes to stderr and is lost. ``cassette report`` puts it on one self-contained HTML page --- no JavaScript, no image files, nothing to install:: $ mhcmatch cassette select --candidates pool.tsv -k 20 --out sel.tsv $ mhcmatch cassette report --cassettes sel.tsv --pool pool.tsv --out design.html Four things, in the order a reviewer asks for them: * **What the design bought, against ranking the same pool.** The top-*k* by score is scored on the same offset and the same axes, so the two columns are levels and their difference is the trade. When the two rules chose the same units the page says so --- on a pool where each unit carries one restriction the couplings often cannot reorder the ranking at all, and that is a finding about the pool rather than a broken run. * **The units**, with the gene and the allotype each was credited to. * **Units per allotype**, drawn against the donor's *whole* genotype, so an allotype holding zero chosen units is visible. Pass ``--universe`` where the pool does not carry every allotype the donor has; without it an allotype that never appears in the pool cannot be drawn as empty. * **The escape routes**: which chosen units fall together on one allotype, and which on one source gene. One loss-of-heterozygosity event or one gene deletion takes a whole group. The page computes nothing of its own --- every number on it is a key :func:`~mhcmatch.cassette.score` already returns. ``escape`` is the one column that is *not* per-unit: :attr:`~mhcmatch.cassette.Cassette.escape` is the mean over the chosen units, written onto every row of that cassette, so the page reports it once and the per-unit table leaves it out. Reading one patient against a cohort ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A ``lam`` on its own is a number in nats. What makes it legible is where it falls among patients whose outcome is known, and ``--reference`` supplies that:: $ mhcmatch cassette report --cassettes sel.tsv --pool pool.tsv --out design.html \ --reference cohort.tsv --reference-outcome responded ``cohort.tsv`` is any TSV with a kill-pressure column (``--reference-column``, default ``lam``) and optionally a binary outcome column. The page then reports the patient's percentile in that cohort, which third of it they fall in, and --- when an outcome column is named --- the observed rate in that third beside the rate across the whole cohort. **The library ships no cohort and no reference table, deliberately.** A distribution baked in here could not be audited by the person relying on it, could not be corrected when the underlying record moved, and would silently apply a lung-cancer cohort to a melanoma patient. Passing the file makes the yardstick explicit. The thirds are thirds and not a fitted cut: a threshold chosen after looking at an outcome is a fitted parameter wearing a threshold's clothes. What the page prints is an observed rate in a cohort, and it says so --- it is not a prediction for the patient in front of you. The calibration offset decides *what is being reported* -------------------------------------------------------- This is the trap, and it is worth a section because it is silent. :func:`~mhcmatch.rank.probability` anchors the mean of **the batch it is handed**. Called once per donor — which is what a per-sample pipeline does without thinking about it — it pins *every* donor's pool mean to the declared prevalence, whatever their pool holds. On 7,261 TCGA donors with pools spanning **1 to 5,221** candidates, every per-donor-anchored pool mean lands on **0.060163** with a standard deviation of **3.37 × 10⁻¹⁷**. Read as a probability, that number is not one, and two donors' numbers are the same number. ================================ ========================== ========================== offset over the batch one offset per donor ================================ ========================== ========================== what ``sum p`` means a **level**: expected an **enrichment**: how far responding units the chosen units sit above that donor's own background pool mean ``p``, range 0.002977 – 0.435027 0.060163 – 0.060163 spread (sd) 2.47 × 10⁻² 3.37 × 10⁻¹⁷ comparable between donors? yes no against an IFN-γ signature ρ = **+0.1261** ρ = **+0.1322** ================================ ========================== ========================== **Neither is wrong and the enrichment is the stronger readout** — on 4,073 TCGA donors across 30 tumour types it correlates better with immune infiltrate on all four independent gene-set constructions. They are two different quantities, and which one you want is a decision. :func:`~mhcmatch.cassette.prob_offset` gives the level, :func:`~mhcmatch.cassette.group_offsets` gives the enrichment for every group at once, and ``mhcmatch cassette score --per-donor-offset`` switches between them at the command line. In the Nextflow module, ``MHCMATCH_CASSETTE_SCORE`` collects every sample before scoring for exactly this reason — it is the one process in that subworkflow that is deliberately not per sample. Python ------ .. code-block:: python import numpy as np from mhcmatch import cassette as CA scores = ... # mhcmatch.rank.aggregate_score over the donor's WHOLE pool peptides = ... # the long window around each mutation, not the minimal epitope alleles = ... # optional; populates the allotype channel of the overlap c = CA.select(scores, peptides, alleles, k=20, tol=3) print(c.k, c.yield_, c.lam, c.channels) s = CA.score(scores, peptides, alleles, chosen=c.index, pool_scores=scores, pool_peptides=peptides, offset=c.offset) print(s["yield"], s["p_at_least"], s["lam"], s["n_effective"]) A pool smaller than *k* returns the whole pool rather than raising: there is nothing to choose, and refusing would delete the donor from a cohort-scale run over a fact the caller can read off ``pool_n``. The next step is assembly — spacers, ordering, junction scanning, back-translation — which is ``mhcmatch cassette build`` and :mod:`mhcmatch.vector`. See :doc:`safety`. API --- Every function above, with its full docstring: :doc:`api` — :mod:`mhcmatch.cassette`.