Store, search and reference data

The entry points, and the deposits everything else reads from.

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: the theory appendix §2-3 (forward per-allele E-value + reverse problem).

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

The ligand lengths tiled per class, and the ONE definition of them. predict.KMER_LENS is this object, not a copy: the package carried four ladders that agreed on class I (8-11) and disagreed on class II in every one – (15,) on the rank path, 13-18 here, 12-15/12-20 in vector, 11-25 in mimics – which is how “a class-II ligand runs to 21” and “we tile 15-mers only” were both true of one library. mimics.CANONICAL_LEN stays separate on purpose: it sizes proteome indexes, where a width costs disk rather than accuracy.

Class I runs to 11 and class II to 21, which is the author’s specification and not a fit to the panel. The panel agrees with it: class-II epitopes run 9-23+ with only 24.5 % at length 15 (84,405 of 343,956), and 12-21 covers 93.5 % of that mix (321,550) against 12-20’s 92.0 %. An earlier revision of this constant stopped at 20 because a benchmark ladder (vector.MHC2_MAP_LENGTHS) did; that was a ladder chosen for a different job, and the stated maximum is 21. The %rank null was always drawn from the full mix, so the single-width tiler was the thing out of step with it, not the wide one. The cost is nothing: 291 -> 2,900 windows on a 305-aa record against a 209 s per-run scorer build.

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. Resolves through fetch_file(), so a local mirror at $MHCMATCH_PMHC_DIR – the dataset root the SLURM profile exports – is used before any download; path= / $MHCMATCH_PMHC, which name the holding directory rather than the dataset root, still win over both.

Parameters:

tier (str)

Return type:

str

mhcmatch.store.fetch_file(relpath)

Download any file of the public HF dataset PMHC_REPO by its repo-relative path.

The escape hatch behind fetch_pmhc() and fetch_proteome(), for the deposits that do not have a named accessor – the immunogenicity corpora, the thymus immunopeptidome, the viral ligandome, the neoantigen screens. It exists so a worked example can run on a whole published deposit rather than a hand-copied excerpt, which is the only version of an example that demonstrates anything about scale.

$MHCMATCH_PMHC_DIR overrides with a local mirror (e.g. ~/hf/pmhc_data), for offline and cluster runs. Each file is fetched once and cached by huggingface_hub.

>>> from mhcmatch import store
>>> store.fetch_file("immunogenicity/chowell_rebuilt.tsv.gz")

Parameters:

relpath (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, vote, enrichment, n_votes, binder, anchor_score=None, rank=None, p_present=None, band=None)

Bases: object

One allele’s presentation call for a peptide – one row of Store.restriction()’s result, ranked by vote (or anchor_score under diffuse=True).

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, tcr_facing, presentation, anchors)

Bases: object

A peptide split into its anchor and TCR-facing halves – Store.decompose()’s result.

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.CORE_LEN: int = 9

Length of the class-II binding core, and the length of a class-I core whenever the peptide is long enough to fill one. Same 9 as mhcmatch.ligand.CORE_LEN.

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

The binding core of peptide and its 0-based offset, NetMHCpan-style.

Returns (core, offset); ("", -1) when the peptide is too short to carry one.

The core is residues, never a padded frame. NetMHCpan defines theirs as “the minimal 9 amino acid binding core directly in contact with the MHC (i.e. excluding potential insertions)”, and the parenthesis is the operative part: where an alignment to a 9-mer motif needs a gap, the inserted position is not part of the core. So this returns CORE_LEN residues whenever the peptide can fill one and the peptide’s own residues when it cannot – 9 for a class-I 9/10/11-mer and for every class-II core, 8 for a class-I 8-mer. A gap character in an amino-acid column would not be neutral anyway: B is Asx in IUPAC, so a reader would take it for a real ambiguity.

Class I is the signed footprint mhcmatch.diffusion.MHC1_CORE – front P1-P5 plus C-terminal P-4..P-1 – resolved by mhc1_positions(), which is the same mapping the scorer uses, so the reported core is the residues the model actually read. Both ends are therefore held and the middle gives way, which is NetMHCpan’s rule: at L=9 the core is the peptide; at L=10 or 11 the central one or two residues drop out, their Gp/Gl deletion. Below 9 the +5 and -4 positions collide, mhc1_positions() yields None for the loser, and that slot is dropped rather than padded – every residue still appears exactly once, so an 8-mer’s core is the 8-mer. The offset is 0: the footprint is anchored at both ends, so there is no N-terminal protrusion and nothing here corresponds to NetMHCpan’s Of > 0.

Class II is peptide[s:s+9] at the register s, and the offset is s – the same quantity NetMHCIIpan reports as Of, “starting position offset of the optimal binding core (starting from 0)”. register_start pins the frame; None falls back to the allele-agnostic heuristic _mhc2_register(). Pass the model register when you have one. The two disagree often on real ligands (see anchor_indices()), which is why every caller that emits a core also emits where the register came from.

>>> binding_core("SIINFEKL", "mhc1")            # 8-mer: its own core, nothing inserted
('SIINFEKL', 0)
>>> binding_core("GILGFVFTL", "mhc1")           # 9-mer: the core is the peptide
('GILGFVFTL', 0)
>>> binding_core("GILGFVFTLA", "mhc1")          # 10-mer: the central residue drops out
('GILGFFTLA', 0)
>>> binding_core("PKYVKQNTLKLAT", "mhc2")       # HA306-318, heuristic register
('YVKQNTLKL', 2)

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.

species: str | None = None

Which species’ panel this store holds, or None for a mixed/unfiltered one.

Set by from_pmhc() from its own species= argument and read by mhcmatch.diffusion.load_vendored_anchor_model(), which otherwise cannot tell a mouse panel from a human one and so cannot pick the right pre-fit model. panel_sha still decides whether a vendored model is valid; this decides which one is worth opening.

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)

Every allele in this class’s panel, as loaded (not filtered by frequency or coverage).

panel_alleles(cls, alleles='all', missing=None)

Public _allele_set(): the queried names in the panel’s spelling, "all" for every one.

A caller outside this class must not reach into self._panel to do this – _allele_set is the one place that folds both of the package’s allele vocabularies together and the one place that reports what it could not match.

Parameters:

missing (list | None)

percent_ranks(peptides, cls=None, alleles='all')

[{allele: %rank}, ...], one dict per peptide – restriction()’s ranking half with the neighbour tally skipped.

Why a second entry point rather than a flag. restriction() opens with panel.tally(peptide) unconditionally, and tally is a scope-widening loop over _SCOPES around KmerIndex.seed_and_gather – the call bench/results/neighbour_search_speed.md measures at 55 queries/s. Its result feeds vote, enrichment, n_votes and the vote half of the binder gate, and a caller that wants only %rank reads none of them. Measured on 400 distinct 9-mers over 6 allotypes: restriction(calibrated=True) costs 0.307 ms/peptide of which the tally is 0.244 ms – 79.5% spent on fields the caller discards. mhcmatch.vector.store_binder() is exactly that caller, and a cassette layout asks it millions of times.

Returning the uncomputed fields as zeros would be worse than not offering the path: vote 0.0 and binder False are readable as answers. So they are absent – this returns ranks and nothing else, and a caller who needs the vote gate calls restriction().

