The recognition axis, reduced#

The recognition score mhcmatch.complement ships as a thirty-column linear head over six feature blocks, backed by 660 fitted log-odds cells (Complementarity: the recognition axis). This page records what happened when its two halves were separated and each was asked to justify its parameters against the neoantigen screens, and what the shipped API now exposes so that the reduction is reproducible from a fresh install. The chemistry half is C_phys, one of the two Complementarity factors the shipped aggregate carries; the other is C_corpus (Corpus complementarity: what the repertoire was shaped by), which is what replaces the identity half argued against below.

The full derivation, with the physics and the literature, is in the manuscript’s Methods and Supplementary Note 1. Every number here is regenerated by a named script in the benchmark repository — see Reproducing all of it.

Note

Nothing on this page changed a default of mhcmatch.complement.score(), which scores exactly as it did, and the score() parameters described below are opt-in. The reduced form itself is not only a measurement: mhcmatch.complement.burial() is the C_phys term of the shipped aggregate (Ranking neoantigens).

The two halves#

The block partition splits cleanly into chemistry and fitted residue identity:

half	blocks	cols	what it is
`C_phys`	`phys`, `role`, `pot`, `motif`	16	property scales read over the anchor and TCR faces. No fitted table.
`C_aa`	`aa`, `kmer`	14	per-role residue log-odds, per length bin and TCR-face third, plus adjacent dipeptides.

Both are exact partial sums of the shipped score, so they add back to it up to the constant that belongs to no block:

from mhcmatch import complement as C
import math, numpy as np

peps = ["GILGFVFTL", "NLVPMVATV", "KRWIILGLNK"]
t = C.table("human", "mhc1")
off = t["logistic"]["intercept"] - math.log(t["prevalence"] / (1 - t["prevalence"]))

phys = C.score(peps, blocks=("phys", "role", "pot", "motif"))
aa = C.score(peps, blocks=("aa", "kmer"))
np.allclose(phys + aa + off, C.score(peps))     # True, to 3.6e-15

The identity is asserted in tests/test_complement.py for both classes, because every comparison between the halves rests on it.

Warning

Whole blocks only. aa_tcr is fitted at -0.8796 while every per-length aa_tcr{8..11} and every aa_tcr_{n,m,c} is positive — the pooled column is collinear with its own decompositions and the fit pushed it negative to compensate. Dropping columns within a block without refitting inverts the dominant term. Dropping a whole block cannot, which is why blocks= takes block names and not column names.

Why the sixteen chemistry columns are not interpretable#

Three properties of the design, none of them properties of the data:

The property PCA re-derives the Kidera basis. Over the twenty residues r(PC1, KF4) = -0.9238 and r(PC2, KF2) = -0.9906; regressed on all ten Kidera factors, PC1 reaches R² = 0.9901 and PC2 R² = 0.9931. So the physical axes carried are hydropathy and side-chain bulk — each present twice. (The amino-acid property basis records the companion fact that the Kidera table is already orthogonal, so PCA on it is degenerate.)
Two columns are exact linear combinations of two others. The roles partition the peptide, so pc1 = pc1_anchor + pc1_tcr identically, and likewise for pc2. Two eigenvalues of the 16×16 correlation matrix are exactly zero and seven components carry 95 % of the variance.
One axis is written three times. On 20,000 corpus peptides, r(pc1_anchor, kf4_anchor) = -0.9322, r(pc1_anchor, mj_anchor) = +0.9409, r(pc1_tcr, kf4_tcr) = -0.9152.

The visible consequence is that the block enters the general model at z = +0.18 and reverses to -0.86 beside the identity half. That is a rank-deficient block whose ridge penalty is choosing how to spread one signal across collinear columns — not an uninformative one. The same block scores 0.6766 AUROC on the GBM cohort alone, above the full shipped score’s 0.6186.

What survives conditional selection#

Selecting features by how well they separate alone is what produced three copies of one axis. The benchmark instead scored all 576 candidate columns — every vendored residue vector over {anchor, TCR} × {sum, mean} — by the BIC change each makes inside the general model, which already carries binding, occupancy, expression, the mimicry channels and the identity half:

step	candidate	χ²	BIC	ΔBIC	outcome
—	baseline	—	4208.7	—	—
1	Rose burial propensity over the TCR face (selected as a sum; shipped as the per-residue mean — see below)	20.0	4201.4	−7.3	admitted
2	SVGER component 3, mean over the TCR face	12.1	4202.0	+0.6	rejected

Exactly one column survives. Step 2 is why the criterion is BIC and not accuracy: admitting it would have raised the within-screen median AUROC from 0.6524 to 0.6909 while worsening BIC — that is selection overfitting on 1,101 positives, caught.

The surviving scale is not a hydrophobicity scale#

