Cassette design

Selecting units, assembling a construct, and what a portfolio is worth.

mhcmatch.cassette module

Choosing the units of a cassette, and scoring one that already exists. The narrative version, with the derivation and the measured numbers, is Cassette design.

Choosing the units of a cassette, and scoring one that already exists.

Two operations, and they are not the same question. select() is handed a donor’s whole candidate pool and returns the k units to manufacture. score() is handed a cassette somebody already built — possibly from another donor, possibly of another size — and returns numbers those cassettes can be compared on. Everything here sits above mhcmatch.rank, which scores one peptide against one allele, and beside mhcmatch.portfolio, which holds the response model and the objective geometry; the assembly step that follows selection is mhcmatch.vector.

The objective is derived from the design goal, not fitted to an outcome cohort. A cassette’s job is that several of its units elicit a detectable response. Writing 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

E[B] = sum_i p_i Var[B] = sum_i s_i^2 + 2 sum_{i<j} rho_ij s_i s_j, s_i = sqrt(p_i (1 - p_i))

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 — per-patient response counts are over-dispersed in every assayed vaccine cohort that has been measured (RHO_ASSAYED). 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 [ p_i - (gamma/2) s_i^2 ] - 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. See goal_energy().

``gamma`` is a statement about a unit, so it is divided by the design effect. For a cassette of k units with mean response pbar and mean pair correlation rho, E[B] = k pbar and Var[B] = k pbar qbar [1 + rho (k - 1)], so

H = k pbar { 1 - (gamma/2) qbar [1 + rho (k - 1)] }

The brace is the certainty-equivalent worth of the average unit, and with gamma held fixed it falls as the cassette grows — the mean of a correlated count is linear in k and its variance is quadratic, so one stated gamma is a stricter trade at every larger size. It reaches zero at

k* = 1 + (2 / (gamma qbar) - 1) / rho

and past k* every unit is a net cost, so the objective prefers a worse unit to a better one and selection inverts. At gamma = 1 and the measured rho = 0.091 that is k* = 16.1 on the TESLA pools and 17.4 on HiTIDE — inside the twenty-unit cassette a trial ships. risk_aversion() divides gamma by the design effect 1 + rho (k - 1), which leaves the brace at 1 - (gamma/2) qbar for every k. One stated preference then means one trade at every size, and no size inverts. It is arithmetic on rho and k, both known before any unit is chosen; no outcome enters it.

Two things that make this different from sorting the candidate list.

Top-m by any pointwise score maximises a modular set function, and H is not modular whenever two units share a mechanism. That is a property of the selection rule, not of the scorer, so no better ranker fixes it — see mhcmatch.portfolio for the measurement.

The calibration offset decides what is being reported. mhcmatch.rank.probability() anchors the mean of the batch it is handed. Handed one donor at a time it pins every donor’s pool mean to the declared prevalence, whatever their pool: 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 2.75e-17. Read as a probability, that number is not one. prob_offset() fits the offset over a batch that does not move, which is what makes two donors comparable; group_offsets() fits one per group, which turns the same sum into how far a donor’s chosen units sit above their own background. Both are useful. They are different quantities, and which one you want is a decision, not a default.

betabinom_rho needs SciPy, which is not a hard dependency; nothing else here does.

mhcmatch.cassette.KMER = 3

k-mer width for the sequence-overlap channel. 3 is what the corpus kernel uses.

mhcmatch.cassette.KAPPA = 7.4635

Mean number of distinct KMER-mers a peptide contributes, measured over the distinct peptides of the neoantigen corpus. Only a scale — halving it doubles the sequence channel and leaves the energy untouched once goal_energy() renormalises — so shared k-mers read as “one peptide’s worth” rather than as a raw count.

mhcmatch.cassette.GAMMA = 1.0

Risk aversion of the objective, in units of “one variance is worth one expected unit”. gamma = 1 says a designer trades one unit of variance in the responding-unit count for one unit of its mean. A stated design preference, not a fitted quantity, and not swept.

It is stated per unit, so select() and score() pass it through risk_aversion() before use — see the module docstring for why a cassette-wide gamma cannot mean the same thing at two sizes. Passing gamma= explicitly bypasses that and uses the number given, which is how the cassette-wide arm stays reproducible.

mhcmatch.cassette.RHO_ASSAYED = 0.091

Default intra-cassette response correlation. Measure your own with mhcmatch.portfolio.betabinom_rho() before relying on this one.

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.45x the independent-Bernoulli variance), 0.091 on IVAC MUTANOME (75 of 125 units, 13 patients, 1.8x), 0.024 on TESLA and 0.010 on HiTIDE. The default is 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 — so it is the only estimate that is not conditioned on which units somebody chose to assay. The dispersion is scale-dependent: a screening pool spanning every allotype shows none, because a pool that wide averages its blocks out. A cassette cannot.

mhcmatch.cassette.MAX_POOL = 2000

Pool size above which select() trims to the top-scoring candidates before building the coupling matrix. J is dense n x n, so an untrimmed 24,366-candidate pool would ask for 4.4 GB to choose twenty units from. Pools are small in practice — median 24 source proteins per TCGA donor, 207 at the 95th percentile — so this bites only in the tail, and the trim keeps the units any objective would have ranked first. select() records trimmed when it fires.

mhcmatch.cassette.tcr_face(peptide, cls='mhc1', register=None)

The peptide with its MHC-facing positions deleted — the residues a TCR can actually read.

Delegates the split to mhcmatch.mimicry.masks(), so the face is the same one every other channel in the package uses: mhcmatch.complement.ANCHORS for class I, and the floating 9-mer core’s P1/P4/P6/P9 for class II. It is deliberately not seqtree’s layout.DEFAULTS["mhc1"], which masks P2 and POmega only.

Deleting the columns is what makes a masked alignment possible at all. No aligner in the stack scores a masked dense matrix – seqtree’s PositionalMatrix reaches only the search engine, only with zero indels, and is silently ignored when its width does not equal the query length. In an ungapped comparison masking a position and deleting it are the same operation, so slicing first and aligning after is exact rather than an approximation.

>>> tcr_face("SIINFEKLL")        # anchors P1-P3 and POmega-1, POmega removed
'NFEK'

Parameters:

• peptide (str)

• cls (str)

• register (int | None)

Return type:

str

mhcmatch.cassette.sequence_overlap(peptides, cls='mhc1', registers=None, mask='face', h=None, threads=1)

BLOSUM-graded pairwise sequence similarity in [0, 1], one (n, n) matrix.

This replaces a channel that was measured to be a duplicate detector. Counting exactly shared 3-mers is zero on almost every real pair: 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. A channel that is zero on 97% of pairs cannot order them, which is why rho_seq does not hold its sign across cassette sizes. 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. Here they score 6 and 19 on the BLOSUM distance.

Built on seqtree.pairwise.dist_matrix(), which returns the sequence-level Gram transform d(a, b) = s(a,a) + s(b,b) - 2 s(a,b) – non-negative, symmetric and zero on the diagonal, so it is a distance rather than a raw score. That symmetry is not decoration: goal_energy() halves the pair sum assuming it, and the obvious alternative (mhcmatch.mimicry.blosum62_kernel() at normalise=True) subtracts the query’s self-score and is therefore asymmetric.

mask="face" compares TCR faces (tcr_face()); "full" compares whole peptides.

The face is the default because the full peptide is confounded by HLA, and by how much is measured rather than argued: on one TESLA donor’s 83 units the correlation between this channel and “these two units share an allele” is +0.0437 on the face against +0.3297 on the whole peptide. Anchor residues are what the allele selects for, so aligning them makes the sequence axis a second reading of the allotype axis — and the whole point of having four axes is that they fail independently. Gaps are handled by the aligner, so units of different length are compared rather than scored zero – which matters, since only 38% of TESLA and 46% of HiTIDE within-donor pairs are of equal length.

h is the bandwidth of sim = exp(-d / h), defaulting to the median off-diagonal distance. It needs no calibration and never did: goal_energy() divides the off-diagonal mean to 1, so any global scale cancels and only the pair ordering survives. That is also why swapping this channel in leaves KAPPA, RHO_ASSAYED and GAMMA untouched.

Speed is not a consideration at cassette scale: 26.1 M pairs/s measured on the HiTIDE pool (1,558 units, 92.9 ms for the whole matrix), in C++ with the GIL released.

A dense matrix is the right shape here, and a thresholded index search is not — measured, because the opposite is the natural guess. Building a seqtree.Index over the faces and batch-querying is 20x faster (4.3 ms against 92.9 ms on that pool), but at max_subs=2 it returns 0.65% of pairs, which is the same near-binary channel this function exists to replace. The reason is the distance distribution: on one donor’s 151 units the off-diagonal BLOSUM distances run min 8, median 65, max 120, so every pair sits within 2x the median and carries similarity above 0.05, with half above exp(-1). A trie prunes nothing against a threshold loose enough to be faithful. The index is the right tool at pool-wide or TCGA scale, where n is 10^5 and the dense matrix cannot be formed at all; it is the wrong one here.

>>> import numpy as np
>>> o = sequence_overlap(["GILGFVFTL", "GILGFVFTL", "NLVPMVATV"])
>>> float(o[0, 1])            # identical peptides
1.0
>>> bool(o[0, 2] < o[0, 1])   # unrelated scores lower
True

Parameters:

• cls (str)

• mask (str)

• h (float | None)

• threads (int)

Return type:

numpy.ndarray

mhcmatch.cassette.not_worse(sel, ref, p, J, exact_max=24)

P(B(S) >= B(R)) — the chance this cassette catches at least as much as the reference.

This is the constraint the v2 objective is posed under. 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. That degeneracy is the design freedom: rather than trading expected responders away for variance, the rule keeps every set that is, with stated probability, no worse than simply sorting, and spends what is left on diversity.

Units in both sets cancel exactly, because they are the same random variable and not merely identically distributed. With A = S \ R and C = R \ S,

D = B(S) - B(R) = B(A) - B(C)

so only the symmetric difference carries any variance at all. That is not an approximation and it is what makes this cheap: the objective shares 17 of 20 slots with the sort on the TESLA pools, so the difference is usually three units against three.

E[D] = sum_A p_i - sum_C p_j and Var[D] = Var[B(A)] + Var[B(C)] - 2 Cov(B(A), B(C)), every covariance read off the same J the coupling already carries — including the exact HLA-loss term when block_live priced one. So the correlation structure still enters, as the scale of the slack rather than as a penalty subtracted from the yield.

Two evaluators, and which one ran is reported by the caller rather than guessed at. Up to exact_max units in the symmetric difference the Poisson-binomial is convolved exactly under an independence approximation within each side; beyond that a normal approximation with a continuity correction is used. Both read the same first two moments.

Returns 1.0 when the two sets are equal — nothing has been given up, so the guarantee is total.

>>> import numpy as np
>>> p = np.array([0.5, 0.5, 0.5, 0.5])
>>> J = np.zeros((4, 4))
>>> float(not_worse([0, 1], [0, 1], p, J))          # the same set
1.0
>>> bool(not_worse([2, 3], [0, 1], p, J) > 0.4)     # a swap between equals
True

Parameters:

exact_max (int)

Return type:

float

mhcmatch.cassette.diversity(sim, sel, how='minmax')

How little a chosen set shares, given per-axis (n, n) similarities. Higher is better.

sim is a mapping {axis name: (n, n) matrix} — allotype, expression, chemistry, sequence. Each axis is reduced to its mean pair similarity over the chosen set, and then:

• "minmax" — 1 - max over axes. A cassette is undone by its worst shared failure mode, not its average one, and averaging is measured to dilute: seven channels (feature-only) lost to two (select-dom) on the same pools, because every axis added makes every existing axis count for less.

• "mean" — 1 - mean over axes. What the v1 overlap did.

Each axis is standardised to unit off-diagonal mean over the pool before the reduction. A max over axes on different scales reads whichever axis has the largest raw spread rather than the one that is actually shared, and the axes here are a 0/1 indicator, a log2 abundance, a Rose propensity and a BLOSUM kernel. "mean" needs it less and gets it too, so the two aggregations are compared on one normalisation rather than on two.

>>> import numpy as np
>>> a = {"x": np.array([[0.0, 1.0], [1.0, 0.0]]), "y": np.zeros((2, 2))}
>>> float(diversity(a, [0, 1], "minmax"))   # the worst axis is fully shared
0.0
>>> float(diversity(a, [0, 1], "mean"))
0.5

Parameters:

how (str)

Return type:

float

mhcmatch.cassette.normalise_axes(sim)

Each axis divided by its own off-diagonal mean, so a max over axes is not a scale contest.

Parameters:

sim (dict)

Return type:

dict

mhcmatch.cassette.swap_for_diversity(sim, p, J, ref, pi=0.5, how='minmax', rounds=8, codes=None, cap=None, must=())

The v2 rule: start from the ranked list, then trade slots for diversity while it stays safe.

ref is the reference cassette — the top-k sort, which maximises the expected number of responding units by construction. Every step swaps one chosen unit for one unchosen one, taking the exchange that raises diversity() most among those keeping not_worse(S, ref) >= pi. It stops when no admissible swap improves diversity.

The sort is therefore the floor, not the rival. The rule cannot wander into a set that is probably worse than sorting, because pi is checked against ref itself at every step rather than against the previous iterate — checking against the iterate would let a chain of individually-safe steps drift arbitrarily far.

The diversity scan is vectorised and exact. For one axis, writing c_i = sum_{j in S} M[ij] (one matrix-vector product per pass), swapping u out for v in changes that axis’s pair sum by c_v - c_u - M[v, u], so the whole k x (n - k) table of deltas is an outer difference rather than k(n-k) submatrix sums.

not_worse is evaluated lazily, in descending order of the diversity it would buy, and the first admissible candidate is taken. Most passes evaluate it once or twice: the probability constraint is slack near the sort and only binds once several slots have moved.

codes / cap / must are the manufacturing constraints greedy() already takes, applied here as a feasibility mask on the same scan.

Parameters:

• sim (dict)

• pi (float)

• how (str)

• rounds (int)

• cap (int | None)

Return type:

list

mhcmatch.cassette.build_axes(peptides, alleles=None, expression=None, physchem=None, coexpr=None, presented=None, cls='mhc1', registers=None, mask='face')

The four ways a pair of units can share a failure, as {axis: (n, n)}, normalised.

One matrix per mechanism, not per column: expression and physchem are (n, d) blocks whose columns are averaged into one axis each, because expr_lvl and expr_norm are two readings of one thing and counting them separately would let the number of columns decide how much abundance matters. That is the dilution diversity() is built to avoid, and it should not be reintroduced one level up.

• allotype — 1[a_i = a_j], or the graded presented-allele overlap when presented is given (allotype_overlap()). The one axis that is discrete, because HLA is.

• expression — (n, d); the source gene’s abundance, its healthy-tissue level and their difference. coexpr is folded in here rather than made its own axis: GTEx tissue-profile similarity is a statement about abundance, and mhcmatch.expression.coexpression() already returns it as a matrix.

• physchem — (n, d); TCR-face burial and charge.

• sequence — sequence_overlap(), BLOSUM-graded over the TCR face.

Every axis is put on unit off-diagonal mean by normalise_axes() before it is returned, so a minmax reduction compares shared-ness rather than scale. An axis the caller cannot supply is simply absent, and which axes were built is part of the result.

Parameters:

• cls (str)

• mask (str)

Return type:

dict

mhcmatch.cassette.prob_offset(scores, prevalence)

The additive offset putting the whole batch’s mean response probability at prevalence.

Same rule and same bisection as mhcmatch.rank.probability() — pick b so that mean_i sigma(s_i + b) = pi — fitted on a batch that does not change, which is the whole difference. The left side is strictly increasing in b from 0 to 1, so the root exists, is unique, and bisection finds it.

Fit this once, over every donor you intend to compare, and then hold it. Re-solving inside each donor is what pins every pool mean to prevalence and makes the resulting sum a statement about nothing. Use group_offsets() when the per-donor quantity is what you actually want.

>>> import numpy as np
>>> b = prob_offset(np.array([3.0, 0.0, -3.0]), 0.25)
>>> round(float((1 / (1 + np.exp(-(np.array([3.0, 0.0, -3.0]) + b)))).mean()), 6)
0.25

Parameters:

prevalence (float)

Return type:

float

mhcmatch.cassette.group_offsets(scores, group, prevalence)

One offset per group, all groups bisected simultaneously.

Same rule as prob_offset() solved for every group at once: b is a vector, the sigmoid is evaluated on the whole column each iteration, and the per-group means come from a single np.bincount. 200 iterations over one array beats one Python-level bisection per group — measured on 7,261 TCGA donors over a 465,343-row column, where the loop form is the bottleneck and this is not.

group is an integer code per row, 0 .. n-1. Returns one offset per code, so scores + offsets[group] is the shifted column.

What this buys is an enrichment, not a level: every group’s mean probability becomes prevalence by construction, so what survives is how far a chosen subset sits above its own group’s background. That is a real and useful quantity — on TCGA it reads out more strongly against immune infiltrate than the pool-anchored level does — but it is not a probability and two groups’ offsets are not comparable.