A NaN rank (MHC-II, an allele with no length-matched background) is left as NaN rather than dropped, so the allele still appears and the caller decides.

Return type:

list

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

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)

Does this one allele present peptide? See is_presented() for “any allele”.

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, *, threads=1)

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', reverse=0.0, rare_max=30, background='ligand', length_prior='score', length_motifs=True, register='marginal', n_motifs=3, pseudocount=0.0, anticore=0.0, pseudo_matrix='blosum62', families=None, route=None, _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 earlier 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(). anticore (MHC-II) weights a pooled flank model of the residues outside the core, so frames are compared on the whole ligand rather than on nine positions each; 0 (default) is off – see mhcmatch.diffusion.AnchorModel._fit_anticore(). route (a dict of parameter overrides, None by default) fits a second model with those overrides and sends alleles at or below rare_max ligands to it – the rare and frequent class-II optima are incompatible in one fit, see mhcmatch.diffusion.RoutedAnchorModel.

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 the theory appendix §5.

class mhcmatch.search.Match(peptide, shared_kmers, score)

Bases: object

One hit from search() – score is the seed-and-gather edit distance under mode (lower is closer), shared_kmers the count of anchor-masked k-mers the pair have in common.

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=1)

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 the proteome – using seqtree.TextIndex. This is a distinct mode from the anchor-masked TCR-facing homology and the presentation-signature searches. See the theory appendix §5 (near-exact source identification).

mhcmatch.proteome.read_fasta(path)

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

mhcmatch.proteome.gene_symbols(path, key='name')

{key: gene} from the UniProt GN= field. key="name" (default) matches read_fasta(); key="accession" matches a bare UniProt accession.

Both keyings exist because two different callers need two different sides of the same header. A SourceHit names its protein as the FASTA’s first whitespace token, sp|Q8WZ42|TITIN_HUMAN, so a proteome scan needs name. The thymic and ligandome deposits record source_protein as a bare accession, Q8WZ42, so mhcmatch.mimicry.safety() needs accession. Neither can reach mhcmatch.expression.safety_profile(), which is keyed on the HGNC symbol TTN, without one of them.

Without it there is no way to ask which tissue a T cell cross-reactive with a given self peptide would attack, and that is the question separating a titin match (Q8WZ42 → TTN → heart left ventricle, 64 TPM) from a testis-restricted one.

The symbol is absent from read_fasta()’s output because that function keeps only the name, and widening its return contract would ripple through every caller. A second pass over the headers is cheap – one second for the human proteome – and additive.

Entries with no GN= map to None rather than being dropped: the 147,506 human records include TrEMBL entries with no assigned symbol, and silently losing them would overstate the coverage of any downstream tissue filter.

Parameters:

key (str)

class mhcmatch.proteome.SourceHit(protein, position, ref_peptide, n_subs, mutations)

Bases: object

One near-exact match of a query peptide against a reference window – the result row of Proteome.find_source() / find_sources() / find_exact_sources().

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, searched by peptide through one seqtree.TextIndex.

classmethod from_fasta(path)

A Proteome over a caller-supplied FASTA, overriding from_hf()’s fetch.

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.

The single-query form of find_sources(). It used to be the wrong entry point for anything but an interactive question, because the index build dominated it; the index is now one shared build of 0.7 s, so asking about one peptide costs one query.

find_sources(peptides, max_subs=1, exclude_exact=False, threads=1, best_only=False)

{peptide: [SourceHit, ...]} for many peptides at once – the batch form of find_source().

One index and one threaded C++ batch query for the whole set, whatever lengths it spans: search_batch releases the GIL; threads=0 opts into allocated cores. The per-length loop this replaced is gone, and with it the advice to ask only for the lengths you need – there is nothing left here that is per length.

best_only=True returns only the nearest non-empty shell. That is what a caller voting on a parent wants (assign_genes()), and asking the search for it is far cheaper than materialising the outer shells in order to discard them: measured on 487,000 peptides at radius 2, 2.2 s for 6,824,481 hits against 21.6 s for 15,822,792.

Duplicate and blank queries are collapsed; the returned dict is keyed by the stripped, upper-cased peptide. A query carrying a residue outside the 20 standard amino acids maps to [], because every window that could have matched one is excluded from the answer anyway.

find_exact_sources(peptides)

{peptide: [SourceHit, ...]} at exactly zero substitutions – find_sources() with max_subs=0.

Same return shape and the same (protein, position, ref_peptide, n_subs, mutations) content; n_subs is 0 and mutations is () for every hit, because that is what an exact match is. Peptides with no source come back with an empty list, so the dict is keyed by every distinct stripped, upper-cased query.

It survives as a name, not as a separate implementation. It existed because the index it would otherwise have queried was a Python loop over every position of every protein (~12.6 GB peak at L = 9), buying an ability an exact question never uses – so it answered membership out of a sorted window array and a pair of np.searchsorted calls instead, in 11.0 s. One TextIndex build is 0.7 s and answers radius 0 as directly as any other radius, so the second construction bought nothing and is gone. The name stays because mhcmatch.vector.self_origin_risk() dispatches on it, and that dispatch is how the code says the safety screen asks an exact question.

windows(L)

Every distinct length-L standard-AA window of the proteome, as a set. Cached, and roughly 1 GB per length for the human proteome – ask for the lengths you need.

window_array(L)

Every distinct length-L standard-AA window, as a sorted |S{L} numpy array.

The vectorized form of windows(), and 2.7x faster on the human proteome: 11.0 s against 30.0 s for 12,073,995 distinct 9-mers, identical output. The loop it replaces ran all(c in _AA for c in w) per window – 12 M windows x 9 residues of Python-level membership tests.

The array form is what a consumer that projects or indexes these wants; windows() still returns a set for the O(1) membership its own callers need, and materialising a 12 M-element Python set is most of the cost that buys.

Packing the residues into uint64 (5 bits each, L <= 12) and sorting integers instead was tried and is 4x slower – 44.5 s – because the shift/or loop costs more than numpy’s fixed-width byte sort saves. Measured, not assumed.

window_genes(peptides, path=None)

{peptide: gene_symbol} for those of peptides that are proteome windows.

The question a neoantigen table has to answer before it can look up expression: a candidate carries a somatic substitution, so it is not itself a window, but its wild type is – and that window names the gene. wildtype() supplies the germline counterpart and this supplies the symbol, which is what GTEx and TCGA are keyed on.

Streams seqs once and keeps only matching windows, rather than indexing the whole proteome and querying it: the query set is known in advance and small (~350k) where the index is ~68M windows per length. path is the FASTA the symbols are read from with gene_symbols(); it defaults to the one this proteome was loaded from.

assign_genes(peptides, max_subs=2, threads=1, path=None)

{peptide: [gene, ...]} – the HGNC symbol(s) of the gene each peptide derives from.

window_genes() answers this for a peptide that is a proteome window. A neoantigen is not: it carries the substitution that made it one, so it has to be found by near-exact search. That is what makes this the entry point a candidate table needs – expression.gene_level and both fitted expression terms are keyed on the symbol, and a row without one contributes a single mean-imputed constant. Measured over the benchmark corpus, 356,387 of 695,811 rows (51.2%) and 5,205 of 5,833 positives (89.2%) carried no deposited symbol; on the VACCIMEL screen expr_norm had standard deviation exactly 0.0000 and AUROC exactly 0.5000 while carrying the largest positive coefficient of the then-shipped EPIC v10 artifact, +0.4950 log-odds per standard deviation. Repairing the symbol is what took that term to +0.2155 in v11, on a measurement rather than a constant.

