API reference

mhcmatch.store module

MHC restriction & presentation from a reference epitope panel.

Productionizes the validated reverse-problem method (seqtree bench/bench_mhc_guess.py): index reference peptides by their anchored presentation signature (seqtree.layout.presentation_features()), widen the search scope around a query until it has enough neighbours, then rank presenting alleles by neighbour vote fraction and score confidence by a binomial-tail enrichment over the panel background. The vote fraction is the ranking statistic (robust to panel skew); the enrichment is the non-binder filter.

Significance theory: appendix/mhcmatch.tex §2-3 (forward per-allele E-value + reverse problem).

mhcmatch.store.fetch_pmhc(tier='full')

Download the pmhc presentation table for tier from the public HF dataset PMHC_REPO and return the local cached path.

Fetches only pmhc/pmhc_<tier>.tsv.gz (~4-12 MB) — never the other dataset directories — and relies on the huggingface_hub cache, so it downloads once and is instant thereafter. This lets a fresh install or a container bootstrap the reference panel with no pre-staged data, which the nextflow/Docker deploy depends on. Override with a local path=/$MHCMATCH_PMHC when present.

Parameters:

tier (str)

Return type:

str

mhcmatch.store.fetch_proteome(name='human')

Download a reference proteome FASTA from the public HF dataset PMHC_REPO (proteome/) and return the local cached path.

name is "human" / "mouse" — the full UniProt proteomes UP000005640 / UP000000589 (for source-protein lookup and peptide-flank extraction) — or a pathogen-proteome stem/filename bundled in the same dataset (e.g. "ecoli_K12_UP000000625", for molecular-mimicry sets). Cached by huggingface_hub, so it downloads once. Feeds mhcmatch.Proteome.from_hf().

Parameters:

name (str)

Return type:

str

mhcmatch.store.infer_class(peptide)

Heuristic class from length: MHC-I if <=11, else MHC-II. Pass cls to override.

Parameters:

peptide (str)

Return type:

str

class mhcmatch.store.Restriction(allele: 'str', vote: 'float', enrichment: 'float', n_votes: 'int', binder: 'bool', anchor_score: 'float | None' = None, rank: 'float | None' = None, p_present: 'float | None' = None, band: 'str | None' = None)

Bases: object

Parameters:

• allele (str)

• vote (float)

• enrichment (float)

• n_votes (int)

• binder (bool)

• anchor_score (float | None)

• rank (float | None)

• p_present (float | None)

• band (str | None)

allele: str

vote: float

enrichment: float

n_votes: int

binder: bool

anchor_score: float | None = None

rank: float | None = None

p_present: float | None = None

band: str | None = None

class mhcmatch.store.Decomposition(peptide: 'str', tcr_facing: 'str', presentation: 'str', anchors: 'tuple')

Bases: object

Parameters:

• peptide (str)

• tcr_facing (str)

• presentation (str)

• anchors (tuple)

peptide: str

tcr_facing: str

presentation: str

anchors: tuple

mhcmatch.store.anchor_indices(peptide, cls, register_start=None)

0-based anchor positions for a peptide: class-I P2/PΩ, class-II core P1/P4/P6/P9.

register_start (class II only) pins the 9-mer core to an explicit frame — e.g. the model’s mhcmatch.diffusion.AnchorModel.best_register(), so a caller that scored with the per-allele register can annotate with the same frame instead of the allele-agnostic heuristic. None keeps the one-pass heuristic register (the default everywhere else).

Parameters:

• peptide (str)

• cls (str)

• register_start (int | None)

Return type:

tuple

mhcmatch.store.resolve_anchor_index(peptide, cls, anchor)

0-based index of a scoring anchor in peptide (or None if out of range).

MHC-I: anchor is a 1-based peptide position (negatives count from the C-terminus). MHC-II: anchor is a 1-based position within the register-anchored 9-mer core (P1..P9).

Parameters:

• peptide (str)

• cls (str)

• anchor (int)

mhcmatch.store.mhc1_positions(length, anchors)

0-based peptide index for each signed MHC-I anchor, with collisions resolved.

Signed anchors collide on short peptides: mhcmatch.diffusion.MHC1_CORE’s +5 and -4 both resolve to index 4 of an 8-mer. Counting that residue twice makes the score an inflated, mis-normalized likelihood ratio (two perfectly-correlated terms), and files the same residue under two positions in Store.anchor_preferences(). Here the first anchor to claim an index keeps it; a losing anchor yields None and contributes nothing.

The return is aligned to ``anchors`` (same length), so callers keep their per-anchor bookkeeping. Returns None if any anchor falls outside the peptide (too short to score).

This is the single mapping shared by the scorer (mhcmatch.diffusion.AnchorModel.score()) and the preference estimator, so training and scoring cannot disagree about which residue sits where.

Parameters:

• length (int)

• anchors (tuple)

Return type:

tuple | None

class mhcmatch.store.Store

Bases: object

Searchable reference panel of presented peptides, partitioned by MHC class.

classmethod from_records(records, impute_alpha=False)

records: dicts with epitope, mhc_a (or mhc), mhc_class; optional weight (default 1.0) confidence-weights the peptide in anchor-preference estimation.

impute_alpha admits class-II records that type only the beta chain, by filling the most likely alpha from mhcmatch.pseudoseq.alpha_prior(); otherwise they are dropped (4,824 human records, 1.5% of the panel, 2,516 of them HLA-DPB1*11:01).

Default off, unlike the lookup path (class2_from_name(), where imputing turns a nan into an answer and is a strict win). Admitting these ligands to the reference panel was measured and it does not help: over the 13 alleles whose reference set grows, held-out AUROC moves -0.0019 and AUPRC -0.0012, and the damage scales with the merge – HLA-DPA10201-DPB11101 gains 2,339 ligands (+89%) and loses 0.0155 AUROC. A study that skipped alpha-typing produced noisier ligand calls too, so the missing alpha is a marker of data quality and not merely of absent metadata. Turn it on only if you want coverage of those ligands more than motif purity.