Parameters:

prevalence (float)

Return type:

numpy.ndarray

mhcmatch.cassette.overlap(peptides, alleles=None, strength=None, features=None, coexpr=None, allotype_graded=None, sequence=None, profile=None, genes=None, kmer=3, combine='mean')

Mechanistic pair overlap in [0, 1]: how much two units share a way of failing.

The mean of whichever channels the caller can populate. Which channels were available is part of the result and should be reported with it — a trial that publishes no per-patient genotype has one fewer than one that does.

• allotype (alleles) — 1 if two units are restricted by the same class-I molecule. Two units on one molecule compete for the same presentation and the same precursor niche, and are lost together if that allele is. Pass allotype_graded — an (n, n) matrix from allotype_overlap() — to use the graded form instead: overlap between the two units’ presented-allele vectors rather than equality of one credited label. It replaces this channel rather than joining it, because they are two readings of one mechanism and averaging them would count presentation twice.

• gene (genes) — 1 if two units come from the same source gene. The antigen-loss twin of the allotype channel: one deletion or one silencing event at a locus removes every unit that locus supplies, exactly as one loss-of-heterozygosity event removes every unit on an allotype. It is the pair term a tumour’s escape route needs and the allotype channel cannot express, two mutations in one gene being restricted by different molecules.

• sequence (always) — shared distinct kmer-mers, in units of KAPPA, clipped at 1. Two units that look alike draw on one repertoire, so the second buys less than its score claims. Pass ``sequence`` — an ``(n, n)`` matrix, normally from :func:`sequence_overlap` — to use the BLOSUM-graded TCR-face form instead. The exact-3-mer count is zero on 97.3% of real within-donor pairs and cannot grade a conservative substitution against a radical one; it is kept as the default only so recorded results reproduce.

• dominance (strength) — closeness on the score axis. This is the score talking to itself: two units are coupled for scoring alike, which is not a mechanism, and the pairwise statistic it corresponds to (rho_dom) fits attractive on the observational arm — where greedy() loses its 1 - 1/e guarantee, which holds only for repulsive couplings. It is kept, and it is now optional rather than unconditional: pass strength=None to drop it.

• features — an (n, d) array of per-unit scalars, one channel per column, each on the same 1 - |f_i - f_j| / span kernel as dominance. This is how chemistry and expression reach the pair term: C_phys_buried, C_phys_charge, expr_lvl, expr_norm and the selectivity delta are all per-unit scalars the ranker already computes and the objective has never seen. A column that is entirely non-finite contributes a zero channel rather than raising.

• profile (profile) — a symmetric (n, n) matrix from profile_overlap(): how much two units owe their scores to the same fitted terms. It is the one channel that reads the model’s own decomposition rather than a proxy for it, and it subsumes dominance — pass strength=None with it, or the same idea enters twice, once as a mechanism and once as the score talking to itself.

• coexpr — a symmetric (n, n) matrix already in [0, 1], averaged in as one further channel. Co-expression is a property of a pair of source genes and cannot be written as |f_i - f_j| on any per-unit scalar, which is why it enters as a matrix; mhcmatch.expression.coexpression() builds one from the GTEx tissue panel.

Vectorised: the sequence channel is one float32 matmul over a k-mer incidence matrix rather than n^2 set intersections.

At ``features=None, coexpr=None, profile=None`` the result is bit-identical to every cassette built before they existed, because the mean is then over exactly the channels it was over before.

>>> import numpy as np
>>> o = overlap(["AAAAAAAAA", "CCCCCCCCC"], features=np.array([[0.0], [1.0]]))
>>> float(o[0, 1])
0.0

Parameters:

• kmer (int)

• combine (str)

Return type:

numpy.ndarray

mhcmatch.cassette.pair_stats(peptides, alleles=None, strength=None, kmer=3)

The three pairwise set statistics of one cassette, each by its exact closed form.

Every one is a sum over pairs divided by C(k, 2), and none is computed pair by pair: a same-category pair count is sum_c C(c, 2) over category occupancies, and the mean absolute difference is a rank-weighted sum of the sorted values. Both are O(k log k), which is what makes them affordable inside a selection loop rather than only after one.

rho_hla — share of pairs on the same allotype. rho_seq — shared k-mer mass per pair, in units of KAPPA. rho_dom — mean absolute strength gap; the sum of |z_i - z_j| over pairs is the Gini numerator, so the entropy/Gini family enters without a surrogate.

>>> s = pair_stats(["AAAAAAAAA", "AAAAAAAAA"], alleles=["A", "A"])
>>> round(s["rho_hla"], 6)
1.0

Parameters:

kmer (int)

Return type:

dict

mhcmatch.cassette.risk_aversion(k, rho=0.091, gamma=1.0)

gamma restated per unit: gamma / (1 + rho (k - 1)).

The design effect 1 + rho (k - 1) is how much a correlated cassette’s count varies above an independent one of the same size — the beta-binomial factor mhcmatch.portfolio. betabinom_rho() estimates rho from. Dividing by it is what makes gamma a preference about a unit rather than about a cassette; the module docstring derives it and gives the size k* at which the undivided form inverts the objective.

>>> round(risk_aversion(1), 4)
1.0
>>> round(risk_aversion(20, rho=0.091), 4)
0.3664

Parameters:

• k (int)

• rho (float)

• gamma (float)

Return type:

float

mhcmatch.cassette.resolve_restriction(alleles, cls='mhc1', collapse=False)

Restriction cells resolved to allotypes: one label per unit, and the presented set.

A screen that did not resolve which of a donor’s alleles restricts a candidate writes the whole genotype into the cell — 'HLA-A*01:01,HLA-A*03:01'. That string is not an allele name, and overlap() couples two units on the allotype axis by string equality of the label it was handed, so such a unit equals no other label in the pool: it couples to nothing, reads as a private allotype of its own, and can hold a slot beside both of its own constituents. Measured on the real pools, the cell is composite on 30 of TESLA’s 736 rows and 53 of HiTIDE’s 1,563, and it put 13 and 12 distinct restriction cells on donors carrying 6 and 5 allotypes.

Five sites keyed that raw string — the similarity channel, pair_stats’ rho_hla, the block index and per-block q, the named-column q lookup, and the coverage floor’s np.unique — so resolving in one of them would have left four wrong. This function is the single boundary they all now sit behind, which is why it lives here rather than in a caller: a pool is clean until somebody hands over a screen that did not resolve its restriction.

mhcmatch.rank.split_alleles() is the splitter and it already existed — it splits on [,;/|], de-duplicates in input order, preserves the spelling as supplied, and drops names the pseudosequence tables do not know. mhcmatch.pseudoseq.normalize_allele() already states the governing rule: a cell naming several alleles is a genotype, not an allele.

Returns

block

one representative label per unit — the first allele the cell resolves to — which stays what the block index, the q lookup, pair_stats and the coverage floor key on.

presented / presented_alleles

an (n, A) 0/1 matrix and its column names, ready for goal_energy(). This is the honest reading of a composite cell: a unit with several routes to the surface is exactly what presented was built to express, and the set form reduces to the single-block form identically, so a pool of singletons reproduces the label reading bit for bit. None under collapse.

resolved

per unit, whether the cell named anything at all.

composite / unresolved

how many cells named several alleles, and how many named none.

A cell that resolves to nothing keeps its unit and loses its allotype. The raw cell stays as the block label, so block_live’s “a label absent from the mapping is never lost” default holds and the unit carries no loss coupling; it gets its own presented column for the same reason. A unit whose restriction is unknown is not a unit that shares a restriction, and the caller is told how many there were — a lookup that returns nothing and a lookup that returns the wrong thing fail the same way when the layer below drops silently.

collapse=True keeps one label per unit and builds no matrix, which is the pre-promiscuity reading and the lossy one. It is here so a recorded result computed that way stays reproducible.

>>> r = resolve_restriction(["HLA-A*01:01,HLA-A*03:01", "HLA-A*01:01"])
>>> r["block"], r["composite"]
(['HLA-A*01:01', 'HLA-A*01:01'], 1)
>>> r["presented_alleles"], r["presented"].tolist()
(['HLA-A*01:01', 'HLA-A*03:01'], [[1.0, 1.0], [1.0, 0.0]])

Parameters:

• cls (str)

• collapse (bool)

Return type:

dict

mhcmatch.cassette.selectivity_delta(expr_lvl, expr_norm)

expr_lvl - expr_norm: tumour-over-normal selectivity in log2-fold, 0 where unknown.

Both terms are log2(1 + TPM/c) on one floor (mhcmatch.rank.expr_level() and mhcmatch.rank.expr_norm_level()), so their difference is a log2 ratio and a unit of it is one doubling of tumour abundance over the same gene’s healthy-tissue median.

A row missing either term takes 0, not nan. nan would propagate into the argmax and silently delete the candidate; 0 is what “no information about this candidate’s selectivity” actually means, and it leaves the unit ranked on everything else. How many rows took it is a number the caller should report – np.isnan(...).sum() on the inputs – rather than absorb.

Return type:

numpy.ndarray

mhcmatch.cassette.goal_energy(p, sim, rho=0.091, gamma=1.0, block=None, block_live=1.0, presented=None, presented_alleles=None)

The mean-variance objective as a field and a coupling: H(S) = sum h - sum_{i<j} J.

See the module docstring for the derivation. h_i = p_i - (gamma/2) s_i^2 and J_ij = gamma rho_ij s_i s_j with s_i = sqrt(p_i (1 - p_i)).

``block_live`` prices HLA loss, and its coupling is derived rather than fitted. Under the response model mhcmatch.portfolio.survival() already uses, a unit responds only if its allotype is live and its own term fires — R_i = B_b eps_i with B_b ~ Bern(q_b), so p_i = q_b r_i. Two units on the same allotype therefore covary by

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 units on different allotypes do not covary at all. So losing an allele contributes exactly gamma (1 - q_b) p_i p_j / q_b to J_ij on same-block pairs and nothing anywhere else. No rho, no overlap heuristic, and no parameter that is not the stated loss rate.

That is worth entering as itself rather than as a channel. overlap() returns the mean of 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. The loss rate is a number a designer states; it enters at full strength or not at all. The sequence and dominance channels keep carrying the residual rho exactly as before.

Promiscuity, when the caller can supply it. The block form above credits each unit to one allotype, so it charges a unit presented by three of a donor’s six the full loss of one. Pass presented — an (n, A) 0/1 matrix of which allotypes present each unit, columns in the order _block_q() returns q — and the same derivation runs over sets. Unit i has a live route iff any allotype in A_i survives, so with L_a ~ Bern(q_a) independent

Q_i = 1 - prod_{a in A_i} (1 - q_a) P(S_i, S_j) = 1 - prod_{A_i}(1-q) - prod_{A_j}(1-q) + prod_{A_i union A_j}(1-q) Cov(R_i,R_j) = (p_i p_j / (Q_i Q_j)) [P(S_i, S_j) - Q_i Q_j]

which is exact, closed-form and still costs no parameter that is not the stated loss rate.

It reduces to the single-block form exactly. With A_i = A_j = {b}: Q = q_b and P(S_i, S_j) = q_b, so Cov = (p_i p_j / q_b^2)(q_b - q_b^2) = (1 - q_b) p_i p_j / q_b — the expression above, term for term. So presented is an extension of the shipped model and not a second one, and a one-hot presented reproduces block.

At ``block_live = 1`` the added term is identically zero, so every cassette built before this existed is reproduced bit for bit — with or without presented, since every Q_i is then 1 and the bracket vanishes.

Where ``rho_ij`` comes from, and why it is not a fit. The scalar rho is the measured mean intra-cassette correlation (RHO_ASSAYED). It is spread over pairs in proportion to sim, then renormalised so the pool’s mean pair correlation is exactly rho. Nothing here is estimated from an outcome: sim is arithmetic on the peptides, rho is one number measured on published per-unit assays, gamma is stated.

Parameters:

• p (per-unit response probability, one entry per candidate.)

• sim (symmetric n x n overlap in [0, 1]; the diagonal is ignored.)

• rho (mean intra-cassette response correlation.)

• block (what a unit is lost with — the allotype, one label per candidate. Required for) – block_live; ignored without it.

• presented (optional (n, A) 0/1 matrix, “does allotype a present unit i”, columns in) – np.unique order over block. Supersedes the one-label-per-unit reading of block for the loss coupling only; block still supplies the labels and the q lookup. A unit whose row is all zero is taken to be presented by its block label alone, which is the pre-promiscuity reading and never a silent zero-survival unit.

• block_live (q, how often each block survives. A scalar, or {label: q} for a per-locus) – loss rate. 1.0 (the default) is “nothing is ever lost”.

• gamma (float)

Returns:

• (h, J) — field of length n, coupling n x n with a zero diagonal, such that

• H(S) = h[S].sum() - J[np.ix_(S, S)].sum() / 2.

• >>> import numpy as np

• >>> h, J = goal_energy([0.5, 0.5], np.array([[0.0, 1.0], [1.0, 0.0]]), rho=0.1)

• >>> float(np.diag(J).sum())

• 0.0

mhcmatch.cassette.greedy(h, J, k, codes=None, cap=None, must=(), weight_coverage=0.0)

Argmax of H over size-k subsets, greedily. Monotone submodular where J >= 0.

One pass per step over a running marginal, so selecting k of n is O(kn) rather than O(C(n, k)) — about 4,000 operations for twenty of two hundred. Ties broken by index, so the result is a function of the data alone and two runs agree.

codes is an integer block index per candidate and turns on the two manufacturing constraints, which are deliberately not objective terms — the coupling already prefers spread, and a second diversity term inside H double-counts unless it is meant (the argument mhcmatch.portfolio.compose() already makes for weight_evenness):

• cap — no block may hold more than this many units. A feasibility mask on the same loop, not a second search.

• must — block codes that each get one unit before the free slots are filled, where the pool supplies one. A block the pool cannot supply is skipped rather than raising: that is a fact about the donor’s candidates, and the caller can read it off the coverage.

weight_coverage is a stated exchange rate, in expected responding units, added to a candidate’s marginal gain where its block holds no unit yet. Coverage is a property of the set and cannot be written as a field term, so it enters here: the bonus is monotone and submodular like the rest of H, which is what keeps the 1 - 1/e bound above. 0.0 — the default — is the loop this always was, exactly. It is the fourth term of its kind beside mhcmatch.portfolio.compose()’s weight_evenness and select’s selectivity and weight_escape, and like them it is charged to the objective and never to p.

Both default off, and with codes=None this is the loop it always was.

Parameters:

• k (int)

• cap (int | None)

• weight_coverage (float)

Return type:

list

mhcmatch.cassette.refine(h, J, sel, rounds=4, codes=None, cap=None, must=(), weight_coverage=0.0)

Improve a chosen set by single swaps until no swap raises H, or rounds are spent.

Greedy is one pass and commits its early slots before it has seen what they cost later. A swap pass is the cheapest repair that cannot make things worse: every accepted move strictly raises H, so the sequence terminates, and rejecting ties keeps it deterministic. Same shape as the bounded 2-opt mhcmatch.vector.order() runs after its greedy layout, and for the same reason.

codes / cap / must are the same manufacturing constraints greedy() takes, and they are enforced here too — a swap pass that quietly violated the floor the greedy respected would be worse than not having one.

O(rounds * k * n). On the pools this is used at that is a few hundred thousand operations.

Parameters:

• rounds (int)

• cap (int | None)

• weight_coverage (float)

Return type:

list

mhcmatch.cassette.log_ek(logw, k)

log e_0 .. log e_k, the elementary symmetric polynomials of exp(logw), in log space.

The textbook recurrence carried in logs, O(n k). This is the exact partition function over every size-k subset of the pool without enumerating any of them, which is what makes lam() computable on a 5,000-candidate pool where C(5000, 20) is not a number anybody is going to sum over.

Parameters:

k (int)

Return type:

numpy.ndarray

mhcmatch.cassette.lam(h, sel, k=None)

How good a cassette is relative to what the donor’s own pool could have given.

\[\lambda(S) = H(S) - \log \sum_{|S'| = k} e^{H(S')} + \log \binom{N}{k}\]

The middle term is the log partition function over every size-k subset of the pool; adding log C(N, k) back turns it into a comparison against the average subset rather than against their sum. So lambda = 0 is a cassette exactly as good as a uniformly random one from the same pool, positive is better, and the units are nats.

This is the quantity that compares cassettes across donors and across sizes, which a raw H or a raw sum p does not: dividing by the donor’s own pool 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.539 nats — below a uniform random subset of the same pool — against +3.417 for the greedy argmax of H, a gain of +4.083 nats.

``lambda`` is computed against whatever field it is handed, so a select run with a selectivity weight reports nats above a uniform subset of the tilted objective, not of the untilted one. That is the coherent reading – both sides move together – but it means two runs at different w are not on one axis, exactly as two runs at different gamma are not.

Exact, with the couplings switched off; h is the field. With couplings the partition function has no closed form and the correction must be estimated, which is a separate job — the coupling-aware estimator moves the score a median 0.666 nats against an inter-quartile spread of 1.914, so the exact field-only value is the one to report unless you have measured otherwise.