Three choices, all load-bearing:

• Only the nearest shell votes. A radius-2 shell is ~85x the size of the radius-1 shell inside it, so pooling them lets a distant coincidence outvote a genuine single-substitution parent. This is best_only=True – a property of the search, not a filter over its output, and 10x cheaper than the filter it used to be (find_sources()).

• Exact matches are excluded. A peptide that is a proteome window is not a neoantigen, and its own gene is not the question being asked.

• Ties come back in full, sorted. Resolving one needs expression data this method does not have – the caller picks among the tied genes (the CLI emits a row each and lets the scorer’s best-per-peptide selection decide). A peptide with no parent, or whose parents carry no GN=, maps to []; neither is an error.

max_subs defaults to 2 because a neoantigen can carry more than one mutation, and the radius is what buys the coverage: at radius 1 VACCIMEL resolves 88.2% of its peptides, at radius 2 96.8% (TESLA 98.5%, ITSNdb 99.5%, GBM 94.0%, Sahin_TNBC 100%). bench/results/gene_resolution.md.

Nothing here is per hit, and at this scale that is the difference that matters. The nearest shell of a 487,000-peptide radius-2 query is 6,824,481 hits; one SourceHit dataclass apiece is ~1.4 GB of Python objects for a function that only ever reads (n_subs, protein). So this reads the flat arrays out of search_batch directly and maps them through a ref_id -> gene integer table built once over the proteome. path is the FASTA the symbols are read from, and defaults to the one this proteome was loaded from, as in window_genes().

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

One peptide at a time. wildtypes() is the batch form and is what a corpus should call: the index is shared either way, and search_batch releases the GIL within the declared native-thread budget.

wildtypes(peptides, max_subs=1, threads=1)

{peptide: wild-type | None} – the batch form of wildtype().

One threaded search_batch for the whole corpus, keyed by the stripped, upper-cased peptide. It replaces two constructions at once: a per-peptide Python loop over the caller’s list, and the L * 19-variant hash-set membership test wildtype() used to run against a _window_set costing ~1 GB per length.

The tie order is part of the contract, and it does not survive the port for free. The hash-set path walked the peptide left to right and _AA_ORDER within each position, so a peptide with two one-substitution parents got the earlier position, and the earlier residue at that position. TextIndex orders hits by (n_subs, ref_id, offset) – a different answer on the same data, feeding agretopicity in a shipped fit. So the radius-1 shell is re-sorted here rather than inherited. Beyond radius 1 there was never a documented order and there is none now: the nearest shell wins, and ties within it go to the index’s own order.

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 the theory appendix §4.

mhcmatch.pseudoseq.trim_allele(a)

IMGT allele name -> its two-field form, dropping any G/P group or expression suffix.

'A*01:01:01G' -> 'A*01:01', 'DRB1*15:01:01' -> 'DRB1*15:01', 'B*44:02:01:02S' -> 'B*44:02'. Names already at two fields, mouse H-2 names and the separator-free pair keys ('HLA-DQA10501-DQB10301') are returned unchanged.

Every HLA typer emits more than two fields. OptiType, kourami, HLA-LA, arcasHLA and HLA-HD all write the G-group form (A*01:01:01G), which is what a donor’s own .alleles.tsv carries – and the pseudosequence tables are keyed at two fields, because that is the depth at which the groove is determined. Without this trim resolve_allele() returns (None, False) for every allele of such a file, and mhcmatch.store.Store._allele_set() drops what it cannot find silently, so the run scores against an empty panel and says nothing. The failure is the one normalize_allele() records for 'H2-Kb', reached from the other side.

A pair name is trimmed chain by chain: 'DQA1*05:01:01-DQB1*03:01:01' -> 'DQA1*05:01-DQB1*03:01', because the pattern matches the digit fields and not the locus.

Parameters:

a (str)

Return type:

str

mhcmatch.pseudoseq.normalize_allele(a)

pmhc allele name -> pseudosequence-FASTA key.

Drops the * ('HLA-A*02:01' -> 'HLA-A02:01') and folds all three mouse H-2 spellings onto one key: pmhc 'H-2Kb', deposit 'H2-Kb' and FASTA 'H-2-Kb' name the same molecule, and mhci_pseudo.fa carries the last two as separate keys on a byte-identical 34-mer. Earlier only the first was folded, so 'H2-Kb' resolved exact=True to a key with zero panel ligands and SIINFEKL scored at presentation %rank 20.19 instead of 0.0040. One molecule, one key – the invariant hla_spellings() already enforces for human class I and class2_from_name() for class II.

Takes one allele name. A cell naming several ('B0801,C0701') is a genotype, not an allele; mhcmatch.rank.split_alleles() is what splits it, and normalising the cell whole is the defect that produced 'HLA-B08:010701' – two names run together into a spelling no table has, which then resolves to nothing.

Parameters:

a (str)

Return type:

str

mhcmatch.pseudoseq.hla_spellings(name)

Both spellings of a class-I HLA name – with the field colon and without it.

The bundled pseudosequence table carries both: HLA-A02:01 and HLA-A0115 are each keys, because the source tables it was built from disagreed. So does every deposited screen, which writes A0201, Cw0401, HLA A0201 or A*02:01 for the same molecule. Offering both spellings is what lets resolve_allele() accept all of them without a caller normalising first – and a benchmark that normalises in its own helper is a second convention nobody else can run. Returns [] for anything that is not a class-I HLA name.

Parameters:

name (str)

Return type:

list

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.REPORT_MODES = ('pair', 'beta', 'isotype')

Class-II reporting granularities, finest first. See class2_report().

mhcmatch.pseudoseq.class2_report(key, mode='pair')

Reduce a class-II key to a reporting granularity.

• "pair" – the key unchanged: 'DRB1_0101', 'HLA-DQA10501-DQB10301'. This is NetMHCIIpan’s own naming and what class2_key() produces, so it is the default and the only mode in which two tools’ outputs are directly comparable as strings.

• "beta" – the beta chain alone, in IMGT form: 'DRB1*01:01', 'DQB1*03:01'.

• "isotype" – 'DR' / 'DP' / 'DQ' (mouse: 'H-2').

Why the coarser modes exist. A class-II key does not lead with the same chain at every isotype: DRA is monomorphic, so DR is keyed by its beta, while DP and DQ keys lead with the alpha. Any comparison that reads the leading gene out of a key is therefore matching DR’s beta against DP/DQ’s alpha – two different genes, and the alpha is the less polymorphic half. It also splits DR against itself, because DRB1 and DRB3 are different leading genes at the same isotype. "beta" and "isotype" both compare like with like; "isotype" is the right granularity for the question “did the two callers even pick the same molecule family”.

Measured on the class-II arm of a two-caller concordance study (10,402 rows where both callers named an allele): leading-gene agreement 0.401, true isotype agreement 0.527. The gap is 1,318 DR-vs-DR pairs differing only in DRB gene.

Parameters:

• key (str)

• mode (str)

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.

Memoised: a miss walks every key and sorts the prefix hits, which is 673 us against 0.3 us for a hit. Callers reach it once per background peptide inside a calibration build, so an unresolvable name used to cost ~6.7 s per allele instead of one lookup.

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.PSEUDO_MATRICES: tuple = ('blosum62', 'blosum45', 'blosum80', 'pam250', 'pam100')

Substitution matrices seqtree carries, and therefore the whole menu for the pseudocount blend. Note what is NOT here. There is no BLOSUM90 or BLOSUM100 to try, and seqtree’s structural is excluded on purpose: its entries are all non-negative (mean +6.6), so it has no positive Karlin-Altschul scale and is a similarity score rather than a log-odds matrix.

The two families run in OPPOSITE directions. BLOSUM-n is built from blocks clustered at >= n% identity, so a HIGHER number means closer relatives and a more conservative substitution model; PAM-n is n accepted point mutations per 100 residues, so a HIGHER number means MORE divergence. Roughly BLOSUM80 ~ PAM120, BLOSUM62 ~ PAM160-200, BLOSUM45 ~ PAM250.

mhcmatch.pseudoseq.substitution_conditional(matrix='blosum62')

{observed: {r: P(r | observed)}} – a substitution conditional from any seqtree matrix.

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.

BLOSUM62_BG is used as p_a for every matrix. That is exact for BLOSUM62 and an approximation for the others, whose own background differs; the identity is P(a|b) ∝ p_a·2^(s_ab/2), so a wrong p tilts the conditional without changing which residues it calls similar. Stated because it bounds what a matrix sweep can conclude.

Parameters:

matrix (str)

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)

Groove similarity in (0, 1], exp(-distance(a, b) / h); 0.0 if either allele’s pseudosequence is unresolved.

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.pseudoseq.blosum62_conditional()

The BLOSUM62 conditional — substitution_conditional() at its default.

Kept as a name because it is what the anchor model has always called; the matrix is now a parameter rather than a constant.

Return type:

dict

mhcmatch.expression module

Reference expression by normal tissue (GTEx) and by tumour type (TCGA), fetched from the public isalgo/pmhc_data dataset. The two are never merged — different measurements, different units.

Reference expression, keyed by normal tissue (GTEx) and by tumour type (TCGA).

A neoantigen ranker needs two different expression questions answered, and they are not the same number, so this module never merges them:

Is the source gene transcribed in the tumour’s lineage? The ranking read. A candidate from a gene that is silent in that tissue is not presented however well it binds. Melanoma is SKCM.

Is it also transcribed in normal tissue? The safety read. A gene expressed everywhere is a toxicity risk, not a target.

Two key types, joined differently, both carrying their own provenance:

• key_type="gene" x a GTEx tissue – per-tissue median TPM over the v11 bulk release. This is what imputes a missing TPM and what answers the safety question.

• key_type="peptide" x a TCGA ``cancer_type`` – the median expressionEB++ over every TCGA sample in which that exact peptide was called. Keyed on the peptide rather than the gene deliberately: the TCGA source table carries ensp and no ENSP->symbol map ships with it, so a gene-level join would be a guess. The peptide-level join is exact, and it answers a question the gene-level number cannot – has this exact neoantigen been seen expressed in this tumour type.

The table is fetched from the public HF dataset mhcmatch.store.PMHC_REPO and cached, like every other reference here.

>>> from mhcmatch import expression
>>> expression.lookup("PMEL", tissue="Skin - Sun Exposed (Lower leg)")
{'median_tpm': 44.1, 'q25_tpm': 21.6, 'q75_tpm': 78.9, 'n': 605, 'source': 'gtex'}

Missing is encoded, never dropped. impute() returns the reference value and a flag saying whether it was observed or imputed, so a caller can carry a missing-indicator column instead of discarding the candidate – which is the standing rule for every partially-covered covariate here.

``species=”mouse”`` reads two deposits where human reads one, and answers the same two questions. Normal tissue is FANTOM5 (REFERENCE_FILE_MOUSE), gene x 35 adult tissues; the tumour side is TUMOR_FILE_MOUSE, gene x 6 syngeneic models. The one structural difference that survives is the key: a human tumor= lookup is peptide-keyed (TCGA has per-peptide rows) and a mouse one is gene-keyed (no mouse deposit has peptide rows at all). The one thing still genuinely absent is the tumour-to-matched-normal map TUMOR_TISSUE, which is keyed by TCGA study code – so tissue_floor() still refuses a mouse tumor= and says to name the tissue instead. Every signature here defaults to "human", so a human call is unchanged.

mhcmatch.expression.REFERENCE_FILE = 'expression/reference_expression.tsv.gz'

Path inside the HF dataset repo.

mhcmatch.expression.REFERENCE_FILE_MOUSE = 'expression/reference_expression_mmu.tsv.gz'

The mouse counterpart, and it is the same table shape rather than an analogue of it. RIKEN FANTOM5 mouse CAGE via EBI Expression Atlas E-MTAB-3579: 659,050 rows, 18,830 gene symbols x 35 adult tissues including thymus, source = "fantom5_mouse", and COLUMNS column for column. That is why species is a file selection here and not a second code path.

Three things it is not, all recorded in the deposit’s own expression/SOURCES.md:

• No peptide rows, so no tumour half. tumor= raises for mouse rather than resolving to anything – there is no mouse TCGA, and a mouse candidate scored against a human tumour’s abundance floor is the failure this parameter exists to prevent.

• ``n`` is 1 on every row. FANTOM5 gives one library per adult tissue, so q25/q75 describe the spread across a gene’s transcripts, not across animals. The human n runs to several hundred donors; the two columns are not the same quantity.

• CAGE tag density reported as TPM is not RNA-seq TPM. It behaves like expression and ranks like it; it is not numerically interchangeable with the human table across a join.

mhcmatch.expression.TUMOR_FILE_MOUSE = 'expression/tumor_expression_mmu.tsv.gz'

The mouse tumour read. The human deposit carries both halves in one file – gene x tissue under key_type="gene" and peptide x TCGA study under key_type="peptide". Mouse has no TCGA and no peptide-keyed tumour measurement at all, so its tumour half is a separate deposit and is gene x syngeneic model: 149,640 rows, 24,940 gene symbols x 6 models (B16F10, CT26, E0771, LLC, MC38, Panc02), three biological replicates each, source = "gse245293_syngeneic", from NCBI GEO GSE245293. COLUMNS column for column, which is why load() can fold it into the same table.

Two properties to know before reading a number off it:

• The units are computed. The deposit is normalised counts with a Length column; length-normalising per sample and rescaling to a million gives a TPM-like quantity comparable within this table and interchangeable with neither FANTOM5 CAGE nor GTEx RSEM TPM. A floor for these contexts must therefore be a quantile of this same column, which is what context_floor() takes.

• ``n`` is 3 on every row, three animals – so unlike REFERENCE_FILE_MOUSE, whose q25/q75 span a gene’s transcripts, these quantiles span replicates.

It covers six models, not every model a mouse experiment names: P815 and the T3/CMS5/Meth A fibrosarcomas are absent and take the pooled floor.

mhcmatch.expression.REFERENCE_TOIL_FILE = 'expression/reference_expression_toil.parquet'