Parameters:

impute_alpha (bool)

classmethod from_pmhc(path=None, tier='full', species=None, classes=('mhc1', 'mhc2'), impute_alpha=False)

Load the isalgo/pmhc_data TSV(.gz). species filters the MHC species ("human" / "mouse"). If path is None it uses $MHCMATCH_PMHC/pmhc_<tier>.tsv.gz when that env var is set, otherwise bootstraps the table from the public HF dataset via fetch_pmhc() (downloads only pmhc/pmhc_<tier>.tsv.gz, cached) — so a fresh install or a container needs no pre-staged data.

Parameters:

impute_alpha (bool)

alleles(cls)

restriction(peptide, cls=None, alleles='all', top=10, alpha=0.05, diffuse=False, calibrated=False)

Rank presenting alleles for peptide (vote fraction), flag binders (enrichment).

alleles: "all", a single allele, or a list. alpha: per-allele significance for the non-binder flag (binder iff binomial-tail p <= alpha and the allele got votes).

calibrated=True (implies diffuse) additionally fills each result’s rank (per-allele %rank vs a random-peptide background, lower = stronger), p_present, and qualitative band (strong/weak/non-binder). The %rank is the cross-allele-comparable score; it also re-ranks the results (ascending %rank).

This rank is on the allele-specificity axis (the model here is background="ligand"): it asks how strongly this allele, versus other alleles, prefers the peptide – so it can band a canonical, widely-shared ligand as “weak” even when the allele is unambiguously the correct restriction. For the presentation axis – is this presented at all, the NetMHCpan %Rank_EL question – score with a proteome null via mhcmatch.predict.predict_windows() / mhcmatch predict (background="proteome").

With diffuse=True the diffusion-shrunk anchor log-odds (mhcmatch.diffusion.AnchorModel) ranks and the neighbour vote/enrichment gates: an allele is a binder if it is vote-significant or the anchors are plausible. On held-out (novel) peptides the anchor log-odds is the far better ranker—the vote method relies on same-allele signature neighbours, which are sparse for a genuinely new peptide, so vote-first ranking buries the true allele; the diffused anchor score scores every allele directly and rescues rare ones. Vote breaks ties. Without diffusion, vote fraction ranks and the call returns [] when there are no neighbours.

“Anchors are plausible” is class-specific, and the difference is load-bearing:

• MHC-II: %rank <= 2 against random peptides of the query’s own length. score is a max over the L-8 register frames, so it climbs with length even on pure noise – the old absolute anchor_score > 0 gate was a length detector (it passed a random 15-mer 85% of the time, a random 21-mer 98%). Scoring the null at the same length puts it through the same frame-max, so the bias cancels. This costs a per-(allele, length) calibration.

• MHC-I: still anchor_score > 0. It is end-anchored – no register search, no max, no length inflation to correct – and its length preference is real modelled biology that a length-conditional null would delete. MHC-I results are unchanged and pay no calibration.

is_binder(peptide, allele, cls=None, alpha=0.05)

is_presented(peptide, cls=None, alpha=0.05)

Overall presentation: does any panel allele present this peptide?

scan_protein(protein, cls='mhc1', alleles='all', lengths=None, alpha=0.05, top=3, correction=None)

Slide all binding-length windows over protein and return presented peptides.

Returns [(position, peptide, [Restriction, ...]), ...] for windows with >=1 binder.

correction controls multiple testing over the (window, allele) presentation tests in the scan (appendix §5): None (default) keeps the per-window per-allele alpha; "bonferroni" controls the family-wise error rate (threshold alpha/m); "bh" controls the Benjamini-Hochberg false-discovery rate. m is the number of voted (window, allele) tests. The vote tail p-value is 10**(-enrichment); corrected calls replace the per-window binder flag.

decompose(peptide, cls=None, allele=None, register_start=None)

Split peptide into anchor and TCR-facing parts, each masked with X.

tcr_facing: anchors -> X (the recognition readout). presentation: TCR-facing -> X (the anchor readout). allele is accepted for forward-compat (allele-specific learned anchors, Phase 1); v0 uses class-default anchor positions.

register_start (class II) pins the 9-mer core frame — pass the model register a caller already scored with (AnchorModel.best_register) so the reported anchors match the scored core; None keeps the allele-agnostic heuristic register (the two systems stay separate, ROADMAP §7).

anchor_model(cls='mhc1', h=2.0, prior_strength=10.0, anchors=None, 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, pseudocount=0.0, _vendored=True, _return_params=False)

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

See mhcmatch.diffusion.AnchorModel. The diffusion rescues rare alleles by borrowing anchor preferences from groove-similar frequent ones, with a bounded prior strength so a large neighbour cannot swamp a rare allele’s own peptides. register_em (MHC-II) runs that many best-frame register-EM passes so training and scoring share the same register. footprint="anchor" (default) scores the primary pockets only; "core" scores the whole binding core (MHC-I P1-P5 + PΩ-3..PΩ, MHC-II 9-mer core) – more discriminative when non-anchor positions carry allele-specific signal. background="ligand" (default) is the allele-specificity null; "proteome" is the presentation null (better for ligand-vs-random screening) – see mhcmatch.diffusion.PROTEOME_AA_FREQ. length_prior="score" (MHC-I) adds the per-allele ligand-length factor the anchor log-odds is blind to – see mhcmatch.diffusion.AnchorModel.length_logodds(). register="marginal" (MHC-II default) integrates the unobserved binding register out under a learned core-offset prior; "max" restores the pre-v0.6 max-over-frames – see mhcmatch.diffusion.AnchorModel.score(). n_motifs (MHC-II) fits that many motif components per allele and scores their mixture; 3 (default) closes ~40% of the frequent-stratum gap to NetMHCIIpan, 1 is the single-PWM model – see mhcmatch.diffusion.AnchorModel._refit_mixture(). pseudocount (β) spreads each anchor’s observed counts onto chemically similar residues with weight β/(n+β); 0 (default) is off – see mhcmatch.diffusion.AnchorModel._add_pseudocounts().