Parameters:

k (int | None)

Return type:

float

class mhcmatch.cassette.Cassette(index=<factory>, p=<factory>, energy=0.0, lam=0.0, offset=0.0, rho=0.091, gamma=1.0, k=0, pool_n=0, trimmed=0, swaps=0, channels=(), block_live=1.0, coverage=<factory>, selectivity=0.0, weight_escape=0.0, escape=0.0, rule='v1', pi=0.0, how='', not_worse=0.0, diversity=0.0, weight_coverage=0.0, n_composite=0, n_unresolved=0)

Bases: object

A chosen set, and every number needed to justify and reproduce the choice.

index indexes the pool as it was passed in — after trimming, if MAX_POOL fired, so trimmed says how many candidates the objective actually saw.

Parameters:

• index (list)

• p (list)

• energy (float)

• lam (float)

• offset (float)

• rho (float)

• gamma (float)

• k (int)

• pool_n (int)

• trimmed (int)

• swaps (int)

• channels (tuple)

• block_live (float | dict)

• coverage (dict)

• selectivity (float)

• weight_escape (float)

• escape (float)

• rule (str)

• pi (float)

• how (str)

• not_worse (float)

• diversity (float)

• weight_coverage (float)

• n_composite (int)

• n_unresolved (int)

index: list

p: list

energy: float = 0.0

lam: float = 0.0

offset: float = 0.0

rho: float = 0.091

gamma: float = 1.0

k: int = 0

pool_n: int = 0

trimmed: int = 0

swaps: int = 0

channels: tuple = ()

block_live: float | dict = 1.0

q, the stated per-block survival probability the objective priced HLA loss at. 1.0 is “nothing is ever lost”, which is what every cassette built before this existed assumed.

coverage: dict

mhcmatch.portfolio.coverage() of the chosen units against universe — so an allotype carrying zero units is visible, which is the whole inequality a floor exists to catch and which a coverage taken over the cassette’s own labels cannot see.

selectivity: float = 0.0

The stated tumour-over-normal exchange rate charged to the field. 0.0 is off.

weight_escape: float = 0.0

The stated exchange rate charged to the field for escape cost — how much expected response the designer will give up per unit of the tumour’s own price for deleting a candidate. 0.0 is off, and at 0.0 every number here reproduces a cassette built before this existed.

escape: float = 0.0

Mean escape cost of the chosen units, on whatever scale the caller supplied. Reported beside yield_ because the pair is the trade, and one of them alone is not.

rule: str = 'v1'

Which selection rule produced this set — "v1" the mean-variance objective, "v2" the degeneracy rule. Recorded because a number cites the rule that produced it.

pi: float = 0.0

the stated floor on P(this set catches at least as much as the sort).

Type:

v2 only

how: str = ''

how diversity was aggregated over axes — "minmax" or "mean".

Type:

v2 only

not_worse: float = 0.0

the realised P(not worse than the sort), which sits at or just above pi.

Type:

v2 only

diversity: float = 0.0

the diversity actually reached, on the same scale how defines.

Type:

v2 only

weight_coverage: float = 0.0

The stated exchange rate charged to the field for allotype coverage — expected responding units the designer will give up to put a unit on an allotype the set does not yet reach. 0.0 is off, and at 0.0 every number here reproduces a cassette built before this existed. Coverage is otherwise a reported measure (see coverage), and making it an objective term is a stated preference rather than a correction.

n_composite: int = 0

How many restriction cells named several alleles — a genotype rather than an allotype — and were therefore resolved to a presented SET by resolve_restriction(). Reported because it changes what the allotype channel and the loss coupling ran on.

n_unresolved: int = 0

How many restriction cells resolved to no allele at all. Those units keep their slot and their score but carry no allotype: they are excluded from coverage, because a unit whose restriction is unknown covers nothing. Reported rather than dropped silently.

property yield_: float

sum p — the expected number of responding units. A level, not a probability.

mhcmatch.cassette.select(scores, peptides, alleles=None, k=20, tol=0, *, prevalence=None, rho=0.091, gamma=None, rounds=4, max_pool=2000, block_live=1.0, universe=None, max_share=None, selectivity=0.0, offset=None, escape=None, weight_escape=0.0, genes=None, expr_lvl=None, expr_norm=None, features=None, feature_names=(), coexpr=None, presented=None, presented_alleles=None, terms=None, terms_cov=None, graded_allotype=False, dominance=False, collapse_allotype=False, overlap_combine='mean', weight_coverage=0.0, rule='v1', pi=0.5, how='minmax', axes=None, reference=None, sequence='kmer', sequence_mask='face')

Choose k units (within tol) from one donor’s candidate pool, maximising H.

scores are aggregate log-odds — what mhcmatch.rank.aggregate_score() returns — for every candidate in the pool, not a pre-selected shortlist. The pool is what defines the background the choice is made against, and handing this function a shortlist that has already been filtered on binding and expression is the one way to make it report nothing: those are the two terms carrying the largest coefficients in the shipped model, and a pool that has been cut on them has no range left along them.

The steps, in order:

offset overrides step 1: pass an offset fitted on a larger pool to make two subsets of one pool comparable. Without it each subset is calibrated to the same declared prevalence and the difference between them is calibrated away.

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

2. rho is the intra-cassette response correlation. The default is a measured background (RHO_ASSAYED); fit your own by maximum likelihood with mhcmatch.portfolio.betabinom_rho() if you have per-patient counts, which is the one parameter here that any assayed readout can improve.

3. gamma defaults to risk_aversion() at the requested k, so the stated preference is per unit and the objective does not invert at large k. Pass gamma= to use a number verbatim; Cassette.gamma records whichever was used.

4. overlap() builds the mechanistic pair similarity, goal_energy() turns it into (h, J).

5. greedy() takes k + tol units, refine() swaps until no single exchange raises H, and the reported size is the one with the largest H in [k - tol, k + tol], the lower end raised to the coverage floor — the number of universe allotypes the pool can supply, since no smaller cassette can hold them.

tol is the manufacturing tolerance: a budget of “twenty units, give or take three” is k=20, tol=3. With tol=0 the size is exactly k.

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 for a fact the caller can read off pool_n.

Four optional parameters, all off by default and all bit-identical at their defaults.

block_live

q, how often each allotype survives — the HLA-loss rate, scalar or {allele: q}. It reaches the objective as the exact covariance a lost allele implies (goal_energy()), and a unit whose marginal p exceeds its own block’s q raises mhcmatch.portfolio.MarginalExceedsBlock rather than being clipped — clipping there would understate the marginal for exactly the strongest units.

universe

the donor’s distinct allotypes. Two jobs: every one of them that the pool can supply gets a unit before the free slots are filled, and Cassette.coverage is computed against it, so an allotype holding zero units is visible. Without it, coverage is taken over the labels the cassette happens to carry and cannot see an allotype it missed.

max_share

no allotype may hold more than this share of the cassette — 0.4 at k = 20 caps each at eight units. A manufacturing constraint, deliberately not an objective term: the loss coupling already prefers spread, and a second diversity term inside H double-counts unless it is meant.

``rule=”v2”`` selects on the degeneracy instead of on a mean-variance trade. 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. v2 takes the top-k sort as its reference, then swaps slots to raise diversity() while not_worse() against that reference stays at or above pi — so the sort is a floor the rule cannot fall below by more than a stated probability, and what it buys is spread across the four axes build_axes() returns. pi = 1.0 returns the sort exactly. how picks the aggregation over axes ("minmax" or "mean"), and axes overrides the built ones. gamma, rho, block_live and the feature channels below all still apply — they build the J whose covariances not_worse reads.

rule="v1" (the default, and what every recorded result was computed under) is the mean-variance objective and is untouched.

Five further optional parameters carry the feature-based couplings, all off by default and all bit-identical when unset.

features / feature_names

an (n, d) array of per-unit scalars, one coupling channel per column, and the names to record on Cassette.channels. This is how chemistry and expression reach the pair term: C_phys_buried, C_phys_charge, expr_lvl, expr_norm and the selectivity delta are per-unit scalars mhcmatch rank already emits and the objective has never seen. Rows are indexed by the same pool order as scores, and are trimmed with it.

coexpr

a symmetric (n, n) pool-order matrix in [0, 1], one further channel. mhcmatch.expression.coexpression() builds one over the GTEx tissue panel, so two units whose source genes are on in the same tissues are coupled — a mechanism no per-unit scalar can express.

presented

an (n, A) 0/1 matrix, “does allotype a present unit i”, columns in np.unique order over alleles. It changes the loss coupling from one-allotype-per-unit to the exact set form (goal_energy()), so a unit with three routes to the surface is not charged the loss of one. Build it from mhcmatch.store.Store.percent_ranks(). presented_alleles names its columns; pass it whenever the donor’s genotype is wider than the allotypes their candidates are credited to, which it normally is.

graded_allotype

True additionally swaps the equality similarity channel for allotype_overlap() on the same presented matrix, so two units sharing their strongest allele and differing on every other stop being scored as fully redundant. Needs presented; the two are separate switches because they answer different questions — one is how much a pair shares, the other is what a pair loses.

dominance

True adds the score-dominance channel to overlap(). Off by default, where it was on before. It is the one channel built from the score rather than from a mechanism, the pairwise statistic it corresponds to fits attractive on the observational arm where greedy() carries no bound, and it never abstains: measured over every within-donor pair of TESLA’s 736 units and HiTIDE’s 1,558, it is zero on 0.03 % and 0.01 % of pairs against 97.46 % and 96.54 % for the exact-3-mer channel, and it supplies 71.4 % and 78.8 % of the total channel mass — so the allotype channel, the only mechanism among the three, entered H at a third weight. Dropping it buys allotype entropy 0.9629 -> 0.9883 of maximum and Gini 0.1745 -> 0.0889 for 4.466 -> 4.368 expected responding units on a 20-unit cassette, over 19 donors.

Pass ``dominance=True`` to restore the three-channel form, and say so beside the number: every published cassette figure was computed with the channel on, and the arm that drops it was published beside it as the mechanism-only rule. rule="v2" has always run without it, so the two rules now agree on channels.

collapse_allotype

True resolves a composite restriction cell to a single allotype instead of to the presented set (resolve_restriction()). The set form is the default because a cell naming several alleles is a unit with several routes to the surface; this is the pre-promiscuity reading, kept so a result recorded that way stays reproducible.

overlap_combine

"worst" reduces overlap()’s channels by a per-axis-normalised max instead of their mean — a cassette is undone by its worst shared failure mode, not its average one, and averaging is measured to dilute (see dominance above and diversity(), which has taken this argument since it was written). "mean" is the default and is bit-identical to every cassette built before this existed.

weight_coverage

a stated exchange rate, in expected responding units, paid to put a unit on an allotype the set 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 greedy() and refine() as a submodular bonus on the marginal gain; the bound above is preserved. 0.0 — the default — is bit-identical, and at 0.0 coverage remains what it has always been here: reported, never optimised. Charged to the objective and never to p, like selectivity and weight_escape, and like them it should be quoted with its counterfactual at zero.

universe / max_share remain hard constraints and are the right tool when a floor must hold. This is the tool for trading spread against expected response continuously, which a constraint cannot express: at k = 20 over four to six allotypes a floor is already met by the unconstrained argmax, so it binds on nothing and changes nothing.

selectivity

a stated exchange rate, in expected responding units per log2-fold of tumour-over-normal abundance, charged to the field as h_i += selectivity * (expr_lvl_i - expr_norm_i). expr_lvl and expr_norm are the two columns mhcmatch rank already emits.

Charged to the objective, never to ``p``. p is a calibrated marginal that mhcmatch.portfolio.survival() reads literally, so discounting it would silently restate the response model as well as the preference — the rule mhcmatch.portfolio.compose() already follows for weight_cost. It is stated rather than fitted for the same reason gamma is: the shipped EPIC model fits both terms positive — v12 puts expr_lvl at +0.5000 and expr_norm at +0.2222 log-odds per standard deviation, and mhcmatch rank --coefficients prints what an install actually carries — because it was fitted on will this respond and a gene transcribed everywhere responds more often. “High in tumour, low in normal” is a different question — a safety preference the designer declares — and imposing it on the fit would assert an answer the data rejects. Both coefficients stay as measured and both terms stay reported.

escape, weight_escape

a per-unit escape cost in [0, 1] and a stated exchange rate on it, charged to the field as h_i += weight_escape * escape_i. The cost is the tumour’s own price for deleting a candidate: a clonal driver in a gene it cannot silence is expensive to lose, a subclonal passenger is free. Selecting against it is the argument goal_energy() already makes for HLA loss, carried to the other escape route — there the objective prices what happens when an allotype goes, here what happens when the antigen does.

Stated, not fitted, for the reason gamma and selectivity are: durability is a designer’s preference over a horizon no response screen observes, and a screen that reads out at one timepoint cannot supply an exchange rate between catching a response now and keeping it later. Charged to the field, never to p.

Report the pair. weight_escape buys escape with yield_, and a cassette that quotes one without the other has not said what it cost — the CLI’s counterfactual, which re-runs at weight_escape = 0, is the form that reports both.

genes

per-unit source-gene labels, adding the gene channel to overlap(): two units from one gene fall together under one deletion. Independent of escape — that is the field, this is the coupling — and both are needed, because a set can be made of expensive units that all sit in one locus.

Parameters:

• k (int)

• tol (int)

• prevalence (float | None)

• rho (float)

• gamma (float | None)

• rounds (int)

• max_pool (int)

• max_share (float | None)

• selectivity (float)

• offset (float | None)

• weight_escape (float)

• graded_allotype (bool)

• dominance (bool)

• collapse_allotype (bool)

• overlap_combine (str)

• weight_coverage (float)

• rule (str)

• pi (float)

• how (str)

• sequence (str)

• sequence_mask (str)

Return type:

Cassette

mhcmatch.cassette.size_for(scores, peptides, alleles=None, *, target=1, confidence=0.9, k_max=40, prevalence=None, rho=0.091, gamma=None, max_pool=2000, block_live=1.0, offset=None, dominance=False, escape=None, weight_escape=0.0, genes=None)

The smallest cassette that reaches P(>= target responses) >= confidence for this donor.

A fixed k asks every donor the same question and gets a different answer. A donor whose best candidates are strong reaches a given confidence in five units; a donor whose pool is shallow, or whose top of the list is not actually that good, does not reach it in twenty — and the honest response to that is a larger cassette, not the same one reported with the same number on it.

The probe walks the objective’s own greedy order and evaluates mhcmatch.portfolio.p_at_least() on each prefix under the block model, the block being the allotype where one is given. It returns the size, not the cassette: hand k back to select(), which re-runs greedy and refine() at that size under its own risk_aversion(). The probe’s own gamma is taken at k_max for the single pass.

k_max is a manufacturing ceiling, not a search bound: when the confidence is unreachable inside it the ceiling is returned with reached = False and p_at_least says how far it got. That is a real answer about the donor and it must not be silently rounded into a smaller cassette that claims the target.

block_live is the HLA-loss rate, and it is what this function was always missing: the block model it evaluates has priced a lost allotype since it was written, and this call site pinned q at 1.0. Below 1 a donor needs more units to reach the same confidence, because some of the ones they have can be lost together — which is the question a designer asking to be protected from losing HLA is actually asking.

escape, weight_escape and genes are passed through so the probe walks the greedy order the caller’s own select() will walk. Without them a weighted selection and its --confidence size are answers to two different questions, and the size comes back for a cassette nobody is going to build.

offset is the same argument select() takes and for the same reason, and it has to be passed whenever select is given one. Sizing a subset of a pool while calibrating it on itself pins its mean to the declared prevalence, so p_at_least is computed on probabilities that no longer say the subset is weaker than its parent, and every group comes back the same size. That is the defect offset exists to prevent, and it reappears here whenever the two are not given the same value.

Returns k · reached · p_at_least at that k · target · confidence · curve, the confidence at every size from 1 to the one returned.

Parameters:

• target (int)

• confidence (float)

• k_max (int)

• prevalence (float | None)

• rho (float)

• gamma (float | None)

• max_pool (int)

• offset (float | None)

• dominance (bool)

• weight_escape (float)

Return type:

dict

mhcmatch.cassette.score(scores, peptides, alleles=None, chosen=None, *, pool_scores=None, pool_peptides=None, offset=None, prevalence=None, rho=0.091, gamma=None, block=None, block_live=1.0, target=1, universe=None)

Score a cassette that already exists, on axes that survive changing donor and changing k.

Two ways to call it. Pass the cassette alone (scores, peptides) with an offset fitted over every cassette being compared, and you get the level — yield is the expected number of responding units, and two cassettes’ levels are comparable because they were calibrated together. Pass the donor’s pool as well (pool_scores, pool_peptides) and you also get lam, which compares the cassette against what that donor’s own pool could have given, and is therefore comparable across donors and across sizes without any shared calibration at all.

Returned keys:

yield sum p, expected responding units · p_mean · p_at_least P(X >= target) under the block model · n_effective how many independent shots the cassette is worth · lam nats above a uniform subset of the donor’s pool, None without a pool · rho_hla / rho_seq / rho_dom the three pairwise statistics · coverage allotype counts, Gini and entropy share, when alleles is given · yield_loh / lost_allotype the expected responding units left after the worst single allotype is lost, and which one that is.