HYDROPHOBICITY["Rose"] is the scale of Rose et al. (Science 1985, doi:10.1126/science.4023714): the mean fraction of a residue’s solvent-accessible surface area that is lost on folding, measured from proteins of known structure. It is an observable of geometry, where Kyte–Doolittle and the octanol scales are transfer free energies from partition experiments. They correlate — r(Rose, KyteDoolittle) = +0.841, r(Rose, KF4) = -0.849 — but they are different measurements.

Summed over the TCR face, it scores the surface area those residues would bury if packed into an interior, which is the currency of an interface that the structural literature describes as a packing problem with poor shape complementarity. There is also direct precedent for a burial propensity beating classical hydrophobicity scales at far lower parameter cost: Zhou & Zhou (Protein Science 2003, doi:10.1110/ps.0305103) matched 100–200-parameter statistical models with a 24-parameter burial scale on transmembrane topology.

import numpy as np
from mhcmatch import complement as C
from mhcmatch.data import aa_tables as T

AA = C.AA
rose = np.array([T.HYDROPHOBICITY["Rose"][a] for a in AA])
_, counts = C.encode(["GILGFVFTL", "NLVPMVATV"], "mhc1")
c_phys = counts["tcr"].astype(float) @ rose        # the selected feature
rose.min(), rose.max()                              # 0.520, 0.910 -- bounded, by construction

It cannot memorise the corpus artefact#

The Chowell construction carries a cysteine gradient — assayed positives against mass-spectrometry-eluted negatives, and MS under-recovers free cysteine — so the residue marks the assay platform. Fitted freely it takes the largest cell in every shipped residue table.

An imported basis cannot express that, and the mechanism is that ρ is a fraction of area and therefore bounded on [0.520, 0.910]: cysteine at the top sits 3.4 % above isoleucine, where a fitted log-odds cell is unbounded and reaches +2.6 nats.

correlation with per-peptide cysteine count#
term	r
shipped `complement` (30 columns)	+0.6881
`C_aa` (identity half)	+0.6809
`C_aa` with `mask_cys=True`	+0.0371
`C_phys_rose` (Rose over the TCR face)	+0.1083

Per residue, not summed — and the difference was 91 % of the variance#

The scale was selected as a sum over the TCR face, and that made it a ruler rather than a chemistry term. The class-I face is L − 5 residues wide and the Rose scale is strictly positive (0.52 to 0.91), so summing it gives roughly 0.75 (L − 5). Measured on 60,000 fit-corpus peptides:

column	Pearson with peptide length
Rose, summed over the face	+0.954
Rose, averaged over the face (shipped)	−0.010
Kidera KF4, summed	+0.052
Kidera KF4, averaged	+0.005

A centred scale like KF4 does not have the problem, which is why the two were not comparable and why the summed Rose column beat everything it was tried against. Dividing by the face width is the same correction mhcmatch.mimicry.corpus_R() makes with its per-window divisor (Corpus complementarity: what the repertoire was shaped by), for the same reason, and it changes what the pair measures: on the summed scale Rose and KF4 correlate −0.20, on the averaged scale −0.836. They were always close to one axis; the length variance was hiding it.

burial(..., per_residue=False) is the summed form, kept so the measurement reproduces.

Two scales, carried together#

The chemistry block is two columns, never one. The first is always burial. The second is the question this section answers, because the obvious choice turns out not to be a second measurement at all.

Burial and hydropathy are one axis#

EPIC v3 shipped C_phys_rose and C_phys_hydrop (Kidera KF4). Rose measures how much surface a residue buries on folding; KF4 measures how it partitions between water and oil. Those are different physical questions, but on real peptides they are not different numbers:

correlation of Rose with Kidera KF4	value	over
per residue, averaged scale	−0.849	the 20 residue types
per peptide, TCR-face mean	−0.837	the 354,909-row fit population of the day

The signature of that in a fit is a Wald test and a likelihood-ratio test disagreeing about the same block: χ² = 3.18, p = 0.204 against LR χ² = 11.0, p = 4.0×10⁻³. The block matters; neither column can be said to be the one that matters. Rotating the pair into principal components does not help — PC2 lands at AUROC 0.5014 / 0.5085, i.e. chance — because a rotation of two collinear columns still spans one direction.

The corroboration is that the fitted basis had already said so. The retired 16-column chemistry head’s own PC1 reproduces KF4 at r = −0.9238: a predictor fitted freely on this data rediscovers hydropathy as its dominant axis, and burial is that axis measured a different way.

No hydropathy scale escapes it, and they carry the cysteine artefact#

Swept over 141 complete residue scales on the HLA-matched Chowell rebuild, every transfer-free- energy scale sits between |r| = 0.74 and 0.95 of Rose. Ranked on raw AUROC the winner is whichever one is closest to burial — Casari at 0.6008, r = +0.9524 with Rose. Ranked on the residual after Rose, which is the quantity that matters, nothing beats KF4. There is no better hydropathy scale; there is no hydropathy scale that is a second axis.