The companion table in which GTEx and TCGA are the same unit, TPM on every row, because both cohorts went through one RSEM pipeline. REFERENCE_FILE cannot answer a question that spans the two – its GTEx half is TPM, its TCGA half is RSEM normalised counts, and both are written into median_tpm – so anything that compares a tumour type with a tissue, or takes a scale from one to divide into the other, reads this one instead. Gene-keyed throughout: 53 GTEx tissues under source="toil_gtex" and 33 TCGA study codes under source="toil_tcga".

Nothing here parses it, and it is not in the bootstrap set. It is 38.6 MB and the record; the scoring path reads MATRIX_FILE and SYNONYMS_FILE, which with FLOORS_FILE are 6.6 MB between them and carry the same numbers. Fetch it with fetch_reference(file=REFERENCE_TOIL_FILE) and read it with polars when an analysis wants the rows themselves.

mhcmatch.expression.MATRIX_FILE = 'expression/toil_matrix.npz'

The whole single-pipeline table as a dense gene x context float32 matrix, which is what every lookup on the scoring path reads. Measured against the same question asked of REFERENCE_TOIL_FILE, which builds a five-million-entry dictionary of dictionaries: 0.05 s against 5.20 s to load, and 29 MB resident against 3,168 MB – 100x the speed on 1/109 of the memory. Both return the same floor to four decimals, which is how the matrix is checked.

mhcmatch.expression.MATRIX_FILE_MOUSE = 'expression/toil_matrix_mmu.npz'

The mouse counterpart of the two files above, and the reason mouse now has a tumour-versus- normal contrast at all. 26,737 genes x 68 contexts: 35 FANTOM5 tissues, 6 GSE245293 syngeneic models and 27 more from the CrownBio panel (prjna1271699+). Same array names and the same source|context key shape as MATRIX_FILE, so one code path reads both.

Harmonised across its three sources, which the two mouse TSV deposits are not. They sat on three scales – pooled mean log2(1 + TPM) of 1.721 (CAGE), 2.513 and 2.059 (two RNA-seq sets) – so before this file a mouse expr_lvl and expr_norm divided by floors from different assays, 0.9964 against 0.8. Here the pooled q25 floor is 0.7174 / 0.7144 / 0.7174. Each source’s non-zero values are quantile-mapped onto the mean of the three quantile functions, per source rather than per context, zero maps to zero, and the map’s grid spans each source’s support – so every within-context rank is preserved exactly, which is what mhcmatch.rank.expr_level() consumes.

Two things it does not fix, both in the deposit’s expression/SOURCES.md: a gene-specific platform effect (Actb reads 3,032 in B16F10 against 10.4 in FANTOM5 thymus, because CAGE at that promoter ranks it mid-distribution), and the detection-rate difference (65.2 / 78.1 / 65.0 % non-zero), which a monotone map respecting zero cannot move.

mhcmatch.expression.FLOORS_FILE = 'expression/toil_floors.tsv'

The same contexts’ abundance floors at three quantiles, human-readable, 88 rows. Nothing on the scoring path parses it – context_floor() computes from the matrix, and this file is what a table or a caption cites.

mhcmatch.expression.FLOORS_FILE_MOUSE = 'expression/toil_floors_mmu.tsv'

The 68 contexts’ floors at q05/q10/q25 plus a pooled row per source, 71 rows, same seven columns as FLOORS_FILE.

mhcmatch.expression.SPECIES = ('human', 'mouse')

The species this module holds a normal-tissue reference for.

mhcmatch.expression.reference_file(species='human')

The reference deposit for species. Raises on anything else, rather than defaulting.

Silently falling back to the human table is the one failure mode worth spending a raise on: a mouse gene symbol misses every HGNC key, so the whole column would impute to the training mean and look like a covariate the model simply did not like.

Parameters:

species (str)

Return type:

str

mhcmatch.expression.matrix_file(species='human')

The dense matrix deposit for species. Raises on anything else, as reference_file() does and for the same reason: silently reading the human matrix for a mouse gene misses every key and imputes the whole column to the training mean.

Parameters:

species (str)

Return type:

str

mhcmatch.expression.SYNONYMS_FILE = 'expression/context_synonyms.tsv'

How a free-text origin becomes a context. 232 rows, long format, one target per row, each carrying where it came from: the TCGA study codes and GTEx tissues present in the matrix, Xena’s _primary_site and detailed_category joined to the study code on the sample barcode, and a short list of curated spelling variants marked curated:spelling. Derived rather than written by hand because an organ maps to more than one study more often than not – Lung is LUAD and LUSC, Kidney is three – and a crosswalk that picks one is silently wrong for the rest. Kept as a deposit rather than a module constant so it can be corrected without a release.

mhcmatch.expression.COLUMNS = ('key', 'key_type', 'source', 'context', 'median_tpm', 'q25_tpm', 'q75_tpm', 'n')

Columns of the reference table, in file order.

mhcmatch.expression.TUMOR_TISSUE: dict[str, tuple[str, ...]] = {'BLCA': ('Bladder',), 'BRCA': ('Breast - Mammary Tissue',), 'CESC': ('Cervix - Ectocervix', 'Cervix - Endocervix'), 'CRC': ('Colon - Transverse', 'Colon - Sigmoid'), 'GBM': ('Brain - Cortex', 'Brain - Frontal Cortex (BA9)'), 'HNSC': ('Minor Salivary Gland', 'Esophagus - Mucosa'), 'KICH': ('Kidney - Cortex',), 'KIRC': ('Kidney - Cortex',), 'KIRP': ('Kidney - Cortex',), 'LIHC': ('Liver',), 'LUAD': ('Lung',), 'LUSC': ('Lung',), 'OV': ('Ovary', 'Fallopian Tube'), 'PAAD': ('Pancreas',), 'PRAD': ('Prostate',), 'SKCM': ('Skin - Sun Exposed (Lower leg)', 'Skin - Not Sun Exposed (Suprapubic)'), 'STAD': ('Stomach',), 'THCA': ('Thyroid',), 'UCEC': ('Uterus',)}

Which normal tissue is a tumour type’s matched normal, so the safety read can be asked without the caller having to know that melanoma pairs with skin.

The two vocabularies are different and neither is clinical, which is worth being explicit about:

• the keys are TCGA study abbreviations (NCI GDC), a research nomenclature. CRC is the one exception – TCGA itself has COAD and READ separately, and the source table merged them.

• the values are GTEx ``SMTSD`` tissue names, GTEx’s own controlled vocabulary.

• neither is ICD-O-3, SNOMED CT or OncoTree. Nothing here maps to a clinical coding system, and a pipeline that needs one has to bring its own crosswalk.

Curated by organ correspondence against the 53 tissues actually present in the reference table, so every value resolves. Ordered best-match first. ``HNSC`` is the weak one and is marked: GTEx has no head-and-neck mucosa, so minor salivary gland and oesophageal mucosa are the nearest epithelia rather than the matched normal. Where a tumour has more than one plausible normal, all of them are listed rather than one being picked silently – SKCM against sun-exposed and sun-protected skin is a different safety question in each.

mhcmatch.expression.TUMOR_TISSUE_APPROXIMATE = ('HNSC',)

Tumour types whose matched normal is an approximation rather than the same organ.