``yield_loh`` is the worst case, not an average, and that is the point. A designer asking to be protected from losing HLA is asking about the bad draw: LOH takes a specific allele, and a cassette whose expected count survives it is a different object from one whose average over losses looks acceptable. It is a level in the same units as yield, so yield_loh / yield reads directly as the share of expected response that does not depend on any one allotype.

Pass ``universe`` — the donor’s distinct allotypes — or coverage is computed over the labels the cassette happens to carry and an allotype holding zero units is invisible, which is exactly the inequality the index exists to report. A patient homozygous at B has five distinct class-I allotypes, not six, so passing it is also what stops a genotype being scored as a design flaw.

``H`` is deliberately not reported here. goal_energy() renormalises the overlap so the set it is handed averages to rho, and the dominance channel of overlap() is scaled by the range of the set it is given — so an H computed on a cassette alone is not the same 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 with overlap() and goal_energy() and evaluate both index sets with energy(); that is five lines and it is exact. lam needs none of that — it is a field-only quantity with a closed form, and it is the axis that already crosses donors and sizes.

block is what a unit’s failures are shared with; the default is the allotype, which is the rule mhcmatch.vector ships. Passing block_live below 1 asserts each block is only live that often, and a unit whose marginal p exceeds its block’s q raises mhcmatch.portfolio.MarginalExceedsBlock rather than being clipped — clipping there would understate the marginal for exactly the strongest units.

Parameters:

• offset (float | None)

• prevalence (float | None)

• rho (float)

• gamma (float | None)

• target (int)

Return type:

dict

mhcmatch.vector module

Assembling the cassette, once selection has chosen the candidates: what to withdraw on safety grounds, how many units each allotype should carry, in what order, and joined by what.

screen() runs first, and it excludes rather than down-ranks: a register that is itself an essential-tissue self peptide, or a target gene transcribed where it was assumed silent, is withdrawn, because the second-best cassette is cheap and myocarditis is not. The two fatal precedents behind that rule, and the measurement that chose self_origin_risk() over a mimicry-similarity screen, are in the function’s own documentation. Competition for a response is then local to the antigen-presenting cell and strongest within an allotype, so expected yield is a sum of independently saturating per-allotype terms rather than one global budget – select() grows each allotype while the next candidate beats that allotype’s own expected yield per slot, and diversification follows from the arithmetic instead of a quota. order() scores every register spanning every junction against the recipient’s own allotypes and picks the spacer and ordering that minimise predicted junctional binding, trying no spacer first.

The two ends of that pipeline are joins to the rest of the library. units_from_context() closes the front: rank_fasta() emits minimal epitopes while a unit is the long window around the mutation, and where that mutation sits is in the FASTA header rather than in the ranking, so the two are combined rather than either being guessed at – and rows are grouped by variant, since twenty registers of one mutation are twenty rows in a ranking and one thing to put in a cassette. back_translate() closes the back, turning sequence into a coding sequence. It is not a codon optimiser: it fixes the two failure modes specific to a concatemer – the m1-pseudouridine +1-frameshift motif (slippery_sites(), whose seams the designer chooses) and synthesis-hostile homopolymers (which spacers like AAA manufacture directly) – and leaves GC content, secondary structure and CpG to a manufacturer’s own tooling. translate() exists so “synonymous” stays checkable rather than asserted.

Assembling a polyepitope vaccine cassette: what to refuse, how many units, in what order, joined by what.

Candidate selection — which mutations are worth targeting — is mhcmatch.rank and the immunogenicity stack behind it. This module is the step after: given ranked candidates with calibrated probabilities, decide what to withdraw on safety grounds, how many to carry, how to lay them out, and what to put between them. Those are four separate questions with four different literatures, and none of them is answered by the candidate score.

Selection and assembly are kept apart because the assembly answer depends on the set, not the candidate: whether to carry a 12th epitope depends on what the first eleven already cover, and the cost of a junction depends on which two units sit either side of it.

Two selection rules ship, and they answer different questions. mhcmatch.cassette.select() maximises the mean–variance energy H at a size you fix: you say twenty units, it returns the best twenty and a score comparable across donors. select() here fixes no size at all — it grows each allotype while the next candidate beats that allotype’s own expected yield per slot, so k is an output and n0, per-allotype capacity, is the input. Use the first when the platform has already decided how many units the construct carries, which is the usual case; use this one when the question is how many are worth carrying at all. The two are not competing estimates of one quantity and neither approximates the other.

What to refuse — an exclusion, and it runs first. screen() withdraws a unit whose own target gene is transcribed in a tissue that must not be attacked, or one of whose registers coincides with a self peptide from an unrelated essential-tissue gene. It excludes rather than down-ranks: the second-best cassette is cheap and myocarditis is not, and capacity spent on a unit that has to be withdrawn is capacity not spent on a safe one. Two patients died of cardiogenic shock and two of necrotising leukoencephalopathy in the trials that define this problem; the references, the shape of each event, and the measurement that chose this screen over the more obvious mimicry-similarity one are in screen() and self_origin_risk().

How many — a per-allotype stopping rule, not a constant.

The clinical numbers (20 for autogene cevumeran across two RNAs, 34 for mRNA-4157, 20 in four pools for NeoVax) derive from no published objective function, and no trial has delivered N versus 2N epitopes at matched dose to measure the cost of the extra ones. What is established is the shape of the competition, and it is not a single global budget:

• Competition is for the antigen-presenting cell, not for MHC. Suppression required co-presentation on the same DC and was reversed by injecting excess pulsed DCs (Kedl et al., J Exp Med 2000, PMID 11034600); the mechanism is CD27 cleavage capturing CD70 on the DC (Burchill et al., Eur J Immunol 2015, PMID 26179759).

• It is strongest within an allotype and can be net-positive across allotypes: with large pre-existing CTL, a multi-epitope vector “failed to prime efficiently new CTL responses that were restricted by the same MHC gene … and vaccine-induced CTL responses restricted by other MHC genes were enhanced” (Sherritt et al., Eur J Immunol 2000, PMID 10671225).

• The total response is not a fixed pool. Deleting four dominant LCMV epitopes gave only “minor response increases … and no new epitopes being recognized” (Kotturi et al., J Immunol 2008, PMID 18641351).

So the expected yield is a sum of independently saturating per-allotype terms. With p_{a,1} >= p_{a,2} >= ... the calibrated probabilities on allotype a and within-allotype competition saturating as n0 / (n0 + n_a):

E = sum_a  n0 * S_a(n_a) / (n0 + n_a),        S_a(n) = sum_{i<=n} p_{a,i}

add the next unit on allotype a   <=>   p_{a,n+1} > S_a(n_a) / (n0 + n_a)

select() implements exactly that line: keep adding to an allotype while the next candidate’s probability beats that allotype’s current expected yield per slot. Because each allotype saturates on its own, the marginal value of a crowded allotype’s next unit falls below an empty allotype’s first one, so diversification across allotypes falls out of the arithmetic instead of being imposed as a quota. n0 is the one free parameter and it means per-allotype capacity; it is not fitted here, because nothing in the public record fits it. Pass the value you can defend and record it — Selection.n0 carries it into the result so a cassette can always name its own assumption.

What between them — scan, do not assume.

The only linker with causal evidence is GPGPG: unspaced concatenation of four HLA-DR epitopes created a high-affinity junctional epitope that suppressed the response to all four, and inserting GPGPG restored all four (Livingston et al., J Immunol 2002, PMID 12023344). Against that, all six orderings of three Fel d I regions produced no detectable junctional responses at all (Rogers et al., Mol Immunol 1994, PMID 7521933). Both results are real, which is why this module measures each junction against the recipient’s own allotypes rather than picking a linker by reputation.

Two mechanistic constraints worth knowing before choosing a spacer:

• AAY ends in tyrosine. ERAP1 prefers hydrophobic C-termini and has low affinity for charged ones (Chang et al., PNAS 2005, PMID 16286653), so a terminal Y genuinely aids processing – while also supplying the C-terminal anchor for A*01:01, A*29:02 and B*35:01. It is a trade-off, not a mistake, and which way it falls is donor-specific.

• KK leaves charged residues at the boundary: the mirror image, poor for ERAP1.

• Gly/Pro-rich spacers use residues disfavoured at MHC-I anchor positions and abundant in the C-terminal regions from which ligands are cleaved (Martin-Galiano & Lopez, PLoS One 2019, PMID 30645615), so they sit in the permissive zone.

SPACERS therefore leads with None. A cassette whose junctions are already clean should carry no spacer at all — every inserted residue is sequence that has to be translated and could itself form a binder.

Order is chosen the way pVACvector chooses it (Hundal et al., Cancer Immunol Res 2020, PMID 31907209): score every register spanning every junction, treat the units as nodes of a complete graph whose edge cost is the strongest predicted binder at that junction, and find a cheap open path. This module uses a deterministic greedy path plus bounded 2-opt rather than simulated annealing — no RNG, so a cassette is reproducible from its inputs.

Scoring is injected, never imported. order() and scan_junctions() take a binder callable, so the layout logic is testable without a Store, a panel or any download, and a caller who wants a different presentation model just passes it:

>>> from mhcmatch import vector
>>> units = [vector.Unit("AAAAAAAAAKAAAAAAAAAAAAAAAAA", 9, "GENE1", "HLA-A*02:01", 0.30),
...          vector.Unit("CCCCCCCCCRCCCCCCCCCCCCCCCCC", 9, "GENE2", "HLA-A*02:01", 0.12)]
>>> sel = vector.select(units, n0=8)
>>> [u.gene for u in sel.units]
['GENE1', 'GENE2']
>>> cas = vector.order(sel.units, binder=lambda peps, alleles: [0.0] * len(peps))
>>> cas.spacer is None and len(cas.units) == 2
True
>>> cas.sequence == units[0].peptide + units[1].peptide
True

A binder returns one number per peptide, higher meaning a stronger predicted binder; -log10(%rank) is the natural choice and is what store_binder() builds.

From the command line, where the whole pipeline is one call:

mhcmatch vector --candidates units.tsv --n0 8 --screen --fasta cassette.fasta

units.tsv carries peptide, gene, allele, p and optionally mutation_index — and peptide is the long window, not rank’s minimal epitope. --screen is opt-in because it builds a whole-proteome index; without it no safety check runs at all and the cassette carries whatever it was handed. mhcmatch deslip <cds> --fix out.fasta is the slippery_sites() half, which takes nucleotides rather than peptides and so is its own command.

mhcmatch.vector.SPACERS: tuple = (None, 'GPGPG', 'GGS', 'AAA', 'AAY', 'GPGPGPG', 'HHAA', 'AAL')

Spacers tried in order, None first. A clean junction needs no spacer, and pVACvector’s default list is tried only when one is needed. See the module docstring for why AAY is a trade-off rather than a default.

class mhcmatch.vector.Linker(name, sequence, family, cls, note)

Bases: object

One named linker preset: its residues, what family it belongs to, and where it comes from.

sequence is the empty string for "none" — direct head-to-tail concatenation — which is a real design choice and not a missing value, so it is a member of the table rather than an absence from it.

Parameters:

• name (str)

• sequence (str)

• family (str)

• cls (str)

• note (str)

name: str

sequence: str

family: str

cls: str

note: str

mhcmatch.vector.LINKERS: dict = {'AAA': Linker(name='AAA', sequence='AAA', family='class-I favouring', cls='mhc1', note="alanine-based, and back-translated it is the construct's main source of synthesis-hostile homopolymer (see MAX_HOMOPOLYMER)"), 'AAY': Linker(name='AAY', sequence='AAY', family='class-I favouring', cls='mhc1', note='the most frequently used proteasome-sensitive spacer; the terminal tyrosine both suits ERAP1 and supplies a C-terminal anchor for three common allotypes'), 'EAAAK': Linker(name='EAAAK', sequence='EAAAK', family='rigid', cls='any', note='an alpha-helical rigid linker'), 'G': Linker(name='G', sequence='G', family='minimal', cls='any', note='a single glycine'), 'G4S': Linker(name='G4S', sequence='GGGGS', family='GS-rich flexible', cls='any', note='the canonical (G4S) flexible linker'), 'G4S2': Linker(name='G4S2', sequence='GGGGSGGGGS', family='GS-rich flexible', cls='any', note='(G4S)2'), 'GPGPG': Linker(name='GPGPG', sequence='GPGPG', family='class-II oriented', cls='mhc2', note='restored all four responses that unspaced concatenation of HLA-DR epitopes had suppressed (PMID 12023344)'), 'GPGPGPG': Linker(name='GPGPGPG', sequence='GPGPGPG', family='class-II oriented', cls='mhc2', note='the seven-residue variant'), 'GS10': Linker(name='GS10', sequence='GGSGGGGSGG', family='GS-rich flexible', cls='any', note="the 10-mer GS linker reported for closely related RNA constructs in the patent literature; the pentatope format's own linker is described, not published"), 'GS10b': Linker(name='GS10b', sequence='GGSGGGSGGS', family='GS-rich flexible', cls='any', note='a compositional variant of the GS 10-mer'), 'GS10c': Linker(name='GS10c', sequence='GGGSSGGGSG', family='GS-rich flexible', cls='any', note='a compositional variant of the GS 10-mer'), 'KK': Linker(name='KK', sequence='KK', family='minimal', cls='mhc1', note='lysine-rich and cleavage-prone; leaves charged residues at the seam, which is the mirror image of AAY and poor for ERAP1'), 'RRRR': Linker(name='RRRR', sequence='RRRR', family='protease-cleavable', cls='any', note='a polybasic furin-type site'), 'furin': Linker(name='furin', sequence='RAKR', family='protease-cleavable', cls='any', note='a furin recognition site'), 'none': Linker(name='none', sequence='', family='minimal', cls='any', note='direct head-to-tail concatenation; every inserted residue is sequence that has to be translated and could itself form a binder')}

Named linker presets, grouped by design intent. Pass a key of this table anywhere a spacer is taken — order(), assemble(), mrna(), --linker on the command line — and it resolves through resolve_linker().

The families, and what each is trying to buy:

GS-rich flexible

Flexible separation with no secondary structure of its own and low intrinsic immunogenicity. This is the family the manufactured pentatope format belongs to: the methodological literature describes those constructs only as joined by “non-immunogenic 10-mer glycine/serine linkers”, and GS10 is the explicit 10-residue sequence reported for closely related RNA constructs in the patent literature. It is a reconstruction, not a sequence read off a published construct, and the compositional variants are here because the format is specified by description rather than by residue.

class-I favouring

Meant to place a preferred cleavage site at the seam so the flanking units are liberated with the exact C-terminus class I needs. AAY ends in tyrosine, which suits ERAP1’s preference for hydrophobic C-termini (Chang et al., PNAS 2005, PMID 16286653) and also supplies the C-terminal anchor for A\*01:01, A\*29:02 and B\*35:01 — a trade-off whose direction is donor-specific, which is why it is scored and not defaulted to.

class-II oriented

GPGPG is the one linker with causal evidence behind it, and the evidence is class II: unspaced concatenation of four HLA-DR epitopes created a junctional epitope that suppressed the response to all four, and inserting GPGPG restored all four (Livingston et al., J Immunol 2002, PMID 12023344).

minimal

The shortest construct and the fewest junctional residues able to form a novel binder. none is a legitimate answer and SPACERS leads with it for that reason.

protease-cleavable

A recognition site for a specific protease, so cleavage is placed rather than predicted.

rigid

An α-helical linker that enforces spatial separation. More usual for multi-domain antigens than for epitope strings; included so a comparison is not confined to flexible linkers.

What this table does not do is rank itself. cls records the class a linker is intended for, which is provenance and not a measurement, and the two mechanisms that would decide the ranking act at different positions: glycine and proline are abundant in the C-terminal regions from which class-I ligands are cleaved (Martin-Galiano & Lopez, PLoS One 2019, PMID 30645615), which is a processing argument for them, while the same residues immediately flanking a class-I epitope inhibit recognition of the epitope on their amino-terminal side (Bergmann et al., J Immunol 1996, PMID 8871618). Choosing between them on either citation alone would be picking a linker by reputation. order() measures each candidate against the recipient’s own allotypes instead, and that measurement — not this table — is what selects.

mhcmatch.vector.resolve_linker(linker)

A preset name, an explicit residue string or None, resolved to residues.

None and "none" both mean no linker and both come back as None, so a caller can write --linker none and get what SPACERS leads with. Anything else is looked up in LINKERS first and taken as a literal sequence only if it is not a key, which is why an unknown name that happens to be a valid peptide — GGS is both — cannot be silently mistaken for a preset it is not. A string that is neither a key nor 20-letter amino acid raises, rather than reaching the construct as an untranslatable residue.

>>> resolve_linker("GS10")
'GGSGGGGSGG'
>>> resolve_linker("none") is None and resolve_linker(None) is None
True
>>> resolve_linker("GGS")
'GGS'

Return type:

str | None