affinity_model(cls='mhc1')

Quantitative IC50 (nM) + neoantigen amplitude/DAI head (mhcmatch.PottsAffinity).

Loads the vendored Potts weights data/affinity_potts_<cls>.npz (fields + peptide×pocket couplings, fit on measured IEDB IC50). For MHC-II it also builds the register oracle (an AnchorModel with the same proteome/core config used at fit time) so the 9-mer core is located consistently. Cached per class. Predict with .predict_ic50(peptide, allele) and the differential .amplitude(wt, mut, allele) / .dai(wt, mut, allele).

binder_score(peptide, alleles='all', cls=None, **kw)

Rank alleles for peptide by the generalized binder score – the geometric mean of the presentation (AnchorModel %rank) and affinity (PottsAffinity %rank), a soft-AND that scores well only when the peptide is both presented and binds. See mhcmatch.predict.binder_score(). Returns list[BinderScore] best-first.

anchor_preferences(cls, anchor, anchors=None, by_length=False)

{allele: Counter(residue)} at a 1-based anchor position (negative from C-term).

anchors (MHC-I): the full footprint. When given, signed-anchor collisions on short peptides are resolved with mhc1_positions() – the same rule the scorer uses – so a residue is filed under exactly one position. Without it an 8-mer’s index-4 residue lands in both +5 and -4, and training would disagree with scoring.

by_length=True returns {peptide_length: {allele: Counter(residue)}} instead. The pooled (default) form mixes every length into one counter, so the motif it yields is really the 9-mer motif (~2/3 of the panel) applied to 8/10/11-mers too – measurably wrong off-9. Splitting by length is what the estimator in mhcmatch.diffusion.AnchorModel._dist_len() backs off from, since per-(allele, length) counts are thin (rare alleles have a median of zero 8-mers).

length_preferences(cls)

{allele: Counter(peptide_length)} over the panel – the per-allele ligand-length distribution, publication-weighted like anchor_preferences().

MHC-I alleles differ strongly here (9-mer share ranges ~0.32-0.96; HLA-B*52:01 is ~65% 8-mers), and the anchor log-odds is blind to it: its term count is length-invariant, so a 9-mer and a 10-mer with the same anchor residues score identically. This feeds mhcmatch.diffusion.AnchorModel._length_logodds(), which restores the missing factor.

logo.motif computes a per-allele length histogram too, but unshrunk and for display only.

mhcmatch.search module

Large-scale peptide similarity search over big peptide sets / proteomes.

Two notions of “similar”, both via the seqtree C++ KmerIndex seed-and-gather:

• mode="tcr" – anchor-masked TCR-facing homology: similar T-cell recognition profile (the basis for cross-reactivity / molecular mimicry).

• mode="mhc" – anchored presentation signature: likely presented by the same MHC.

For neoantigen mimicry with per-allele presentation-aware E-values, use find_mimics() (re-exported from seqtree). See appendix/mhcmatch.tex §5.

class mhcmatch.search.Match(peptide: 'str', shared_kmers: 'int', score: 'int')

Bases: object

Parameters:

• peptide (str)

• shared_kmers (int)

• score (int)

peptide: str

shared_kmers: int

score: int

mhcmatch.search.search(query, peptides, mode='tcr', cls='mhc1', k=4, max_subs=1, min_shared=1, exclude_self=True, threads=0)

Peptides in peptides similar to query under mode ("tcr" or "mhc").

mhcmatch.proteome module

Near-exact source-peptide lookup against a reference proteome.

Given a query peptide (e.g. a neoantigen), find the nearly-exact self peptide it derives from and its parent protein / position via full-sequence (unmasked) <= max_subs search over all windows of the proteome of the query’s length – using the seqtree Hamming fast path. This is a distinct mode from the anchor-masked TCR-facing homology and the presentation-signature searches. See appendix/mhcmatch.tex §5 (near-exact source identification).

mhcmatch.proteome.read_fasta(path)

{name: sequence} from a (optionally gzipped) FASTA; name = first whitespace token.

class mhcmatch.proteome.SourceHit(protein: 'str', position: 'int', ref_peptide: 'str', n_subs: 'int', mutations: 'tuple')

Bases: object

Parameters:

• protein (str)

• position (int)

• ref_peptide (str)

• n_subs (int)

• mutations (tuple)

protein: str

position: int

ref_peptide: str

n_subs: int

mutations: tuple

class mhcmatch.proteome.Proteome(seqs)

Bases: object

A reference proteome with lazily-built per-length window indices.

classmethod from_fasta(path)

classmethod from_hf(name='human')

Load a reference proteome by name, auto-fetched from the public HF dataset (no manual download). name = "human" / "mouse" (UP000005640 / UP000000589) or a pathogen stem; see mhcmatch.store.fetch_proteome().

find_source(peptide, max_subs=1, exclude_exact=False)

Self peptides within max_subs substitutions of peptide, nearest first.

Returns [SourceHit, ...]. exclude_exact=True drops perfect (0-mismatch) matches – useful to find the wild-type a mutated neoantigen derives from when the query is itself self.

wildtype(peptide, max_subs=1)

The wild-type self peptide a mutated peptide derives from, or None.

A self peptide exactly one substitution away (its point-mutation origin) – the position-aligned WT counterpart needed for agretopicity / DAI when the caller has no WT window (e.g. a bare neoantigen list like TESLA). None when nothing is one sub away (indel / spliced / non-self, or the peptide is itself an exact self peptide with no mutated origin). Ties resolve to the first variant found (position, then residue order).

For max_subs=1 this uses a hash-set fast path (generate the L*19 single-sub variants and test proteome membership – microseconds/peptide, so it scales to large corpora); larger max_subs falls back to the general find_source() fuzzy search.

mhcmatch.pseudoseq module

MHC pseudosequence allele-similarity & cross-allele diffusion.