mhcmatch.expression.C_MIN = 0.05

The floor is clamped to this range, in TPM. Measured over all 86 human contexts it runs 0.1000 (whole blood) to 0.4000 (testis), so the clamp is far wider than anything real and exists only so that a degenerate input – an empty context, a table read wrong, a caller passing a filter of zero – cannot produce a floor of 0 and a division that is not one.

The mouse range is not the human one, and it reaches the upper bound. Measured over all 35 FANTOM5 contexts it runs 0.60 (aorta) to 2.00 (pancreas, cecum, os femoris, stomach, testis) – three to five times the human floors, because CAGE tag density is a different measurement reported in the same units. Nothing is clamped today (2.00 is the bound, not past it), but five contexts sit exactly on it, so a redeposit that moves any of them upward would be silently truncated. Compare floors within a species, never across one.

mhcmatch.expression.C_MAX = 2.0

The floor is clamped to this range, in TPM. Measured over all 86 human contexts it runs 0.1000 (whole blood) to 0.4000 (testis), so the clamp is far wider than anything real and exists only so that a degenerate input – an empty context, a table read wrong, a caller passing a filter of zero – cannot produce a floor of 0 and a division that is not one.

The mouse range is not the human one, and it reaches the upper bound. Measured over all 35 FANTOM5 contexts it runs 0.60 (aorta) to 2.00 (pancreas, cecum, os femoris, stomach, testis) – three to five times the human floors, because CAGE tag density is a different measurement reported in the same units. Nothing is clamped today (2.00 is the bound, not past it), but five contexts sit exactly on it, so a redeposit that moves any of them upward would be silently truncated. Compare floors within a species, never across one.

mhcmatch.expression.MIN_SHARED = 1000

A batch scale needs an unconditioned view of the transcriptome, and these two guards are what enforce it. Absolute floor first: fewer shared genes than this is not an estimate.

mhcmatch.expression.MIN_COVERAGE = 0.5

…and then the one that actually matters – the shared genes as a fraction of the genes the reference context has switched on. A TCGA context expresses 23,508 to 30,561 genes (median 26,518), a whole-transcriptome profile covers 0.6 to 0.9 of one, and no candidate list can reach 0.5: the largest screen in the fitting corpus manages 4,772 shared genes, or 0.18. Measured on the three screens carrying their own RNA-seq, a scale taken from their candidate lists comes out at 1.78, 2.18 and 3.15 – all above 1, none of them a unit difference. A mutation reaches a candidate list only if it was seen in RNA, so a candidate’s abundance is conditioned on having been detected, and the ratio measures that conditioning instead of the library. Counting more candidates cannot fix it, which is why the guard is coverage.

mhcmatch.expression.GAMMA_MIN = 1e-06

raw counts against TPM differ by three orders of magnitude or more, and absorbing exactly that is what the estimate is for. These bounds catch a non-finite or structurally broken result, not an unusual unit – a clamp tight enough to call 1000x implausible would reject the commonest case it exists to handle.

Type:

The magnitude a scale is allowed to reach. Wide on purpose

mhcmatch.expression.GAMMA_MAX = 1000000.0

raw counts against TPM differ by three orders of magnitude or more, and absorbing exactly that is what the estimate is for. These bounds catch a non-finite or structurally broken result, not an unusual unit – a clamp tight enough to call 1000x implausible would reject the commonest case it exists to handle.

Type:

The magnitude a scale is allowed to reach. Wide on purpose

mhcmatch.expression.matched_tissues(tumor)

The GTEx tissue(s) that are tumor’s matched normal, best match first.

() for a tumour type with no entry, which is the honest answer – a wrong matched normal turns the safety read into a confident wrong one. See TUMOR_TISSUE for what the two vocabularies are and what they are not.

Parameters:

tumor (str)

Return type:

tuple[str, …]

mhcmatch.expression.resolve_context(text, path=None, approximate=True, detail=False, species='human')

(TCGA study codes, GTEx tissue names) for a free-text origin.

For species="mouse" the first element is always empty – there are no mouse study codes – and the second is the FANTOM5 tissue, matched case- and separator-insensitively. There is no organ back-off there: the 35 contexts are already organ names, so a miss is a miss.

Accepts what a submission actually carries: a study code (SKCM), an organ (liver, Lung), or a GTEx tissue verbatim (Skin - Sun Exposed (Lower Leg)). Matching is case-insensitive throughout, and an organ that is more than one study returns all of them rather than one of them.

An unrecognised string raises. A tumour type that silently became the pooled reference would return a plausible number computed from the wrong distribution, and nothing downstream could tell. An origin that is genuinely unknown is expressed by passing nothing, not by passing a word that does not resolve.

>>> resolve_context("liver")
(('LIHC',), ('Liver',))
>>> resolve_context("lung")
(('LUAD', 'LUSC'), ('Lung',))

Parameters:

• text (str)

• path (str | None)

• approximate (bool)

• detail (bool)

• species (str)

mhcmatch.expression.fetch_reference(path=None, file=None, species='human')

Local path to a reference table, downloading it from HF on first use.

file selects which one – None means species’ own reference (REFERENCE_FILE or REFERENCE_FILE_MOUSE), or name REFERENCE_TOIL_FILE / MATRIX_FILE / FLOORS_FILE / SYNONYMS_FILE directly. Those four are human-only deposits and take no species.

$MHCMATCH_EXPRESSION overrides, for offline and cluster runs, matching how $MHCMATCH_PMHC overrides the presentation panel. Pointed at a directory it resolves file’s basename inside it, so one setting serves every deposit here.

Pointed at a file it overrides REFERENCE_FILE and that alone – whatever the file is called, which is the long-standing contract. The other deposits resolve beside it and fall through to the download when they are not there. One path cannot stand in for four different files, and letting it try is how a caller ends up handing a gzipped TSV to np.load.

Parameters:

• path (str | None)

• file (str | None)

• species (str)

Return type:

str

mhcmatch.expression.fetch_matrix(path=None, species='human')

Local path to species’ dense matrix, downloading it from HF on first use.

Parameters:

• path (str | None)

• species (str)

Return type:

str

mhcmatch.expression.fetch_synonyms(path=None)

Local path to SYNONYMS_FILE, downloading it from HF on first use.

Parameters:

path (str | None)

Return type:

str

mhcmatch.expression.load(path=None, species='human')

{(key_type, key, context): {median_tpm, q25_tpm, q75_tpm, n, source}}, read once and cached.

Read with the stdlib rather than a dataframe library: the package’s runtime dependencies are numpy and huggingface_hub, and a dict keyed by the join tuple is what every caller here wants anyway. species selects the deposit; the parsing is identical because the schema is.

Mouse reads two files into this one table. Its tumour half is a separate deposit (TUMOR_FILE_MOUSE) and is gene-keyed rather than peptide-keyed, so it is folded in under key_type="tumor". Stamping the rows is what keeps every kt == "gene" filter in this module – tissues(), _by_gene(), _tissue_quantile(), and through them safety_profile() – excluding a tumour model by construction rather than by a rule each would have to repeat and one of them would eventually miss. An explicit path reads exactly that file and nothing beside it.

Parameters:

• path (str | None)

• species (str)

Return type:

dict

