Presentation — the P of EPIC

Per-allele %rank, affinity and the calibration that makes them comparable.

mhcmatch.diffusion module

Anchor-factored presentation scoring with cross-allele kernel-shrinkage diffusion.

A per-allele anchor log-odds predictor – a small PWM over the anchor positions (MHC-I N-pocket + C-pocket, MHC1_ANCHORS) – whose per-allele anchor residue distributions are smoothed toward groove-similar alleles via mhcmatch.Pseudoseq. With raw=True (or bandwidth h -> 0) there is no borrowing and a rare allele scores off its own few peptides; with diffusion on, it borrows from frequent groove-neighbours, rescuing rare alleles. This is the forward per-allele E-value’s data-rescued null of the theory appendix §4.

mhcmatch.diffusion.load_markov1()

Order-1 human-proteome transition matrix {prev_residue: {residue: P(residue|prev)}} for background="markov" – a context-conditional presentation null. Vendored from UP000005640 (data/proteome_markov1.tsv). Opt-in and not the default: measured against the order-0 proteome null it is slightly worse on MHC-I rare-allele screening (AUPRC 0.820 vs 0.839, −0.019; AUROC −0.006; PPV −0.020 – compare_mhc1_human_random_{markov,proteome}bg.md) and neutral on medium/frequent. Kept for the adjacent-position covariance it injects, which may help elsewhere; it is not a win on the axis measured so far.

class mhcmatch.diffusion.AnchorModel(store, cls='mhc1', anchors=None, h=2.0, prior_strength=10.0, anticore=0.0, learn_weights=True, prune_dpi=False, weights='learned', register_em=2, footprint='anchor', rare_max=30, background='ligand', length_prior='score', length_motifs=True, register='marginal', n_motifs=3, families=None, reverse=0.0, pseudocount=0.0, pseudo_matrix='blosum62')

Bases: object

Per-allele anchor presentation model with optional cross-allele diffusion.

Built from a mhcmatch.Store. anchors are 1-based positions (default MHC-I MHC1_ANCHORS = N-pocket P1/P2/P3 + C-pocket PΩ-1/PΩ; MHC-II MHC2_ANCHORS = P1/P4/P6/P9). Per-anchor groove-position weights are learned by mutual information unless learn_weights is False; the kernel bandwidth h controls how much rare alleles borrow.

weights: "learned" (per-anchor MI over the panel, default) or "uniform". learn_weights=False forces uniform.

register_em (MHC-II only): number of GibbsCluster-style register EM passes. The anchor preferences are first estimated on the one-pass heuristic register; each pass then re-assigns every training peptide to the frame its own model scores best and re-estimates the preferences, so training and scoring use the same (best-frame) register. The default 2 lifts held-out binder-vs-decoy AUC across rare/medium/frequent MHC-II alleles (frequent +0.10); 0 keeps the one-pass heuristic register. Ignored for MHC-I (end-anchored). "converge" runs each allele to its own fixed point instead of a shared count – see _converge_registers(). No global pass count is right for every allele: HLA-DP is still improving at 32 passes while the rare stratum is done by 8, so 2 is an early stop that flatters rare rather than a correct value.

"converge-frequent" is "converge" plus an allele-frequency gate on the mixture, and the gate is there because that is where convergence actually reaches a thin allele. Measured on the 65-allele class-II shortlist: under "converge" a rare allele’s own register never moves (0 of 187 rare training frames, 0 of 225 raw count cells change against 2 – thin alleles reach their fixed point inside the first two passes and freeze), yet 567 of 675 of their mixture cells do, because _refit_mixture() re-searches every peptide’s per-component frame afterwards using the converged model. So the frequency gate holds alleles at or below rare_max out of the component refit; _dist()’s n_k = 0 backoff then returns the pooled single-PWM motif for them identically, with no new code path. Fitting n_motifs components on 30 ligands is the same overfit the adaptive footprint already guards against for MHC-I.

length_prior (MHC-I only) adds the per-allele ligand-length factor the anchor log-odds is structurally blind to – see length_logodds(). "score" (default) folds it into score(), so %rank and everything downstream inherit it; "post" only exposes length_logodds() for a caller that composes it itself; False is the earlier length-blind behaviour.

length_motifs (MHC-I only) estimates the residue distributions per peptide length instead of pooling every length into one counter – see _dist_len(). Complementary to length_prior: the prior is over L, the motifs are over residues given L.

register (MHC-II only) decides how the unobserved binding register enters score(): "marginal" (default) integrates it out under a learned core-offset prior; "max" is the earlier max-over-frames. See score().

n_motifs (MHC-II only) fits that many motif components per allele by EM and scores their mixture – see _refit_mixture(). 3 (default, human MHC-II) closes ~40% of the frequent-stratum AUPRC gap to NetMHCIIpan-4.3i; 1 is the single-PWM model (bit-identical to the pre-mixture code – it never enters the mixture path). Measured on human MHC-II only; thin alleles back off to the single PWM regardless of K.

families (MHC-II) gives each component its own gap placement: a list of (positions, k), meaning k components that score only positions. Positions outside a component’s family contribute log-odds 0 – that component asserts the allele is indistinguishable from background there, which is a statement about the likelihood, unlike a soft per-position weight (measured worse on every cell that moved: bench/results/anchor_position_weights.md). n_motifs becomes the sum of the k. self.anchors is untouched – families are subsets of it – so anchor_terms, score_sd, best_register and the learned groove weights are unaffected, and families=[(tuple(anchors), K)] is bit-identical to n_motifs=K.

The motivation is per-locus: MHC2_ANCHORS is DR’s pocket set applied to every locus (bench/results/anchor_position_weights.md), and the crystals put P7 above P4 in mouse H-2 I-A (bench/results/mhc2_gap_families.md). Pinning component k to family F_k also fixes the label-switching that stops components being kernel-shrunk across alleles – see _dist(). Ships off: the default is one implicit family over the whole footprint.

reverse_ctx = None

{context key: {residue: log-odds}} over the six intra-ligand context positions, learned by _fit_reverse_context() (reverse="auto+ctx"). Shifts the reverse prior per peptide on top of the per-allele one. None is off; class-level for the same unpickle reason.

anticore_w = 0.0

anticore = None

Class-level so an AnchorModel unpickled from a vendored model serialised before the anticore existed still answers “off” – __init__ does not run on unpickle, and _register_logprior reads this on every class-II score. Without it every shipped anchor_model_mhc2_*.pkl.gz raises AttributeError on first use.

reverse = 0.0

Prior mass on the reverse (C-to-N) reading of a class-II peptide; 0.0 is off and bit-identical. Class-level because score reads it and __init__ does not run on unpickle.

reverse_by_allele = None

{allele: p_a} learned by _fit_reverse() (reverse="auto"), overriding the scalar reverse per allele. None is “no per-allele prior”, which is every model built before it existed – class-level for the same unpickle reason as reverse itself.

length_logodds(length, allele, eps=0.001)

log P(L | allele) - log P_bg(L) – the ligand-length factor, in nats. MHC-I only; 0.0 when the model was built without length_prior.

The anchor log-odds is structurally length-blind: it sums a length-invariant number of per-position terms, so a 9-mer and a 10-mer with the same anchor residues score identically. But MHC-I length preference is strong and allele-specific (9-mer share ~0.32-0.96), and a screen tiles every length, so P_bg is uniform. The exact factorization

log P(pep|ligand,a)/P(pep|decoy) = [log P(L|a) - log P_bg(L)] + [log P(res|L,a)/P(res|L,decoy)]