mhcmatch.vector.JUNCTION_LENGTHS: tuple = (8, 9, 10, 11)

MHC-I register lengths scanned across a junction. Class II cores are 9-mers read out of a longer span, so a class-II junction scan wants MHC2_JUNCTION_LENGTHS instead.

mhcmatch.vector.MHC2_JUNCTION_LENGTHS: tuple = (12, 13, 14, 15)

Class-II junction registers. The bound core is 9 residues but the presented span is longer, so the scan needs the window a core could be read from.

mhcmatch.vector.flank_identity(a, ai, b, bi, length, k=10)

Fraction of matching residues in the k positions either side of a shared register.

a/b are the two contexts the register was found in (a unit’s 27-mer and a reference protein), ai/bi the register’s 0-based start in each, length its length. Positions truncated by either sequence’s end are not compared, and a register flush against both ends scores 0.0 — nothing was compared, so nothing supports homology.

This is what separates a paralog from a coincidence, and the separation is the whole reason the report tier is usable. Two proteins from one locus share their flanks as well as the register, so a d = 1 match between them is descent, not mimicry, and the T cell that sees one is a T cell tolerance already had to deal with. A match whose flanks are unrelated is an independent occurrence: the same nine residues arrived twice by chance, and that is the object the titin and MAGE-A12 deaths are drawn from. Measured on 178 validated immunogenic neoantigens, the cut at 0.5 removes 156 of 230 different-gene d = 1 hits at L=9 — 130 of them at 90% or better, which is one locus under two symbols rather than a coincidence at all.

The comparison is bounded by whatever context the caller gives, and a 27-mer unit gives ±9-10 residues. That is enough to separate loci, not superfamilies: NRAS → KRAS scores 0.23 and is reported, which is the wanted behaviour — a T cell raised on an NRAS Q61 neoantigen that cross-reacts to wild-type KRAS is a real on-target/off-tumour concern, and KRAS is transcribed everywhere.

Parameters:

• a (str)

• ai (int)

• b (str)

• bi (int)

• length (int)

• k (int)

Return type:

float

mhcmatch.vector.ESSENTIAL_TISSUES: tuple = ('Heart', 'Brain', 'Nerve', 'Lung', 'Liver', 'Kidney', 'Adrenal Gland', 'Pituitary', 'Muscle - Skeletal', 'heart muscle', 'kidney', 'liver', 'lung', 'adrenal gland', 'pituitary gland', 'skeletal muscle', 'smooth muscle', 'cerebellum', 'cerebral cortex', 'midbrain', 'basal ganglia', 'hippocampal formation', 'hypothalamus', 'amygdala', 'spinal cord', 'choroid plexus', 'retina')

GTEx SMTSD tissue prefixes whose destruction is not survivable or not repairable, so a candidate whose self-mimic is transcribed there is excluded rather than ranked down. Matched by prefix because GTEx splits organs into regions – Brain alone covers twelve of the 53 tissue names, Heart two.

The list is a clinical judgement, not a fitted parameter, and both fatal precedents behind it are in screen()’s docstring. Override it: a cassette for a patient who has already lost an organ is a different calculation, and so is one for a tissue the protocol accepts damaging. The nine organs whose damage is not survivable or not recoverable. Both naming schemes. The expression table carries 123 distinct context names from two sources – GTEx-style ("Brain - Caudate (basal ganglia)", "Muscle - Skeletal") and HPA-style lowercase ("basal ganglia", "skeletal muscle") – and a nine-entry Title-Case tuple matched by str.startswith sees only 22 of the 123. Thirteen essential organs were invisible: heart muscle, kidney, liver, lung, adrenal gland, pituitary gland, cerebellum, cerebral cortex, midbrain, hippocampal formation, spinal cord, skeletal muscle, smooth muscle.

Measured cost of that omission, on the shipped table: of the 7,527 genes reaching 50 TPM in an essential tissue, 1,517 (20.2 %) were invisible to the screen – among them CEACAM5 (4.65 seen vs 28.50 actual; Parkhurst 2011 colitis, 3 of 3 patients), CDH13 (7.09 vs 64.20) and albumin (26,217 vs 198,524). This is a false negative in the fatal direction, so the list is not a scope choice: it is the same nine organs, spelled the way the data spells them.

class mhcmatch.vector.Unit(peptide, mutation_index, gene, allele, p, cls='mhc1', kind='missense')

Bases: object

One vaccine unit: the long peptide carrying one mutation, plus what it is worth.

peptide is placed into the cassette verbatim, so build it with unit() rather than slicing by hand — the mutation has to sit far enough from both ends that every register containing it is generated (see unit()).

allele is the restriction the unit is credited to for the per-allotype budget in select(). A long peptide usually presents on several allotypes; credit it to the one whose probability p refers to, and carry a second Unit for the same mutation only if a different allotype genuinely needs a different window.

p is a calibrated probability of eliciting a detectable response, on the operating prior the caller intends to deploy at — not a corpus-prevalence posterior, which overstates it. Nothing here checks that; it is the caller’s contract.

Parameters:

• peptide (str)

• mutation_index (int)

• gene (str)

• allele (str)

• p (float)

• cls (str)

• kind (str)

peptide: str

mutation_index: int

gene: str

allele: str

p: float

cls: str = 'mhc1'

kind: str = 'missense'

What kind of variant produced the neoepitope. "missense" is a single substitution against a self protein; anything else – frameshift, fusion, splice, retained_intron, ORF, editing – is non-conventional and is charged to its own arm by mhcmatch.portfolio.compose(). The distinction earns a quota of its own because a non-conventional product is foreign over a stretch rather than at one position, so it fails differently from a missense: whatever makes the missense arm miss (a wrong wild type, a tolerised residue) does not make this arm miss.

class mhcmatch.vector.Selection(units=<factory>, dropped=<factory>, n0=0.0, trace=<factory>, keys=<factory>)

Bases: object

What select() kept, what it dropped, and the rule that decided.

trace is one row per considered candidate in the order the rule saw it, carrying the threshold it was compared against. A cassette that cannot explain why its 12th unit is in and the 13th is out is not auditable, and the threshold is cheap to record.

Parameters:

• units (list)

• dropped (list)

• n0 (float)

• trace (list)

• keys (list)

units: list

dropped: list

n0: float = 0.0

trace: list

keys: list

property expected_yield: float

sum_b n0 * S_b / (n0 + n_b) — expected responses under the saturation model.

Grouped by whatever select() blocked on, not unconditionally by allotype: the yield has to be computed against the same partition the stopping rule spent its budget on, or it describes a cassette that was never built.

per_allele()

{allele: (n_units, summed p, saturated yield)} — where the budget actually went.

Return type:

dict

per_block()

{block key: (n_units, summed p, saturated yield)} for the rule’s own partition.

Return type:

dict

class mhcmatch.vector.Cassette(units, spacer, sequence, boundaries=<factory>, junctions=<factory>, cost=0.0)

Bases: object

An ordered, spaced cassette and the junction evidence behind its layout.

sequence is the epitope cassette only — no start codon, no stop, no leader and no trafficking domain. Those belong to the vector backbone, and a cassette that silently included them could not be cloned into one that already had them.

Parameters:

• units (list)

• spacer (str | None)

• sequence (str)

• boundaries (list)

• junctions (list)

• cost (float)

units: list

spacer: str | None

sequence: str

boundaries: list

junctions: list

cost: float = 0.0

property worst_junction: float

Strongest predicted junctional binder anywhere in the cassette.

mhcmatch.vector.unit(context, mutation_offset, length=27, **kw)

Centre mutation_offset of context in a length-mer window, and wrap it as a Unit.

length defaults to 27 with the mutation at position 14 — the configuration BioNTech’s backbone carries (Kreiter et al., Nature 2015, PMID 25901682) and the one that guarantees every 8-to-14-residue register containing the mutation is present, so all lengths and all allotypes are covered by a single unit and duplication buys nothing.

Centring is also why a minimal epitope should never be a unit. A minimal peptide can load directly onto MHC-I of any cell, including non-professional APCs with no costimulation, and does: injected alone it “transiently activated CD8+ effector T cells, which eventually failed to undergo secondary expansion or to kill target cells”, while simply extending it to 30 residues restored both, “independent of T cell help, because the longer CTL peptide was predominantly presented in the locally inflamed draining lymph node” (Bijker et al., J Immunol 2007, PMID 17911588). Short units are not merely less efficient, they are the tolerising configuration.

Near a protein terminus the window is clamped and the mutation sits off-centre; that is unavoidable and Unit.mutation_index records where it actually landed.

Parameters:

• context (str)

• mutation_offset (int)

• length (int)

Return type:

Unit

mhcmatch.vector.units_from_context(rows, records, length=27, cls='mhc1')

[Unit] from ranked minimal epitopes plus the window FASTA they were called on.

This is the join between mhcmatch.rank.rank_fasta() and this module. rank emits minimal epitopes and a score; a unit is the long window around the mutation, and where that mutation sits is in the FASTA header rather than in rank’s output – so neither side alone can build one. rows are dicts carrying peptide (the minimal epitope), gene, allele and p; records is mhcmatch.predict.parse_fasta()’s output for the very FASTA rank was pointed at.

One unit per variant, not per epitope. Twenty registers of one mutation are twenty rows in rank and one thing to put in a cassette, and select() spends capacity per unit – so rows are grouped by their source window and the group’s best-scoring row supplies the allotype and the score. Anything whose epitope matches no window is returned to the caller’s attention by being absent; counts are the caller’s to report.

All four header families are admitted, not only Somatic:. Each carries its novel residue or span somewhere different, and _centred_context() reads all four; skipping the other three discarded 317 of the 489 non-missense records in one real cohort and left the nonconventional quota arm with nothing to fill itself from – which is the one arm whose whole point is that it fails differently from the missense arm.

Parameters:

• length (int)

• cls (str)

Return type:

list

mhcmatch.vector.junction_windows(left, right, spacer=None, lengths=(8, 9, 10, 11))

Every lengths-mer that spans the left/right boundary, as (peptide, offset).

A window spans the boundary when it contains at least one residue from each side, so windows lying wholly inside either unit are excluded — those are the intended epitopes, not artefacts of concatenation. With a spacer, a window containing only spacer plus one side still counts: it did not exist in either unit and did not exist in any genome.

offset is the window’s start relative to the start of left, so a caller can point at the exact residues to change.

Parameters:

• left (str)

• right (str)

• spacer (str | None)

Return type:

list

mhcmatch.vector.scan_junctions(units, binder, spacer=None, lengths=(8, 9, 10, 11), alleles=None, binder_threshold=None)

Score every junction of units laid out in the given order.

Returns one dict per junction: left, right (unit indices), score (the strongest predicted binder spanning it), peptide and offset of that worst window, and n_windows.

The strongest rather than the mean, because a junction is a hazard if it forms one good binder; averaging it against the many bad windows around it hides exactly the case being looked for. This is pVACvector’s “lowest binding score” convention read on a scale where higher is a better binder.

Parameters:

• spacer (str | None)

• binder_threshold (float | None)

Return type:

list

mhcmatch.vector.screen(units, risk, lengths=(8, 9, 10, 11), notes=None)

(kept, rejected) — drop units carrying a register that mimics an essential-tissue self peptide. Run this before :func:`select`: capacity spent on a unit that has to be withdrawn is capacity not spent on a safe one.

The precedent is not hypothetical and it is not a binding-prediction failure. An affinity-enhanced TCR against the HLA-A*01:01-restricted MAGE-A3 epitope EVDPIGHLY killed the first two patients infused, by cardiogenic shock within days; autopsy found T-cell infiltrate and myocardial damage with no MAGE-A3 expressed in heart at all, and the off-target was ESDPIVAQY from titin (Linette et al., Blood 2013;122(6):863-71, PMID 23770775; Cameron et al., Sci Transl Med 2013;5(197):197ra103, PMID 23926201). Separately, a TCR recognising MAGE-A3/A9/A12 caused necrotising leukoencephalopathy and two deaths, because MAGE-A12 turned out to be transcribed in human brain (Morgan et al., J Immunother 2013;36(2):133-51, PMID 23377668).

Read the two together and they give this function its shape:

• The hazard is in the source gene’s tissue, not the candidate’s score. Both events were invisible to binding prediction and visible in expression: MAGE-A12 is transcribed in brain, and titin in heart. So the check joins a peptide to a protein to a tissue.

• Two different questions, and the unit answers both. Is the unit’s own target gene transcribed somewhere it must not be attacked (MAGE-A12)? And does any register of it coincide with a self peptide from an unrelated essential-tissue gene?

• Exclusion, not down-ranking. The second-best cassette is cheap; myocarditis is not.

The unit’s own gene has to be excluded from the register test, or the screen rejects everything. A 27-mer is native context by design: its flanking registers are self peptides from its own parent protein, and the mutated register sits one substitution from that protein’s wild type. Screened naively at any useful radius, every unit of every cassette fires. Those matches are also the ones tolerance already covers — the flanks are presented in normal tissue daily. What is not covered is a register that coincides with a different protein, which is why risk is handed the unit and not just its registers.

risk(unit, registers) -> [reason, ...] returns zero or more reasons, empty meaning safe. Each reason is a dict; a "register" key naming which register triggered it is carried into the record, and its absence means the reason is unit-level. self_origin_risk() builds one from the human proteome and reference expression. Injected for the same reason binder is — the policy here is testable with no panel, no proteome and no download, and a site with its own toxicity list substitutes it wholesale.

rejected is [(unit, register, reason)], register being None for a unit-level reason. A withdrawn candidate has to say what withdrew it: “the screen dropped 3 of 40” is not a safety argument, and the reason is what a clinician overrides or accepts.

A reason may decline to withdraw. r["veto"] = False marks a finding that is recorded but does not exclude – the graded mode of self_origin_risk(), where a hit below veto_tpm is a cost to composition rather than a refusal. The key’s absence means True, so a risk callable that never sets it behaves exactly as it always did. Pass notes=[] to collect those non-vetoing findings; they arrive in the same (unit, register, reason) shape as rejected and are what offtarget_cost() reads.

One batch query for the whole candidate list, not one per unit. Where risk exposes a prepare(registers) (self_origin_risk() does), it is handed the deduplicated registers of every unit before any unit is judged. A 27-mer carries ~70 registers and a real cohort carries thousands of units, so the per-unit call pattern made ~19,000 proteome queries where 1 suffices, and units share registers heavily.

A junction can manufacture a self-mimic just as it can manufacture a binder, but that check needs the layout, so it belongs after order() rather than here.

Return type:

tuple

mhcmatch.vector.offtarget_cost(findings)

{unit: cost} from screen()’s notes (or its rejected) – the size of a unit’s off-target fingerprint, being the number of distinct (clause, gene) pairs it reaches.

Distinct genes, not reasons: a gene is transcribed in many tissues and self_origin_risk() reports one finding per tissue, so counting reasons would charge a unit for the breadth of GTEx rather than the breadth of its off-targets. Units absent from findings are absent from the dict; read it with .get(u, 0.0).

This is the number mhcmatch.portfolio.compose() subtracts under weight_cost. It is a count and not a probability on purpose – there is no calibration behind “how much worse is two off-target genes than one”, so the weight is the caller’s to set and to record.

Return type:

dict

mhcmatch.vector.presented(findings, binder, threshold=-1.47712, alleles=None)

Keep the near-identity findings whose off-target variant is actually presented.

A d = 1 coincidence is only a hazard if a T cell can see it, and seeing it means the off-target’s own sequence — not the unit’s register — is presented on the allotype the unit was selected for. A variant that no allotype presents is a sequence coincidence and nothing more, so it is dropped from the fingerprint rather than reported as a safety consideration. Findings with no "variant" key (clauses 1 and 2, and any sub-veto finding under graded) pass through untouched: presentation is not what they are about.

binder(peptides, alleles) -> [score] is store_binder()’s contract, score being -log10(%rank) so higher is stronger. Every finding on one allotype goes in one call — the alternative is a restriction() per finding, and a cassette’s report tier carries thousands. alleles overrides the per-unit allotype for callers whose units carry none.

``threshold`` defaults to ``-log10(30)``, and the conventional 2% rank would be wrong here. This gate is a safety read-out, so the expensive error is missing a hazard, and the cut belongs where the positives are rather than at a number borrowed from a different scorer. Measured on the 176 assayed immunogenic neoantigen/allotype pairs in isalgo/pmhc_data, scored by this scorer on their own allotype, the median sits at 0.69% rank and the 5th percentile at 14.3%:

%rank cut	assayed immunogenic peptides kept	units carrying a report
none	100.0%	27 of 174 (15.5%)
30	97.2%	14 of 174 (8.0%)
20	96.0%	12 of 174 (6.9%)
15	94.9%	9 of 174 (5.2%)
5	88.6%	4 of 174 (2.3%)
2	70.5%	0 of 174 (0.0%)

At 2% the gate discards three in ten genuinely immunogenic peptides, which on a safety question is the error that costs something. 30% still halves the tier, 27 units to 14.

Each finding gains a "variant_binder" key with its score, kept even when it fails, because “we looked and it is not presented” is the part of a safety argument that a bare absence cannot make.