There is also a reason not to want one here. The Chowell family carries a 12.5× cysteine enrichment among immunogenic peptides which does not appear on the neoantigen screens, so a scale that loads on cysteine imports it. Every candidate scale reports its cysteine loading, and the spread is the point: KF4 loads −0.1033, Casari +0.1861.

The second axis is charge, and what it buys is burial’s stability#

The column being fixed here is burial, not the one being swapped in. Beside a collinear partner burial’s coefficient is inflated and its standard error with it; beside an orthogonal one both come down and the evidence goes up. AF5’s own coefficient in the model is small and its own p is not the statistic to read on it.

What AF5 has to be is two things, and it is both. Orthogonal to burial — r = +0.008 per peptide, where every one of the 39 transfer scales swept sits at 0.74 to 0.95. And carrying its own signal, so that the second slot holds a measurement rather than noise: on the selection corpora it is above chance everywhere (0.528 / 0.550 / 0.536 across the Chowell arms) and it is the strongest of the three chemistry columns on kesmir_vanilla mouse at 0.576, against Rose’s 0.519 and KF4’s 0.481 — the corpus where the contrast is not confounded by self versus non-self. An orthogonal axis that were merely noise would stabilise burial too, and would deserve none of the slot.

C_phys_charge is Atchley AF5, the fifth factor of Atchley et al. (PNAS 2005, doi:10.1073/pnas.0408677102) — electrostatic charge. It is not a hydropathy scale, which is exactly why it works:

	Rose + Kidera KF4	Rose + Atchley AF5
r between the two columns, per peptide	−0.837	+0.008
cysteine loading of the second scale	−0.1033	−0.0028 (lowest of the 141 swept)
burial’s coefficient	+0.1491	+0.1146
burial’s bootstrap sd	0.0874	0.0487
burial’s z / p	+1.71 / 0.088	+2.34 / 0.020
burial’s sign stability	96.5 %	100 %
second column’s sign stability	69.5 %	90 %
BIC	4181.7	4180.3
leave-one-screen-out mean / median	0.6661 / 0.6418	0.6688 / 0.6497
CV, grouped on peptide / twin	0.6469 / 0.6418	0.6566 / 0.6497
held-out verdict	0i/6t/1r	2i/5t/0r

Burial’s coefficient gets smaller and its evidence gets stronger, which is what removing a collinear partner does — the standard error halves, and the collinearity is gone. The pair is orthogonal without any rotation, so the two columns keep their physical names.

Rose + AF5 is what ships. The arm above was measured on the 354,909 rows / 958 positives over nine screens the corpus held at the time, 400 cluster bootstraps; the shipped fit has moved since and mhcmatch.rank.PHYS_COLUMNS names the two scales it currently carries. KF4 itself is not retired — complement.burial(peptides, scale="KIDERA:KF4") computes it, as does any of the other vendored residue vectors, so every number recorded under the KF4 arm keeps its meaning and the comparison stays runnable. What changed is which two the model fits.

Charge is also the axis a structure argues for. In the MART-1 decamer ELAGIGILTV the P1 glutamate contacts germline TCR α CDR1 in eight of nine native complexes — a positive Arg partner in three at close range, down to 2.71 Å in 3QDM, and a polar Gln partner in six (3.05–4.25 Å). Several other charged epitopes in the same set behave the same way. A charged TCR-facing residue is read as charge, by a germline element, and burial has nothing to say about it.

What it costs, recorded#

Charge transfers across species worse than burial does: human→mouse 0.554 and mouse→human 0.539, against Rose’s 0.581 and 0.579. And on the Chowell family it is the weaker of the two marginals (0.528 / 0.550 / 0.536 against Rose’s 0.641 / 0.579 / 0.618), which is what a minor axis looks like next to the major one. Recorded as measured.

Which ships#

The vendored artifact carries Rose + AF5 (C_phys_buried and C_phys_charge). The pair was selected on BIC 4173.3 → 4172.4, leave-one-screen-out mean 0.6583 → 0.6602 and median 0.6309 → 0.6385, and CV 0.6307 → 0.6386 grouped on peptide and 0.6309 → 0.6385 on twin group, against the KF4 arm on identical rows at equal parameter count — and, the reason for the swap, on burial going from p = 0.064 at 98 % sign stability to p = 0.0202 at 100 %.

The library computes the two keys the shipped model fits — C_phys_buried and C_phys_charge — from mhcmatch.rank.PHYS_COLUMNS, and aggregate_score reads only the names the artifact’s own features list asks for. The v3 pair (C_phys_rose, C_phys_hydrop) is gone: nothing shipped read it, and computing two matrix products per peptide for a model that is not shipped is not free.