is over two different variables, so this term adds to the anchor sum and cannot double-count it. Weight is fixed at 1 – it is a log-likelihood ratio, not a tunable feature.

P(L|a) is the panel’s per-allele length histogram, kernel-shrunk toward groove-similar alleles by the same bounded-prior estimator used for residues (Pseudoseq.shrink(), which is generic over the key type), so a rare allele borrows a length profile instead of trusting a handful of ligands. anchor=None gives uniform groove weights – correct here, since length preference is whole-groove (A/B/F pocket geometry), not a single pocket’s property.

register_entropy(allele, length=15)

Normalised entropy H/Hmax of the fitted core-offset prior, in [0, 1].

A pure function of the fitted model – no rival, no labels – and a usable health check: it tracks agreement with NetMHCIIpan’s own Core at Spearman -0.885 and the class-II AUROC gap at -0.703. Alleles below 0.85 average +0.0094 of AUROC against NetMHCIIpan and those at or above it -0.1208 (bench/results/mhc2_register_deficit.md). A near-uniform value means the register EM never locked on for this allele and its core placement should not be trusted, whatever the score says.

Returns 0.0 for MHC-I, which is end-anchored and has no register to be uncertain about.

panel_key(allele)

allele in this model’s key space – the spelling its panel was fit under.

The panel is keyed on the raw corpus string ('HLA-A*02:01', 'H-2Kb') while every lookup path resolves through normalize_allele() ('HLA-A02:01', 'H-2-Kb'), so one molecule had two key spaces and which one you got depended on how the caller typed the name. For human that is lossy – the kernel fallback over ~200 alleles is nearly right. For mouse it is fatal: 'H-2-Kb' and 'H2-Kb' are both keys in mhci_pseudo.fa on a byte-identical 34-mer, so resolve_allele called them exact while the panel held 35,037 ligands under 'H-2Kb' and none under either. SIINFEKL scored presentation %rank 0.0040 under the panel spelling and 20.19 under the resolver’s own.

Canonicalising to the panel’s spelling rather than the FASTA’s is deliberate: it is what every output string, result table and calibrator key already says, so this repairs the lookup without renaming 'HLA-A*02:01' in a user’s output. Unknown names pass through, so an allele the model has never seen still falls out downstream rather than raising.

best_register(peptide, allele, raw=False, eps=0.001)

Best-scoring binding register of peptide for allele, as (start, score).

For MHC-II every 9-mer core frame is scored and the winning one is returned (NNAlign/GibbsCluster-style, per allele) – start is its 0-based offset in peptide. MHC-I anchors are peptide-end-relative, so there is no register search and start is 0. Returns (-1, -inf) when the peptide is too short for the anchors. Ties are broken leftmost.

This is the register the model infers. It is not the allele-agnostic heuristic register (mhcmatch.store._mhc2_register) used for signatures, decompose and logos; on real ligands the two disagree often. Both are kept on purpose – see ROADMAP.

Parameters:

• peptide – the ligand (MHC-II) or peptide (MHC-I).

• allele – panel allele key.

• raw – score off the allele’s own anchor frequencies, without cross-allele borrowing.

• eps – log-odds regularizer.

Returns:

(start, score).

score(peptide, allele, raw=False, eps=0.001)

Anchor log-odds of peptide for allele vs the panel background.

raw=True uses the allele’s own anchor frequencies (no borrowing); the default diffuses over groove-similar alleles. Returns -inf if the peptide is too short for the anchors.

MHC-I anchors are peptide-end-relative, so there is no register search. For MHC-II the binding register is unobserved and register decides how it is handled:

"marginal" (default) integrates it out – log Σ_r P(r | L, allele) · exp(s_r) over frames r (see _offset_logprior()). The offset prior is real signal, not bookkeeping: a decoy’s best frame lands at a low-prior offset about as often as not, while a real ligand’s lands at the peaked one, and because the prior is normalized within a length the term still separates length-matched candidates.

"max" is the earlier behaviour, max_r s_r – a max over L-8 frames, which grows with peptide length even under the null (bench/results/binder_gate_length_bias.md).

With n_motifs > 1 the motif mixture wraps that marginal – log Σ_k π_k Σ_r P(r | L, allele) · exp(s_{k,r}) – one log Σ exp per latent, register inside, component outside (_refit_mixture()). The two compose because the background is common to every component and every frame, so it factors out of both sums.

Neither mode is comparable across peptide lengths. "marginal" normalizes the frame count away and roughly halves the inflation, but a Jensen residual remains (measured on random peptides, DRB1_1501, 9mer -> 21mer: +4.44 nats under "max", +2.28 under "marginal") – it saturates towards log E[e^s] rather than growing like ln n, but it is not zero. So an absolute binder call needs a length-matched %rank, not this score, and candidate ligand spans must be ranked by mhcmatch.ligand’s flank model – ranking them here would still just prefer the longest span.

score_many(peptides, allele, raw=False, eps=0.001, *, batch_bytes=8388608)

Score an iterable for one allele, in order, without retaining peptide intermediates.

Class-II canonical peptides use NumPy table lookups over bounded batches of core frames. Anchor additions and register/mixture reductions keep the scalar order, so results equal score() bit for bit. Markov, reverse and anticore models and noncanonical residues use that same scalar scorer. No worker pool or BLAS operation is used; callers may parallelize independent alleles with one kernel thread per worker.

batch_bytes bounds the estimated temporary array working set (8 MiB by default), excluding input/output and immutable model tables. A single peptide is the minimum unit of work. The returned list and fitted-model lookup tables remain caller-owned.

anchor_terms(peptide, allele, raw=False, eps=0.001)

Per-position log-odds components at the best register, one per self.anchors position (the full footprint, ignoring the rare-allele mask), or None if the peptide is too short.

Unlike score() (their sum) this exposes the vector, so a downstream regressor can weight positions differently – e.g. the affinity head (mhcmatch.affinity) learns pocket weights for binding energy rather than presentation specificity.

The width is always len(self.anchors): a signed-anchor collision (an 8-mer’s +5/-4) contributes 0.0 rather than dropping a column, so the vector stays a fixed-width feature.

score_sd(peptide, allele, raw=False, eps=0.001)

Posterior SD of score(), in nats. nan if the peptide is too short to score.

How much of a score is this allele’s own ligands and how much is borrowed. shrink with a prior_strength is a Dirichlet posterior – concentration alpha_r = own[r] + tau * nbr[r]/m, total alpha0 = n_own + tau – so the uncertainty is already determined and needs no sampling:

Var(theta_r) = theta_r (1 - theta_r) / (alpha0 + 1) Var(log theta_r) ~= Var(theta_r) / theta_r^2 = (1 - theta_r) / (theta_r (alpha0 + 1))

by the delta method, summed over the scored anchors – they are independent under the factorised model, and the background is a constant that drops out of a variance. Only the positions score() actually sums are counted, so a rare allele under the adaptive footprint reports the SD of the 5 anchors it was scored on, not of 9 it was not.

Two things it separates that a ligand count alone does not. Across alleles it tracks support (Spearman -0.945 against log n_ligands over 107 human class-I alleles, SD 0.032 nats at A*02:01’s 115,408 ligands to 6.07 at n=2). Within one allele, where alpha0 is fixed, all the variation is residue-driven – a peptide putting a rare residue on an anchor is uncertain even on a well-sampled allele – and it still predicts error: binned by SD quartile against length-matched decoys, the top quartile loses 0.09-0.18 AUROC against the best bin at A*02:01, B*07:02 and B*27:05.