Parameters:

threshold (float)

Return type:

list

mhcmatch.vector.self_origin_risk(proteome, symbols, tissues=('Heart', 'Brain', 'Nerve', 'Lung', 'Liver', 'Kidney', 'Adrenal Gland', 'Pituitary', 'Muscle - Skeletal', 'heart muscle', 'kidney', 'liver', 'lung', 'adrenal gland', 'pituitary gland', 'skeletal muscle', 'smooth muscle', 'cerebellum', 'cerebral cortex', 'midbrain', 'basal ganglia', 'hippocampal formation', 'hypothalamus', 'amygdala', 'spinal cord', 'choroid plexus', 'retina'), min_tpm=0.25, max_subs=0, *, novel_kinds=frozenset({'frameshift', 'fusion', 'inframe_deletion', 'inframe_insertion', 'missense', 'protein_altering', 'start_lost', 'stop_lost'}), veto_tpm=5.0, graded=False, report_subs=0, report_identity=0.5, report_flank=10, report_min_length=9)

A risk callable for screen(): near-exact self origin, joined to tissue.

A register is risky when mhcmatch.Proteome.find_source() places it within max_subs of a human protein whose gene is transcribed above min_tpm in a tissue named by ESSENTIAL_TISSUES. Reasons are [{"protein", "gene", "subs", "position", "tissue", "tpm"}].

This is a near-identity test, not a similarity test, and that distinction is the whole design. The obvious alternative — score each register with mhcmatch.mimicry and flag the ones resembling a tolerance-side reference — was built and measured against this one on 1,000 viral epitopes (which cannot be self, so every firing is a false positive) and 1,000 thymic peptides from essential-tissue genes (bench/results/vector_safety_screen.md):

route	false pos.	true pos.
mimicry, masked	0.693	0.944
self origin, this one	0.020	0.940

Equal sensitivity, 35× the false positives. The reason is the one mimicry_collinear.md already records: anchor-channel similarity to a presented reference is presentation, not recognition, so an anchor-masked match fires for every peptide sharing the allele’s motif — the influenza epitope GILGFVFTL draws 14 essential-tissue hits. Nobody withdraws two-thirds of a cassette, so that route excludes nothing in practice and is not offered here.

find_source separates instead: ESDPIVAQY resolves to sp|Q8WZ42|TITIN_HUMAN at 0 substitutions, EVDPIGHLY to sp|P43357|MAGA3_HUMAN at 0 (and MAGE-A6 at 1), and GILGFVFTL to nothing at all. The 2 % that remains is not obviously noise — a viral 9-mer within one substitution of a human protein is what molecular mimicry means — and it is the floor on how specific this can get.

Two clauses, and the reason carries which one fired. "target gene" — the unit’s own gene is transcribed in an essential tissue, the MAGE-A12 case, and no register search is needed to see it. "unrelated self origin" — a register coincides within max_subs of a protein that is not the unit’s own parent. Hits to the parent are dropped, because a long peptide is native context by design and tolerance already covers it; without that exclusion the screen rejects every unit of every cassette.

Clause 1 is skipped for a product whose sequence is not in the normal proteome, and that is a category error being corrected rather than a threshold being relaxed. MAGE-A12 is a cancer-testis antigen: a shared, unmutated self protein, so its 0.33 TPM in brain caudate is the hazard exactly because the construct encodes a sequence brain tissue also presents. A somatic neoantigen is a different object — a missense, a frameshift, an inframe indel, a fusion junction all encode a sequence that is absent from normal tissue by construction — so the parent gene’s expression is not that hazard. What is a hazard for it is clause 2, and clause 2 tests it for every kind, unchanged. Measured on a 37-donor cohort, clause 1 as an unconditional rule withdrew a candidate for the fact that its parent gene exists: 10 of 37 donors lost every unit they had, and one lost 1,098 of 1,618 to clause 1 alone.

novel_kinds is mhcmatch.predict.NOVEL_PRODUCTS and is matched against Unit.kind. An isoform, a cnv locus, a wild-type or overexpressed target is the MAGE-A12 case and keeps clause 1.

Clause 2 is asked only of the registers that carry novel sequence, and for the same reason. A 27-mer unit is thirteen-fourteenths wild type by construction, and the unrestricted clause read that design as the hazard. Measured on 178 experimentally immunogenic somatic neoantigens from isalgo/pmhc_data, rebuilt as the 27-mer units they would enter a cassette as: 178 of 178 (100 %) trip clause 2, at a median of 36 self registers each. 36 is 12 + 10 + 8 + 6 — exactly the count of 8/9/10/11-mer windows of a 27-mer that cannot contain a centred mutation — and the measured self fraction tracks that geometry at every length (L=8 60.02 % against 60.0 predicted, L=9 52.6 / 52.6, L=10 44.4 / 44.4, L=11 35.2 / 35.3; 6,350 hits, 99.1 % of the geometric ceiling). At the minimal-epitope level the clause is clean: 0 of 178 mutant epitopes are in the proteome and 178 of 178 wild types are. There were essentially no genuine coincidences to find — the veto was arithmetic, not evidence.

So a window that does not contain novel sequence is structurally exempt: it is wild type, it was always going to be in the proteome, and no cassette avoids it short of not using long units. Which windows those are depends on the product, and mhcmatch.predict.TRACT_PRODUCTS is the split: a frameshift or fusion is novel from Unit.mutation_index to the end of the unit, everything else in novel_kinds at that one index. n_registers_spanning and n_hit_spanning ride on every clause-2 reason so the exemption is auditable rather than silent.

The exemption is gated on the same novel_kinds as clause 1 — one list, two rules that cannot disagree. For an isoform, a cnv or an unannotated unit every register is judged, as before, because for those the self-ness of the sequence is the finding. One case it reads generously: an in-frame fusion’s downstream tract is genuine second-parent sequence, and it is exempted with the rest of the tract rather than charged to the second gene.

An unknown or empty kind keeps clause 1 — fail closed. The screen may not exempt a unit because nobody annotated it. Every clause-1 reason therefore carries "kind", so a rejection can be read as this is a shared self antigen or as nothing said what this was. Note the one thing this cannot see: Unit.kind defaults to "missense" in the dataclass, so a Unit constructed in Python with no kind is indistinguishable from one annotated as a somatic missense and is exempted. Annotating the unit is the caller’s contract — units_from_context() and the CLI’s unit table both fill it from the pipeline header.

``veto_tpm = 5.0`` separates a veto from a cost, and it is not the same line as ``min_tpm``. min_tpm = 0.25 stays what it always was: the reporting floor, set under MAGE-A12’s 0.33 TPM so the fatal case is always visible. What 0.25 cannot also be is the exclusion line — at that level nearly every human gene is “detectable somewhere”, which is what made the screen withdraw almost everything. veto_tpm is the conventional 5 TPM “is it expressed” cut, and with graded=True a finding below it is reported with "veto": False: screen() keeps the unit, offtarget_cost() turns the finding into a per-unit cost, and mhcmatch.portfolio.compose() prices it against the response model instead of any one register vetoing a 27-mer. graded=False is the default and is the shipped veto behaviour.

``report_subs=1`` adds a third clause that reports and never withdraws. The two deaths this screen is shaped around were both near-identity, not identity: titin’s ESDPIVAQY differs from MAGE-A3’s EVDPIGHLY at four positions, and MAGE-A12 is a different gene altogether. So the exact clause 2 cannot be the whole answer — but neither can a d=1 veto, and the reason is measured rather than argued. On 178 validated immunogenic somatic neoantigens the exact clause withdraws 2 units (1.1%), while d=1 to any different expressed gene reaches 125 units (70.2%): a veto there costs two thirds of every cassette to buy a hazard the exact clause has largely already taken. Clause 3 therefore emits "veto": False unconditionally — independently of graded — so screen() keeps the unit and notes carries the finding.

Three filters keep that annotation readable, and the first carries most of it. report_min_length = 9 excludes 8-mers, because at d = 1 an 8-mer’s ball is mostly chance: 152 neighbours against 68,398,087 proteome windows in a space of 20**8 expects 0.41 coincidences per register, where a 9-mer’s 171 neighbours in 20**9 expect 0.023 – 18x fewer. Exact matching is unaffected and keeps its 8-mers, a d=0 8-mer expecting 0.0027 hits, which is why max_subs=0 can scan a length report_subs=1 must not. Then report_identity = 0.5 drops hits whose flanks are homologous to the unit’s own context (flank_identity()), because a match to a paralog is descent rather than mimicry; and the off-target gene must clear min_tpm in an essential tissue, since a hazard needs something to be expressed. report_flank is how far either side the identity is read. Feed the survivors to presented() for the fourth and last filter. report_subs is 0 by default, which is exactly the two-clause screen as previously shipped.

``d=2`` is refused, not merely discouraged. At radius 2 every expression floor from 0 TPM to 100 TPM flags 178 of 178 units, with a median of 20 off-target genes each. The hazard genuinely does live out there — EPS8L2 at d=2, titin at d=4 — and no parameter reaches it without taking the entire cassette with it. That boundary is the finding, and it is why the screen stops at 1 and hands the rest to composition.

``max_subs=0`` — exact coincidence — because the decision is per unit while the search is per register, and that multiplies. A 27-mer carries ~70 class-I registers and is withdrawn if any one of them fires, so a per-register false-positive rate that reads as small is not the rate a cassette experiences. Measured on six random 27-mers that carry no hazard, plus one burying the real titin epitope (bench/results/vector_screen_radius.md) — units falsely withdrawn, out of the six:

max_subs	9-mers	9-11	8-11
0	0	0	0
1	1	1	4

Radius 0 is clean at every length set. Radius 1 is not clean anywhere, and collapses once 8-mers enter: an 8-mer plus its 152 one-substitution neighbours is ~153 of 20⁸ sequences against the proteome’s ~68 M windows, so a chance hit per register is expected and the ~20 8-mer registers in a unit make it near-certain. Every setting still catches the titin unit, so what radius 1 buys is nothing and what it costs is most of the cassette.

``min_tpm`` defaults to 0.25 because the two precedents disagree by two orders of magnitude and the lower one is what has to be caught. Titin is 64.4 TPM in heart left ventricle and 351.4 in skeletal muscle, so any sane floor finds it. MAGE-A12 is 0.33 TPM in brain caudate and 0.31 in putamen — expression that killed two patients, and that a conventional 5-TPM “is it expressed” cut would have waved through. The floor sits just under the fatal case, not at the conventional line. It still separates: MAGE-A3’s own non-testis medians are 0.00.

What this does not catch, stated because a safety screen that oversells itself is worse than none. It would not have caught the titin event as it happened. There the cassette contained MAGE-A3, whose profile is clean — 13.4 TPM in testis, 0.00 elsewhere — and the cross-reactive titin peptide was never in the construct. Four TCR-facing substitutions separate the two, so no distance threshold reaches it from the candidate, and the affinity-enhanced TCR that bridged them was the actual cause. What this catches is the adjacent and commoner failure: a register that is a self peptide from an essential-tissue gene, and a target gene like MAGE-A12 that is transcribed where it was assumed silent.

symbols is {accession: gene} from mhcmatch.proteome.gene_symbols(path, key="accession")(); the search names proteins as sp|P43357|MAGA3_HUMAN and mhcmatch.expression.safety_profile() is keyed on MAGEA3. It is required rather than defaulted because a missing map resolves nothing and so returns “no risk” for every peptide — the one wrong answer this must never give quietly.

proteome needs find_sources(), the batch form, and at max_subs=0 find_exact_sources() is used where the object has it. The dispatch is no longer about cost – one seqtree.TextIndex answers every radius off the same 0.7 s build – but it is kept because it is how the code says, in one line, that this screen asks an exact question. The returned callable carries a prepare(registers) that screen() hands every register of every unit, so the whole candidate list resolves in one query rather than one per unit — screen everything in one process, never one unit per invocation.

Parameters:

• min_tpm (float)

• max_subs (int)

• veto_tpm (float)

• graded (bool)

• report_subs (int)

• report_identity (float)

• report_flank (int)

• report_min_length (int)

mhcmatch.vector.select(candidates, n0, cls=None, block=None)

Apply the per-allotype stopping rule (module docstring) to ranked candidates.

Candidates are grouped by Unit.allele and sorted by Unit.p descending within each group, then each group grows while p_next > S_a / (n0 + n_a). The first unit on an allotype is always taken: S_a = 0 makes the threshold 0, and any candidate with p > 0 clears it.

n0 is per-allotype capacity and must be positive. There is no default, deliberately — the literature does not fix it (see the module docstring), so a caller who has not decided what to assume has not finished designing the cassette.

block chooses what the budget saturates against. The default is the allotype, which is the rule as shipped and as described in the module docstring. Passing a callable Unit -> hashable blocks on something else, and the intended use is a pair: allotype together with the mechanism a unit was selected on, e.g. block=lambda u: (u.allele, corner[u.peptide]). The arithmetic does not care what the key means; what it assumes is that two units sharing a key share a way of failing. Diversification across whatever that is falls out of the saturation, exactly as it does across allotypes — see mhcmatch.portfolio for the response model this is the greedy rule for, and for the measured intra-patient correlation that motivates blocking on more than the allotype.

Parameters:

• n0 (float)

• cls (str | None)

Return type:

Selection

mhcmatch.vector.assemble(units, linker=None)

Lay units out in the order given, joined by linker, with no prediction of any kind.

This is the assembly step on its own: the caller has already decided which units, in what order, and what to put between them, and wants the construct. order() is the same assembly with the ordering and the linker chosen by measurement against the recipient’s allotypes, and it needs a binder callable to do that; this needs nothing.

linker is a preset name from LINKERS, an explicit residue string, or None. The returned Cassette carries no junction evidence, because none was computed — junctions is empty and cost is 0.0, which is unknown and not clean. Run scan_junctions() on it if the junctions matter, or use order() from the start.

boundaries tiles the sequence exactly, so the linker spans are the gaps between them.

>>> us = [Unit("AAAAAAAAAKAAAAAAAAAAAAAAAAA", 9, "G1", "HLA-A*02:01", 0.3),
...       Unit("CCCCCCCCCRCCCCCCCCCCCCCCCCC", 9, "G2", "HLA-A*02:01", 0.1)]
>>> cas = assemble(us, "GS10")
>>> cas.spacer
'GGSGGGGSGG'
>>> cas.sequence[cas.boundaries[0][1]:cas.boundaries[1][0]]
'GGSGGGGSGG'

Return type:

Cassette

mhcmatch.vector.order(units, binder, spacers=(None, 'GPGPG', 'GGS', 'AAA', 'AAY', 'GPGPGPG', 'HHAA', 'AAL'), lengths=(8, 9, 10, 11), alleles=None, threshold=None, objective='sum', binder_threshold=None, linker=None)

Choose a spacer and an ordering that minimise predicted junctional binding.

``objective`` matters and the two choices disagree, so it is explicit.

"sum" total of the strongest predicted binder at each junction. A junction is a hazard if

one good binder forms there, which is pVACvector’s logic (PMID 31907209). The junction count is n-1 whatever the spacer, but a longer spacer creates more registers per junction and so a stochastically larger maximum — this objective therefore has a real bias toward the shortest spacer, up to and including none.

"rate" predicted binders per register, needing binder_threshold. Length-neutral, and it

is the metric a junction sweep naturally reports.

On one measured payload the two picked different spacers — "sum" chose no spacer where a rate sweep put AAA ahead of it — so a caller who has not chosen has not finished designing.

Spacers are tried in spacers order and the first one whose worst junction falls at or below threshold wins; with threshold=None every spacer is tried and the one with the lowest total junction cost wins. Because SPACERS leads with None, a cassette that needs no spacer gets none — which is the right default, since every inserted residue is translated sequence that could itself form a binder.

For each spacer the layout is an open path over a complete graph whose edge i -> j costs the strongest predicted binder spanning that junction, solved by _greedy_2opt().

linker= pins one linker instead of sweeping — a preset name from LINKERS, an explicit residue string, or None/"none" for direct concatenation. Every entry of spacers is resolved the same way, so spacers=("none", "AAY", "GS10") is a sweep over three presets. Pinning is what a caller does when the construct format is already decided and the question is only what order the units go in; sweeping is what a caller does when it is not. Use assemble() when the order is decided too and no prediction is wanted at all.

binder(peptides, alleles) -> [float], higher meaning a stronger binder. Use store_binder() to build one from a Store.

Parameters:

• threshold (float | None)

• objective (str)

• binder_threshold (float | None)

Return type:

Cassette

mhcmatch.vector.slippery_sites(cds)

Codon positions where N1-methylpseudouridine drives +1 ribosomal frameshifting.

Returns [{codon_index, nt_offset, codon, next_codon}, ...].

m1Ψ is not translationally neutral. Mulroney et al. (Nature 2024, PMID 38057663) measured +1 frameshifting at ~8% of the in-frame product in m1Ψ mRNA, localised to a slippery motif — m1Ψ m1Ψ m1Ψ X with X = m1Ψ or C at the first position of the following codon, i.e. a TTT codon followed by a codon starting T or C. Six such sites sit in the BNT162b2 spike coding sequence, and BNT162b2-vaccinated humans mounted a significantly higher IFN-γ response against the +1 frameshifted product than controls.