Each allele is a 34-residue groove pseudosequence (NetMHCpan-style; vendored in data/{mhci,mhcii}_pseudo.fa). Allele similarity is an anchor-factored kernel over these positions: K_j(a,b) = exp(-d_j(a,b)/h) where d_j is a position-weighted Hamming distance and the per-anchor weights w_j say which groove residues govern peptide anchor j (e.g. MHC-I P2 vs PΩ). learn_anchor_weights() learns w_j from data (mutual information between a groove position and the allele’s anchor-residue choice) – the “feature importance” of each pocket.

Kernel-weighted shrinkage (Pseudoseq.shrink()) borrows presented-peptide statistics from similar alleles to rescue rare ones, lifting the seqtree limitation “distinct alleles are distinct nulls”. See appendix/mhcmatch.tex §4.

mhcmatch.pseudoseq.normalize_allele(a)

pmhc allele name -> pseudosequence-FASTA key.

Drops the * ('HLA-A*02:01' -> 'HLA-A02:01') and repairs the mouse H-2 dash (pmhc 'H-2Kb' -> FASTA 'H-2-Kb').

Parameters:

a (str)

Return type:

str

mhcmatch.pseudoseq.alpha_prior()

DP/DQ beta chain -> most likely alpha chain, for typings that omit the alpha.

Learned from the IEDB-derived panel and vendored (data/mhc2_alpha_prior.tsv); a beta is listed only when its 34-mer groove is >=95% determined over >=50 fully-typed ligands. See class2_key().

Return type:

dict

mhcmatch.pseudoseq.class2_key(mhc_a, mhc_b='', impute_alpha=True)

pmhc class-II allele -> pseudosequence-FASTA key (locus-aware).

DR (the DRA chain is monomorphic) is keyed by the beta chain alone, e.g. 'HLA-DRB1*01:01' -> 'DRB1_0101'. DP/DQ are keyed by the alpha-beta pair, e.g. ('HLA-DPA1*01:03', 'HLA-DPB1*04:01') -> 'HLA-DPA10103-DPB10401'. With no beta chain the input is returned unchanged (mouse H-2 and fallbacks).

impute_alpha (default on) fills a missing DP/DQ alpha from alpha_prior(), so a beta-only typing resolves to a real groove instead of the unscorable '-DPB11101'. This is the polymorphic-locus analogue of what DR already gets for free from monomorphic DRA. It fires only where the panel pins the groove to >=95% over >=50 ligands – DQA1’s polymorphism sits in the alpha1 domain the pseudosequence samples, so a name- or 2-digit-group-level rule is not a substitute: DQA1*01:02 and DQA1*01:05 share the group DQA1*01 but not the 34-mer, which reads as 100% certain while the sequence is a 58/42 coin flip. Rare DQ betas are left unresolved on purpose – a wrong groove scores silently, which is worse than not scoring.

Parameters:

• mhc_a (str)

• mhc_b (str)

• impute_alpha (bool)

Return type:

str

mhcmatch.pseudoseq.class2_from_name(name, impute_alpha=True)

Class-II allele name (user- or IEDB-typed) -> mhc2 pseudoseq key, locus-aware.

Handles DR (beta-only 'HLA-DRB1*15:01' -> 'DRB1_1501'), the DP/DQ alpha-beta pair given as 'HLA-DQA1*05:01/DQB1*03:01', a DP/DQ beta given alone ('HLA-DPB1*11:01' -> 'HLA-DPA10201-DPB11101', the alpha imputed via alpha_prior() – see class2_key()), and mouse ('H2-IAb' / 'I-Ab' -> 'H-2-IAb'). Falls back to normalize_allele() for anything already in key form.

Parameters:

• name (str)

• impute_alpha (bool)

Return type:

str

mhcmatch.pseudoseq.resolve_allele(name, cls)

Resolve a user-typed allele name to a pseudosequence key for cls.

Returns (key, exact). exact=True when name (after normalize_allele(), or the locus-aware class2_from_name() for cls=="mhc2") is a known key; otherwise the closest key by name—a missing HLA- prefix is repaired and a too-short (e.g. two-field 'HLA-A02:01') name is completed by prefix to its first matching key—with exact=False; (None, False) if nothing matches. Serotype names ('HLA-A2') are not expanded. Lets callers accept messy input ('A*02:01', 'HLA-A0201') and report when a requested allele is unknown rather than silently dropping it.

Parameters:

• name (str)

• cls (str)

mhcmatch.pseudoseq.load_pseudo(cls)

allele-id -> 34-mer for the bundled pseudosequence FASTA of a class.

Alleles sharing a 34-mer are collapsed to one FASTA record whose header lists every such allele (>A B C|n=3), so all of them are keys here. Listing only the first would silently make the rest unscorable – they are not rare variants: 8,854 of the source table’s 12,997 alleles (68%) are non-representatives, among them HLA-B*14:02, B*18:05 and C*03:04.

Parameters:

cls (str)

Return type:

dict

mhcmatch.pseudoseq.BLOSUM62_BG = {'A': 0.0742, 'C': 0.0247, 'D': 0.0536, 'E': 0.0543, 'F': 0.0474, 'G': 0.0741, 'H': 0.0262, 'I': 0.0679, 'K': 0.0582, 'L': 0.0989, 'M': 0.025, 'N': 0.0446, 'P': 0.0385, 'Q': 0.0343, 'R': 0.0516, 'S': 0.0572, 'T': 0.0509, 'V': 0.0729, 'W': 0.013, 'Y': 0.0323}

BLOSUM62’s own background – the Blocks pair marginals p(i,*) of Henikoff & Henikoff’s blosum62.qij (PMID 8743679). The matrix’s lambda and this background are jointly determined: s_ab = nint(2·log2(q_ab / (p_a·p_b))) holds only with these frequencies. Deliberately not mhcmatch.diffusion.PROTEOME_AA_FREQ, which answers a different question (the scoring null).