This is the SD of the estimator, not of the biology, and not of discriminability. It says how well the panel pins this score down. It does not know that an allele’s ligands were all measured in one assay, and it cannot see model misspecification.

So do not select on it. Measured as a coverage curve on per-pMHC held-out ligands vs 19:1 length-matched decoys (bench/results/sd_coverage.md), keeping the lowest-SD fraction makes AUROC worse, monotonically: pooled 0.9921 at full coverage to 0.9847 at 25%, and the same direction in all three rarity strata. The mechanism is that SD is low wherever theta is high at the peptide’s anchor residues – which is equally true of a decoy carrying canonical anchor residues, i.e. of the hardest negatives. Low SD selects confidently-estimated hard cases, not easy ones.

What it is for is reporting: an allele-level and query-level statement of how much of this score is borrowed, so a rare-allele call can be shown as provisional. Never fit a weight on it – the moment it acquires a coefficient it stops being a posterior and becomes a hyperparameter.

For MHC-II this is evaluated at the best register rather than marginalised over registers, matching anchor_terms() and not score().

mhcmatch.diffusion.vendored_models(species='human')

The vendored-model registry for species. A store with no species declared reads human.

None is a mixed panel – Store.from_pmhc(species=None) loads both – and no vendored model can match one, so which registry it reads changes nothing except that the human names are the ones tried and missed.

Parameters:

species (str | None)

Return type:

dict

mhcmatch.diffusion.panel_sha(store, cls)

Content hash of the cls panel rows (epitope + allele, stored/build order). Cached on the store so the vendored-model guard is a one-off ~50 ms, not a per-call cost.

Return type:

str

class mhcmatch.diffusion.RoutedAnchorModel(frequent, rare, counts, rare_max)

Bases: object

Two fitted AnchorModel instances, one per allele-frequency stratum, dispatched per query.

bench/results/mhc2_register_frequency_gate.md §5 established that the class-II rare and frequent optima are incompatible inside one fitted model. Gating everything a rare allele borrows – its register, its mixture, its pooled null, its donor table and its tau – recovers about half the rare loss and stops there, and the reason is not a missing channel: tau="auto" at register_em=2 reaches rare screening AUPRC 0.689 while the same tau on a fully gated converge reaches 0.639. The rare optimum needs the donors themselves to be under-converged, which no gate on the borrower can supply – so the two have to be two fits.

The cut is counts[a] <= rare_max, which introduces no new number: 30 is already both the library’s rarity threshold (AnchorModel._rare_max, the adaptive-footprint switch) and the benchmark’s rare stratum boundary.

Only the four methods the library actually calls on an anchor model are dispatched; everything else – .ps, .cls, .anchors, .panel_key and friends – delegates to the frequent model, which is the one whose panel-wide quantities are shared.

ponytail: raw score values from two different fits are not comparable across alleles, only within one. Every shipped cross-allele path is already safe – RankCalibrator is per allele and restriction/vote ship calibrated=True – so this is a documented ceiling, not a bug to design around. If a caller ever needs raw cross-allele comparability, calibrate.

score(peptide, allele, raw=False, eps=0.001)

Dispatches to AnchorModel.score() on the frequent or rare fit, by allele’s count.

score_many(peptides, allele, raw=False, eps=0.001, *, batch_bytes=8388608)

Dispatches the whole batch to the same fit as score().

best_register(peptide, allele, raw=False, eps=0.001)

Dispatches to AnchorModel.best_register() on the frequent or rare fit.

score_sd(peptide, allele, raw=False, eps=0.001)

Dispatches to AnchorModel.score_sd() on the frequent or rare fit.

anchor_terms(peptide, allele, raw=False, eps=0.001)

Dispatches to AnchorModel.anchor_terms() on the frequent or rare fit.

mhcmatch.diffusion.load_vendored_anchor_model(store, cls, params)

The pre-fit AnchorModel for (cls, footprint, background) when one is shipped and the mhcmatch version, panel hash and full params all match; else None (caller builds).

Which file is tried comes from store.species; whether it is used still comes from panel_sha. Keeping those two separate is what makes a wrong species a cache miss rather than a wrong answer – and it is why the guard was already correct before mouse models existed, just slow.

mhcmatch.diffusion.save_vendored_anchor_model(store, cls, path, **kw)

Build the cls model (kw overrides, e.g. footprint= / background=; the rest are Store.anchor_model() defaults) and serialize it, gzipped, with a version / panel / params guard, to path. The release-time regenerator (mhcmatch build anchor).

mhcmatch.calibrate module

Per-allele score calibration: turn the allele-incomparable anchor log-odds into a cross-allele-comparable %rank (NetMHCpan %Rank_EL analogue) plus a calibrated presentation probability and a qualitative binding band.

The raw mhcmatch.AnchorModel.score() is a log-odds with a per-allele offset, so scores are not comparable across alleles. %rank fixes that: it is the percentile of a query score in the allele’s own random-peptide background (lower = stronger, exactly NetMHCpan’s definition), which is scale/offset-free and therefore comparable across alleles and directly usable as a binder threshold.

mhcmatch.calibrate.CACHE_ENV = 'MHCMATCH_CALIBRATION_CACHE'

Override for the on-disk per-allele calibration cache. Point it at a shared path to let a SLURM array or a Nextflow run reuse each other’s work; set it to "0", "off", "none" or "false" to disable caching entirely.

mhcmatch.calibrate.CACHE_DEFAULT = 'mhcmatch/calibration'

Where the cache lives when nothing overrides it – $XDG_CACHE_HOME if set, else ~/.cache, which is also what macOS users get and is fine there.

mhcmatch.calibrate.cache_dir()

The calibration cache directory, or None if caching is off. Created on first use.

On by default. A per-allele background is a random-peptide draw scored under one allele’s model: ~0.95 s to build, and a pure function of (allele, model, background, footprint, seed, library version) – every one of which is in the cache key, so a stale entry cannot be served across a refit or a version bump. Before this it was opt-in through CACHE_ENV and essentially nothing set it, so every process rebuilt every allele it touched on every run. On the neoantigen feature build, 2,093 distinct alleles in fourteen workers, that was the entire cost of the stage.

One file per (scorer fingerprint, allele), ~250 kB each – a 10,000-score background and an isotonic fit. A run that touches the whole 2,093-allele neoantigen panel through all three calibrators leaves on the order of 1 GB behind. Deleting the directory is always safe: it is derived data and rebuilds on demand.

Return type:

str | None

mhcmatch.calibrate.corpus_stats(peptides)

(aa_freq: Counter, length_dist: Counter) over an iterable of peptides.

mhcmatch.calibrate.random_peptides(aa, lens, n, rng, length_bg='corpus')

n random peptides with residue ~ aa frequency and length ~ lens distribution.

length_bg selects the length composition of the null:

• "corpus" (default): length ~ lens, i.e. the reference ligands’ own distribution (~9-mer heavy). Kept for MHC-II and for backwards compatibility.

• "uniform": equal numbers of each length in lens. This is what a screen actually sees – scan_protein/predict_windows tile every length, and a proteome yields ~n-L+1 windows per length (uniform to <1% for n >> L). It is also the convention of the %rank-style predictors mhcmatch is compared against. Use it for MHC-I, where the length preference is real biology that the score must be allowed to express against a length-neutral null.

