Cassette composition#

Read Designing a cassette: what goes in, and what it is worth first — it is the practical page, and this one is the geometry underneath it: the block response model, the measured over-dispersion, coverage, redundancy, and what no weighted score can ever select.

The objective a cassette is actually judged on is whether at least one unit elicits a response — better, at least \(k\). Sorting by a score and keeping the top \(m\) maximises \(\sum_{i \in S} p_i\), the expected number of responding units. Those two objectives agree only if the units respond independently.

mhcmatch.portfolio is what the difference costs. It computes nothing new about a peptide — it takes the scores the rest of the library produces and says what a proposed set of them is worth. mhcmatch.vector.select() is the rule and Safety, prior evidence, and what goes in the cassette derives it, expected-yield formula included; this module is the diagnostics.

Why the naive objective is wrong#

Write a unit’s response as a conjunction of shared and private events:

\[y_i = B_{a(i)} \cdot G_{c(i)} \cdot \varepsilon_i, \qquad B_a \sim \mathrm{Bern}(\beta_a),\ G_c \sim \mathrm{Bern}(\gamma_c),\ \varepsilon_i \sim \mathrm{Bern}(\eta_i)\]

with \(B_a\) the event that allotype \(a\) is live in this donor and \(G_c\) the event that the mechanism the unit was selected on is live. Marginally \(p_i\) is what the ranker reports; the joint is not a product.

Saturation. If every unit shares one block pair, \(\Pr(\ge 1) \le \beta_a \gamma_c\) for every \(m\). The scalar objective grows without bound while the one a vaccine needs has a ceiling no further unit passes.

This is not a hypothetical. On the adjuvant TNBC mRNA vaccine trial of Sahin et al. (Nature 2026;651:1088–1096) — 13 patients, 216 assayed units between them — the intra-patient correlation is \(\rho = 0.124\) at \(p = 1.0 \times 10^{-3}\), 3.45x the binomial variance. Measure it on your own readout before assuming a value:

from mhcmatch import portfolio

m = [20, 5, 20, 20, 20, 20, 1, 10, 20, 20, 20, 20, 20]    # units assayed per patient
k = [8, 5, 2, 8, 6, 2, 0, 0, 3, 1, 2, 2, 2]               # of which positive

portfolio.dispersion(m, k)["ratio"]      # observed / independent-Bernoulli variance
portfolio.betabinom_rho(m, k)            # {'rho': ..., 'D': ..., 'p_value': ...}

Keep the zero-response patients. They carry most of the information about dispersion, and a minimum-pool-size filter deletes exactly them.

The observational read: a set nobody chose#

Everything above is about a cassette some rule selected. mhcmatch.portfolio.visibility() answers the same question for the set a tumour happens to present:

\[V(S) \;=\; 1 - \prod_{i \in S}\bigl(1 - p_i\bigr)\]

over every presented unit — weight one each, no live probability \(q\), no coupling, and neither \(\gamma\) nor \(\rho\). It agrees with p_at_least() at \(q = 1\), which is the no-loss case \(V\) already assumes.

The missing parameters are the modelling decision, not a gap. mhcmatch.cassette designs, and a design carries \(\gamma\), \(\rho\), the pairwise couplings and a greedy search over subsets. A tumour selects nothing in that sense: it is under selection towards escape, and that dynamics is not represented anywhere in this library, so there is no designer’s risk to price and nothing to maximise over. The command line keeps the two apart for the same reason — mhcmatch visibility is its own command and does not accept --gamma or --rho.

\(V\) saturates, which is what makes it a whole-tumour read rather than a count: once a few strong antigens are present, further units move it very little, so it does not inherit mutational burden the way \(\sum_i p_i\) does.

Read a cohort on the log-complement, not on \(V\). Saturation is not an edge case: over 7,261 recorded tumours \(V\) runs median 0.7976, q95 0.99999929 and q99 = max = 1.00000000, so at six decimals the whole top percentile prints as one indistinguishable block of 1.000000. All of that information sits in the complement, so visibility returns \(S = \sum_i \log(1 - p_i)\) as log1m — in nats, negative, linear in the units and unbounded below — and mhcmatch visibility emits it as log1m_vis beside vis_escape. Rank donors on that.

HLA loss of heterozygosity comes out in the same units. Pass block — one allotype label per unit — and losing allotype \(a\) deletes its terms. With \(S = \sum_i \log(1 - p_i)\) and \(S_a\) that sum over the allotype’s own units, visibility after the loss is \(1 - e^{S - S_a}\); the worst single loss is the allotype whose \(S_a\) is most negative, and loh_cost is how much visibility the tumour’s most protective allotype carries.