This matters far more for a designed polyepitope than for a natural ORF, for two reasons. A concatemer has many more codon-boundary junctions per kilobase, and the residues at those junctions are the designer’s choice — glycine/serine linkers are encoded by GGN/AGY/TCN, which is exactly how U-runs end up at seams. And the consequence is worse: a frameshift inside a polyepitope does not merely lose protein, it translates an entire downstream out-of-frame cassette that is itself presented, so the construct delivers a second, unintended and unscreened antigen payload.

Scanning is therefore mandatory for any m1Ψ construct and pointless for an unmodified-uridine one — BioNTech’s cancer platform deliberately uses unmodified uridine, so which applies depends on the platform, not the sequence.

Only the codon-aligned motif as published is reported. Whether non-codon-aligned U-runs also induce slippage was not characterised in that work, so it is not guessed at here.

Parameters:

cds (str)

Return type:

list

mhcmatch.vector.deslip(cds)

(cds, n_fixed) — remove every slippery_sites() motif synonymously.

TTT and TTC both encode phenylalanine, so rewriting the upstream codon breaks the U-run without touching the protein and without disturbing the downstream codon’s own optimisation. This is the fix Mulroney et al. validated: single U*187C / U*208C substitutions strongly reduced frameshifting and the double mutant produced none detectable.

Idempotent — a second call finds nothing to do.

Parameters:

cds (str)

Return type:

tuple

mhcmatch.vector.CODON_USAGE_HUMAN: dict = {'AAA': ('K', 24.4), 'AAC': ('N', 19.1), 'AAG': ('K', 31.9), 'AAT': ('N', 17.0), 'ACA': ('T', 15.1), 'ACC': ('T', 18.9), 'ACG': ('T', 6.1), 'ACT': ('T', 13.1), 'AGA': ('R', 12.2), 'AGC': ('S', 19.5), 'AGG': ('R', 12.0), 'AGT': ('S', 12.1), 'ATA': ('I', 7.5), 'ATC': ('I', 20.8), 'ATG': ('M', 22.0), 'ATT': ('I', 16.0), 'CAA': ('Q', 12.3), 'CAC': ('H', 15.1), 'CAG': ('Q', 34.2), 'CAT': ('H', 10.9), 'CCA': ('P', 16.9), 'CCC': ('P', 19.8), 'CCG': ('P', 6.9), 'CCT': ('P', 17.5), 'CGA': ('R', 6.2), 'CGC': ('R', 10.4), 'CGG': ('R', 11.4), 'CGT': ('R', 4.5), 'CTA': ('L', 7.2), 'CTC': ('L', 19.6), 'CTG': ('L', 39.6), 'CTT': ('L', 13.2), 'GAA': ('E', 29.0), 'GAC': ('D', 25.1), 'GAG': ('E', 39.6), 'GAT': ('D', 21.8), 'GCA': ('A', 15.8), 'GCC': ('A', 27.7), 'GCG': ('A', 7.4), 'GCT': ('A', 18.4), 'GGA': ('G', 16.5), 'GGC': ('G', 22.2), 'GGG': ('G', 16.5), 'GGT': ('G', 10.8), 'GTA': ('V', 7.1), 'GTC': ('V', 14.5), 'GTG': ('V', 28.1), 'GTT': ('V', 11.0), 'TAA': ('*', 1.0), 'TAC': ('Y', 15.3), 'TAG': ('*', 0.8), 'TAT': ('Y', 12.2), 'TCA': ('S', 12.2), 'TCC': ('S', 17.7), 'TCG': ('S', 4.4), 'TCT': ('S', 15.2), 'TGA': ('*', 1.6), 'TGC': ('C', 12.6), 'TGG': ('W', 13.2), 'TGT': ('C', 10.6), 'TTA': ('L', 7.7), 'TTC': ('F', 20.3), 'TTG': ('L', 12.9), 'TTT': ('F', 17.6)}