Note "uniform" is not the same as a length-conditional (per-length) background: that would normalize each length to its own null and delete the length signal, which is wanted for the MHC-II register-max gate but is exactly wrong for MHC-I.

Parameters:

• aa (Counter)

• lens (Counter)

• n (int)

• length_bg (str)

class mhcmatch.calibrate.RankCalibrator(model, alleles, corpus, n=10000, seed=0, positives=None, length_bg='corpus', fingerprint=None)

Bases: object

Per-allele %rank (and optional calibrated P(present)) from a random-peptide background.

model is an mhcmatch.AnchorModel; alleles the panel to calibrate; corpus an iterable of reference peptides (for the background AA/length distribution). If positives (a {allele: [peptides]} map of known ligands) is given, a monotone isotonic P(present) is fit per allele from those positives vs the background. length_bg – see random_peptides(); "uniform" is the right null for MHC-I once the score carries a length prior.

Parameters:

• n (int)

• seed (int)

• length_bg (str)

• fingerprint (str | None)

clear()

Release loaded calibration distributions and isotonic fits, retaining the seeded null.

The next query recomputes or reloads them. This does not mutate the fitted scorer or delete disk entries. Call between independent work units, when no queries are active.

percent_rank(allele, score, length=None)

Percentile of score in the allele’s background: % of random peptides scoring higher (lower = stronger binder). nan if the allele has no background.

length conditions the null on that peptide length (_ensure_len()) instead of marginalising over the corpus length mix – required for any absolute threshold (a binder gate), since the raw score is length-inflated. Leave it None to rank peptides of a single length against each other, where the marginal null is what preserves MHC-I’s real length preference.

Above the background maximum the rank is extrapolated, not floored at zero. An empirical rank over n draws resolves 100/n, so with the default n=10,000 every peptide beating all 10,000 returned exactly 0.0 and they all tied. On the NCI exome scan that was 79 of 420,786 rows – and 6 of the 104 assayed-immunogenic candidates, a 307x enrichment, because the strongest binders are exactly where the answer is. In -log10 terms the emitted column could take values in [-2, 2] u {4} with literally nothing in between, and crossing that empty two-log-unit gap was worth +2.94 log-odds of EPIC score under the shipped standardizer for no information at all. That figure is coef * 2 / sigma and so moves with the fit — it was +2.05 under v9, and v11’s heavier binder makes the artefact worse, not better, which is the reason this path is not optional.

The tail of a sum of per-position terms is close to exponential, so a mean-excess fit to the top percentile of the background extends the rank smoothly past its last draw. It costs no extra peptide scoring – the background is already sorted – and it is monotone, so it can only break ties that were previously exact. It cannot move an AUROC; it moves a shortlist’s ordering and any threshold that sits inside the saturated block.

Parameters:

• allele (str)

• score (float)

• length (int | None)

Return type:

float

p_present(allele, score)

Isotonic-calibrated P(present | score) if positives were supplied, else a rank-derived fallback 1 - %rank/100.

Parameters:

• allele (str)

• score (float)

Return type:

float

mhcmatch.calibrate.band(percent_rank, strong=0.5, weak=2.0)

Qualitative binding band from %rank (NetMHCpan class-I thresholds): strong/weak/non-binder.

Parameters:

• percent_rank (float)

• strong (float)

• weak (float)

Return type:

str

mhcmatch.affinity module

Quantitative binding-affinity head: turn the presentation anchor log-odds into a calibrated IC50 (nM) and the neoantigen quantities that need it.

mhcmatch’s AnchorModel.score() is a presentation/specificity log-odds with a per-allele offset. Here we (1) center it against the allele’s own random-peptide background to make it cross-allele comparable, then (2) map that to the measured-affinity scale y = 1 - log(IC50)/log(50000) with a small ridge fitted offline on IEDB competition-binding IC50 (bench/affinity/train.py; coefficients vendored in data/affinity_<cls>.json). Predict back IC50 = 50000^(1-y) nM.

The headline use is the differential for neoantigen fitness – for a single-mutation WT/MT pair on the same allele the per-allele offset and systematic biases cancel, so the ratio is far more robust than either absolute nM:

• AffinityModel.amplitude() – Łuksza’s A = Kd_WT / Kd_MT with the 500 nM-cutoff correction (Łuksza et al. 2017 Nature, eq. 7/9), the amplitude of the neoantigen fitness model.

• AffinityModel.dai() – the differential agretopicity index (Duan 2014; Ghorani 2018), log10(Kd_WT / Kd_MT).

mhcmatch.affinity.ic50_to_y(nm)

Measured IC50 (nM) -> the NetMHC log50k regression target 1 - log(IC50)/log(50000) in [0,1].

Parameters:

nm (float)

Return type:

float

mhcmatch.affinity.y_to_ic50(y)

Inverse of ic50_to_y(): log50k score -> IC50 (nM), clamped to [0,1] first.

Parameters:

y (float)

Return type:

float

mhcmatch.affinity.fit_ridge(X, y, lam=1.0)

Closed-form ridge weights (XᵀX + λI)⁻¹ Xᵀy (numpy). Intercept column must be in X.

Parameters:

lam (float)

class mhcmatch.affinity.PottsAffinity(cls_name='mhc1', anchor_model=None)

Bases: object

Shipped affinity predictor: a Potts / direct-coupling energy model mapped to IC50 (nM).

The binding energy is E = Σ_i h_i(core_i) + Σ_j g_j(pocket_j) + Σ_{i,j} J_{ij}(core_i, pocket_j) – single-site fields on the 9-mer peptide core and the 34-mer MHC pseudosequence, plus pairwise couplings between every core and pocket position (the peptide×pocket interaction a purely additive model cannot represent). Weights are vendored (data/affinity_potts_<cls>.npz, fit on measured IEDB IC50 by bench/affinity/fit_potts.py), so prediction is a one-hot sparse dot product – numpy only, no sklearn at runtime.

MHC-I is end-anchored (core = the peptide, N5+C4). MHC-II’s 9-mer core is located by anchor_model.best_register (register-EM on presentation data). Same differential API as AffinityModel. Benchmark (per-allele held-out median Spearman ρ vs measured log-IC50): MHC-I common 0.70 / rare 0.49; MHC-II human 0.53 / mouse 0.51 (NetMHCpan/IIpan lead, but with IEDB train/test overlap). Build via mhcmatch.Store.affinity_model().

Parameters:

cls_name (str)

predict_y(peptide, allele)

log50k score (higher = stronger binder), or nan if the allele can’t be resolved.

Return type:

float

predict_y_batch(peptides, allele)

predict_y() over an iterable, as a numpy array – one gather instead of a loop.

Peptides the model cannot score (unresolvable allele, a core shorter than the footprint) come back nan in place, so the result lines up with the input and a caller never has to reconcile two differently-filtered lists.

predict_ic50(peptide, allele)

Predicted IC50 in nM (nan if the allele is unknown).

Return type:

float

amplitude(wt, mut, allele)

Łuksza amplitude A = Kd_WT/Kd_MT · 1/(1 + Kd_WT·ε/[L]) (eq. 9).

Return type:

float

dai(wt, mut, allele)

Differential agretopicity index log10(Kd_WT/Kd_MT) (>0 when the mutant binds better).

Return type:

float

class mhcmatch.affinity.AffinityModel(anchor_model, corpus, coef=None, n_bg=2000, seed=0)

Bases: object

Predict IC50 (nM) and neoantigen amplitude/DAI from an mhcmatch.AnchorModel.