from mhcmatch import portfolio

p = [0.6, 0.1, 0.1]
portfolio.visibility(p)                            # {'visibility': 0.676, 'log1m': -1.127,
                                                   #  'n_units': 3}
portfolio.visibility(p, ["A", "B", "B"])           # + visibility_loh 0.19, lost_allotype 'A'

Selecting against blocks#

mhcmatch.vector.select() saturates a budget per block. The shipped default blocks on the allotype; pass block a key that pairs it with the mechanism a unit was selected on:

from mhcmatch import portfolio, vector
import numpy as np

# Z: one row per candidate, one column per objective, HIGHER IS BETTER on every column.
# A %rank is lower-is-better and has to enter as -log10(rank).
corner = portfolio.corner(Z, groups={0: "presentation", 1: "presentation",
                                     2: "recognition", 3: "abundance"})

sel = vector.select(units, n0=20.0,
                    block=lambda u: (u.allele, corner_of[u.peptide]))
sel.expected_yield          # computed against the partition the rule actually used
sel.per_block()             # where the budget went

n0 is per-block capacity and has no default; Safety, prior evidence, and what goes in the cassette owns that argument. Sweeping it retrospectively on 178 validated-immunogenic neoantigens puts the selection-layer optimum near 20, which is a starting point for the dose-matched trial and not a substitute for it.

Read the cassette back#

p     = np.array([...])            # calibrated per-unit probabilities, in cassette order
block = np.array([...])            # block index per unit
q     = 0.5                        # probability a block is live in this donor

ge1 = portfolio.p_at_least(p, block, q, k=1)
portfolio.survival(p, block, q)    # the whole tail: element k is P(X >= k)
portfolio.n_effective(p, ge1)      # how many INDEPENDENT units this cassette is worth

p_at_least() refuses a marginal it cannot represent: a unit cannot respond more often than its own block is live, so p_i > q raises rather than silently clipping.

Both are exact, and cheaply so. \(X = \sum_b B_b S_b\) with \(S_b\) a Poisson binomial over block b’s units, and \(B_b S_b\) has pmf \((1-q_b)\delta_0 + q_b\,\mathrm{pmf}(S_b)\). The blocks are independent, so the pmf of X is the convolution of those — \(O(B m^2)\), with no \(2^B\) enumeration over live sets and no Monte Carlo. The 200,000-draw sampler this replaced agreed to its own noise; the convolution has none.

Composing to quotas#

A cassette is usually specified as quotas, not as a top-m: eight class-I slots of which at least two should respond, four class-II of which one, three non-conventional of which one. That is what compose() fills.

comp = portfolio.compose(
    units,
    {"mhc1": (8, 2), "mhc2": (4, 1), "nonconventional": (3, 1)},
    q=0.5,                          # P(a block is live) in this donor
    universe=donor_allotypes,       # the donor's DISTINCT allotypes -- see homozygosity below
    weight_evenness=0.0)

comp.arms["mhc1"]["p_at_least"]     # attained P(>= 2) on the class-I arm
comp.joint                          # every quota met at once
comp.coverage                       # Gini and H/Hmax over class-I allotypes
comp.trace                          # one row per greedy step, with the gain it bought

Or from the command line, which emits both — the composed cassette and the same slot budgets filled by score alone — so the comparison is laid out on your own candidates rather than asserted:

$ mhcmatch cassette build --candidates units.tsv --n0 20 \
      --quota 'mhc1=8:2,mhc2=4:1,nonconventional=3:1' --block-live 0.5 \
      --alleles "$(cat donor.hla)" \
      --fasta cassette.faa --fasta-nt cassette.fna --map cassette.map.tsv

With a quota, --fasta and --fasta-nt carry two records, cassette_composed and cassette_topk; --map describes the composed one. Without a quota each carries the single cassette record it always did.

Note

--block-live is a ceiling on every unit’s own p. A block is an allotype, so a unit cannot respond more often than its allotype is live; a candidate whose p exceeds q makes the marginal unrepresentable and survival() refuses rather than silently clipping. Feeding rank’s p_response at a pool prevalence of a few per cent leaves plenty of headroom under the default q = 0.5; feeding a raw sigmoid of the log-odds does not.

Note