mhcmatch.expression.lookup(key, tissue=None, tumor=None, path=None, species='human')

Reference expression for a gene in a normal tissue or a peptide in a tumor type.

Exactly one of tissue / tumor is required – they are different measurements in different units from different studies, and silently falling back from one to the other is how a tumour abundance ends up reported as a normal-tissue TPM.

``key`` is a peptide for a human ``tumor=`` and a gene for a mouse one, because the two tumour deposits are keyed differently: TCGA is peptide x study code, and the mouse syngeneic models are gene x model with no peptide rows anywhere. Same rung, different key. An unrecognised mouse model raises rather than returning None, so “not in this deposit” cannot be read as “not expressed in this tumour”.

Parameters:

• key (str)

• tissue (str | None)

• tumor (str | None)

• path (str | None)

• species (str)

Return type:

dict | None

mhcmatch.expression.impute(key, observed=None, tissue=None, tumor=None, path=None, species='human')

(value, was_imputed) – the observed value if there is one, else the reference median.

Never returns “drop this row”. A candidate with no expression measurement still has a source gene whose typical expression in that tissue is known, and the flag lets a model carry a missing-indicator rather than losing the sample.

Parameters:

• key (str)

• observed (float | None)

• tissue (str | None)

• tumor (str | None)

• path (str | None)

• species (str)

Return type:

tuple[float | None, bool]

mhcmatch.expression.tissues(path=None, species='human')

Every normal tissue in species’ reference table.

Parameters:

• path (str | None)

• species (str)

Return type:

list[str]

mhcmatch.expression.tumor_types(path=None, species='human')

Every tumour context in species’ reference table.

Human: the 33 TCGA study codes of the peptide-keyed half (SKCM is melanoma). Mouse: the six syngeneic models of TUMOR_FILE_MOUSE, which are gene-keyed – a different key type answering the same question, which is why both are read out here rather than only the one the human file happens to use.

Parameters:

• path (str | None)

• species (str)

Return type:

list[str]

mhcmatch.expression.safety_profile(gene, top=10, path=None, species='human')

[(tissue, median_tpm)] for a gene across normal tissues, highest first.

The safety read: a target expressed only in the tumour’s lineage is very different from one expressed in heart and lung too, and the ranking score alone does not show that.

Parameters:

• gene (str)

• top (int)

• path (str | None)

• species (str)

Return type:

list[tuple[str, float]]

mhcmatch.expression.context_floor(tumor=None, tissue=None, q=0.25, path=None, clamp=True, prefilter=0.0, detail=False, species='human')

The abundance floor c for mhcmatch.rank.expr_level(), in TPM.

The q-th percentile of non-zero median abundance over every gene in a context, taken from the table where GTEx and TCGA share a pipeline, so a tumour type and a tissue are the same unit and the two can be mixed without a conversion.

Resolution order, and it prefers the tumour deliberately:

tumor  -> that TCGA study's own transcriptome        SKCM 0.1600, LUAD 0.2000 TPM
tissue -> that GTEx tissue's transcriptome           Lung 0.3500, Liver 0.1800 TPM
neither-> the pooled TCGA reference                  0.1800 TPM

A tumour’s floor is not its matched normal’s. Measured across the pairs, a study sits at roughly half the tissue it arose in – SKCM 0.1600 against skin 0.3050, BLCA 0.1700 against bladder 0.3600, LUAD 0.2000 against lung 0.3500. Scoring a tumour candidate against a normal floor puts the whole term about one unit low.

tumor may be a study code or an organ; resolve_context() decides, and an organ that is several studies pools them. prefilter is an expression cut the candidates already passed, in TPM, and raises the floor to meet it – a filter removes the range this term resolves. clamp holds the result inside C_MIN–C_MAX.

With detail=True returns {"floor", "contexts", "n_genes", "pooled", "clamped"}.

``species=”mouse”`` walks the same three rungs off its own two deposits: tumor= is one of the six syngeneic models of TUMOR_FILE_MOUSE, tissue= one of FANTOM5’s 35 adult tissues, and neither pools over the normal reference. contexts comes back as bare context names rather than the human table’s source|context keys, because each rung has one source.

The mouse floors are an order of magnitude above the human ones – 0.60 TPM (aorta) to 2.00 (pancreas, stomach, testis) against 0.10–0.40 – because CAGE tag density concentrates on fewer genes than RSEM TPM does. Compare floors within a species, never across one.

Parameters:

• tumor (str | None)

• tissue (str | None)

• q (float)

• path (str | None)

• clamp (bool)

• prefilter (float)

• detail (bool)

• species (str)

mhcmatch.expression.gene_level(gene, tumor=None, tissue=None, path=None, species='human')

{"tumor", "normal", "pan", "found"} – one gene’s level in TPM, three ways.

All three from the single-pipeline table, so they are directly comparable and their differences mean something:

tumor   the gene's median in that TCGA study, pooled over studies where the origin is
        an organ. ``None`` when no tumour type is given
normal  the gene's median in the matched normal tissue(s). ``None`` when neither a tissue
        nor a tumour with a matched normal is given
pan     the gene's median across the normal tissues where it is transcribed at all, and
        0.0 for a gene silent everywhere. **Always defined**, so a resolution chain has a
        last step that cannot fail

found is whether the gene is in the reference at all. A gene that is not is not a gene with a level of zero, and the two must not be collapsed – one is silence and the other is ignorance.

For species="mouse" all three resolve, off the two deposits rather than one table: tumor from TUMOR_FILE_MOUSE, normal and pan from FANTOM5. They are not in the same unit – one is length-normalised RNA-seq and the other is CAGE tag density – so a mouse tumor/normal difference is not the like-for-like contrast the human single-pipeline table gives. Read them as two terms on their own floors, which is how mhcmatch.rank.expr_level() and mhcmatch.rank.expr_norm_level() use them.

Parameters:

• gene (str)

• tumor (str | None)

• tissue (str | None)

• path (str | None)

• species (str)

Return type:

dict

mhcmatch.expression.batch_scale(values, genes, tumor=None, path=None, detail=False)

(scale, n shared genes, fell back) – what to divide a submitted column by to reach TPM.

A median of ratios against the reference:

scale = median over g of  x(g) / r(g)

over the genes carrying a positive value in both, where r is the tumour type’s own reference level, or the pooled reference where no tumour type is given. It is the standard robust per-sample size factor, and it does what a unit declaration cannot: it absorbs FPKM against TPM, raw counts against TPM, and one pipeline’s TPM against another’s, without the caller having to know which of those they have.

Only genes positive in both count. A candidate at zero is a measurement of silence, not of scale, and letting it in drags the median toward zero in proportion to how many genes the tumour happens to have switched off. A negative value raises. Zero is a measurement and is skipped for a stated reason; a negative one is not a measurement at all, and skipping it too would estimate the scale from whatever part of the column happened to be valid.

Pass a whole-transcriptome profile, not a candidate list, and the guards enforce it: the shared genes must clear MIN_SHARED and cover MIN_COVERAGE of the genes the reference context has switched on. A candidate list cannot, by two-and-a-half fold, and it must not – a mutation is only called where it was seen in RNA, so candidate abundances are conditioned on detection and their ratio to the reference measures that conditioning. On the three screens carrying their own RNA-seq it returns 1.78, 2.18 and 3.15, all above 1, none of them a unit.