anchor_model supplies the presentation log-odds; corpus an iterable of reference peptides for the per-allele random background (same idea as mhcmatch.calibrate.RankCalibrator). coef is the vendored fit {"b": [...], "lengths": [...]}; pass None to fit one with fit().

Parameters:

• n_bg (int)

• seed (int)

features(peptide, allele)

Feature row [1, <per-position z>..., <length one-hot>] or None if the peptide can’t be scored. Each z_i = the position-i log-odds centered by the allele’s background.

predict_y(peptide, allele)

The fitted linear model’s raw output – y_to_ic50() is what turns this into nM.

nan when features() cannot build a vector (peptide too short for the allele’s anchors, or no background z-scores yet computed).

Return type:

float

predict_ic50(peptide, allele)

Predicted IC50 in nM (nan if the peptide is too short for the allele’s anchors).

Return type:

float

amplitude(wt, mut, allele)

Łuksza amplitude A = Kd_WT/Kd_MT · 1/(1 + Kd_WT·ε/[L]) (eq. 9). A>1 when the mutation improves binding relative to self – the neoantigen-fitness amplitude.

Return type:

float

dai(wt, mut, allele)

Differential agretopicity index log10(Kd_WT/Kd_MT) (>0 when the mutant binds better).

Return type:

float

fit(rows, lam=1.0, lengths=(8, 9, 10, 11))

Fit the ridge on rows = iterable of (peptide, allele, ic50_nm); sets self.coef. Returns the number of usable training points.

Parameters:

lam (float)

classmethod load(anchor_model, corpus, cls_name='mhc1', path=None, **kw)

Load vendored coefficients data/affinity_<cls_name>.json (or an explicit path).

mhcmatch.predict module

Predict presented epitopes from a variant peptide-window FASTA.

Scores every binding-length k-mer of each window for a patient’s HLA alleles and emits two views. The input is one FASTA record per variant with the mutated window as the sequence and the variant annotation on the header – what pvacseq generate_protein_fasta writes after VEP, and what any producer can write, since parse_variant_header() reads a self-describing key=value header.

• native (write_native()) – one row per predicted binder with presentation %rank, P(present), band, IC50 (nM), the wild-type counterpart + agretopicity / amplitude / DAI, the synthesise / model peptides, and the anchor / TCR-facing decomposition. This is the output to read; everything mhcmatch computes is in it.

• scored-csv (write_scored_csv()) – the same calls in a fixed 57-column wide CSV (SCORED_COLUMNS). A compatibility export, kept so a caller whose downstream already reads that schema can swap its class-I/class-II binding predictor without touching anything else. It is not part of the generic path and it carries columns this package never fills.

mhcmatch scores per-allele presentation %rank / P(present) / band (mhcmatch.calibrate.RankCalibrator, the NetMHCpan %Rank_EL analogue) and quantitative IC50 (nM) via the Potts affinity head (mhcmatch.PottsAffinity). The export fills affinity (nM), affinity_percentile (%rank), and – for k-mers that span the somatic mutation – agretopicity (Kd_MT/Kd_WT vs the position-aligned wild-type peptide); expression / immunogenicity / composite-score columns are left to their own modules.

Alleles are used in whatever form the pipeline supplies (class I HLA-A*02:01; class II DRB1_1301 / HLA-DPA10103-DPB10401): built with Store.from_pmhc(), the panel keys match, and AnchorModel.score() normalizes internally for pseudosequence diffusion, so panel-absent alleles (e.g. HLA-B*15:07) are still scored zero-shot.

mhcmatch.predict.KMER_LENS = {'mhc1': (8, 9, 10, 11), 'mhc2': (12, 13, 14, 15, 16, 17, 18, 19, 20, 21)}

Binding-length k-mers tiled per class (pipeline params.mhcI_epit_len / mhcII_epit_len). An alias, not a copy – see mhcmatch.store.LIGAND_LENGTHS for why there is one ladder.

mhcmatch.predict.RANK_STRONG: dict = {'mhc1': 0.5, 'mhc2': 2.0}

%rank cut-offs, per class, in the NetMHCpan vocabulary. Strong and weak binder are not the same number in the two classes, and one threshold for both is the mistake this pair exists to stop: NetMHCpan calls class I strong at %rank <= 0.5 and weak at <= 2.0; NetMHCIIpan calls class II strong at <= 2.0 and weak at <= 10.0. A single 2.0 is therefore the weak cut for class I and the strong cut for class II.

They live here, not in mhcmatch.vector, because vector imports from this module and the cut belongs to whatever applies it. vector re-exports both names, so a caller that learned them there keeps working.

mhcmatch.predict.RANK_NONE: float = 100.0

every scored pair is emitted. A %rank is a percentile, so 100 is “no cut” exactly.

Type:

none

mhcmatch.predict.RANK_DEFAULT_TIER: str = 'none'

What predict/rank keep unless told otherwise. No cut. A tier here does not report, it filters – a dropped row is gone before ranking, and a caller cannot tell an empty table from a donor with nothing to offer. Measured on one class-II window pair against DRB1\*15:01: the old flat 2.0 default kept 0 of 56 scored pairs, discarding a best window at %rank 2.364 – an ordinary weak binder by the published convention, and the de novo arm returned an empty table with returncode 0. Choose the cut deliberately; wb is the conventional one.

mhcmatch.predict.RANK_TIERS: dict = {'sb': {'mhc1': 0.5, 'mhc2': 2.0}, 'strong': {'mhc1': 0.5, 'mhc2': 2.0}, 'wb': {'mhc1': 2.0, 'mhc2': 10.0}, 'weak': {'mhc1': 2.0, 'mhc2': 10.0}}

Spellings accepted by resolve_rank_threshold(), NetMHCpan’s on the left.

mhcmatch.predict.band_for(percent_rank, cls='mhc1')

strong/weak/non-binder at this class’s published cut-offs.

mhcmatch.calibrate.band() defaults to the class-I pair (0.5 / 2.0) and every call site here used to take the default, so a class-II ligand at %rank 5.0 – a textbook weak binder – came back labelled non-binder. The label is the thing a reader trusts without checking the number beside it, so getting it wrong is worse than not printing it.

Parameters:

• percent_rank (float)

• cls (str)

Return type:

str

mhcmatch.predict.KEEP_COLUMN: str = 'keep'

Name of the flag column a whitelist writes. 1 where a rule matched, 0 everywhere else.

mhcmatch.predict.KEEP_REASON_COLUMN: str = 'keep_reason'

Name of the column saying which rule matched. Two whitelists make two different claims about a row – “this gene is a driver” and “this peptide is one substitution from a validated immunogenic neoantigen” – and a single 1 cannot tell them apart. A reader who cannot see which rule fired will read a gene hit as evidence about the epitope, which it is not.

mhcmatch.predict.KEEP_BUILTIN: str = 'builtin'

The name that selects a shipped corpus instead of a caller-supplied list.

mhcmatch.predict.KEEP_INDEX_FILE: str = 'known_neoantigens.idx'

The shipped 1-mismatch index of validated immunogenic neoantigens, and its version sidecar. A seqtree.Index is opaque binary and cannot carry a version, so the stamp lives beside it.

mhcmatch.predict.KEEP_REASONS: tuple = ('epitope', 'epitope~1', 'gene')

Report order when more than one rule fires. An exact hit against a validated epitope is direct evidence about this peptide; a one-substitution hit is an inference from a neighbour; a gene symbol says nothing about the peptide at all. Same principle as mhcmatch.known.lookup().