The same q is what cassette select --block-live prices HLA loss with, and it is documented once, here. What differs is where it lands: compose() uses it inside P(X >= target), while goal_energy() uses it to add the covariance a lost allele implies, \(\gamma (1 - q_b) p_i p_j / q_b\), to same-allotype pairs of the coupling. Both read the same block model; neither fits anything. Designing a cassette: what goes in, and what it is worth has the derivation and what it is worth measured.

The arms are disjoint on purpose. A frameshift neoepitope is presented on MHC-I, so if it counted toward both budgets, “at least one non-conventional epitope responds” could be satisfied for free by the class-I arm and would never change a cassette. Charged to its own arm (default_arm() reads Unit.kind), it has to earn a slot. It earns one because it fails differently: a non-conventional product is foreign over a stretch rather than at one position, so whatever makes the missense arm miss — a wrong wild type, a tolerised residue — does not make this arm miss.

Why this is not top-m, shown rather than asserted#

Nine candidates: five strong ones all restricted to A*02:01, four weaker ones spread over four other allotypes. Four slots, target “at least one responds”, \(q = 0.5\).

rule

P(≥ 1)

Gini

H/Hmax

allotypes taken

top-4 by score

0.4806

0.800

0.000

A*02:01 ×4

compose

0.6550

0.200

0.861

four distinct

Nothing was told to diversify. The spread falls out of the objective, because a block that is already represented contributes less than a fresh one.

And the instinct is wrong when the target is \(k \ge 2\). Two units in one block need that one block live; two units in two blocks need both live, which at \(q = 0.5\) costs a factor of two. On the same pool at target 2, compose concentrates — P(≥ 2) 0.3829 concentrated against 0.2929 spread — and it is right to. “Diversify” is a heuristic for \(\Pr(\ge 1)\); the tail probability is the thing, and it does not always agree.

Coverage evenness, and homozygosity#

When spread matters for reasons the response model does not price — manufacturing risk, an uncertain genotype, provisional typing — weight_evenness adds \(w\,\Delta(H/H_{\max})\) to the objective. It costs, and the cost is reported: on the pool above at target 2, weight_evenness=0.2 buys H/Hmax 0.000 → 0.646 for P(≥ 2) 0.3829 → 0.2929.

coverage() takes universe — the donor’s distinct allotypes — and this is the whole point when the donor is homozygous. A patient homozygous at B has five distinct class-I allotypes, not six, so a cassette spread evenly over five is perfectly even; scoring it against a denominator of six would report a genotype as a design flaw.

What a scalar score cannot select#

For any \(\beta \ge 0\), top-\(m\) by \(\beta^\top z\) selects only candidates on the upper convex hull of the objective cloud. Pareto-efficiency is necessary for reachability but not sufficient — measured on 178 validated-immunogenic neoantigens, 45 of the 161 Pareto-efficient ones are ranked first by no non-negative weighting at all.

front = portfolio.pareto_front(Z)          # non-dominated
portfolio.linearly_supported(Z, i)         # exact, by LP: is it on the hull?
portfolio.crowding_distance(Z[front])      # NSGA-II tie-break within a front

That limit belongs to the weighted sum, not to scalarization. chebyshev_score() reaches the whole front, and the optimal weights for a given candidate are closed form:

d = (Z.max(0) + 1e-6) - Z[i]
w = (1.0 / d) / (1.0 / d).sum()            # equalises the weighted shortfalls
portfolio.chebyshev_score(Z, w).argmax()   # == i, even inside the hull

A gradient-boosted score is in the same position, and on real candidate pools it is worth considerably more than either: a boosted classifier on the same eleven objectives captured 136 of 178 validated neoantigens at a 30-unit budget against a fitted linear score’s 113.

What none of them escapes. Top-\(m\) by any pointwise score maximises \(\sum_{i \in S} s_i\), a modular set function, while \(\Pr(\ge k \mid S)\) is submodular whenever two units share a block. The limitation is a property of the selection rule, not of the scorer, so it cannot be fitted away at any model capacity.

Limits#

  • The mechanism corner from corner() is a proxy for a latent variable: it says which axis a candidate stands out on, not why it works.

  • \(\Pr(\ge k)\) treats a response as binary at the assay’s threshold; magnitude is discarded.

  • An absolute \(\Pr(\ge k)\) inherits the calibration of its inputs, and corpus prevalence varies by four orders of magnitude across published screens. The ordering results do not.

  • linearly_supported and betabinom_rho need SciPy, which is not a hard dependency.

API#

Full signatures: mhcmatch.portfolio in API reference.