mhcmatch.pseudoseq.blosum62_conditional()

{observed: {r: P(r | observed)}} – the BLOSUM62 substitution conditional.

The q(a|b) of Nielsen et al. 2004 (PMID 14962912), used to spread an anchor’s observed residue counts onto chemically similar residues (see mhcmatch.diffusion.AnchorModel._add_pseudocounts()).

No q_ij table and no new dependency are needed. BLOSUM half-bits are s_ab = 2·log2(q_ab / (p_a·p_b)), so q_ab = p_a·p_b·2^(s_ab/2) and

P(a|b) = q_ab / p_b = p_a · 2^(s_ab/2) (normalized over a)

– only the 20 background frequencies survive. Reads .similarity() (the raw signed half-bits); .penalty() is the Gram form s_aa + s_bb - 2·s_ab, which forces the diagonal to zero and so cannot recover the log-odds.

Return type:

dict

mhcmatch.pseudoseq.mutual_information(xs, ys)

MI(X;Y) in bits for two aligned categorical sequences.

Return type:

float

mhcmatch.pseudoseq.learn_anchor_weights(pseudo_seqs, anchor_residue, prune_dpi=False, tol=0.0)

Per-position relevance w[p] = MI(groove position p residue ; anchor residue) across alleles, normalized to mean 1. anchor_residue: {allele: residue} (e.g. the modal residue at one peptide anchor for that allele). Positions that discriminate the anchor get more weight.

Raw MI is inflated by linkage between groove positions (they co-vary across alleles), so many positions look relevant and the per-pocket profile is smeared. With prune_dpi=True an ARACNE data-processing-inequality prune removes indirect links: position p’s edge to the pocket is dropped if some other position q is more informative about the pocket and about p (I(p;pocket) <= min(I(q;pocket), I(p;q))), leaving the direct pocket positions sparse and distinct.

Parameters:

• pseudo_seqs (dict)

• anchor_residue (dict)

• prune_dpi (bool)

• tol (float)

Return type:

list

class mhcmatch.pseudoseq.Pseudoseq(cls, h=2.0, weights=None, metric='blosum')

Bases: object

Allele-similarity kernel and diffusion over groove pseudosequences for one MHC class.

h: kernel bandwidth. weights: per-position list (one kernel) or {anchor: [34 weights]} (anchor-factored, from learn_anchor_weights()). metric: "blosum" (default) scores each position by the BLOSUM62 Gram distance (conservative substitutions cost less); "identity" counts plain mismatches.

kernel(a, b, anchor=None)

Return type:

float

neighbors(allele, candidates=None, anchor=None, top=10, min_k=0.0)

[(allele, kernel), ...] most groove-similar to allele (self excluded).

cluster(alleles, anchor=None, threshold=0.5)

Single-linkage clusters: merge alleles with kernel >= threshold. O(n^2); use on a panel (~hundreds of alleles), not the full 4k-allele set.

shrink(prefs, allele, anchor=None, candidates=None, prior_strength=None)

Kernel-weighted empirical-Bayes pooling of a per-anchor residue distribution.

prefs: {allele: Counter(residue -> count)} for one anchor. Returns the shrunk probability dict for allele.

With prior_strength=None (default) this is the counts-weighted form (n_a π_a + Σ_b K_ab n_b π_b) / (n_a + Σ_b K_ab n_b) with limits h -> 0 (raw per-allele) and h -> ∞ (global pool). With prior_strength=τ it uses the fixed-concentration form (n_a π_a + τ m_a) / (n_a + τ) where m_a is the kernel-weighted neighbour mean – a bounded prior that prevents one large neighbour from swamping a rare allele’s own peptides and self-adapts to n_a (appendix §4, Prop. on bias–variance). The latter is the recommended default for the forward scorer.

Return type:

dict

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 appendix/mhcmatch.tex §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, 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, pseudocount=0.0)

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.

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 length-blind v0.4 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 pre-v0.6 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.

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.

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 pre-v0.6 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.

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.

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

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

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 (tools/build_anchor_models.py).

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

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)

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.

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_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)

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

Structure-based pMHC Miyazawa–Jernigan contact energy and WT/MT ΔΔG (optional tcren dep).

Threads a peptide onto a template pMHC crystal groove and sums the MJ contact potential (all shipped by tcren). For a query allele with no own template we borrow the groove-closest template (the same pseudosequence kernel the diffusion uses). For a single-mutation WT/MT pair the ΔΔG = MJ(mut) − MJ(wt) on one backbone is a physics-based differential affinity estimator (adaptive double threading, Jojic et al. 2006) – the structural complement to the sequence-based mhcmatch.affinity.

On measured HLA-A*02:01 the MJ energy tracks log-IC50 at Spearman ≈ 0.55 (see bench/affinity/bench_structure.py), at ~0.02 ms/peptide after a one-time template build.

Needs the [structure] extra:

pip install 'mhcmatch[structure]'      # pulls tcren

Template structures are not vendored (they live in tcren’s Canonical2026 set, which the tcren wheel deliberately does not ship). Resolution order for the templates: the structure_dir argument, then $MHCMATCH_STRUCTURES, then tcren’s own data dir ($TCREN_DATA_DIR, or an editable tcren checkout’s data/) under Canonical2026.

class mhcmatch.structure.StructureScorer(structure_dir=None, templates=None, pseudoseq=None, cutoff=5.0)

Bases: object

MJ contact-energy scorer over template pMHC structures.

pseudoseq (a mhcmatch.Pseudoseq) enables borrowing the groove-closest template for alleles without their own; omit it to restrict to exact-allele templates.

template_for(allele, length)

(pdb, chains) of the best template for allele at peptide length: exact allele if present, else the groove-closest templated allele (needs pseudoseq). None if none.

mj_energies(peptides, allele)

{peptide: MJ energy} for equal-length peptides on allele’s template (lower = stronger binding). Empty dict if no matching-length template. One batch swap.

mj_energy(peptide, allele)