mhcmatch.predict.keep_genes(spec)

Gene symbols that survive any %rank cut, from a comma list or a file.

Case is folded, so tp53 and TP53 are one entry. None/""/"none" is an empty set, which every caller reads as “no gene whitelist”.

"builtin" is not available: no driver-gene list ships yet. It raises rather than resolving to nothing, because a silent empty whitelist is indistinguishable from one that matched no row – the caller would see a table with every driver dropped and no way to tell why.

Return type:

frozenset

mhcmatch.predict.keep_epitope_index(spec, mismatch=0, quiet=False)

A seqtree.Index over the epitope whitelist, or None for no whitelist.

Three sources, and one match path for all three:

none / None

no epitope whitelist.

builtin

the shipped index of validated immunogenic neoantigens – every peptide that mhcmatch.known collects into its neoantigen set, i.e. an assay called it positive. Loaded from a pre-built file in ~1 ms and never rebuilt at run time: a thousand-sample Nextflow run would otherwise pay the build a thousand times and race on any cache it wrote.

a file or comma list

the caller’s own peptides, indexed here.

mismatch is the Hamming radius: 0 is exact, 1 also matches one substitution. Insertions and deletions are always off, so a hit is an equal-length peptide – a 9-mer query never matches a 20-mer known epitope by containment, which is a different question.

The index is C++ (seqtree), built once and queried concurrently; there is no Python dictionary anywhere on this path and no per-row rebuild.

Parameters:

• mismatch (int)

• quiet (bool)

class mhcmatch.predict.Keep(genes=None, epitopes=None, mismatch=0, quiet=False, *, threads=1)

Bases: object

Two independent whitelists – gene symbols and epitope sequences – and a Hamming radius.

Independent on purpose: a gene symbol keeps every candidate in a driver gene, and an epitope sequence keeps the ones with a validated response. They answer different questions, so one list matched against both fields (which is what --keep did) cannot say which claim a surviving row rests on.

Parameters:

• mismatch (int)

• quiet (bool)

threads

genes

mismatch

index

params

reasons(peptides, genes=(), *, threads=None)

Why each row is kept, "" where it is not – one batched C++ call for the table.

search_batch releases the GIL within the declared thread budget, so the epitope side costs one call for N rows rather than N calls. Returns a list as long as peptides.

Return type:

list

reason(peptide='', gene='')

Why this one row is kept, "" where it is not. Prefer reasons() for a table.

Parameters:

• peptide (str)

• gene (str)

Return type:

str

mhcmatch.predict.as_keep(spec)

Whatever a caller passed, as a Keep or None – built once, never per row.

A Keep passes through; None stays None; anything else is the deprecated flat --keep list and becomes one Keep with the same entries on both sides, which is exactly what that flag did with it.

mhcmatch.predict.keep_set(spec)

Deprecated – the flat --keep list, matched against gene and peptide alike.

Kept so a command line written against that flag still runs. Use Keep with separate genes= and epitopes=: one list matched both ways cannot report which claim kept a row, and it cannot do the 1-substitution epitope match at all.

Return type:

frozenset

mhcmatch.predict.is_kept(keep, peptide='', gene='')

Does this row match the whitelist? Accepts a Keep or a deprecated flat set.

Parameters:

• peptide (str)

• gene (str)

Return type:

bool

mhcmatch.predict.resolve_rank_threshold(spec, cls='mhc1')

A tier name or a bare percentage to the %rank cut for this class.

sb/strong and wb/weak resolve per class off RANK_STRONG and RANK_WEAK; none/all is RANK_NONE; anything numeric is taken as a percentage and used as given, so 25 means %rank <= 25 in either class.

The point of naming the tiers is that a number cannot be class-aware and a name can. A caller who writes 2.0 gets 2.0 in both classes, which is the weak cut in one and the strong cut in the other; a caller who writes wb gets 2.0 and 10.0 respectively, which is what they meant. Numbers stay honoured because “top 25%” is a real request that no tier expresses.

Parameters:

cls (str)

Return type:

float

mhcmatch.predict.SCORED_COLUMNS = ['type', 'subtype', 'chrom', 'pos', 'gene_name', 'gene_id', 'transcript_id', 'uniprot_id', 'tpm', 'ffpm', 'epitope', 'epitope_context', 'cluster_consensus', 'group', 'best_allele', 'agretopicity', 'affinity', 'affinity_percentile', 'CDR3', 'TCR-score', 'cellular_prevalence', 'rna_alts', 'rna_cov', 'ref_seq', 'seq', 'junction_reads', 'spanning_frags', 'isoform', 'orf_len', 'cov', 'fpkm', 'sv_len', 'cnv_score', 'paired_ref', 'paired_alt', 'single_ref', 'single_alt', 'ref', 'alt', 'd_signature', 'scaled_tpm', 'scaled_ffpm', 'score_expr_gene', 'score_expr_local_total', 'score_expr_local_ratio', 'score_expr_local', 'score_agretopicity', 'score_affinity', 'score_affinity_percentile', 'score_agretopicity_scaled', 'score_expr_gene_scaled', 'score_expr_local_scaled', 'score_affinity_percentile_scaled', 'score_signature', 'score', 'is_driver', 'driver_class']

The pipeline’s .epitopes.scored.csv header (57 columns, exact order). mhcmatch fills the variant-annotation and presentation columns; the rest are left empty for downstream modules.

mhcmatch.predict.NATIVE_COLUMNS = ('source', 'type', 'variant_type', 'gene_name', 'chrom', 'pos', 'ref', 'alt', 'peptide', 'offset', 'best_allele', 'cls', 'percent_rank', 'p_present', 'band', 'affinity_nm', 'affinity_rank', 'binder_rank', 'binder_band', 'wt_peptide', 'wt_affinity_nm', 'agretopicity', 'amplitude', 'dai', 'synth_peptide', 'model_peptide', 'anchors', 'tcr_facing', 'keep', 'keep_reason')

type is the header’s provenance (Somatic / Fusion / Isoform / CNV); variant_type beside it is the product class variant_product() derives, the same value rank emits under that name, so the two tables join on it. Appended to the group it belongs with rather than at the end, because this is mhcmatch’s own native format and not the fixed 57-column pipeline contract (SCORED_COLUMNS), which is unchanged.

mhcmatch.predict.CORE_COLUMNS = ('core', 'core_offset', 'core_source')

the 57-column SCORED_COLUMNS is a pipeline contract, and write_scored_csv’s extrasaction="ignore" would silently drop these rather than fail, so the opt-in is the flag and not the schema.

Type:

Appended by --core to every output that carries it. Never in a default header

class mhcmatch.predict.Prediction(source, peptide, allele, offset, cls, percent_rank, p_present, band, anchors, tcr_facing, affinity_nm=nan, wt_peptide='', wt_affinity_nm=nan, agretopicity=nan, amplitude=nan, dai=nan, affinity_rank=nan, binder_rank=nan, binder_band='', n_alleles_presenting=0, alleles_presenting='', keep=0, keep_reason='', core='', core_offset=-1, core_source='', synth_peptide='', model_peptide='', var=<factory>)

Bases: object

One predicted epitope: a window k-mer, its best-presenting allele, and its annotations.

Parameters:

• source (str)

• peptide (str)

• allele (str)

• offset (int)

• cls (str)

• percent_rank (float)

• p_present (float)

• band (str)

• anchors (tuple)

