The recognition axis, reduced ============================= The recognition score :mod:`mhcmatch.complement` ships as a thirty-column linear head over six feature blocks, backed by 660 fitted log-odds cells (:doc:`complementarity`). 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`` (:doc:`corpus`), 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 :ref:`burial-repro`. .. note:: Nothing on this page changed a default of :func:`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: :func:`mhcmatch.complement.burial` is the ``C_phys`` term of the shipped aggregate (:doc:`neoantigen`). The two halves -------------- The block partition splits cleanly into chemistry and fitted residue identity: .. list-table:: :header-rows: 1 :widths: 12 40 8 40 * - 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: 1. **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. (:doc:`property_basis` records the companion fact that the Kidera table is already orthogonal, so PCA *on it* is degenerate.) 2. **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. 3. **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: .. list-table:: :header-rows: 1 :widths: 6 44 12 12 12 14 * - 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. .. code-block:: python 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. .. list-table:: correlation with per-peptide cysteine count :header-rows: 1 :widths: 55 20 * - 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: .. list-table:: :header-rows: 1 :widths: 40 30 * - 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 :func:`mhcmatch.mimicry.corpus_R` makes with its per-window divisor (:doc:`corpus`), 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: .. list-table:: :header-rows: 1 :widths: 46 27 27 * - 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: .. list-table:: :header-rows: 1 :widths: 40 30 30 * - - 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 :data:`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 :data:`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. :func:`mhcmatch.complement.score` grew ``mask_cys=`` for the same reason — :mod:`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: .. list-table:: :header-rows: 1 :widths: 8 34 16 12 * - 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 :data:`mhcmatch.complement.BLOCKS`. Exact partial sum, no refit. ``score(..., mask_cys=True)`` Zero cysteine in every fitted log-odds table, as :mod:`mhcmatch.posbayes` does. ``score(..., positions="profile")`` Read the TCR-facing chemistry over the **positional contact profile** (:func:`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 :data:`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. .. _burial-repro: Reproducing all of it --------------------- Scripts live in ``2026-mhcmatch-code`` (private; released with the manuscript) and write to ``bench/results/``: .. list-table:: :header-rows: 1 :widths: 46 54 * - 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 -------------------------------------- :func:`mhcmatch.complement.burial` returns this column directly: .. code-block:: python 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 :data:`mhcmatch.data.aa_tables.HYDROPHOBICITY`, a ``"FAMILY:COMPONENT"`` key into :data:`~mhcmatch.data.aa_tables.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.