MJ contact energy of peptide on allele’s template (nan if untemplated).

ddg(wt, mut, allele)

MJ(mut) − MJ(wt) on the shared backbone (the WT/MT differential; <0 = mutant binds better). wt and mut must be the same length. nan if untemplated.

mhcmatch.ligand module

Full ligand spans: extend a binding core to the peptide that is actually presented.

Given a 9-mer binding core located in its source protein, presented_span() returns the most likely observed eluted-ligand span around it – the peptide a wet lab would synthesise, rather than the bare core. Three tiers of evidence, weakest last:

• observed – a reference ligand in the panel that contains the core and occurs in the protein. A real eluted span: the gold standard when it exists.

• modeled – the highest-scoring feasible span under SpanModel, a flank/context model fit to mass-spectrometry ligandome data.

• fixed – caller-specified flank sizes, clipped at the protein termini (fixed_span()).

This is not a cleavage predictor, and not an immunogenicity predictor. MHC-II peptides are generated bind-first-trim-later: the groove protects the core while exopeptidases erode the flanks, so there is no strong sequence-specific endoprotease step to simulate (Paul et al. 2018, PMID 30127785 – a dedicated MHC-II cleavage motif reaches AUC 0.767 on ligands and has zero predictive power on CD4 epitopes). What this models is P(observed ligand span | source protein), a convolution of protease specificity, HLA-DM editing, binding, stability and mass-spectrometry detection bias. Context/flank models are known to improve ligand prediction while degrading CD4 T-cell epitope benchmarks (Reynisson et al. 2020, PMID 32406916). Use this to enumerate and choose ligands to synthesise or model structurally – never to rank epitopes by immunogenicity.

For MHC-I the peptide is the ligand: there is nothing to extend, so there is no span function. processing_score() instead scores an 8-11mer’s source-protein context, the shape MHCflurry-2.0 uses for antigen processing (PMID 32711842). Class I and class II are deliberately different entry points – a 9-mer class-II core is always <=11 residues and would silently misroute through any length-based class inference.

mhcmatch.ligand.PAD = '-'

Out-of-protein context position (the span abuts a protein terminus). Modelled, not dropped: a ligand ending exactly at the protein’s C-terminus is evidence about where spans end.

mhcmatch.ligand.LIGAND_KEYS = ('ligN+1', 'ligN+2', 'ligN+3', 'ligC-3', 'ligC-2', 'ligC-1')

The ligand’s own terminal residues – inside the detected peptide, so subject to MS bias.

mhcmatch.ligand.FLANK_KEYS = ('flankN-3', 'flankN-2', 'flankN-1', 'flankC+1', 'flankC+2', 'flankC+3')

The residues flanking the ligand in the source protein – never in the detected peptide.

mhcmatch.ligand.CTX_KEYS = ('flankN-3', 'flankN-2', 'flankN-1', 'ligN+1', 'ligN+2', 'ligN+3', 'ligC-3', 'ligC-2', 'ligC-1', 'flankC+1', 'flankC+2', 'flankC+3')

3 upstream + 3 ligand-N + 3 ligand-C + 3 downstream. This is the NetMHCIIpan -context window (PMID 30446001); half the signal sits inside the ligand.

Type:

All 12 context positions

mhcmatch.ligand.STRUCTURE_FLANK = 2

across the 93 pMHC-II crystals of the Canonical2026 set the resolved peptide has a median length of 13 with a median of 2 flanking residues on each side, and only 13% resolve <=11 residues. So the core ± 1 (11mer) that TCRmodel2 and the fine-tuned-AlphaFold pipelines feed their networks is an input convention, not a statement about what is ordered – it discards real density in most structures. Reproduce with bench/pdb_flanks.py.

Type:

Flank size for structure prediction (core ± 2 = a 13mer). Measured, not guessed

mhcmatch.ligand.ASSAY_FLANK = 6

Flank size for a synthesised assay peptide (core ± 6 = a 21mer). What matters for a CD4 assay is not hitting the eluted boundaries exactly – the APC re-trims whatever you give it – but that the peptide contains the natural ligand. Measured on held-out eluted ligands, the fraction of cores whose full observed ligand is contained in the emitted peptide: 13mer 11%, 15mer 31%, 17mer 52%, 19mer 67%, 21mer 80%. Longer also tracks the MHC-II affinity optimum of ~18-20 aa (O’Brien et al. 2008, PMID 19036163). The conventional 15mer covers only 31%.

class mhcmatch.ligand.Span(peptide, start, end, core, core_start, source, score=0.0, n_alternatives=0, clipped=(0, 0), support=0)

Bases: object

A ligand span located in its source protein.

Parameters:

• peptide (str)

• start (int)

• end (int)

• core (str)

• core_start (int)

• source (str)

• score (float)

• n_alternatives (int)

• clipped (tuple)

• support (int)

peptide: str

start: int

end: int

core: str

core_start: int

source: str

score: float = 0.0

n_alternatives: int = 0

clipped: tuple = (0, 0)

support: int = 0

property flanks

(n_left, n_right) residues flanking the core within this span.

class mhcmatch.ligand.SpanModel(ctx, lens, padbg=0.02, background='markov', _m1=None)

Bases: object

Ligandome-fit flank/context model: P(observed ligand span | source protein).

ctx is allele-agnostic by construction: exopeptidase trimming is a property of the proteolytic machinery, not of the groove. That is measured, not assumed – per-allele context PWMs sit within JSD 0.003-0.010 of the pooled one for MHC-II – and pooling also unlocks the ~70% of class-II eluted-ligand records whose restriction is only a placeholder.

lens is a ligand-length prior, not a core-relative flank-length prior: defining an N-/C- flank length requires a binding core, and the allele-agnostic register is tied across >=2 frames on ~66% of real ligands, so such a prior would encode a tie-breaking rule rather than biology. N/C asymmetry is instead carried by the context positions, which are fit independently per side.