Homo sapiens codon usage as {codon: (amino acid, occurrences per thousand codons)}. Kazusa Codon Usage Database (https://www.kazusa.or.jp/codon/), Homo sapiens [gbpri] – 93,487 CDSs, 40,662,582 codons; retrieved 2026-09-20.

Stored as the measured frequencies, not as a pre-reduced “best codon per residue” map, so back_translate()’s choice is derived from data under a stated rule and can be re-derived under a different one. Substituting a host’s own table is then a one-argument change.

mhcmatch.vector.MAX_HOMOPOLYMER: int = 4

Run of one nucleotide past which back_translate() reaches for a rarer synonymous codon. Homopolymers are a synthesis constraint, not a translation one: vendors reject or mis-assemble long single-base runs, and a spacered concatemer manufactures them directly (AAA and GPGPG are both in SPACERS). Four is the conventional screening floor.

It is a target, not a guarantee, because the choice is greedy and per-codon. Measured over 5,000 random 20-60mers under the default table: longest run 6, and 84% of sequences at or below 4; the same peptides back-translated by most-frequent-codon alone reach 13, with 6% above 6.

mhcmatch.vector.translate(cds)

Amino-acid sequence of a coding sequence, * for a stop.

Exists so “synonymous” is checkable rather than asserted – deslip() and back_translate() both claim it, and a caller who supplies their own codon table needs the same check. U is read as T; a trailing partial codon is ignored.

Parameters:

cds (str)

Return type:

str

mhcmatch.vector.back_translate(peptide, usage=None, *, avoid_slip=True, max_run=4)

Coding sequence for peptide – the cassette’s nucleotide half.

Highest-usage synonymous codon per residue, backing off to the next one whenever the first would extend a single-nucleotide run past max_run, then deslip() to remove the m1-psi +1-frameshift motif. Deterministic: the same peptide and table always give the same CDS.

The backoff is greedy, so max_run is a target rather than a bound – see MAX_HOMOPOLYMER for what it is measured to be worth.

What this is not is a codon optimiser. It fixes the two things that make a polyepitope construct fail where a natural ORF would not – the frameshift motif (slippery_sites(), which a concatemer hits far more often because the designer chooses the seam residues) and synthesis-hostile homopolymers (which spacers like AAA manufacture directly). It does not touch GC content, secondary structure, splice sites or CpG, and a manufacturer’s own optimiser should be preferred where one is available; this exists so a cassette ships with a usable CDS rather than none.

Emits the epitope cassette only – no start codon, no stop, no leader, no trafficking domain – matching Cassette.sequence, because those flanks are the vector’s, not the payload’s.

>>> cds = back_translate("SIINFEKL")
>>> translate(cds)
'SIINFEKL'
>>> slippery_sites(cds)
[]

Parameters:

• peptide (str)

• usage (dict)

• avoid_slip (bool)

• max_run (int)

Return type:

str

class mhcmatch.vector.MRNA(sequence, protein, cds, parts=<factory>, linker=None)

Bases: object

An assembled mRNA construct, its parts, and the checks that were actually run on it.

sequence is the whole molecule 5’→3’ in DNA letters — every element the caller supplied, in order, with the coding sequence this library generated in the middle. as_rna() gives the same molecule in RNA letters.

parts tiles sequence exactly: consecutive, non-overlapping, 0-based half-open, and the concatenation of every part’s slice is the whole molecule. Every part carries the same keys — kind, name, start, end, aa_start, aa_end — with the two amino-acid offsets None on a part that is not translated, so a reader never has to know which kind it is holding before it can index one. That is the property that makes the record auditable — an element that silently went missing shows up as a gap, and one that was added twice shows up as an overlap, neither of which a length check would catch.

Parameters:

• sequence (str)

• protein (str)

• cds (str)

• parts (list)

• linker (str | None)

sequence: str

protein: str

cds: str

parts: list

linker: str | None = None

as_rna()

The construct in RNA letters. Nothing else changes — this is a transliteration.

Return type:

str

of_kind(kind)

Every part of one kind, in order: utr5, start, leader, unit, linker, trailer, stop, utr3, poly_a.

Parameters:

kind (str)

Return type:

list

slice_of(part)

The nucleotides one part occupies. Exists so a caller never re-derives an offset.

Parameters:

part (dict)

Return type:

str

property checks: dict

What was verified about the construct, as numbers rather than as a pass/fail.

translates is the one that must hold: the coding sequence read back in the frame the construct sets it in gives exactly protein. It is the check that catches a frame broken by a supplied 5’ element, which is otherwise invisible.

The sequence composition figures are over the coding sequence only — a poly(A) tail is a homopolymer by construction and would swamp both. length_nt is the whole molecule.

mhcmatch.vector.mrna(cassette, linker=None, *, leader='', trailer='', utr5='', utr3='', poly_a=0, start=True, stop='TGA', usage=None, avoid_slip=True, max_run=4)

Build the mRNA for a cassette, joined by a chosen linker. The assembly step, end to end.

cassette is a Cassette, an iterable of Unit, an iterable of peptide strings, or a single peptide string. linker is a preset name from LINKERS, explicit residues, or None; given a Cassette and no linker, the cassette’s own spacer is kept, and given both, the units are re-joined by the linker asked for.

What the library supplies and what it refuses to. The coding sequence is generated here: the whole open reading frame is back-translated in one pass, so homopolymer avoidance and the m1Ψ frameshift repair (back_translate(), deslip()) work across the seams the designer created rather than within each unit separately. The backbone — 5’ and 3’ UTRs, a signal peptide or a trafficking domain, the tail — is not supplied and defaults to nothing, because those sequences belong to a particular vector and inventing a plausible one is worse than emitting none. Pass the ones your construct uses and they are placed, checked and mapped:

cas  = vector.order(units, binder=b, linker="GS10")
m    = vector.mrna(cas, leader=SEC, trailer=MITD, utr5=U5, utr3=U3, poly_a=100)
m.checks["translates"]                      # the frame survived every element
[p["name"] for p in m.of_kind("unit")]      # what is in it, in the order it is in

leader and trailer are amino acids and are translated in frame with the payload — a secretory signal and an MHC trafficking domain are the two that matter for a class-I cassette, and they are the reason translates is worth checking rather than assuming. utr5/utr3 are nucleotides and are not translated.

start prepends an initiator methionine unless the open reading frame already begins with one, so the codon is real and MRNA.protein shows the residue it encodes rather than hiding it. stop is one codon and is checked to be one; pass stop="" when the stop lives in the vector rather than in the payload. poly_a is a count of adenosines; a template-encoded tail is a property of the vector and its length is the caller’s to set.

>>> m = mrna(["SIINFEKLA", "KVAELVHFL"], "GS10")
>>> m.protein
'MSIINFEKLAGGSGGGGSGGKVAELVHFL'
>>> m.checks["translates"], m.checks["n_units"]
(True, 2)
>>> [(p["kind"], p["name"]) for p in m.parts]
[('start', 'start'), ('unit', 'u1'), ('linker', 'GGSGGGGSGG'), ('unit', 'u2'), ('stop', 'TGA')]

Parameters:

• leader (str)

• trailer (str)

• utr5 (str)

• utr3 (str)

• poly_a (int)

• start (bool)

• stop (str)

• usage (dict)

• avoid_slip (bool)

• max_run (int)

Return type:

MRNA

mhcmatch.vector.store_binder(store, alleles, cls='mhc1')

A binder callable over a Store: -log10(%rank) of the best allele, so higher means a stronger predicted binder.

Kept as a thin adapter, and out of order(), so the layout logic stays testable with no panel, no download and no allele list.

Parameters:

cls (str)

mhcmatch.vector.rebuild(cassette, **kw)

Re-lay an existing cassette’s units under different settings, keeping the units fixed.

The point of comparison is the incumbent: re-order and re-space the same payload to separate “these units were the wrong choice” from “these units were laid out badly”.

Parameters:

cassette (Cassette)

Return type:

Cassette

mhcmatch.vector.from_sequence(sequence, spacer, lengths=(8, 9, 10, 11))

Split an existing cassette on a known spacer into pseudo-units, for auditing.

Only usable when the spacer is unambiguous — a cassette joined by alternating tokens has to be split by the caller, who knows the grammar. mutation_index is unknown from sequence alone and is set to the window centre; p is 0.0. These units are for junction scanning, not selection.

Parameters:

• sequence (str)

• spacer (str)

Return type:

list

mhcmatch.vector.MHC2_MAP_LENGTHS: tuple = (12, 13, 14, 15, 16, 17, 18, 19, 20)

Class-II ligand lengths scanned when mapping a cassette. Wider than MHC2_JUNCTION_LENGTHS, which exists to score junctions and only needs the shortest windows a core can be read from: a map is annotating what a cassette actually presents, and a class-II ligand runs to 25.

class mhcmatch.vector.Feature(id, kind, start, end, seq, cls='', allele='', rank=nan, unit=0, gene='', core_start=0, core_end=0, core='', overlaps=())

Bases: object

One annotated span of an assembled cassette, in 1-based inclusive amino-acid coordinates.

kind is unit (a vaccine unit), linker (the spacer between two of them) or epitope (a predicted binder). Units and linkers tile the cassette exactly; epitopes overlay it and may span a junction, which is the case unit = 0 marks.

Parameters:

• id (str)

• kind (str)

• start (int)

• end (int)

• seq (str)

• cls (str)

• allele (str)

• rank (float)

• unit (int)

• gene (str)

• core_start (int)

• core_end (int)

• core (str)

• overlaps (tuple)

id: str

kind: str

start: int

end: int

seq: str

cls: str = ''

allele: str = ''

rank: float = nan

unit: int = 0

gene: str = ''

core_start: int = 0

core_end: int = 0

core: str = ''

The binding core itself (mhcmatch.store.binding_core()), 9 residues, both classes. core_start/core_end stay class-II-only because a class-I core is not a contiguous span – it drops the bulge – so it has no cassette coordinates to give.

overlaps: tuple = ()

property length: int

Inclusive span, end - start + 1 – residue count, not a half-open Python length.

mhcmatch.vector.store_ranker(store, alleles, cls='mhc1', calibrated=True)

A ranker callable over a Store: [(allele, %rank), ...] per peptide, one entry per allele that presents it.

Distinct from store_binder(), which collapses to the best allele because a layout cost only needs to know whether some binder forms. A map needs the allele: a heterozygote presents the same peptide on two molecules, and that is two facts about the cassette rather than one.

Parameters:

• cls (str)

• calibrated (bool)

mhcmatch.vector.RANK_DEFAULT_TIER: str = 'weak'

What the cassette map annotates unless asked otherwise. Weak, because the map is a report – it selects nothing and removes nothing – and under-reporting help a construct genuinely carries is the more costly error here.

Not the same thing as mhcmatch.predict.RANK_DEFAULT_TIER, which is "none". Same name, two jobs: this one chooses a label, that one chooses what to delete. A report should default to the generous convention and a filter should default to off, so the two differ on purpose – see --block-live in CLAUDE.md for the last time one name covered two knobs.

mhcmatch.vector.rank_cutoffs(tier='weak')

{"mhc1": f, "mhc2": f} for "strong" or "weak"; raises on anything else.

Parameters:

tier (str)

Return type:

dict

mhcmatch.vector.epitope_map(cassette, ranker1=None, ranker2=None, threshold=2.0, lengths1=(8, 9, 10, 11), lengths2=(12, 13, 14, 15, 16, 17, 18, 19, 20), threshold2=None, stats=None)

Annotate an assembled cassette: units, linkers, predicted epitopes, and which class-I and class-II epitopes overlap each other.

Why the overlap is the point and not a decoration. A cassette that carries a CD8 epitope and borrows its CD4 help from an unrelated universal helper (PADRE, HBVcore) raises no T-cell response against the tumour antigen on the class-II side. Kissick et al. built one 27-mer around the HLA-A*02:01 SIM2_237-245 epitope so that a class-II epitope from the same protein overlapped it, and it replaced the exogenous HBVcore helper outright: the long peptide alone raised both the CD8 IFN-γ recall response to the 9-mer and a CD4 IL-2 response to SIM2_240-254, and 137 class-II binders were predicted across DR/DP/DQ from that one 27-mer (PLoS One 2014;9(4):e93231, PMID 24690990, doi:10.1371/journal.pone.0093231). A unit whose class-I epitope has no overlapping class-II epitope is the configuration that needed the borrowed helper, and this map is what says which units those are.

ranker1 / ranker2 are ranker(peptides) -> [[(allele, %rank), ...], ...], per class; store_ranker() builds one from a Store. Either may be None, which simply omits that class. Injected for the same reason order()’s binder is — the whole map is testable with no panel, no download and no calibration.

Every (peptide, allele) at or below threshold %rank is its own Feature, so a peptide presented by two of the patient’s alleles appears twice. That is deliberate: at a heterozygous locus the two molecules are two independent presentation events, they are what the per-allotype capacity in select() is spent on, and collapsing them would under-count the cassette’s coverage of exactly the patients it was personalised for.

Coordinates are 1-based inclusive over Cassette.sequence, which is the epitope cassette only — no start codon, no leader, no tag. An mRNA construct that adds those must offset.

Parameters:

• cassette (Cassette)

• threshold (float)

• threshold2 (float | None)

• stats (dict | None)

Return type:

list

mhcmatch.vector.MAP_COLUMNS: tuple = ('id', 'kind', 'start', 'end', 'length', 'seq', 'cls', 'allele', 'rank', 'unit', 'gene', 'core', 'core_start', 'core_end', 'overlaps', 'n_overlaps')

Column order of the cassette map, one source of truth for the TSV and the JSON.

mhcmatch.vector.map_rows(features)

The map as plain dicts in MAP_COLUMNS order — the payload both writers share.

Return type:

list

mhcmatch.vector.map_summary(cassette, features)

Per-unit coverage: how many class-I and class-II epitopes, on how many allotypes, and whether the unit’s class-I epitopes have class-II help from within the same unit.

The last column is the one a reviewer reads first. A unit with class-I epitopes and no overlapping class-II epitope is the configuration that needed a borrowed universal helper.

Parameters:

cassette (Cassette)

Return type:

dict

mhcmatch.vector.write_map(cassette, features, tsv_path=None, json_path=None)

Write the cassette map as TSV and/or JSON, and return the summary dict.

The TSV is the flat table — one row per feature, one value per cell, so it sorts and joins. The JSON carries the same rows plus the per-unit summary and the cassette sequence, which is what a viewer needs to draw the thing without recomputing anything.

Parameters:

• cassette (Cassette)

• tsv_path (str | None)

• json_path (str | None)

Return type:

dict

mhcmatch.portfolio module

The composition layer above select(): objective-space geometry (Pareto front, crowding, hull membership, Chebyshev scalarization) and the block response model that says what a proposed cassette is worth. Narrative and worked examples: Cassette composition.

Cassette composition: the objective geometry and the response model behind vector.select.

A cassette is a set, and the quantity that decides whether it works is not how good its units are on average but whether at least one of them elicits a response — better, at least k. Sorting by a score and keeping the top m answers that question correctly only if the units respond independently. They do not, and this module is what the difference costs.

Two facts, both measured rather than assumed, and both recorded in the benchmark repository:

Response counts within a patient are over-dispersed. On the adjuvant TNBC mRNA vaccine trial of Sahin et al. (Nature 2026;651:1088-1096, PMID 41708868) — 13 patients, 20 assayed units each, a pooled per-unit response rate of 19.0% — the intra-patient correlation is rho = 0.124 at p = 1.0e-3, 3.45x the binomial variance. TESLA gives rho = 0.024 and HiTIDE rho = 0.010. The dispersion is scale-dependent: a 4,967-candidate screening pool spanning every allotype shows none, because a pool that wide averages its blocks out. A cassette cannot, and a shortlist ranked on one score averages them least of all, since ranking is what concentrates the units in the first place. Use betabinom_rho() to measure it on your own readout before assuming a value.

A weighted sum cannot select part of what the objectives describe. For any beta >= 0, top-m by beta @ z selects only candidates on the upper convex hull of the objective cloud, so 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 whatsoever. That limit belongs to the weighted sum, not to scalarization: chebyshev_score() reaches the whole front, and so, in principle, does any sufficiently rich nonlinear model. What none of them escapes is separability — top-m by any pointwise score maximises a modular set function, and P(>= k | S) is not modular whenever two units share a block. That is a property of the selection rule, not of the scorer, so it cannot be fitted away.

The rule this module supports is already in mhcmatch.vector.select(); pass its block argument a key that pairs the allotype with the mechanism a unit was selected on. Everything here is diagnostic: it says what a proposed cassette is worth and why, and fits nothing.

linearly_supported and betabinom_rho need SciPy, which is not a hard dependency; both import it lazily and say so if it is missing.

mhcmatch.portfolio.pareto_front(Z, *, batch_bytes=8388608)

Boolean mask of non-dominated rows of Z (n x K), higher is better on every column.

Orient every column that way before calling: a %rank is lower-is-better and has to enter as -log10(rank) or the front is the wrong end of the cloud. Dominance comparisons use bounded array blocks (estimated batch_bytes working set), with one reference row as the minimum block. No candidate-by-candidate matrix is retained.

Parameters:

batch_bytes (int)

Return type:

numpy.ndarray

mhcmatch.portfolio.nondominated_rank(Z, max_fronts=3)

Front index per row, 0 = first front. Rows beyond max_fronts share the last index.

Parameters:

max_fronts (int)

Return type:

numpy.ndarray

mhcmatch.portfolio.crowding_distance(Z)

NSGA-II crowding distance: per-objective normalised gap to the two flanking neighbours.

Boundary rows on any objective get inf, which is what keeps the extremes of the front from being pruned. Use it to break ties within a front, never across fronts.

Return type:

numpy.ndarray

mhcmatch.portfolio.linearly_supported(Z, i)

Is row i ranked first by some beta >= 0? Exact, by linear-programming feasibility.

True exactly when Z[i] lies on the upper convex hull. A Pareto-efficient row inside the hull returns False: no weighting of the objectives ever puts it on top, and tuning them is wasted effort. Use chebyshev_score() for those. The small feasibility problem uses one HiGHS thread; solver failures raise instead of reporting an unsupported point. An incompatible pre-existing HiGHS scheduler is a failure too; this function does not reset process-global state owned by another caller.

Parameters:

i (int)

Return type:

bool

mhcmatch.portfolio.chebyshev_score(Z, weights, ideal=None, aug=0.001)

Augmented weighted Chebyshev score, higher is better.

s(z) = -(max_k w_k (z*_k - z_k) + aug * sum_k w_k (z*_k - z_k)) with z* the ideal point (per-column max, plus a nudge, unless given). Unlike a weighted sum this reaches every Pareto-efficient point for some w — the classical guarantee (Bowman 1976; Steuer and Choo 1983) — including the concave stretches of the front a linear score cannot support. aug breaks ties towards properly efficient points and should stay small.

Parameters:

aug (float)

Return type:

numpy.ndarray

mhcmatch.portfolio.corner(Z, groups=None)

Assign each row the objective it is relatively strongest on — its mechanism corner.

Ranks within the pool per column and takes the arg-max, so the answer is scale-free and does not move when one objective is rescaled. groups maps column index -> label to pool related columns (three agretopicity parameterisations are one mechanism, not three).

This is a proxy for a latent variable, not the variable: it says which axis a candidate stands out on, which is a defensible stand-in for why it might work and nothing more.

Return type:

numpy.ndarray

mhcmatch.portfolio.survival(p, block, q)

P(X >= j) for every j, under y_i = B_{block(i)} * eps_i, B_b ~ Bern(q_b).

p is the marginal per-unit probability a ranker reports, so the unit-specific term is eps_i ~ Bern(p_i / q_{block(i)}). That requires p_i <= q_{block(i)}: a unit cannot respond more often than its own block is live. Clipping there would understate the marginal for exactly the strongest units, so this raises instead.

Exact, and cheap. X = sum_b B_b S_b with S_b a Poisson binomial over the block’s units, and B_b S_b has pmf (1 - q_b) delta_0 + q_b pmf(S_b). The blocks are independent, so the pmf of X is the convolution of those – no 2^B enumeration over live sets and no Monte Carlo. Returns an array of length len(p) + 1; element k is P(X >= k), so survival(...)[0] == 1.

With every unit in one block the tail is capped at q however large the cassette grows, which is the whole reason to block on more than the allotype.

>>> float(survival([0.5, 0.5], [0, 1], 1.0)[1])
0.75

Return type:

numpy.ndarray

mhcmatch.portfolio.p_at_least(p, block, q, k=1)

P(at least k responses) under the block model. See survival(), which it reads.

>>> round(p_at_least([0.5, 0.5], [0, 0], 1.0), 4)
0.75

Parameters:

k (int)

Return type:

float

mhcmatch.portfolio.visibility(p, block=None)

V = 1 - prod_i (1 - p_i) over a set nobody chose, and what losing one allotype costs.

The observational read, and it carries no design parameter. survival() and p_at_least() answer this question for a cassette some rule selected, under a block model with a live probability q. This answers it for the set a tumour happens to present: every unit, weight one, no q, no coupling, no gamma and no rho. A tumour is under selection too, but towards escape, and that dynamics is not modelled here — so there is no designer’s risk to price and nothing to maximise over. Keeping the two apart is the point: see mhcmatch.cassette.goal_energy() for the design side, which is where those parameters belong.

V is the probability that at least one presented unit elicits a response, so it saturates: once a few strong antigens are present, further units move it very little. That is what makes it a whole-tumour read rather than a count, and why it decouples from mutational burden where a sum does not.

HLA loss of heterozygosity, 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 same sum over the allotype’s own units, visibility after the loss is 1 - exp(S - S_a); the worst single loss is the allotype minimising it, which is the one whose S_a is most negative. loh_cost is how much visibility the tumour’s most protective allotype carries — the McGranahan escape route, priced on the same scale as V.

Returns visibility, log1m and n_units, and with block also visibility_loh, loh_cost, lost_allotype and n_allotypes. log1m is S in nats — negative, linear in the units and unbounded below, so it separates the saturated top of a cohort that visibility itself cannot.

Why this is a closed form and not a call into p_at_least(). At q = 1 that function returns exactly this number whatever the blocking, so the design record asked for the call rather than a second spelling of one quantity — normally the right instinct, and this repository names that hazard in five other places. It does not hold here: survival() convolves one Poisson binomial per block, which is quadratic in the unit count, and V is defined over a whole presented set — hundreds to thousands of units, on every donor of a cohort. The closed form is linear. The two spellings are held together by a test that asserts they agree rather than by one calling the other (test_visibility_agrees_with_p_at_least_where_the_two_overlap).

>>> round(visibility([0.5, 0.5])["visibility"], 4)
0.75
>>> v = visibility([0.5, 0.5], ["A", "B"])
>>> round(v["visibility_loh"], 4), v["lost_allotype"]
(0.5, 'A')

Return type:

dict

mhcmatch.portfolio.n_effective(p, p_ge1)

Independent Bernoulli(mean p) units that would buy the same P(>= 1).

n_eff <= len(p), with equality only under independence. The bound is imposed rather than reported: block correlation can only lower P(>= 1), so exceeding it is Monte-Carlo noise.

Parameters:

p_ge1 (float)

Return type:

float

mhcmatch.portfolio.coverage(labels, universe=None)

How evenly a cassette spreads over allotypes: counts, Gini, and share of maximum entropy.

universe is the donor’s distinct allotypes, and passing it is the whole point when the donor is homozygous. A patient homozygous at B has five distinct class-I allotypes, not six, so an even cassette over five is perfectly even; scoring it against a denominator of six would report a genotype as a design flaw. Allotypes in universe with no unit are counted as zeros, which is exactly the inequality the index should see.

Returns gini (0 = every allotype equally covered, -> 1 = all units on one) and entropy_ratio = H / log(|universe|), the “% of maximum entropy” reading. With one allotype both are defined as perfectly even, because there is nothing to be uneven about.

>>> c = coverage(["A", "A", "B", "B"])
>>> round(c["gini"], 6), round(c["entropy_ratio"], 6)
(0.0, 1.0)
>>> round(coverage(["A", "A", "A", "B"])["entropy_ratio"], 4)
0.8113

Return type:

dict

mhcmatch.portfolio.dispersion(m, k)

Observed vs independent-Bernoulli variance of the per-patient response rate.

One (m, k) per patient: m units assayed, k positive. Descriptive only — the ratio is inflated when patients differ widely in m, because a patient with m = 1 contributes k/m in {0, 1} whatever the biology. Use betabinom_rho() to test.

Return type:

dict

mhcmatch.portfolio.betabinom_rho(m, k, profile=True)

Intra-patient correlation rho and a likelihood-ratio test against the binomial.

The null rho = 0 sits on the boundary of the parameter space, so the reference is the 50:50 mixture of chi2_0 and chi2_1 (Self and Liang 1987) — p = P(chi2_1 > D) / 2. Simulation under the null puts the realised type-I error below nominal at the cohort sizes this is used at (0.022 at alpha = 0.05 for 13 patients x 20 units), so the p-value is conservative.

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

profile=True (the default) holds p at the pooled rate and profiles the likelihood over rho alone; profile=False maximises over (p, rho) jointly. The two agree to about a thousandth on the cohorts this has been run on — the pooled rate is very nearly the joint maximiser — and the profile form is the default because it is a one-dimensional bounded search with no starting point to get wrong. The joint form exists so a caller who needs the fitted p reported beside rho does not have to write a second estimator, which is how this function acquired a duplicate in the first place.

Parameters:

profile (bool)

Return type:

dict

class mhcmatch.portfolio.Composition(units=<factory>, arms=<factory>, trace=<factory>, coverage=<factory>)

Bases: object

A cassette built to a set of quotas, and the arithmetic that justifies each slot.

arms carries one entry per arm: the units chosen, the slot budget, the response target, and the attained P(X >= target). trace is one row per greedy step with the marginal gain the step bought, so a cassette that cannot explain why its 12th unit is in and the 13th is out is not shipped.

Parameters:

• units (list)

• arms (dict)

• trace (list)

• coverage (dict)

units: list

arms: dict

trace: list

coverage: dict

property joint: float

prod_arm P(X_arm >= target_arm) – every quota met at once, arms independent.

The independence is across arms, not across units: within an arm the block model is carried in full. Two arms sharing a patient are not independent in truth, so read this as the product of three separately-meaningful numbers rather than a calibrated joint.

mhcmatch.portfolio.compose(candidates, quotas, q, block=None, arm=None, weight_evenness=0.0, universe=None, cost=None, weight_cost=0.0)

Fill each arm’s slots to maximise P(at least target responses), not the mean score.

quotas is {arm: (slots, target)} – e.g. {"mhc1": (8, 2), "mhc2": (4, 1), "nonconventional": (3, 1)} reads eight class-I slots, of which at least two should respond. q is the per-block live probability of the response model (survival()); block a callable Unit -> hashable (default the allotype); arm a callable Unit -> str (default default_arm()).

This is not top-m by score, and the difference is the point. P(X >= k) is not a modular set function whenever two units share a block, so no pointwise score – however well fitted – can be sorted to maximise it. The greedy step here takes the unit with the largest gain in P(X >= target), and because a block that is already represented contributes less than a fresh one, diversification across allotypes and mechanisms falls out of the objective rather than being bolted on as a rule. A cassette of eight units all restricted to the same allotype is capped at q for that block no matter how good the eight are.

weight_evenness adds w * delta(H / H_max) over the arm’s allotypes (coverage()), for when spreading matters beyond what the response model already pays for – manufacturing risk, an uncertain genotype, a donor whose typing is provisional. Pass universe (the donor’s distinct allotypes) so homozygosity is not scored as a design flaw. Default 0: the block model already prefers spread, and stacking a second diversity term on top of it double-counts unless you mean it.

cost is a callable Unit -> float and weight_cost the price the objective pays for it: the greedy value becomes P(X >= target) - weight_cost * sum(cost(u)). The intended supply is mhcmatch.vector.offtarget_cost(), the size of a unit’s off-target fingerprint under the graded safety screen, so a unit with a sub-veto essential-tissue hit is priced rather than withdrawn. The cost is charged to the objective, never to Unit.p: p is a calibrated marginal that survival() reads literally, and discounting it would silently restate the response model as well as the preference. weight_cost = 0.0, the default, leaves the composition bit-identical to one computed without a cost at all.

Every arm is filled independently, which is exact here because default_arm() makes them disjoint, so the objective separates.

Parameters:

• weight_evenness (float)

• weight_cost (float)

Return type:

Composition

exception mhcmatch.portfolio.MarginalExceedsBlock(n_over, n_total, worst_p, worst_q, arm=None)

Bases: ValueError

A unit’s marginal p exceeds its block’s live probability q, so p / q > 1 is not a probability and survival() cannot represent it.

A ValueError subclass, so except ValueError still catches it, but a named one carrying arm, n_over, n_total and the worst offending pair — because the bare raise it replaces told a donor’s operator only that something was 4-digit-too-large, and the actionable facts are which arm, how many of its units, and how far --block-live has to move.

arm is filled in by compose(), which is the only caller that knows it; the message is built from the attributes on each str() so setting it after the fact is enough.

mhcmatch.portfolio.VISIBILITY_COLUMNS: tuple = ('donor', 'pool_n', 'k', 'offset', 'field', 'vis_escape', 'log1m_vis', 'vis_escape_drv', 'vis_escape_pas', 'at_least_k', 'vis_at_least', 'vis_loh', 'loh_cost', 'lost_allotype', 'n_allotypes', 'foot_lam', 'lam_degenerate', 'foot_energy', 'foot_gap')

The ordered header mhcmatch visibility emits, exported for the same reason mhcmatch.cassette.SELECT_COLUMNS is: a workflow stub has to type this header and must not type it by hand. Two stubs had already drifted before those tuples existed.

field names which field foot_lam was read on and is always score, the raw aggregate log-odds. That is not a placeholder for a future option: the design’s lam is read on the objective field p - (gamma/2) p (1 - p) inside cassette score, and offering that reading here would put gamma back into an arm that has none. The column exists so a reader of the table never has to guess which of the two lam scales it is holding.