Magnitude is deliberately not policed beyond GAMMA_MIN–GAMMA_MAX, which only catches a non-finite or structurally broken result: raw counts sit three or more orders of magnitude from TPM, and refusing that would reject the commonest case this exists to handle.

detail=True adds the spread of the underlying ratios, (q75 - q25) / median, and the coverage the estimate rests on. A pure rescale has spread 0; the measured screens run 2.1 to 3.0. Tumour biology moves individual genes on its own, so spread is reported for a caller to judge rather than gated on a threshold – coverage is the gate.

>>> # a batch that is the reference times 7 recovers 7
>>> batch_scale([7 * v for v in ref], names, tumor="SKCM")
(7.0, 412, False)

Parameters:

• tumor (str | None)

• path (str | None)

• detail (bool)

mhcmatch.expression.coexpression(genes, path=None, absolute=False)

Pairwise similarity of GTEx tissue profile, (n, n) in [0, 1], one row per gene given.

What this is, exactly. Each gene is taken as its vector of median TPM over the 53 GTEx tissues of the Toil recompute (_matrix()), put on log2(1 + TPM), centred and scaled, and correlated with every other. It is therefore similarity of tissue-specificity profile across tissue medians — two genes that are on in the same organs and off in the same organs score high. It is not co-regulation within one tumour, which would need per-sample expression the deposit does not carry, and it should not be described as such.

What it is for. Two cassette units whose source genes are on together are lost together when that transcriptional programme is silenced, and that is a way of failing that no per-unit scalar can express — which is why mhcmatch.cassette.overlap() takes it as a matrix rather than as a feature column.

A gene the matrix does not carry, or one flat across every tissue, gets a zero row: no information about its pairs, which leaves the unit coupled to nothing and ranked on everything else. It is never nan, which would propagate into an argmax and delete the candidate. genes may repeat and may be empty strings; the result is positionally aligned with it.

absolute takes |r|, so a gene pair that is reciprocally regulated counts as related. The default keeps the sign and clips at 0, because the mechanism claimed is on together.

>>> import numpy as np
>>> c = coexpression(["TP53", "MDM2", "TP53"])
>>> float(c[0, 2])                                         # same gene twice
1.0

Parameters:

• path (str | None)

• absolute (bool)

mhcmatch.known module

Built-in known-epitope reference sets for exact-match lookup, assembled from the public deposits: confirmed tumour neoantigens, peptides the screens tested and found negative, IEDB-immunogenic epitopes, the thymic self-immunopeptidome and the viral ligandome. An exact match is stronger evidence than any model output, so mhcmatch.rank reports it as a flag and never folds it into the score.

Known-epitope reference sets, assembled from the public deposits, for exact-match lookup.

An exact match is stronger evidence than any model output, so it is a flag and never a score. mhcmatch.rank reports which set a candidate was found in and floats those candidates into a tier of their own, with the model score still shown beside them. Burying “this peptide is a confirmed NCI neoantigen” inside a weighted sum lets a mediocre model score dilute the one piece of direct evidence in the row.

Five sets, each answering a different question about a candidate:

set	a hit means
neoantigen	a confirmed immunogenic tumour neoantigen – the NCI/Gartner deconvolved minimal peptides, the epitope-resolution screens (NCI, HiTIDE, TESLA) and the aggregated cohorts, all restricted to rows the assay called positive. This candidate has already been shown to work.
neoantigen_neg	screened and found non-immunogenic. The same deposits, negative rows. Not evidence of nothing: it is the one label that says this exact peptide was tested and did not respond.
immunogenic	a positive T-cell assay in IEDB against any source – viral, bacterial, tumour. Immunogenic in some context, not necessarily this one.
self	present in the thymic self-immunopeptidome (HLA Ligand Atlas). The tolerance argument, and a cross-reactivity/autoimmunity flag for a vaccine: reactive T cells were plausibly deleted.
viral	a presented pathogen-derived peptide. A pre-existing anti-pathogen repertoire may cross-react.

neoantigen and neoantigen_neg are disjoint by construction: a peptide reported positive anywhere is a positive, because a negative call in one assay does not overturn a positive one in another. That rule is applied after pooling, in load().

Every set is fetched from the public HF dataset (mhcmatch.store.PMHC_REPO) and cached, so a fresh install needs no pre-staged data. MHCMATCH_PMHC_DIR points at a local mirror instead.

>>> from mhcmatch import known
>>> refs = known.load()
>>> "GILGFVFTL" in refs["viral"]
True

mhcmatch.known.SOURCES: dict[str, list[tuple]] = {'immunogenic': [('immunogenicity/chowell_rebuilt.tsv.gz', 'label', {'1'})], 'neoantigen': [('neoantigens/nci_gartner_mmp.tsv.gz', 'immunogenicity', {'1', 'cd8', 'positive'}), ('neoantigens/neoantigens_tested_peptides.tsv.gz', 'immunogenicity', {'1', 'cd8', 'positive'}), ('neoantigens/neoag_tested.tsv.gz', 'immunogenicity', {'1', 'cd8', 'positive'}), ('neoantigens/neoag_tested_hsa.tsv.gz', 'immunogenicity', {'1', 'cd8', 'positive'}), ('neoantigens/neoag_tested_mmu.tsv.gz', 'immunogenicity', {'1', 'cd8', 'positive'})], 'neoantigen_neg': [('neoantigens/neoantigens_tested_peptides.tsv.gz', 'immunogenicity', {'0', 'negative', 'non-immunogenic'}), ('neoantigens/neoag_tested.tsv.gz', 'immunogenicity', {'0', 'negative'}), ('neoantigens/neoag_tested_hsa.tsv.gz', 'immunogenicity', {'0', 'negative'}), ('neoantigens/neoag_tested_mmu.tsv.gz', 'immunogenicity', {'0', 'negative'})], 'self': [('thymus/thymus_immunopeptidome.tsv.gz', None, None)], 'viral': [('ligandome/viral_foreign_iedb.tsv.gz', None, None)]}

set -> [(repo-relative file, label column or None, values that count as a hit)]. A None label column means every row of the file belongs to the set.

mhcmatch.known.SET_NAMES = ('neoantigen', 'neoantigen_neg', 'immunogenic', 'self', 'viral')

strongest direct evidence first, so lookup() names the best one.

Type:

Report order

mhcmatch.known.load(names=None)

{set name: frozenset of peptides}, downloaded on first use and cached.

names restricts which sets are built – each one costs a download and a full-file scan, so ask for what you need. Default is all of SET_NAMES.

A peptide reported immunogenic anywhere is removed from neoantigen_neg: one assay calling it negative does not overturn another calling it positive, and leaving it in both would make the flag report whichever set happened to be checked first.

Parameters:

names (tuple[str, ...] | None)

Return type:

dict[str, frozenset]

mhcmatch.known.lookup(peptide, refs=None)

The first set in SET_NAMES containing this peptide exactly, else "".

Order is deliberate: a peptide that is both a confirmed neoantigen and a thymic self-peptide should report as the neoantigen, because that is the stronger and more actionable evidence – the self hit is still visible by checking refs directly.

Parameters:

• peptide (str)

• refs (dict | None)

Return type:

str