The span score is a plain log-likelihood, log P(L) + context log-odds, with no free parameters. A tuned weight on the length prior was tried – it looked better on the training fold and did not transfer (held-out set-recall 0.155 vs 0.158 unweighted, within noise), so it was dropped rather than shipped.

Parameters:

• ctx (dict)

• lens (dict)

• padbg (float)

• background (str)

• _m1 (dict)

ctx: dict

lens: dict

padbg: float = 0.02

background: str = 'markov'

context_score(protein, start, end, flank_only=False)

Log-odds of the context around span [start, end) vs the proteome null.

Parameters:

• protein – the source protein sequence.

• start – 0-based span start.

• end – 0-based exclusive span end.

• flank_only – score only the 6 FLANK_KEYS. The 6 ligand-internal positions carry the peptide’s own anchor signal, so including them partly measures binding; the flank-only score is the honest processing signal.

Returns:

Summed log-odds. Positions outside the protein score against padbg.

best_span(protein, core_start, core_len=9, delta=1.0)

Highest-scoring feasible span containing the core, as (start, end, score, n_alt).

Only spans that are real substrings of protein are enumerated, so the result never runs off a terminus. n_alt counts other spans within delta log-odds of the best – nested sets mean several spans are often legitimately correct.

The binding term is identical for every span sharing this core, so it cancels in the argmax and is omitted: ranking is driven purely by the length prior and the flank context. (Do not substitute AnchorModel.score() here – it is a max over register frames and so grows with peptide length, which would just select the longest span.)

mhcmatch.ligand.load_span_model(cls='mhc2', background='markov')

Load the vendored SpanModel for cls ("mhc1" | "mhc2").

Fit from IEDB mass-spectrometry eluted ligands against UniProt reference proteomes; see src/mhcmatch/data/PROVENANCE.md and bench/train_spans.py.

mhcmatch.ligand.observed_spans(core, protein, corpus, core_start=None)

Reference ligands that contain core and occur in protein, best first.

corpus is any iterable of peptide strings (e.g. store._panel['mhc2'].epitopes). Because the caller supplies the protein, ligand in protein is the provenance check – no source accession is needed. Occurrences that do not bracket the core are rejected: a ligand may also appear elsewhere in the protein.

This is a lookup, not a prediction. Never fold its hit rate into a prediction metric, and never report a training ligand as a novel result – check Span.source.

mhcmatch.ligand.fixed_span(core, protein, left, right, strict=False, core_start=None)

Extend core by left/right residues, clipped at the protein termini.

A requested flank that runs off the protein is reported, not silently shortened: the shortfall lands in Span.clipped.

Parameters:

• core – the binding core (must occur in protein).

• protein – the source protein sequence.

• left – residues requested upstream of the core.

• right – residues requested downstream of the core.

• strict – raise instead of clipping when the flank does not fit.

• core_start – 0-based start of core; defaults to its first occurrence.

Raises:

ValueError – if core is not in protein, or strict and the flank does not fit.

mhcmatch.ligand.presented_span(core, protein, model=None, corpus=None, mode='auto', flanks=(3, 3), core_start=None)

The most likely presented ligand span around an MHC-II binding core.

Parameters:

• core – the 9-mer binding core, located in protein.

• protein – the source protein sequence.

• model – a SpanModel; defaults to the vendored MHC-II model.

• corpus – reference ligands for the observed tier (e.g. panel epitopes). Optional.

• mode – "auto" (observed -> modeled), or force one of "observed" | "modeled" | "fixed". Benchmarks must pass "modeled": leaving observed on turns the metric into a coverage statistic.

• flanks – (left, right) for mode="fixed".

• core_start – 0-based start of core; defaults to its first occurrence. Pass it explicitly when the core repeats in the protein.

Returns:

A Span, or None when mode="observed" and no reference ligand contains the core – itself informative: the core has never been eluted.

Raises:

ValueError – if core is not a substring of protein, or is not 9 residues.

Warning

Do not pick a peptide to synthesise from the ``modeled`` span alone. Held-out, it places both boundaries within +-1 residue only 31% of the time and within +-2 only 47% – barely better than simply centring a 15mer on the core (28% / 50%), which actually wins at +-3 (79% vs 62%). Its real edge is the exact-span hit rate (0.158 vs 0.069), i.e. the question “what was actually eluted?” – not “what should I make?”. For the latter use fixed_span() with ASSAY_FLANK (a 21mer, which contains the true ligand 80% of the time) or STRUCTURE_FLANK (a 13mer, the median resolved crystal). See bench/results/spans_mhc2_human.md.

mhcmatch.ligand.processing_score(peptide, protein, model=None, flank_only=False, start=None)

Source-protein context log-odds of an MHC-I peptide – a score, never a span.

For MHC-I the peptide is the ligand, so there is nothing to extend. This scores how ligand-like its context looks (the antigen-processing signal MHCflurry-2.0 models, PMID 32711842) and composes into ranking, not into an emitted peptide.

Parameters:

• peptide – the 8-11mer, located in protein.

• protein – the source protein sequence.

• model – a SpanModel; defaults to the vendored MHC-I model.

• flank_only – score only the 6 flanking positions. The ligand’s own termini carry its anchor signal, so the full 12-position score partly measures binding rather than processing.

• start – 0-based start of peptide; defaults to its first occurrence.

Raises:

ValueError – if peptide is not in protein, or is not 8-11 residues.

mhcmatch.logo module

Per-allele motif logos (information content) + length distributions.

motif() returns the numeric logo (PWM, per-position bits, length histogram) – pure-Python, always available. render() draws it with logomaker (optional [logo] extra). MHC-I logos use peptides of a fixed length (default the modal length); MHC-II uses register-anchored 9-mer cores. See appendix/mhcmatch.tex §6.

mhcmatch.logo.motif(store, allele, cls, length=None)

Logo data for allele’s presented peptides.

Returns {allele, cls, width, n, pwm, bits, length_hist} where pwm[i] is a residue->freq dict (sums to 1), bits[i] the information content (log2(20) - entropy) in [0, log2(20)], and length_hist a length->count dict over all the allele’s peptides.