• tcr_facing (str)

• affinity_nm (float)

• wt_peptide (str)

• wt_affinity_nm (float)

• agretopicity (float)

• amplitude (float)

• dai (float)

• affinity_rank (float)

• binder_rank (float)

• binder_band (str)

• n_alleles_presenting (int)

• alleles_presenting (str)

• keep (int)

• keep_reason (str)

• core (str)

• core_offset (int)

• core_source (str)

• synth_peptide (str)

• model_peptide (str)

• var (dict)

source: str

peptide: str

allele: str

offset: int

cls: str

percent_rank: float

p_present: float

band: str

anchors: tuple

tcr_facing: str

affinity_nm: float = nan

wt_peptide: str = ''

wt_affinity_nm: float = nan

agretopicity: float = nan

amplitude: float = nan

dai: float = nan

affinity_rank: float = nan

binder_rank: float = nan

binder_band: str = ''

n_alleles_presenting: int = 0

How many of the queried allotypes present this peptide, and which, banded on the presentation %rank. The scoring loop already computes that rank for every allele and keeps only the best, so counting the rest is free; banding on binder_rank instead would mean an affinity call and a Fisher combine per allele, which is not. A peptide presented by three of a donor’s six class-I allotypes is a different bet from one presented by one: it spans three blocks of the response model in mhcmatch.portfolio.

alleles_presenting: str = ''

keep: int = 0

The 9-residue binding core, its 0-based offset, and which register produced it – the NetMHCpan core/Of pair. core_source is footprint for class I (both ends anchored, no register to choose), model when the class-II register came from mhcmatch.diffusion.AnchorModel.best_register(), and heuristic when it came from the allele-agnostic one-pass scan. The provenance is a column and not a docs sentence because the two registers disagree often on real ligands, and a core nobody can attribute is not evidence. 1 when a --keep whitelist named this peptide or its gene. Such a row is never dropped by a --rank-threshold, however strict, and says so in its own column rather than only by surviving – a reader cannot tell “kept because whitelisted” from “kept because it scored well” without one.

keep_reason: str = ''

epitope (exact hit against a validated immunogenic neoantigen), epitope~1 (one substitution from one), gene (the gene symbol is whitelisted), or "". A gene hit is not evidence about the peptide, so the two claims are reported apart rather than collapsed into the keep flag.

Type:

Which whitelist rule kept the row

core: str = ''

core_offset: int = -1

core_source: str = ''

synth_peptide: str = ''

model_peptide: str = ''

var: dict

mhcmatch.predict.parse_fasta(path)

[(header, sequence)] from a .peptide.fasta (header without the leading >).

Parameters:

path (str)

Return type:

list

mhcmatch.predict.NOVEL_PRODUCTS = frozenset({'frameshift', 'fusion', 'inframe_deletion', 'inframe_insertion', 'missense', 'protein_altering', 'start_lost', 'stop_lost'})

Product classes whose encoded sequence is absent from the normal proteome by construction – every somatic consequence _PRODUCT maps, plus fusion, whose novelty is a junction rather than a coding change and so has no consequence term to be mapped from.

Derived from _PRODUCT.values() rather than written out again, so a consequence added there cannot silently fail to be novel here. The union is {missense, frameshift, inframe_deletion, inframe_insertion, fusion, stop_lost, start_lost, protein_altering}.

What it is for. mhcmatch.vector.self_origin_risk()’s first clause – is the unit’s own gene transcribed in an essential tissue – is the MAGE-A12 hazard, and MAGE-A12 is a cancer-testis antigen: a shared, unmutated self protein whose transcription in brain is exactly what the construct would teach a T cell to attack. A product in this set is not that object. Its sequence is not in normal tissue, so its parent gene’s expression is not the hazard; what remains is the unrelated self-origin clause, which is tested separately and for every kind. An isoform, a wild-type or an overexpressed target is the MAGE-A12 case and is deliberately absent here.

mhcmatch.predict.TRACT_PRODUCTS = frozenset({'frameshift', 'fusion'})

a frameshift reads out of frame from its variant offset to the end of the product, and a fusion reads across a junction, so everything C-terminal of the offset is novel rather than one residue. Everything else in NOVEL_PRODUCTS alters a single position, and a window that does not contain it is wild-type sequence.

The distinction is what lets mhcmatch.vector.self_origin_risk() ask its second clause only of the registers that carry novel sequence. See that function for the measurement.

Type:

The subset of NOVEL_PRODUCTS whose novelty is a tract, not a position

mhcmatch.predict.parse_variant_header(header)

Parse a variant window header into variant-annotation fields.

The generic form is ``key=value``, separated by semicolons or whitespace and in any order:

>id=S1_1;gene_name=BRAF;chrom=chr7;pos=140753336;ref=A;alt=T;
 subtype=missense_variant;tpm=47.9;wt_window=...;mut_window=...

A producer emits only the fields it has, the reader needs no agreed field order, and adding a field cannot shift the meaning of the ones beside it – which is the one thing a positional schema cannot promise. Unknown keys are ignored, and _KV_ALIASES accepts the obvious synonyms (gene, consequence, chromosome, expression).

Four positional families are also read, because they are what the pipelines this module was first written against emit: Somatic:, Fusion: and CNV: are colon-delimited and Isoform: is pipe-delimited after its type. All five key into the same field names, so a consumer reads one dict. Best-effort and never raising – a field the header does not carry comes back empty.

The non-``Somatic`` families are parsed rather than skipped because they are the non-conventional neoepitopes, the ones a cassette holds a quota for, and dropping their gene and expression on the floor imputes the model’s largest coefficient on exactly the candidates that most need it.

Parameters:

header (str)

Return type:

dict

mhcmatch.predict.variant_product(var)

The product class of a parsed header – what kind of neoepitope this is.

Somatic is only the provenance of a variant; the product is its consequence, which lives in subtype. Returning the former is what let every candidate be charged to the non-conventional arm, since mhcmatch.portfolio.default_arm() asks only whether the kind is "missense".

missense / frameshift / inframe_deletion / … for a Somatic: window; the lower-cased type (fusion, isoform, cnv) otherwise, because a fusion is non-conventional whether its junction is in frame or not. Empty when the header says nothing – the caller’s own default then applies, rather than a guess made here.

Parameters:

var (dict)

Return type:

str

mhcmatch.predict.tile(seq, lengths)

[(kmer, offset)] for every standard-AA window of a length in lengths.

Parameters:

seq (str)

Return type:

list

mhcmatch.predict.SCORER_EPOCH = 9

Bumped by hand when the scoring code changes what a head returns, independently of any version bump or artifact rebuild. This exists because everything else in _fingerprint is data, and no hash over data can see a code change: the corpus length prior in PottsAffinity.predict_y and an extrapolated upper tail in calibrate.percent_rank both arrived within one released version, so a background cached before those and one cached after shared a key and the stale one was served. Same discipline as the EPIC model version – an int, moved deliberately.

1 = the original heads. 2 = length-aware Potts + extrapolated %rank tail. 3 = canonical allele keys (H-2Kb / H2-Kb / H-2-Kb collapsed to one molecule, so a cached background is no longer keyed on which spelling the caller typed). 4 = the background=”ligand” null leaves the queried allele out. 5 = the expression reference is keyed by species, so expr_lvl and expr_norm on a mouse row read FANTOM5 instead of missing GTEx and imputing to the training mean.

