From 79479108a215348855a6895059499e58bcdb5900 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Fri, 10 Jul 2026 10:36:44 -0400 Subject: [PATCH] Consolidate predictor wrappers + BigMHC consistency fixes + docstring cleanup Follow-up to the code-review correctness PR (#247). Structural + cosmetic: net -129 lines in the wrapper files with no behavior change (except the BigMHC fixes below). Consolidation (new mhctools/wrapper_base.py): - NewModelPredictorMixin holds the members that were copy-pasted across the standalone wrappers: predict_dataframe, supported_kinds, __repr__, and peptide normalization (single-string + strip/upper). - AlleleFreePredictor adds the allele-independent kind_support shape. - Calis / DeepTAP / ERAMER now subclass AlleleFreePredictor; BigMHC uses NewModelPredictorMixin. Each keeps only its constructor, _default_pred_kind, tool-specific validation (bounds/messages differ), and scoring. BigMHC consistency fixes (were flagged in review): - Build predictions from the *normalized* peptide, not the raw input, so "siinfekl " is no longer scored as SIINFEKL but returned verbatim. - Rename _pred_kind -> _default_pred_kind to match every other wrapper (generic consumers can now call it uniformly). Minor code: - eramer: early-return on empty input (was loading the whole PWM first); drop the unreachable score-is-None branch; delete the dead _predictor_name. - calis/deeptap: delete the dead _predictor_name. - pred: add reduce_op(kind, field) and use it in PeptideResult.best_by and AnnotationSpec.direction_op (the max/min mapping was written twice). - mixmhc2pred: resolve the %Rank_/Score_ column positions once, not per row. Docs / LLMism cleanup: - Trim the copy-pasted "Note on interpretation" caveats to one sentence each (calis/prime/eramer/deeptap); drop PRIME's benchmark editorializing and the marketing/hedging phrasing; remove **bold** from reST docstrings. - Fix stale docstrings: calis position_weights ("leading slice"), eramer's class docstring (add $ERAMER_PWM to the resolution order). - Trim the bloated best_direction docstring, filler class docstrings, an over-defensive annotate comment, and the verbose nettcr suppression docstring. Tests: new tests/test_wrapper_base.py covers the shared base and asserts the wrappers inherit it (and BigMHC standardizes on _default_pred_kind). Full public subset: 496 passed. Version 3.28.0 -> 3.29.0. Claude-Session: https://claude.ai/code/session_01LZahFhBSCiehXTESCYQ7wG --- .github/workflows/tests.yml | 1 + mhctools/__init__.py | 2 +- mhctools/annotate.py | 15 +++---- mhctools/bigmhc.py | 27 +++--------- mhctools/calis.py | 85 +++++++++++------------------------- mhctools/deeptap.py | 56 ++++++------------------ mhctools/eramer.py | 86 +++++++++++++------------------------ mhctools/mixmhc2pred.py | 24 ++++++----- mhctools/nettcr.py | 21 +++------ mhctools/pred.py | 42 ++++++------------ mhctools/prime.py | 12 ++---- mhctools/wrapper_base.py | 76 ++++++++++++++++++++++++++++++++ tests/test_wrapper_base.py | 68 +++++++++++++++++++++++++++++ 13 files changed, 265 insertions(+), 250 deletions(-) create mode 100644 mhctools/wrapper_base.py create mode 100644 tests/test_wrapper_base.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 9790a99..c1750ac 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -58,6 +58,7 @@ jobs: tests/test_deeptap.py \ tests/test_calis.py \ tests/test_eramer.py \ + tests/test_wrapper_base.py \ tests/test_processing_predictor.py \ tests/test_pepsickle.py \ tests/test_bigmhc.py \ diff --git a/mhctools/__init__.py b/mhctools/__init__.py index bf76f56..bd84e3e 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -85,7 +85,7 @@ def __getattr__(name): raise AttributeError( "module %r has no attribute %r" % (__name__, name)) -__version__ = "3.28.0" +__version__ = "3.29.0" __all__ = [ "Prediction", diff --git a/mhctools/annotate.py b/mhctools/annotate.py index c2875f6..8236a34 100644 --- a/mhctools/annotate.py +++ b/mhctools/annotate.py @@ -48,7 +48,7 @@ import pandas as pd from .allele_normalization import normalize_allele_name_or_raw -from .pred import Kind, best_direction +from .pred import Kind, reduce_op # Whitespace, comma, or semicolon separate multiple alleles packed into one @@ -140,8 +140,7 @@ def build_predictor(self, alleles): def direction_op(self): """``max`` or ``min`` callable for reducing candidate predictions.""" - return max if best_direction(self.kind, self.prediction_field) == "max" \ - else min + return reduce_op(self.kind, self.prediction_field) def parse_annotation_spec(token): @@ -338,7 +337,7 @@ def annotate_table( predictor = spec.build_predictor(union_alleles or None) results = predictor.predict(union_peptides) by_pair, by_peptide = _build_lookup(results, spec) - reduce_op = spec.direction_op() + reducer = spec.direction_op() field = spec.prediction_field values = [] @@ -347,15 +346,13 @@ def annotate_table( candidates = [by_pair[(peptide, a)] for a in alleles if (peptide, a) in by_pair] if not candidates: - # Allele-free predictors (e.g. processing) emit allele-less - # predictions indexed by peptide only. Only such predictors - # populate ``by_peptide``, so this fallback never masks a - # genuine unsupported-allele miss from a binding predictor. + # Allele-free predictors (e.g. processing) index by peptide + # only; fall back to that lookup. allele_free = by_peptide.get(peptide) if allele_free is not None: candidates = [allele_free] if candidates: - best = reduce_op(candidates, key=lambda p: getattr(p, field)) + best = reducer(candidates, key=lambda p: getattr(p, field)) values.append(getattr(best, field)) best_alleles.append(best.allele or None) else: diff --git a/mhctools/bigmhc.py b/mhctools/bigmhc.py index cdfe027..7d2cf1f 100644 --- a/mhctools/bigmhc.py +++ b/mhctools/bigmhc.py @@ -32,11 +32,11 @@ import torch from .pred import ( - COLUMNS, Kind, PeptideResult, Prediction, ) +from .wrapper_base import NewModelPredictorMixin def _find_bigmhc_dir(bigmhc_path=None): @@ -79,7 +79,7 @@ def _import_bigmhc_modules(bigmhc_dir): sys.path.remove(src_dir) -class BigMHC: +class BigMHC(NewModelPredictorMixin): """Wrapper for BigMHC presentation and immunogenicity predictions. Models are loaded lazily on the first call to :meth:`predict` and @@ -130,26 +130,19 @@ def __str__(self): return "BigMHC(mode=%s, alleles=%s, %s)" % ( self.mode, self.alleles, loaded) - def __repr__(self): - return str(self) - - def _pred_kind(self): + def _default_pred_kind(self): if self.mode == "im": return Kind.immunogenicity return Kind.pMHC_presentation def kind_support(self): return { - self._pred_kind(): { + self._default_pred_kind(): { "mhc_dependence": "single_allele", "mhc_class": "I", } } - @property - def supported_kinds(self): - return tuple(self.kind_support()) - def _predictor_name(self): return "bigmhc_%s" % self.mode @@ -256,8 +249,7 @@ def predict(self, peptides): list of PeptideResult One entry per peptide; each contains one Prediction per allele. """ - if isinstance(peptides, str): - peptides = [peptides] + peptides = self._normalize_peptides(peptides) # Build all (peptide, allele) combinations all_peptides = [] @@ -269,7 +261,7 @@ def predict(self, peptides): scores = self._predict_raw(all_peptides, all_alleles) - kind = self._pred_kind() + kind = self._default_pred_kind() name = self._predictor_name() idx = 0 @@ -288,13 +280,6 @@ def predict(self, peptides): results.append(PeptideResult(preds=tuple(preds))) return results - def predict_dataframe(self, peptides, sample_name=""): - """``predict()`` flattened to a DataFrame.""" - dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] - if not dfs: - return pd.DataFrame(columns=COLUMNS) - return pd.concat(dfs, ignore_index=True) - class BigMHC_EL(BigMHC): """BigMHC in presentation (eluted-ligand) mode.""" diff --git a/mhctools/calis.py b/mhctools/calis.py index db2ed0a..3f2fc84 100644 --- a/mhctools/calis.py +++ b/mhctools/calis.py @@ -19,28 +19,19 @@ importance, with the anchor positions (P1, P2, and the C-terminus) masked out because their identity is driven by MHC binding rather than TCR recognition. -Unlike every other predictor in mhctools, this needs **no external install and -no downloaded weights** — the ~30 published parameters are reproduced here -directly (they come from the open-access CC-BY paper, Calis et al., -*PLoS Comput. Biol.* 2013, ``10.1371/journal.pcbi.1003266``). It is a fast, -dependency-free baseline that sits alongside the learned immunogenicity -predictors (``BigMHC_IM``, ``PRIME``). - -The model is **allele-independent** here: it always masks P1/P2/C-terminus, the -default the IEDB tool uses when no allele is supplied. (The IEDB tool can swap -in allele-specific anchor masks for a fixed ~50-allele catalogue, but those -anchor sets almost all reduce to the P1/P2/C-terminus default anyway.) - -Note on interpretation: like every current CD8 immunogenicity predictor, Calis -is useful for ranking but generalizes weakly — independent benchmarks put the -whole field near AUC 0.5-0.65 on unseen neoepitopes. Treat scores as a -prioritization aid, not ground truth. A score > 0 leans immunogenic, < 0 leans -non-immunogenic. +The ~30 parameters come from the open-access paper (Calis et al., +PLoS Comput. Biol. 2013, 10.1371/journal.pcbi.1003266) and are reproduced here, +so the predictor is self-contained. It is allele-independent: it always masks +P1/P2/C-terminus, the default the IEDB tool uses when no allele is given. A +score > 0 leans immunogenic, < 0 non-immunogenic. + +Like every current CD8 immunogenicity predictor it ranks better than it +generalizes (independent benchmarks put the field near AUC 0.5-0.65 on unseen +neoepitopes) — a prioritization aid, not ground truth. """ -import pandas as pd - -from .pred import COLUMNS, Kind, PeptideResult, Prediction +from .pred import Kind, PeptideResult, Prediction +from .wrapper_base import AlleleFreePredictor # Per-amino-acid log-enrichment score, log(freq(aa | immunogenic) / # freq(aa | non-immunogenic)) from the Calis 2013 training set. Higher = more @@ -67,7 +58,8 @@ def position_weights(peptide_length): For 9-mers this is :data:`IMMUNOWEIGHT` verbatim. For longer peptides the Calis rule pads the interior (after P5) with ``0.30`` weights, matching the - IEDB tool. Shorter peptides use the leading slice of :data:`IMMUNOWEIGHT`. + IEDB tool. For shorter peptides it returns ``IMMUNOWEIGHT`` unchanged and the + caller uses only the leading ``peptide_length`` entries. """ if peptide_length > 9: return ( @@ -104,12 +96,12 @@ def immunogenicity_score(peptide): return round(score, 5) -class Calis: +class Calis(AlleleFreePredictor): """The Calis (IEDB) class-I immunogenicity predictor. - Self-contained (no external tool, no downloaded weights). Allele-independent: - ``predict()`` returns one :class:`~mhctools.pred.Prediction` per peptide with - an empty ``allele`` and ``kind == Kind.immunogenicity``. + Allele-independent: ``predict()`` returns one + :class:`~mhctools.pred.Prediction` per peptide with an empty ``allele`` and + ``kind == Kind.immunogenicity``. Parameters ---------- @@ -121,6 +113,8 @@ class Calis: but are outside its intended scope). """ + mhc_class = "I" + def __init__(self, min_peptide_length=8, max_peptide_length=11): if min_peptide_length < 3: # With P1/P2/C-term masked, peptides shorter than this have no @@ -137,27 +131,9 @@ def __str__(self): return "Calis(min_peptide_length=%d, max_peptide_length=%d)" % ( self.min_peptide_length, self.max_peptide_length) - def __repr__(self): - return str(self) - - def _predictor_name(self): - return "calis" - def _default_pred_kind(self): return Kind.immunogenicity - def kind_support(self): - return { - Kind.immunogenicity: { - "mhc_dependence": "none", - "mhc_class": "I", - }, - } - - @property - def supported_kinds(self): - return tuple(self.kind_support()) - def _check_peptides(self, peptides): for peptide in peptides: n = len(peptide) @@ -182,24 +158,13 @@ def predict(self, peptides): ``Kind.immunogenicity`` prediction (empty ``allele``; ``score`` is the Calis score, higher = more immunogenic). """ - if isinstance(peptides, str): - peptides = [peptides] - peptide_list = [str(p).strip().upper() for p in peptides] + peptide_list = self._normalize_peptides(peptides) self._check_peptides(peptide_list) - - results = [] - for peptide in peptide_list: - pred = Prediction( + return [ + PeptideResult(preds=(Prediction( kind=Kind.immunogenicity, score=immunogenicity_score(peptide), peptide=peptide, - predictor_name="calis") - results.append(PeptideResult(preds=(pred,))) - return results - - def predict_dataframe(self, peptides, sample_name=""): - """``predict()`` flattened to a DataFrame.""" - dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] - if not dfs: - return pd.DataFrame(columns=COLUMNS) - return pd.concat(dfs, ignore_index=True) + predictor_name="calis"),)) + for peptide in peptide_list + ] diff --git a/mhctools/deeptap.py b/mhctools/deeptap.py index f5356e6..4d83ca0 100644 --- a/mhctools/deeptap.py +++ b/mhctools/deeptap.py @@ -20,25 +20,23 @@ mhctools otherwise has no TAP predictor; DeepTAP fills that gap and emits ``Kind.tap_transport``. -Like the cleavage predictors, DeepTAP is **allele-independent**: it scores each +Like the cleavage predictors, DeepTAP is allele-independent: it scores each peptide once, so ``predict()`` returns one prediction per peptide (with an empty ``allele``). -DeepTAP ships its pretrained weights in-repo and is Apache-2.0 licensed, but -pins an old scientific-Python stack (``pytorch-lightning==1.9.2``, torch). To -keep that stack out of the mhctools environment, this wrapper shells out to -DeepTAP's own ``deeptap.py`` CLI in a user-provided checkout, run by a -user-provided interpreter (the Tulip pattern). The checkpoints load fine under -modern Lightning too, so any interpreter with torch + pytorch-lightning works. +DeepTAP ships its pretrained weights in-repo (Apache-2.0) but pins an old stack +(pytorch-lightning==1.9.2, torch). To keep that out of the mhctools environment, +this wrapper shells out to DeepTAP's own ``deeptap.py`` CLI in a user-provided +checkout (``DEEPTAP_HOME``), run by a user-provided interpreter +(``DEEPTAP_PYTHON``, default the current one); the checkpoints also load under +modern Lightning. Upstream: https://github.com/zjupgx/DeepTAP -Cite: Chen et al., *Comput. Biol. Med.* 2023 — "DeepTAP: An RNN-based method of +Cite: Chen et al., Comput. Biol. Med. 2023 — "DeepTAP: An RNN-based method of TAP-binding peptide prediction in the selection of tumor neoantigens". -Note on interpretation: DeepTAP's evaluation is self-reported, and no -independent TAP benchmark exists for any tool (true of the whole TAP field). -Treat the score as a useful pathway signal for prioritization, not a validated -oracle. +DeepTAP's evaluation is self-reported and no independent TAP benchmark exists — +a useful pathway signal, not a validated oracle. """ import os @@ -49,8 +47,9 @@ import pandas as pd from .cleanup_context import CleanupFiles -from .pred import COLUMNS, Kind, PeptideResult, Prediction +from .pred import Kind, PeptideResult, Prediction from .process_helpers import run_command +from .wrapper_base import AlleleFreePredictor # DeepTAP one-hot encodes the 20 standard amino acids plus "X" (unknown / # padding) and pads every peptide to a fixed width of 17; longer peptides @@ -83,7 +82,7 @@ def _find_deeptap_home(deeptap_home=None): return candidate -class DeepTAP: +class DeepTAP(AlleleFreePredictor): """Wrapper for the DeepTAP TAP-transport predictor. Parameters @@ -131,27 +130,9 @@ def __str__(self): return "DeepTAP(task_type=%r, deeptap_home=%r)" % ( self.task_type, self.deeptap_home) - def __repr__(self): - return str(self) - - def _predictor_name(self): - return "deeptap_%s" % self.task_type - def _default_pred_kind(self): return Kind.tap_transport - def kind_support(self): - return { - Kind.tap_transport: { - "mhc_dependence": "none", - "mhc_class": "none", - }, - } - - @property - def supported_kinds(self): - return tuple(self.kind_support()) - def _check_peptides(self, peptides): for peptide in peptides: if not peptide: @@ -178,9 +159,7 @@ def predict(self, peptides): ``"reg"`` mode the prediction's ``value`` is the predicted affinity in nM. """ - if isinstance(peptides, str): - peptides = [peptides] - peptide_list = [str(p).strip().upper() for p in peptides] + peptide_list = self._normalize_peptides(peptides) self._check_peptides(peptide_list) if not peptide_list: return [] @@ -220,13 +199,6 @@ def predict(self, peptides): for peptide in peptide_list ] - def predict_dataframe(self, peptides, sample_name=""): - """``predict()`` flattened to a DataFrame.""" - dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] - if not dfs: - return pd.DataFrame(columns=COLUMNS) - return pd.concat(dfs, ignore_index=True) - def parse_deeptap_results(filename, task_type="cla"): """Parse a DeepTAP prediction CSV into ``{peptide: [Prediction]}``. diff --git a/mhctools/eramer.py b/mhctools/eramer.py index e4b709f..f852b65 100644 --- a/mhctools/eramer.py +++ b/mhctools/eramer.py @@ -21,31 +21,27 @@ position-weight matrix and scores a precursor by averaging the PWM specificity over each residue trimmed off as it is cut down to a target epitope length. -Licensing / why nothing is vendored ------------------------------------ -ERAMER is **GPLv3** and its PWM ships in a GPL-licensed ``PWM.xlsx``; mhctools is -Apache-2.0 and vendors neither. Instead this is a clean-room Python-3 -reimplementation of the (simple, factual) trimming-cascade average — the -upstream tool is Python 2.7 — that loads the PWM from a user-provided ERAMER -checkout at runtime, exactly as the netMHC / MixMHCpred wrappers read -user-provided model files. Point at the checkout with ``ERAMER_HOME`` (or pass -``eramer_home=`` / ``pwm_path=``). +Licensing: ERAMER is GPLv3 and its PWM ships in a GPL-licensed ``PWM.xlsx``; +mhctools is Apache-2.0 and vendors neither. This is a clean-room Python-3 +reimplementation of the trimming-cascade average (the upstream tool is Python +2.7) that loads the PWM from a user-provided ERAMER checkout at runtime, as the +netMHC / MixMHCpred wrappers read user-provided model files. Point at the +checkout with ``ERAMER_HOME`` (or pass ``eramer_home=`` / ``pwm_path=``). Upstream: https://github.com/aalokaily/ERAMER -Cite: Al-okaily et al., *Comput. Biol. Med.* 2024 — "ERAMER: A novel in silico +Cite: Al-okaily et al., Comput. Biol. Med. 2024 — "ERAMER: A novel in silico tool for prediction of ERAP1 enzyme trimming". -Note on interpretation: ERAP1 trimming is a genuine but noisy processing signal -and ERAMER's evaluation is self-reported; treat the score (roughly -1..1, higher -= more likely trimmed) as a pathway prior, not a validated oracle. +ERAMER's evaluation is self-reported and ERAP1 trimming is intrinsically noisy — +treat the score (roughly -1..1, higher = more likely trimmed) as a pathway +prior, not a validated oracle. """ import os from os.path import isdir, isfile, join -import pandas as pd - -from .pred import COLUMNS, Kind, PeptideResult, Prediction +from .pred import Kind, PeptideResult, Prediction +from .wrapper_base import AlleleFreePredictor # ERAP1 processes precursors in this length range; ERAMER ships one PWM sheet # per precursor length in [9, 16]. @@ -165,7 +161,7 @@ def eramer_score(precursor, epitope_length, weights_by_length): return sum(scores) / len(scores) -class ERAMER: +class ERAMER(AlleleFreePredictor): """ERAMER ERAP1-trimming predictor (clean-room reimplementation). Allele-independent: ``predict()`` returns one ``Kind.erap_trimming`` @@ -179,12 +175,15 @@ class ERAMER: from the precursor length down to ``epitope_length + 1``. Default 8 (so every valid 9-16mer precursor yields a score). Must be in 8-15. eramer_home : str, optional - Path to an ERAMER checkout (containing ``PWM.xlsx``). Resolved from the - argument, then ``$ERAMER_HOME``, then ``~/ERAMER``. + Path to an ERAMER checkout (containing ``PWM.xlsx``). If omitted, the PWM + is resolved from ``$ERAMER_PWM``, then ``$ERAMER_HOME``, then + ``~/ERAMER``. pwm_path : str, optional Direct path to ``PWM.xlsx`` (overrides *eramer_home*). """ + mhc_class = "I" + def __init__(self, epitope_length=8, eramer_home=None, pwm_path=None): if not (8 <= epitope_length <= ERAMER_MAX_PRECURSOR_LENGTH - 1): raise ValueError( @@ -199,27 +198,9 @@ def __str__(self): return "ERAMER(epitope_length=%d, pwm_path=%r)" % ( self.epitope_length, self.pwm_path) - def __repr__(self): - return str(self) - - def _predictor_name(self): - return "eramer" - def _default_pred_kind(self): return Kind.erap_trimming - def kind_support(self): - return { - Kind.erap_trimming: { - "mhc_dependence": "none", - "mhc_class": "I", - }, - } - - @property - def supported_kinds(self): - return tuple(self.kind_support()) - def _ensure_loaded(self): if self._weights_by_length is None: self._weights_by_length = load_pwm(self.pwm_path) @@ -255,27 +236,20 @@ def predict(self, peptides): ``Kind.erap_trimming`` prediction (empty ``allele``; ``score`` is the ERAMER trimming score, higher = more likely trimmed). """ - if isinstance(peptides, str): - peptides = [peptides] - peptide_list = [str(p).strip().upper() for p in peptides] + peptide_list = self._normalize_peptides(peptides) self._check_peptides(peptide_list) + if not peptide_list: + return [] weights_by_length = self._ensure_loaded() - results = [] - for peptide in peptide_list: - score = eramer_score( - peptide, self.epitope_length, weights_by_length) - preds = () if score is None else (Prediction( + # _check_peptides guarantees each precursor is longer than + # epitope_length, so eramer_score always yields a value here. + return [ + PeptideResult(preds=(Prediction( kind=Kind.erap_trimming, - score=float(score), + score=float(eramer_score( + peptide, self.epitope_length, weights_by_length)), peptide=peptide, - predictor_name="eramer"),) - results.append(PeptideResult(preds=preds)) - return results - - def predict_dataframe(self, peptides, sample_name=""): - """``predict()`` flattened to a DataFrame.""" - dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] - if not dfs: - return pd.DataFrame(columns=COLUMNS) - return pd.concat(dfs, ignore_index=True) + predictor_name="eramer"),)) + for peptide in peptide_list + ] diff --git a/mhctools/mixmhc2pred.py b/mhctools/mixmhc2pred.py index b23a812..97cb453 100644 --- a/mhctools/mixmhc2pred.py +++ b/mhctools/mixmhc2pred.py @@ -211,21 +211,23 @@ def parse_mixmhc2pred_results(filename, alleles, cli_names): "MixMHC2pred output missing column %r (columns: %s)" % (column, list(df.columns))) + # itertuples mangles '%Rank_...' names, so index by position. Resolve the + # column positions once, not per row. + peptide_idx = df.columns.get_loc("Peptide") + rank_idx = [df.columns.get_loc("%Rank_" + name) for name in cli_names] + score_idx = [df.columns.get_loc("Score_" + name) for name in cli_names] + results = {} for row in df.itertuples(index=False): - peptide = getattr(row, "Peptide") - preds = [] - for allele, cli_name in zip(alleles, cli_names): - # itertuples mangles '%Rank_...' names, so index by position via - # the column list captured from df.columns. - rank = float(row[df.columns.get_loc("%Rank_" + cli_name)]) - score = float(row[df.columns.get_loc("Score_" + cli_name)]) - preds.append(Prediction( + peptide = row[peptide_idx] + preds = [ + Prediction( kind=Kind.pMHC_presentation, - score=score, + score=float(row[score_idx[i]]), peptide=peptide, allele=allele, - percentile_rank=rank, - predictor_name="mixmhc2pred")) + percentile_rank=float(row[rank_idx[i]]), + predictor_name="mixmhc2pred") + for i, allele in enumerate(alleles)] results[peptide] = preds return results diff --git a/mhctools/nettcr.py b/mhctools/nettcr.py index 9be5de2..1fad974 100644 --- a/mhctools/nettcr.py +++ b/mhctools/nettcr.py @@ -26,9 +26,8 @@ (``ai-edge-litert``, ``tflite-runtime``, or ``tensorflow``). The BLOSUM50 encoding, per-feature padding lengths, and name-based input -tensor assignment here reproduce ``src/predict.py`` from NetTCR-2.2 exactly -(verified to ~1e-7 against its own output). Scores are the mean over the -pan cross-validation ensemble, matching the published usage. +tensor assignment here reproduce ``src/predict.py`` from NetTCR-2.2. Scores are +the mean over the pan cross-validation ensemble, matching the published usage. """ from __future__ import annotations @@ -47,19 +46,13 @@ @contextlib.contextmanager def _suppress_native_stderr(): - """Silence C-level stderr for the duration of the block. + """Silence native (C-level) writes to fd 2 for the duration of the block. The TFLite runtimes print ``INFO: Created TensorFlow Lite XNNPACK delegate - for CPU.`` from native code straight to file descriptor 2 (during - ``allocate_tensors``), bypassing Python's ``logging``. Redirect fd 2 to - ``os.devnull`` so that chatter is dropped; Python-level exceptions still - propagate normally (only native writes to stderr are hidden). Restores the - original fd afterward, so ordinary stderr keeps working. - - Note: fd 2 is process-global, so for the (brief) duration of the block a - concurrent thread's stderr and any non-fatal native warning are also - dropped. It wraps only interpreter setup (construction / allocation), never - the prediction loop, so the window is small. + for CPU.`` straight to fd 2 during interpreter setup, bypassing Python's + ``logging``. Redirecting fd 2 is process-global, so this wraps only + construction / allocation, never the prediction loop; Python exceptions + still propagate. """ sys.stderr.flush() # Acquire both descriptors inside the try so a failure partway through diff --git a/mhctools/pred.py b/mhctools/pred.py index 668b2f0..a1add29 100644 --- a/mhctools/pred.py +++ b/mhctools/pred.py @@ -80,30 +80,12 @@ class Kind: def best_direction(kind, field) -> str: - """Canonical "best" direction for a ``(kind, field)`` pair. - - Returns ``"max"`` or ``"min"``. - - Parameters - ---------- - kind : str or Kind - Prediction kind (e.g. ``Kind.pMHC_affinity``). - field : str - Column name within the kind's predictions — - ``"score"``, ``"percentile_rank"``, or ``"value"``. - - Raises - ------ - ValueError - If ``field`` is unknown, or if ``field == "value"`` for a kind - without a registered direction in :data:`VALUE_BEST_DIRECTIONS`. - - Notes - ----- - ``score`` and ``percentile_rank`` directions are uniform across all - kinds. ``value`` is kind-dependent: e.g. ``pMHC_affinity`` reports - IC50 in nM (lower better) while ``pMHC_stability`` reports half-life - (higher better). + """Canonical "best" direction (``"max"`` or ``"min"``) for ``(kind, field)``. + + ``score`` (max) and ``percentile_rank`` (min) are uniform across kinds; + ``value`` is kind-dependent (see :data:`VALUE_BEST_DIRECTIONS`). Raises + ``ValueError`` for an unknown ``field``, or for ``value`` on a kind with no + registered direction. """ direction = FIELD_BEST_DIRECTIONS.get(field) if direction is not None: @@ -122,6 +104,11 @@ def best_direction(kind, field) -> str: ) +def reduce_op(kind, field): + """``max`` or ``min`` — the reducer that picks the best ``(kind, field)``.""" + return max if best_direction(kind, field) == "max" else min + + COLUMNS = ( "sample_name", "peptide", @@ -142,7 +129,7 @@ def best_direction(kind, field) -> str: @dataclass(frozen=True, repr=False) class Prediction: - """Single prediction from one model on one peptide. Self-contained.""" + """Single prediction from one model on one peptide.""" kind: str score: float peptide: str = "" @@ -206,7 +193,7 @@ def from_dict(cls, d): @dataclass(repr=False) class PeptideResult: - """All predictions for one peptide. Contains a tuple of Prediction objects.""" + """All predictions for one peptide (a tuple of ``Prediction`` objects).""" preds: tuple[Prediction, ...] = () def __repr__(self): @@ -371,8 +358,7 @@ def best_by(self, kind, field) -> Optional[Prediction]: allele, falls back to allele-less predictions (e.g. processing predictors that emit allele-independent scores). """ - direction = best_direction(kind, field) - op = max if direction == "max" else min + op = reduce_op(kind, field) def has_value(p): return getattr(p, field) is not None diff --git a/mhctools/prime.py b/mhctools/prime.py index 02e4cc2..3d492ea 100644 --- a/mhctools/prime.py +++ b/mhctools/prime.py @@ -24,16 +24,12 @@ ``mixmhcpred_path`` if it is not on ``PATH``. Upstream: https://github.com/GfellerLab/PRIME -Cite: Gfeller et al., *Cell Systems* 2023 — "Improved predictions of antigen +Cite: Gfeller et al., Cell Systems 2023 — "Improved predictions of antigen presentation and TCR recognition with MixMHCpred2.2 and PRIME2.0". -Note on interpretation: like every current CD8 immunogenicity predictor, PRIME -is useful for ranking in the well-characterized regime but generalizes poorly -to truly novel neoepitopes (independent benchmarks put the field near -AUC 0.5-0.65). PRIME's training positives are dominated by viral / cancer-testis -antigens, so it does relatively better on infectious-disease epitopes and, in -the one neutral head-to-head (NeoaPred, Bioinformatics 2024), trails BigMHC_IM -on cancer neoepitopes. Treat scores as a prioritization aid, not ground truth. +Like every current CD8 immunogenicity predictor, PRIME ranks better than it +generalizes to novel neoepitopes (independent benchmarks put the field near +AUC 0.5-0.65) — a prioritization aid, not ground truth. """ from os import remove diff --git a/mhctools/wrapper_base.py b/mhctools/wrapper_base.py new file mode 100644 index 0000000..511ee12 --- /dev/null +++ b/mhctools/wrapper_base.py @@ -0,0 +1,76 @@ +# Copyright (c) 2016. Mount Sinai School of Medicine +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Shared helpers for predictors that emit the new ``PeptideResult`` model +directly, rather than going through ``BasePredictor``/``BindingPrediction``. + +``NewModelPredictorMixin`` holds the members that were otherwise copy-pasted +across the standalone wrappers (``predict_dataframe``, ``supported_kinds``, +``__repr__``, peptide normalization). ``AlleleFreePredictor`` adds the +allele-independent ``kind_support`` shape used by Calis, DeepTAP, and ERAMER. +""" + +import pandas as pd + +from .pred import COLUMNS + + +class NewModelPredictorMixin: + """Members shared by every new-model wrapper (allele-free or not).""" + + def __repr__(self): + return str(self) + + @property + def supported_kinds(self): + """Prediction kind strings this predictor can emit.""" + return tuple(self.kind_support()) + + @staticmethod + def _normalize_peptides(peptides): + """Accept a single string or an iterable; strip + upper-case each.""" + if isinstance(peptides, str): + peptides = [peptides] + return [str(p).strip().upper() for p in peptides] + + def predict_dataframe(self, peptides, sample_name="", n_flanks=None, + c_flanks=None): + """``predict()`` flattened to a DataFrame. + + ``n_flanks``/``c_flanks`` are accepted for a uniform flank-aware API; + allele-free wrappers do not use flanking context and ignore them. + """ + dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] + if not dfs: + return pd.DataFrame(columns=COLUMNS) + return pd.concat(dfs, ignore_index=True) + + +class AlleleFreePredictor(NewModelPredictorMixin): + """Base for allele-independent new-model predictors. + + Subclasses set ``mhc_class`` and implement ``_default_pred_kind()`` plus the + scoring in ``predict()``; peptide length/character validation stays with the + subclass since each tool has its own bounds and messages. + """ + + mhc_class = "none" + + def kind_support(self): + return { + self._default_pred_kind(): { + "mhc_dependence": "none", + "mhc_class": self.mhc_class, + }, + } diff --git a/tests/test_wrapper_base.py b/tests/test_wrapper_base.py new file mode 100644 index 0000000..def8599 --- /dev/null +++ b/tests/test_wrapper_base.py @@ -0,0 +1,68 @@ +# Copyright (c) 2016. Mount Sinai School of Medicine +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for the shared new-model wrapper base classes.""" + +from mhctools.pred import COLUMNS, Kind, Prediction, PeptideResult +from mhctools.wrapper_base import AlleleFreePredictor, NewModelPredictorMixin + + +def test_normalize_peptides_single_string_and_strip_upper(): + norm = NewModelPredictorMixin._normalize_peptides + assert norm("siinfekl") == ["SIINFEKL"] # single string -> list + assert norm([" gilgfvftl ", "NlV"]) == ["GILGFVFTL", "NLV"] + + +class _FakeAlleleFree(AlleleFreePredictor): + mhc_class = "I" + + def _default_pred_kind(self): + return Kind.immunogenicity + + def predict(self, peptides): + return [ + PeptideResult(preds=(Prediction( + kind=Kind.immunogenicity, score=1.0, peptide=p, + predictor_name="fake"),)) + for p in self._normalize_peptides(peptides)] + + +def test_allele_free_kind_support_and_supported_kinds(): + p = _FakeAlleleFree() + support = p.kind_support()[Kind.immunogenicity] + assert support == {"mhc_dependence": "none", "mhc_class": "I"} + assert p.supported_kinds == (Kind.immunogenicity,) + + +def test_inherited_predict_dataframe(): + p = _FakeAlleleFree() + df = p.predict_dataframe(["siinfekl"], sample_name="s") + assert list(df["peptide"]) == ["SIINFEKL"] + assert df.iloc[0]["kind"] == Kind.immunogenicity + # empty input -> empty frame with the canonical columns + empty = p.predict_dataframe([]) + assert list(empty.columns) == list(COLUMNS) + assert len(empty) == 0 + + +def test_wrappers_use_the_shared_base(): + # The standalone wrappers inherit the shared helpers rather than + # re-implementing them. + from mhctools.calis import Calis + from mhctools.bigmhc import BigMHC + assert issubclass(Calis, AlleleFreePredictor) + assert issubclass(BigMHC, NewModelPredictorMixin) + # BigMHC standardizes on _default_pred_kind (not the old _pred_kind). + assert not hasattr(BigMHC, "_pred_kind") + assert hasattr(BigMHC, "_default_pred_kind")