mhcmatch.logo.render(m, ax=None)

Render a motif() result as an information-content sequence logo (needs [logo]).

mhcmatch.predict module

Predict presented epitopes from a variant peptide-window FASTA.

Scores every binding-length k-mer of each window (the Gamaleya nextflow_vaccine pipeline’s .peptide.fasta) for a patient’s HLA alleles and emits two views:

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

• scored-csv (write_scored_csv()) – the same calls in the pipeline’s 57-column .epitopes.scored.csv schema, so mhcmatch can stand in for the MHCflurry/TLimmuno2 predictors.

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': (15,)}

Binding-length k-mers tiled per class (pipeline params.mhcI_epit_len / mhcII_epit_len).

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.

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='', 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)

• 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 = ''

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.parse_variant_header(header)

Parse a pipeline window header into variant-annotation fields.

Somatic: headers follow the fixed colon schema. Fusion: / CNV: use different internal delimiters, so only their type (and any of the shared trailing fields that line up) is extracted – best-effort, never raising: unknown fields come back empty.

Parameters:

header (str)

Return type:

dict

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

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)

peptide: str

allele: str

cls: str

presentation_rank: float

affinity_nm: float

affinity_rank: float

binder_rank: float

band: str

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.predict_windows(store, cls, records, alleles, rank_threshold=2.0, top=None, background='proteome', footprint='adaptive', seed=0)

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 (non-binders). Each kept 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)

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

Parameters:

path (str)

Return type:

None

mhcmatch.predict.write_scored_csv(preds, path)

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.

Parameters:

path (str)

Return type:

None

mhcmatch.mimics module

Molecular-mimicry annotation for strong binders.

For each strong-binding neoantigen, search reference peptide sets for mimics — near-identical presented peptides — and report the presentation-aware E-value (mhcmatch.search.find_mimics(), lower = more significant mimicry) per category:

• thymus — the thymic self-immunopeptidome (HLA Ligand Atlas). A significant thymic mimic means the neoantigen resembles a self-peptide presented during negative selection: reactive T cells were likely deleted (reduced immunogenicity) and it flags cross-reactivity / autoimmune risk for a cancer vaccine.

• viral / bacterial — foreign presented peptides / pathogen proteomes. A foreign mimic can raise immunogenicity (a pre-existing anti-pathogen repertoire cross-reacts) — molecular mimicry.

• neoag — the tested-neoantigen database: has this (or a near-identical) neoantigen been reported.

This scores cross-reactivity, not presentation or immunogenicity directly; compose it with the presentation / affinity scores from mhcmatch.predict. Reference data: the isalgo/pmhc_data compendium (thymus/, ligandome/, immunogenicity/, proteome/).

mhcmatch.mimics.DEFAULT_REFS = {'neoag': ('immunogenicity/neoag_tested.tsv.gz', 'database'), 'thymus': ('thymus/thymus_immunopeptidome.tsv.gz', 'self'), 'viral': ('ligandome/viral_foreign_iedb.tsv.gz', 'foreign')}

(folder/file under pmhc_data, kind). self is the tolerance reference passed as find_mimics’ self_set; the rest are foreign/database sets.

Type:

Default reference categories

class mhcmatch.mimics.MimicResult(binder, allele, category, n_exact, n_near, top_mimic, top_subs, e_value, n_hits, significant)

Bases: object

Per-(binder, category) mimicry summary.

A mimic is a reference peptide of the same length within near_subs substitutions of the binder (T cells cross-react across a few substitutions). n_exact / n_near count identical and near-identical mimics; top_mimic / top_subs are the closest one. e_value / n_hits are the raw presentation-aware search stats, kept for reference.

Parameters:

• binder (str)

• allele (str)

• category (str)

• n_exact (int)

• n_near (int)

• top_mimic (str)

• top_subs (int)

• e_value (float)

• n_hits (int)

• significant (bool)

binder: str

allele: str

category: str

n_exact: int

n_near: int

top_mimic: str

top_subs: int

e_value: float

n_hits: int

significant: bool

mhcmatch.mimics.load_peptides(pmhc_dir, rel_path, cls, species='human')

The peptide column of a compendium TSV, filtered to cls / species and plausible presented lengths. Rows without a class/species field are kept (some sets are unlabelled).

Parameters:

• pmhc_dir (str)

• rel_path (str)

• cls (str)

• species (str)

Return type:

list

mhcmatch.mimics.load_reference_sets(pmhc_dir, cls, species='human', refs=None)

(self_set, foreign_sets) for scan(). self_set is the single tolerance reference (the self-kind entry, thymus by default); foreign_sets is {name: [peptides]} for the rest. refs overrides DEFAULT_REFS.

Parameters:

• pmhc_dir (str)

• cls (str)

• species (str)

Return type:

tuple

mhcmatch.mimics.scan(binders, self_set, foreign_sets, cls='mhc1', max_subs=2, near_subs=2, self_name='thymus')

Mimic-scan an iterable of (peptide, allele) binders. Returns list[MimicResult] (one per binder × category with >=1 same-length reference peptide within near_subs substitutions).

self_set is the tolerance reference (category self_name); foreign_sets is {name: [peptides]}. max_subs is the fuzzy-search radius. find_mimics() excludes the exact query (a neoantigen’s identical peptide is its source, not a mimic), so n_exact is a direct set-membership check and n_near counts same-length reference peptides 1..``near_subs`` substitutions away (from the fuzzy hits, by exact Hamming distance). One find_mimics() call per binder scores every category at once.

mhcmatch.mimics.patient_summary(results, binders)

Aggregate scan() output into patient-level counts for a dashboard row.

binders is the full strong-binder list (so “0 mimics” binders are counted too).

Return type:

dict

mhcmatch.mimics.write_table(results, path)

Write per-(binder, category) mimic results as a TSV (one row per category with a near mimic).

Parameters:

path (str)

Return type:

None