Both columns of the pair in force are emitted, so the block is auditable rather than a single number with a footnote.

mhcmatch.complement.score() grew mask_cys= for the same reason — mhcmatch.posbayes masks cysteine by construction and this module never did. It is off by default, so no recorded number moves:

C.score("GILGFCFTL")[0] - C.score("GILGFVFTL")[0]                  # +2.8971
C.score("GILGFCFTL", mask_cys=True)[0] - C.score("GILGFVFTL", mask_cys=True)[0]   # -0.4718

Masking reproduces the posbayes construction rather than approximating it: over the 19 non-cysteine residues the masked tables differ from the independently deposited posbayes ones by mean 0.00000, sd 1.3e-04 (human).

Where the halves rank#

Fitted together in the general model, standardized coefficients:

rank	term	coefficient	z
1	expression	+0.3368	+6.20
2	C_phys	+0.3091	+4.57
3	self mimicry (TCR face)	+0.2656	+2.85
5	binder score	+0.1448	+4.09
7	C_aa	+0.0994	+1.97

They are not competing for the same variance — r(C_phys, C_aa) = +0.1430 — so the identity half loses its slot to the parameter penalty, not to collinearity. On BIC the chemistry term alone is the best design tested (4195.5 against the incumbent’s 4209.8, within-screen median 0.6511).

Other parameters this work added#

All opt-in, all defaulting to the shipped behaviour:

score(..., blocks=(...)): Restrict to a subset of mhcmatch.complement.BLOCKS. Exact partial sum, no refit.
score(..., mask_cys=True): Zero cysteine in every fitted log-odds table, as mhcmatch.posbayes does.
score(..., positions="profile"): Read the TCR-facing chemistry over the positional contact profile (mhcmatch.immuno.contact_profile(), per-position TCR↔peptide contact frequency from 8,062 contacts over 370 crystals) instead of the binary anchor mask. Touches role/pot/motif only; aa/kmer come back bit-identical.
score(..., paratope="contact"): Read the pot block off mhcmatch.complement.PARATOPE_CONTACT — the TCRen potential marginalised over the receptor residues that actually contact peptide, P(aa | contact = 1), rather than over the flat composition of the whole CDR3 loop. Only 35.4 % of TRB loop residues ever contact, and the germline flanks carrying most of the flat mass sit at P(contact) = 0, so conditioning on contact is the flank trim.

positions and paratope are independent: one chooses which peptide positions are read, the other which receptor residues the potential was averaged over.

Reproducing all of it#

Scripts live in 2026-mhcmatch-code (private; released with the manuscript) and write to bench/results/:

claim	script → result
PCA re-derives Kidera; rank deficiency	`bench/immuno/basis_redundancy.py` → `basis_redundancy.md`
576-candidate conditional selection	`bench/immuno/physchem_select.py` → `physchem_select.md`
identity-half structure ablation	`bench/immuno/aa_structure.py` → `aa_structure.md`
EM mixture over residue tables	`bench/immuno/aa_mixture.py` → `aa_mixture.md`
final ranks, BIC, cysteine loadings	`bench/immuno/reduced_recognition.py` → `reduced_recognition.md`
per-arm standalone AUROC, 21 arms	`bench/immuno/recognition_axes.py` → `recognition_axes.md`
the 84-cell basis × scheme sweep	`bench/immuno/physchem_families.py` → `physchem_families.md`
cysteine in every vendored table	`bench/immuno/artifact_composition_audit.py` → `artifact_composition.md`

Scope#

The selection ran on MHC class I, human; the mouse arm was used only for the identity-half ablation and the reduced form is not yet measured on class II. The chemistry column was chosen from 576 candidates on 1,101 positives — BIC is the guard, but a prospective cohort should confirm it. And the identity half is rejected by the parameter penalty in this model, not shown uninformative: it holds the third-largest coefficient when fitted without the chemistry term, and it carries the Kešmir-construction arms where the corpora invert.

Recomputing it, and swapping the basis#

mhcmatch.complement.burial() returns this column directly:

from mhcmatch import complement
complement.burial(["GILGFVFTL"])                       # the shipped Rose basis
complement.burial(["GILGFVFTL"], scale="KIDERA:KF4")   # exploration only

scale= accepts any of the 45 keys of mhcmatch.data.aa_tables.HYDROPHOBICITY, a "FAMILY:COMPONENT" key into DESCRIPTORS ("KIDERA:KF4", "CRUCIANI:PP1", …), or a dict over the twenty residues, and raises on an unknown name rather than guessing.

It exists for comparison, not for scoring. "Rose" was selected out of 576 candidates by the BIC change it produced inside the general model; a column computed on another basis re-parameterises that result and must be reported as a comparison. The Kidera factors are the ones worth quoting, since the Chowell-family literature is usually written against them — and they lose to "Rose" on the neoantigen corpus.