5 moves no human number and is bumped anyway, because this int is load-bearing in two repos: the benchmark’s feature frame keys its freshness guard on it, and that frame carries mouse rows whose expression column does change. A frame built under epoch 4 accepted under epoch 5 would fit mouse coefficients on a human-imputed column.

6 = the mouse expression halves are read off one harmonised matrix (toil_matrix_mmu.npz), so expr_lvl and expr_norm stop dividing by floors taken from two different assays – RNA-seq TPM at 0.9964 against CAGE tag density at 0.8000 – and every mouse value in both terms moves.

6 -> 7: rank._expression_for now ends its chain at the gene’s pan-tissue median instead of nan, so expr_lvl moves on every row that names a gene and no tissue – 485 of 968 mouse class-I rows and 289 of 522 class-II. A frame built under epoch 6 carries the imputed column.

mhcmatch.predict.build_scorer(store, cls, background='proteome', footprint='adaptive', seed=0, n_bg=10000)

(model, calibrator, affinity) for cls: an AnchorModel, a per-allele %rank calibrator, and the quantitative IC50 head (PottsAffinity), or None if unavailable.

background="proteome" puts the presentation score on the presentation axis (ligand-vs- proteome), matching NetMHCpan’s %Rank_EL; "ligand" measures allele-specificity instead.

Memoised on store: the result depends only on the panel, never on the query alleles, so scoring many samples against one store reuses a single build. The two costly MHC-II AnchorModel EM builds (this scorer + the affinity register oracle) are served from the vendored pre-fit models when the panel matches (see Store.anchor_model()), so the pipeline’s one-process-per-sample pattern pays no rebuild; RankCalibrator fills its per-allele background lazily.

class mhcmatch.predict.BinderScore(peptide, allele, cls, presentation_rank, affinity_nm, affinity_rank, binder_rank, band, p_binder=nan, presentation_sd=nan)

Bases: object

Generalized binder score for one (peptide, allele): a calibrated combined %rank that fuses presentation and affinity – a soft-AND scoring well only when the peptide is both presented and binds. It is the per-allele %rank of Fisher’s combined statistic -(ln p_pres + ln p_aff) against a random-peptide background, so binder_rank is itself a true %rank (lower = stronger, correctly banded) and is cross-allele comparable with no candidate pool. (Fisher’s statistic is monotone with the geometric mean of the two %ranks, so it induces the same ranking; calibration is what makes it a proper %rank and absorbs the presentation<->affinity correlation.)

Parameters:

• peptide (str)

• allele (str)

• cls (str)

• presentation_rank (float)

• affinity_nm (float)

• affinity_rank (float)

• binder_rank (float)

• band (str)

• p_binder (float)

• presentation_sd (float)

peptide: str

allele: str

cls: str

presentation_rank: float

affinity_nm: float

affinity_rank: float

binder_rank: float

band: str

p_binder: float = nan

presentation_sd: float = nan

Posterior SD of the presentation score, in nats (AnchorModel.score_sd()).

How much of presentation_rank is this allele’s own ligands and how much is borrowed from groove-similar neighbours. It is the SD of the estimator, not of the biology: it says how well the panel pins this score down, and it cannot see model misspecification or an allele whose ligands all came from one assay.

Report it; do not select on it. Across 107 human class-I alleles it runs 0.032 nats (A*02:01, 115,408 ligands) to 6.07 (n=2), Spearman -0.945 against log ligand count – so it says cleanly how provisional a rare-allele call is. But keeping the lowest-SD fraction makes AUROC worse (bench/results/sd_coverage.md), because low SD also picks out decoys with canonical anchor residues. See AnchorModel.score_sd().

mhcmatch.predict.binder_score(store, peptide, alleles='all', cls=None, background='proteome', footprint='adaptive', seed=0)

Rank alleles for peptide by the generalized binder score (presentation x affinity).

Motivation: the presentation head (AnchorModel %rank) and the affinity head (PottsAffinity) disagree along the binding-strength axis – presentation rescues weak-but-well-presented ligands, affinity rescues strong-but-atypical binders – so their geometric-mean %rank is a more robust binder index than either alone (measured: on the diverse NCI-423k neoantigen set the combined immunogenicity AUROC 0.965 beats presentation 0.945 and affinity 0.925; on affinity-labelled TESLA the affinity head alone is marginally better).

Returns list[BinderScore] sorted by binder_rank ascending (best first).

mhcmatch.predict.binder_ranks(store, peptides, allele, cls=None, background='proteome', footprint='adaptive', seed=0)

The transpose of binder_score(): one allele, many peptides, one call.

binder_score takes one peptide and ranks alleles for it. Scoring a benchmark is the other way round – a corpus of peptides against a known allele – so this is the natural call shape there, and it hoists the class inference, the build_scorer() entry and the three memo lookups out of the loop.

It is not a speed fix, and it was measured before being described as one. On a warm allele it runs 5,000 peptides at 82,241/s against 72,966/s for the per-peptide loop – 1.13x. The cost in a real feature build is the cold per-allele calibrator background, ~0.95 s the first time an allele is touched in a process; over the neoantigen corpus’s 2,093 distinct alleles, most of which carry a single peptide, that is the whole bill and no amount of batching the peptide loop reaches it. What would: persisting the per-allele background, which is a pure function of (allele, model, background, footprint, seed).

Returns four float arrays aligned with peptides – (presentation_rank, affinity_rank, binder_rank, ic50_nm) – with nan wherever binder_score would have dropped the row (an allele with no background, or a peptide the models cannot score).

Score-identical to binder_score by construction: same calibrators, same combine, same rounding. tests/test_predict.py pins that over the shipped panel.

mhcmatch.predict.predict_windows(store, cls, records, alleles, rank_threshold=None, top=None, background='proteome', footprint='adaptive', seed=0, keep=None)

Predict presented epitopes over records ([(header, sequence)]) for alleles.

For each window k-mer the best-presenting allele is chosen (lowest %rank); k-mers whose best %rank is above rank_threshold are dropped – and nothing is dropped by default, because rank_threshold=None resolves to RANK_DEFAULT_TIER. Pass "wb"/"sb" for the published class-aware cut, a number for an arbitrary percentile, or a keep whitelist of gene symbols and peptides that survive any cut and are flagged in the keep column. Each binder is annotated with its IC50 (nM), the wild-type counterpart’s IC50 + agretopicity / Luksza amplitude / DAI (when the k-mer spans the mutation), and the synthesise / model peptides. top optionally caps binders per window (strongest first). Returns list[Prediction].

mhcmatch.predict.predict_fasta(store, cls, fasta_path, alleles, **kw)

Convenience: parse_fasta() then predict_windows().

mhcmatch.predict.write_native(preds, path, core=False)

Write predictions as a native TSV (one row per predicted binder).

core=True appends CORE_COLUMNS; the default header is unchanged.

Parameters:

• path (str)

• core (bool)

Return type:

None

mhcmatch.predict.write_scored_csv(preds, path, core=False)

Write predictions in the pipeline’s 57-column .epitopes.scored.csv schema.

mhcmatch fills the variant-annotation columns (from the header) and the binding columns: best_allele, affinity (IC50 nM), affinity_percentile (%rank), and agretopicity (Kd_MT/Kd_WT for mutation-spanning k-mers). The expression / immunogenicity / composite-score columns are left empty for their own pipeline modules to populate.

core=True appends CORE_COLUMNS. The 57 are a contract with those modules, so the default stays byte-identical; widening it is the caller’s explicit choice.

Parameters:

• path (str)

• core (bool)

Return type:

None