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:

mhcmatch cassette select --candidates pool.tsv -k 20 --tol 3 --out cassette.tsv
mhcmatch cassette score  --cassettes cassette.tsv --pool pool.tsv

Safety, prior evidence, and what goes in the cassette is the step before — which units to withdraw before capacity is spent on them. Cassette composition is the geometry underneath: the response model, the Pareto/reachability results, and the measured over-dispersion. This page is the operational middle.

The objective, and where it comes from

It is derived from the design goal, not fitted to an outcome cohort. Write \(R_i\) for unit i’s response indicator and \(p_i = E[R_i]\) for its calibrated probability. The breadth of a cassette \(S\) is \(B(S) = \sum_{i \in S} R_i\), with

\[\begin{split}E[B] &= \sum_i p_i \\ Var[B] &= \sum_i s_i^2 + 2 \sum_{i<j} \rho_{ij} s_i s_j, \qquad s_i = \sqrt{p_i (1 - p_i)}\end{split}\]

On the mean alone the optimiser is a sort and there is no design problem. The design problem is entirely in the variance, and the variance is not the independent one. A designer who wants at least m units to respond is worse off with a positively correlated portfolio of the same mean, because that is the one with the fatter lower tail. So the objective is mean–variance:

\[H(S) = \sum_i \left[ p_i - \tfrac{\gamma}{2} s_i^2 \right] - \gamma \sum_{i<j} \rho_{ij} s_i s_j\]

which is exactly \(\sum h_i - \sum J_{ij}\). The Potts form is not imposed on the goal; it falls out of it.

Three inputs, and not one of them is an outcome cohort:

p_i

the calibrated response probability from probability(), fitted on immunogenicity screens.

rho

one number, the mean intra-cassette response correlation, measured once on published per-unit assays. Four cohorts have been measured and they do not agree: 0.124 on the Sahin TNBC mRNA trial (41 of 216 assayed units, 13 patients, 3.45× the independent-Bernoulli variance), 0.091 on IVAC MUTANOME (75 of 125 units, 13 patients, 1.8×), 0.024 on TESLA, 0.010 on HiTIDE. RHO_ASSAYED defaults to IVAC’s, because it is the only one of the four whose corpus carries a measured label on every manufactured unit — its 50 non-responding units are measured negatives rather than units nobody looked at. betabinom_rho() fits your own by maximum likelihood; it needs SciPy, so pip install 'mhcmatch[stats]'.

gamma

a stated design preference, 1.0: one unit of variance in the responding-unit count is worth one expected unit. Not fitted, not swept — but stated per unit, so it is divided by the design effect before use. For k units of mean response \(\bar p\) and mean pair correlation \(\rho\), \(H = k\bar p\,\{1 - \tfrac{\gamma}{2}\bar q\,[1 + \rho(k-1)]\}\); the brace is the worth of the average unit, and because a correlated count’s mean is linear in k while its variance is quadratic, an undivided gamma makes it fall with k and turn negative at \(k^\star = 1 + (2/\gamma\bar q - 1)/\rho\) — past which the optimiser prefers a worse unit to a better one. risk_aversion() divides by \(1 + \rho(k-1)\), which holds the brace at \(1 - \tfrac{\gamma}{2}\bar q\) at every size. rho is measured and k is given by the design, so no outcome enters the correction; passing gamma= uses the number given instead.

\(\rho_{ij}\) spreads rho over pairs in proportion to overlap() — the mechanistic similarity of the pair — then renormalises so the pool’s mean pair correlation is exactly rho. The overlap combines whichever channels the data supports — by default their mean, or a per-axis-normalised maximum under overlap_combine="worst", on the argument diversity() has carried since it was written: a cassette is undone by its worst shared failure mode, not its average one, and every axis added makes every existing axis count for less. The normalisation is not optional there — the channels run at off-diagonal means 0.011, 0.781 and 0.303 on the real pools, so an unnormalised maximum would read whichever axis has the largest raw spread every time. "mean" is the default and the channels are:

channel	what it says two units share
allotype	the same class-I molecule, so the same presentation and the same precursor niche, and they are lost together if that allele is. Graded rather than binary when presented is supplied (allotype_overlap())
sequence	BLOSUM-graded similarity of the TCR face (sequence_overlap()), so a conservative substitution reads as more shared than a radical one
physchem	closeness on TCR-face burial and charge, passed as features
expression	closeness on the source gene’s abundance, and GTEx tissue-profile similarity through coexpr (mhcmatch.expression.coexpression())
profile	how much two units owe their scores to the same terms, from the fitted model’s own decomposition (aggregate_terms(), profile_overlap()). Pass terms and terms_cov to select(), with dominance left off, which is the default
dominance	closeness on the score axis. Off by default, and always off in v2 — it is the one channel built from the score rather than from a mechanism, its pairwise statistic fits attractive on the observational arm where greedy() carries no bound, and it never abstains: it is zero on 0.03% of within-donor pairs against 97.3% for the 3-mer channel, so it supplied 71-79% of the total channel mass and the allotype channel — the only mechanism of the three — entered \(H\) at a third weight. Pass dominance=True to restore the three-channel form, and say so beside the number

Note

The profile channel is what dominance was reaching for. Dominance couples two units for scoring alike; this one couples them for scoring alike because of the same thing. A row of aggregate_terms() is a unit’s score broken into one contribution per fitted term, so two rows pointing the same way name the same failure mode — both carried by presentation, both carried by abundance — and the cosine between them is non-negative by construction, which keeps greedy() inside its 1 - 1/e bound.

Whiten against the cohort, never against the pool. Whitening n points against a covariance estimated from those same n points sends them to the vertices of a regular simplex, where every pairwise cosine is exactly -1/(n-1) whatever the data said — the coupling then carries no information and carries it silently. epic_axes() raises below SELF_COV_MIN rows per column rather than let that happen, and select() refuses a pool too small to estimate its own.

On the two labelled donor pools the channel is informative — off-diagonal mean 0.136, sd 0.20 over TESLA’s eight donors — and it does not out-catch the arms already there. It ships because it is the coupling the objective’s derivation asks for and because a cohort with more donors can test it; it is not claimed to catch more units.

Which channels were available is part of the result. A trial that published no per-patient genotype has one fewer, and Cassette.channels 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. sequence_overlap() scores those two pairs 6 and 19 on the BLOSUM distance. The k-mer form is kept, as overlap(..., features=None) with KMER, so results recorded under it reproduce.

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.

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

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 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. goal_energy() turns (p, overlap, rho, gamma) into (h, J).

4. greedy() takes k + tol units in \(O(kN)\) — about 4,000 operations for twenty of two hundred — 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 \(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 \(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. size_for() returns the smallest cassette reaching \(\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.

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 \(\sigma(s + b)\) with a single additive offset: the slope is measured and is right — \(\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. 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

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 \(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

\[\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 (log_ek()) — so it never enumerates a subset, and \(\binom{5000}{20}\) is not a number anybody was going to sum over. Adding \(\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 \(\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. 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 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, 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.

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.

survival() has modelled this since it was written: a unit responds only if its block is live and its own term fires, \(R_i = B_b \varepsilon_i\) with \(B_b \sim \mathrm{Bern}(q_b)\). --block-live is that \(q\). What it buys the objective is not a heuristic but a covariance — for two units on one allotype,

\[\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 \(\gamma (1 - q_b) p_i p_j / q_b\) to \(J_{ij}\) and nothing anywhere else. No rho, no overlap channel, and no parameter that is not the loss rate the designer stated. At \(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. overlap() returns the mean of its two or three channels, so with all three populated a same-allotype pair reaches \(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 \(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 \(H\): the loss coupling already prefers spread, and stacking a second diversity term inside the objective double-counts unless you mean it — the argument 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 \(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:

\[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 \(\log_2(1 + \mathrm{TPM}/c)\) on one floor, so their difference is a log2 ratio.

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 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. 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.

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, C_SUBCLONAL for detected-but-not-clonal, 0 for neither. A caller holding a CCF rather than a call can threshold it at 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: D_PASSENGER where neither fires, 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.

$ 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 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.

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 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: 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 score() already returns. escape is the one column that is not per-unit: 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.

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.

what sum p means	a level: expected responding units	an enrichment: how far 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.

prob_offset() gives the level, 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

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 mhcmatch.vector. See Safety, prior evidence, and what goes in the cassette.

API

Every function above, with its full docstring: API reference — mhcmatch.cassette.