From 50ed087327a356c0a6fc0722bd86a86e047f5738 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 20:32:08 -0400 Subject: [PATCH 1/3] Add NetCleave predictor for MHC-I/II C-terminal cleavage (#213) Adds NetCleave, filling the MHC-II antigen-processing gap: NetChop and pepsickle only cover MHC-I proteasomal cleavage, and no class-II (endolysosomal) option existed. - New Kind.endolysosomal_cleavage (MHC-II, cathepsin) alongside the existing proteasome_cleavage (MHC-I, cytosolic); PeptideResult.endolysosomal_cleavage accessor. Naming per the decision recorded on #213. - NetCleave predictor (+ NetCleave_I / NetCleave_II subclasses). Unlike NetChop/pepsickle it emits one C-terminal score per peptide, so it shells out to NetCleave.py (NetChop-style subprocess) rather than subclassing ProcessingPredictor. Class I -> proteasome_cleavage, class II -> endolysosomal_cleavage. - predict(peptides, c_flanks=...) scores peptides given their downstream residues; predict_proteins() scans proteins so peptides are scored in real context. DataFrame variants included. NetCleave ships its Keras .h5 weights in-repo (git clone, no download); the R dependency in its README is only for training-data generation, not prediction (confirmed: not referenced by any Python in the predict path). Tests reproduce NetCleave's own CLI output exactly (abs 1e-3) on both the class-I and class-II pan models, from upstream `NetCleave.py --predict` (non-circular). Model-free constructor/validation tests run in CI without an install; model tests skip-gate on NETCLEAVE_DIR. A class-I > class-II signal check reflects the paper's AUC 0.91 vs 0.66. Version 3.15.0 -> 3.16.0. Claude-Session: https://claude.ai/code/session_01LZahFhBSCiehXTESCYQ7wG --- README.md | 37 +++- mhctools/__init__.py | 6 +- mhctools/netcleave.py | 416 ++++++++++++++++++++++++++++++++++++++++ mhctools/pred.py | 6 + tests/test_netcleave.py | 206 ++++++++++++++++++++ 5 files changed, 666 insertions(+), 5 deletions(-) create mode 100644 mhctools/netcleave.py create mode 100644 tests/test_netcleave.py diff --git a/README.md b/README.md index d03fe40..73f2873 100644 --- a/README.md +++ b/README.md @@ -166,7 +166,8 @@ The canonical prediction kind strings are defined in `mhctools.pred.Kind`. | `pMHC_TCR_binding` | TCR recognition of a peptide-MHC (pMHC:TCR binding) | | `immunogenicity` | T-cell immunogenicity | | `antigen_processing` | Combined processing score | -| `proteasome_cleavage` | Proteasomal cleavage score | +| `proteasome_cleavage` | Proteasomal (MHC-I, cytosolic) C-terminal cleavage score | +| `endolysosomal_cleavage` | Endolysosomal (MHC-II, cathepsin) C-terminal cleavage score | | `tap_transport` | TAP transport score (reserved, not yet used) | | `erap_trimming` | ERAP trimming score (reserved, not yet used) | @@ -204,6 +205,8 @@ Examples: | `MHCflurry` haplotype mode | `pMHC_presentation` | `haplotype` | `I` | | `MHCflurry` per-allele panel mode | `pMHC_presentation` | `single_allele` | `I` | | `Pepsickle` | `proteasome_cleavage` | `none` | `none` | +| `NetCleave_I` | `proteasome_cleavage` | `none` | `I` | +| `NetCleave_II` | `endolysosomal_cleavage` | `none` | `II` | | `NetTCR` | `pMHC_TCR_binding` | `none` | `I` | For MHCflurry presentation, `presentation_allele_mode="haplotype"` treats the @@ -265,10 +268,36 @@ affinity, hours for stability). `percentile_rank` is always optional, |---|---|---| | `Pepsickle` | proteasome cleavage | `pip install pepsickle` ([paper](https://doi.org/10.1093/bioinformatics/btab628)) | | `NetChop` | proteasome cleavage | [NetChop](https://services.healthtech.dtu.dk/services/NetChop-3.1/) | +| `NetCleave_I` / `NetCleave_II` | proteasomal (I) / endolysosomal (II) C-terminal cleavage | [NetCleave](https://github.com/BSC-CNS-EAPM/NetCleave) clone (set `NETCLEAVE_DIR`) | -Processing predictors use configurable scoring to aggregate per-position -cleavage probabilities into peptide-level scores. See `ProcessingPredictor` -and `ProteasomePredictor` for details. +`Pepsickle` and `NetChop` use configurable scoring to aggregate per-position +cleavage probabilities into peptide-level scores (see `ProcessingPredictor` +and `ProteasomePredictor`). + +`NetCleave` is different: it emits a **single C-terminal cleavage score per +peptide** and covers **both** the MHC-I proteasomal (`NetCleave_I` → +`proteasome_cleavage`) and MHC-II endolysosomal (`NetCleave_II` → +`endolysosomal_cleavage`) pathways — MHC-II processing is otherwise a gap in +the predictor set. It needs the residues downstream of the peptide to build +the cleavage site, so pass `c_flanks` (or scan proteins). Its weights ship in +the git repo; the R dependency in NetCleave's README is only for its training +pipeline, not prediction. + +```python +from mhctools import NetCleave_II + +predictor = NetCleave_II() # resolves NETCLEAVE_DIR / ~/NetCleave +# score peptides with their C-terminal flanking residues (>= 3) +results = predictor.predict(["SIINFEKL"], c_flanks=["DGH"]) +results[0].endolysosomal_cleavage.score + +# or scan a protein so each peptide is scored in real context +by_protein = predictor.predict_proteins({"TP53": "MEEPQ..."}, peptide_lengths=[15]) +``` + +> ⚠️ NetCleave's own paper reports class-II C-terminal cleavage is a much +> weaker signal than class I (AUC ~0.66 vs ~0.91). Treat +> `endolysosomal_cleavage` scores accordingly. ### TCR specificity diff --git a/mhctools/__init__.py b/mhctools/__init__.py index af85cf4..38ab889 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -36,6 +36,7 @@ ) from .proteasome_predictor import ProteasomePredictor from .netchop import NetChop +from .netcleave import NetCleave, NetCleave_I, NetCleave_II from .pepsickle import Pepsickle from .netmhc import NetMHC from .netmhc3 import NetMHC3 @@ -77,7 +78,7 @@ def __getattr__(name): raise AttributeError( "module %r has no attribute %r" % (__name__, name)) -__version__ = "3.15.0" +__version__ = "3.16.0" __all__ = [ "Prediction", @@ -114,6 +115,9 @@ def __getattr__(name): "score_nterm_cterm_anti_max_internal", "score_nterm_cterm_anti_mean_internal", "NetChop", + "NetCleave", + "NetCleave_I", + "NetCleave_II", "Pepsickle", "NetMHC", "NetMHC3", diff --git a/mhctools/netcleave.py b/mhctools/netcleave.py new file mode 100644 index 0000000..7aa6dd7 --- /dev/null +++ b/mhctools/netcleave.py @@ -0,0 +1,416 @@ +# 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. + +"""Wrapper for NetCleave (https://github.com/BSC-CNS-EAPM/NetCleave). + +NetCleave predicts **C-terminal antigen processing** for both the MHC-I +(proteasomal) and MHC-II (endolysosomal) pathways. It scores the cleavage +site (4 residues of the peptide C-terminus + 3 downstream residues) using a +neural network over QSAR descriptors, and emits **one score per peptide** +(unlike NetChop/pepsickle, which emit a per-position profile). + +Because the output granularity is per-peptide rather than per-position, this +wrapper does not subclass :class:`ProcessingPredictor`; it shells out to +NetCleave's ``NetCleave.py --predict`` CLI (the NetChop-style subprocess +model) and parses the resulting CSV. The class-I and class-II models emit +different kinds: + +* class I -> :attr:`~mhctools.pred.Kind.proteasome_cleavage` +* class II -> :attr:`~mhctools.pred.Kind.endolysosomal_cleavage` + +NetCleave ships its pretrained weights (small Keras ``.h5`` files) directly +in its git repository, so only a ``git clone`` is required — no separate +download. Its ``--predict`` path needs Python (tensorflow/keras, scikit-learn, +biopython, pandas, numpy); the R dependency in NetCleave's README is only for +its training-data generation and is **not** used for prediction. + +.. warning:: + NetCleave's own paper reports class-II C-terminal cleavage is a much + weaker signal than class I (AUC ~0.66 vs ~0.91): the class-II peptide + C-terminus lies outside the binding groove and cathepsin specificity is + diffuse. Treat ``endolysosomal_cleavage`` scores accordingly. +""" + +import logging +import os +import subprocess +import sys +import tempfile + +import pandas as pd + +from .base_predictor import ( + _check_flank_inputs, + _normalize_sequence_dict, + _peptide_contexts, +) +from .pred import COLUMNS, Kind, PeptideResult, Prediction + +logger = logging.getLogger(__name__) + +NETCLEAVE_TIMEOUT_SECONDS = 600 + +# NetCleave needs 3 residues downstream of the peptide C-terminus to build +# the 4+3 cleavage site; without them it cannot score the site. +_C_FLANK_REQUIRED = 3 + + +def _find_netcleave_dir(netcleave_path=None): + """Resolve the NetCleave installation directory. + + Checks, in order: the *netcleave_path* argument, the ``NETCLEAVE_DIR`` + environment variable, then ``~/NetCleave`` and ``~/code/NetCleave``. An + explicitly-provided path is validated up front. + """ + clone_hint = "Clone from https://github.com/BSC-CNS-EAPM/NetCleave" + for source, path in ( + ("netcleave_path argument", netcleave_path), + ("NETCLEAVE_DIR", os.environ.get("NETCLEAVE_DIR"))): + if path: + if not os.path.isdir(path): + raise FileNotFoundError( + "NetCleave directory from %s does not exist: %s. %s" + % (source, path, clone_hint)) + return path + home = os.path.expanduser("~") + for candidate in ( + os.path.join(home, "NetCleave"), + os.path.join(home, "code", "NetCleave")): + if os.path.isdir(candidate): + return candidate + raise FileNotFoundError( + "NetCleave not found. Set NETCLEAVE_DIR or pass netcleave_path= to " + "the constructor. %s" % clone_hint) + + +class NetCleave(object): + """Wrapper for NetCleave C-terminal cleavage predictions. + + Parameters + ---------- + mhc_class : str + ``"I"`` (proteasomal) or ``"II"`` (endolysosomal). Selects the + default model and the emitted :class:`~mhctools.pred.Kind`. + mhc_allele : str + Allele label used to pick the pretrained model directory + (``data/models/{class}_mass-spectrometry_{allele}``). Default + ``"HLA"`` (the pan-allele model for the class). + netcleave_path : str, optional + Path to the cloned NetCleave repository. Resolved from + ``NETCLEAVE_DIR`` / ``~/NetCleave`` when omitted. + model_path : str, optional + Full path to a specific model directory, overriding + ``mhc_class`` / ``mhc_allele`` selection. + python_executable : str, optional + Interpreter used to run ``NetCleave.py``. Defaults to the current + interpreter; override if NetCleave's dependencies live in a separate + environment. + subprocess_timeout : int + Timeout (seconds) for a single NetCleave invocation. + + Notes + ----- + NetCleave is distributed under GPL-v2; this wrapper only *runs* a + user-provided installation and vendors none of it. + """ + + VALID_CLASSES = ("I", "II") + + def __init__(self, mhc_class="I", mhc_allele="HLA", netcleave_path=None, + model_path=None, python_executable=None, + subprocess_timeout=NETCLEAVE_TIMEOUT_SECONDS): + if mhc_class not in self.VALID_CLASSES: + raise ValueError( + "mhc_class must be one of %s, got %r" + % (self.VALID_CLASSES, mhc_class)) + self.mhc_class = mhc_class + self.mhc_allele = mhc_allele + self.netcleave_dir = _find_netcleave_dir(netcleave_path) + self.python_executable = python_executable or sys.executable + self.subprocess_timeout = subprocess_timeout + + self._script = os.path.join(self.netcleave_dir, "NetCleave.py") + if not os.path.isfile(self._script): + raise FileNotFoundError( + "NetCleave.py not found in %s" % self.netcleave_dir) + + if model_path is None: + model_path = os.path.join( + self.netcleave_dir, "data", "models", + "%s_mass-spectrometry_%s" % (mhc_class, mhc_allele)) + self.model_path = model_path + weights = os.path.join( + model_path, "%s_model.h5" % os.path.basename(model_path)) + if not os.path.isfile(weights): + raise FileNotFoundError( + "NetCleave model weights not found: %s. Use --mhc_options to " + "list available models, or pass model_path=." % weights) + + self._call_counter = 0 + + def __str__(self): + return "NetCleave(mhc_class=%s, model=%s)" % ( + self.mhc_class, os.path.basename(self.model_path)) + + def __repr__(self): + return str(self) + + def _predictor_name(self): + return "netcleave" + + def _pred_kind(self): + return (Kind.proteasome_cleavage if self.mhc_class == "I" + else Kind.endolysosomal_cleavage) + + def kind_support(self): + return { + self._pred_kind(): { + # Cleavage is MHC-independent (no allele in the prediction). + "mhc_dependence": "none", + "mhc_class": self.mhc_class, + } + } + + @property + def supported_kinds(self): + return tuple(self.kind_support()) + + # ------------------------------------------------------------------ + # Subprocess + # ------------------------------------------------------------------ + + def _run_netcleave(self, epitopes, protein_seqs): + """Run NetCleave (pred_input 3) on parallel epitope/protein lists. + + Returns a list of scores (float or None) aligned to the inputs; + None means NetCleave could not build a cleavage site (e.g. missing + downstream context). + """ + self._call_counter += 1 + # Basename must contain no '.' before the extension — NetCleave + # derives the output filename via ``basename.split('.')[0]``. + basename = "mhctools_netcleave_%d_%d" % (os.getpid(), self._call_counter) + tmp_dir = tempfile.mkdtemp(prefix="mhctools_netcleave_") + input_csv = os.path.join(tmp_dir, basename + ".csv") + output_csv = os.path.join( + self.netcleave_dir, "output", basename + "_NetCleave.csv") + + pd.DataFrame({ + "epitope": [e.upper() for e in epitopes], + "protein_seq": [s.upper() for s in protein_seqs], + "protein_name": [str(i) for i in range(len(epitopes))], + }).to_csv(input_csv, index=False) + + try: + try: + result = subprocess.run( + [self.python_executable, self._script, + "--predict", input_csv, "--pred_input", "3", + "--model_path", self.model_path, + "--mhc_class", self.mhc_class], + cwd=self.netcleave_dir, + capture_output=True, + timeout=self.subprocess_timeout, + ) + except subprocess.TimeoutExpired as e: + raise RuntimeError( + "NetCleave timed out after %d seconds on %d epitopes" + % (self.subprocess_timeout, len(epitopes))) from e + except OSError as e: + raise RuntimeError( + "Could not run NetCleave with %r: %s" + % (self.python_executable, e)) from e + + stderr_text = result.stderr.decode("utf-8", errors="replace").strip() + if result.returncode != 0: + raise RuntimeError( + "NetCleave exited with code %d.\nstdout: %s\nstderr: %s" + % (result.returncode, + result.stdout.decode("utf-8", errors="replace").strip(), + stderr_text)) + if not os.path.isfile(output_csv): + raise RuntimeError( + "NetCleave produced no output file at %s.\nstderr: %s" + % (output_csv, stderr_text)) + + out = pd.read_csv(output_csv) + if len(out) != len(epitopes): + raise RuntimeError( + "NetCleave returned %d rows for %d epitopes" + % (len(out), len(epitopes))) + scores = [] + for value in out["prediction"].tolist(): + # NaN (no cleavage site) is not equal to itself. + scores.append(None if value != value else float(value)) + return scores + finally: + for path in (input_csv, output_csv): + try: + os.remove(path) + except OSError: + pass + try: + os.rmdir(tmp_dir) + except OSError: + pass + + # ------------------------------------------------------------------ + # Public API + # ------------------------------------------------------------------ + + def predict(self, peptides, n_flanks=None, c_flanks=None): + """Predict C-terminal cleavage scores for peptides. + + NetCleave needs the residues **downstream** of the peptide to build + the cleavage site, so ``c_flanks`` (>= 3 residues each) is required. + ``n_flanks`` is accepted for interface parity and recorded on each + :class:`Prediction`, but does not affect the C-terminal score. + + Parameters + ---------- + peptides : list of str + n_flanks : list of str, optional + c_flanks : list of str + C-terminal flanks (>= 3 residues) — one per peptide. + + Returns + ------- + list of PeptideResult + One per peptide; empty (no preds) for peptides whose cleavage + site could not be built. + """ + peptide_list, n_flank_list, c_flank_list = _check_flank_inputs( + peptides, n_flanks, c_flanks) + if c_flank_list is None: + raise ValueError( + "NetCleave.predict requires c_flanks (>= %d residues per " + "peptide) to build the C-terminal cleavage site. To score " + "peptides in a protein context, use predict_proteins()." + % _C_FLANK_REQUIRED) + + epitopes, proteins, meta = [], [], [] + for i, peptide in enumerate(peptide_list): + n_flank = n_flank_list[i] if n_flank_list is not None else "" + c_flank = c_flank_list[i] + if len(c_flank) < _C_FLANK_REQUIRED: + logger.warning( + "NetCleave: peptide %r has c_flank %r shorter than %d " + "residues; cannot score its cleavage site", + peptide, c_flank, _C_FLANK_REQUIRED) + meta.append(None) + continue + meta.append((peptide, n_flank, c_flank)) + epitopes.append(peptide) + proteins.append(n_flank + peptide + c_flank) + + scores = self._run_netcleave(epitopes, proteins) if epitopes else [] + + results, score_idx = [], 0 + for entry in meta: + if entry is None: + results.append(PeptideResult(preds=())) + continue + peptide, n_flank, c_flank = entry + score = scores[score_idx] + score_idx += 1 + results.append(self._make_result( + peptide, score, n_flank=n_flank, c_flank=c_flank)) + return results + + def predict_proteins(self, sequence_dict, peptide_lengths=None, + flank_length=3): + """Score peptides scanned from full protein sequences. + + Each peptide is scored in its real protein context, so NetCleave sees + the true downstream residues. This is the preferred entry-point. + + Parameters + ---------- + sequence_dict : dict or str + peptide_lengths : list of int, optional + Default: ``[9]`` for class I, ``[15]`` for class II. + flank_length : int + Residues of flank recorded on each Prediction (default 3). + + Returns + ------- + dict mapping sequence_name -> list of PeptideResult + """ + sequence_dict = _normalize_sequence_dict(sequence_dict) + if peptide_lengths is None: + peptide_lengths = [9] if self.mhc_class == "I" else [15] + if isinstance(peptide_lengths, int): + peptide_lengths = [peptide_lengths] + + contexts = _peptide_contexts( + sequence_dict, peptide_lengths, flank_length) + epitopes = [c.peptide for c in contexts] + proteins = [sequence_dict[c.source_sequence_name] for c in contexts] + scores = self._run_netcleave(epitopes, proteins) if epitopes else [] + + from collections import defaultdict + results = defaultdict(list) + for context, score in zip(contexts, scores): + results[context.source_sequence_name].append(self._make_result( + context.peptide, score, + n_flank=context.n_flank, c_flank=context.c_flank, + source_sequence_name=context.source_sequence_name, + offset=context.offset)) + return dict(results) + + def _make_result(self, peptide, score, n_flank="", c_flank="", + source_sequence_name=None, offset=0): + if score is None: + return PeptideResult(preds=()) + return PeptideResult(preds=(Prediction( + kind=self._pred_kind(), + score=score, + peptide=peptide, + n_flank=n_flank, + c_flank=c_flank, + source_sequence_name=source_sequence_name, + offset=offset, + predictor_name=self._predictor_name(), + ),)) + + def predict_dataframe(self, peptides, n_flanks=None, c_flanks=None, + sample_name=""): + """``predict()`` flattened to a DataFrame.""" + dfs = [pp.to_dataframe(sample_name) + for pp in self.predict(peptides, n_flanks, c_flanks)] + if not dfs: + return pd.DataFrame(columns=COLUMNS) + return pd.concat(dfs, ignore_index=True) + + def predict_proteins_dataframe(self, sequence_dict, peptide_lengths=None, + flank_length=3, sample_name=""): + """``predict_proteins()`` flattened to a DataFrame.""" + dfs = [] + for pp_list in self.predict_proteins( + sequence_dict, peptide_lengths, flank_length).values(): + for pp in pp_list: + dfs.append(pp.to_dataframe(sample_name)) + if not dfs: + return pd.DataFrame(columns=COLUMNS) + return pd.concat(dfs, ignore_index=True) + + +class NetCleave_I(NetCleave): + """NetCleave for the MHC-I (proteasomal) pathway.""" + def __init__(self, mhc_allele="HLA", **kwargs): + NetCleave.__init__(self, mhc_class="I", mhc_allele=mhc_allele, **kwargs) + + +class NetCleave_II(NetCleave): + """NetCleave for the MHC-II (endolysosomal) pathway.""" + def __init__(self, mhc_allele="HLA", **kwargs): + NetCleave.__init__(self, mhc_class="II", mhc_allele=mhc_allele, **kwargs) diff --git a/mhctools/pred.py b/mhctools/pred.py index 679afaa..4dca7bd 100644 --- a/mhctools/pred.py +++ b/mhctools/pred.py @@ -50,6 +50,7 @@ class Kind: immunogenicity = "immunogenicity" antigen_processing = "antigen_processing" proteasome_cleavage = "proteasome_cleavage" + endolysosomal_cleavage = "endolysosomal_cleavage" tap_transport = "tap_transport" erap_trimming = "erap_trimming" @@ -276,6 +277,11 @@ def cleavage(self) -> Optional[Prediction]: """Best proteasomal cleavage prediction, or None.""" return self.best_by_score(Kind.proteasome_cleavage) + @property + def endolysosomal_cleavage(self) -> Optional[Prediction]: + """Best endolysosomal (MHC-II) cleavage prediction, or None.""" + return self.best_by_score(Kind.endolysosomal_cleavage) + @property def tcr_binding(self) -> Optional[Prediction]: """Best pMHC:TCR binding prediction, or None.""" diff --git a/tests/test_netcleave.py b/tests/test_netcleave.py new file mode 100644 index 0000000..c595ff8 --- /dev/null +++ b/tests/test_netcleave.py @@ -0,0 +1,206 @@ +# 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. + +import os + +import pytest + +from mhctools.netcleave import NetCleave, NetCleave_I, NetCleave_II +from mhctools.pred import COLUMNS, Kind + + +# --------------------------------------------------------------------------- +# Constructor / validation — no NetCleave install required. +# --------------------------------------------------------------------------- + +def test_init_invalid_class(): + with pytest.raises(ValueError, match="mhc_class"): + NetCleave(mhc_class="bad") + + +def test_init_missing_path_raises(): + with pytest.raises(FileNotFoundError, match="does not exist"): + NetCleave(netcleave_path="/nonexistent/netcleave") + + +# --------------------------------------------------------------------------- +# Model tests — require a cloned NetCleave (weights) and a Python with +# NetCleave's own deps (tensorflow/keras, scikit-learn, biopython). +# --------------------------------------------------------------------------- + +NETCLEAVE_DIR = None +for candidate in [ + os.environ.get("NETCLEAVE_DIR", ""), + os.path.join(os.path.expanduser("~"), "NetCleave"), + os.path.join(os.path.expanduser("~"), "code", "NetCleave"), +]: + if candidate and os.path.isfile(os.path.join(candidate, "NetCleave.py")) \ + and os.path.isfile(os.path.join( + candidate, "data", "models", "I_mass-spectrometry_HLA", + "I_mass-spectrometry_HLA_model.h5")): + NETCLEAVE_DIR = candidate + break + +requires_netcleave = pytest.mark.skipif( + NETCLEAVE_DIR is None, + reason="NetCleave not installed (set NETCLEAVE_DIR or clone to ~/NetCleave)") + + +# Reference C-terminal cleavage scores from NetCleave's OWN CLI +# (`NetCleave.py --predict ... --pred_input 3`) on the pan-allele models, +# with protein_seq = peptide + c_flank. Upstream code, not this wrapper. +# Keys: (peptide, c_flank). +CLASS_I_REF = { + ("SIINFEKL", "DGH"): 0.842449, + ("GILGFVFTL", "QRS"): 0.908024, + ("NLVPMVATV", "KKK"): 0.582904, + ("LLWNGPMAV", "QRS"): 0.696255, +} +CLASS_II_REF = { + ("SIINFEKL", "DGH"): 0.360996, + ("GILGFVFTL", "QRS"): 0.260713, + ("NLVPMVATV", "KKK"): 0.245104, + ("LLWNGPMAV", "QRS"): 0.227258, +} + + +@pytest.fixture(scope="module") +def predictor_I(): + return NetCleave_I(netcleave_path=NETCLEAVE_DIR) + + +@pytest.fixture(scope="module") +def predictor_II(): + return NetCleave_II(netcleave_path=NETCLEAVE_DIR) + + +def _predict_ref(predictor, ref): + pairs = list(ref) + peptides = [p for p, _ in pairs] + c_flanks = [c for _, c in pairs] + return pairs, predictor.predict(peptides, c_flanks=c_flanks) + + +@requires_netcleave +def test_reproduces_cli_class_I(predictor_I): + pairs, results = _predict_ref(predictor_I, CLASS_I_REF) + for (pep, cflank), pp in zip(pairs, results): + got = pp.cleavage.score + assert got == pytest.approx(CLASS_I_REF[(pep, cflank)], abs=1e-3) + + +@requires_netcleave +def test_reproduces_cli_class_II(predictor_II): + pairs, results = _predict_ref(predictor_II, CLASS_II_REF) + for (pep, cflank), pp in zip(pairs, results): + got = pp.endolysosomal_cleavage.score + assert got == pytest.approx(CLASS_II_REF[(pep, cflank)], abs=1e-3) + + +@requires_netcleave +def test_class_I_emits_proteasome_kind(predictor_I): + pp = predictor_I.predict(["SIINFEKL"], c_flanks=["DGH"])[0] + pred = pp.preds[0] + assert pred.kind == Kind.proteasome_cleavage + assert pred.predictor_name == "netcleave" + assert pred.c_flank == "DGH" + assert pred.allele == "" + assert 0.0 <= pred.score <= 1.0 + + +@requires_netcleave +def test_class_II_emits_endolysosomal_kind(predictor_II): + pp = predictor_II.predict(["SIINFEKL"], c_flanks=["DGH"])[0] + assert pp.preds[0].kind == Kind.endolysosomal_cleavage + + +@requires_netcleave +def test_kind_support_metadata(predictor_I, predictor_II): + assert predictor_I.kind_support() == { + Kind.proteasome_cleavage: {"mhc_dependence": "none", "mhc_class": "I"}} + assert predictor_II.kind_support() == { + Kind.endolysosomal_cleavage: {"mhc_dependence": "none", "mhc_class": "II"}} + + +@requires_netcleave +def test_class_I_stronger_than_class_II(predictor_I, predictor_II): + """NetCleave's class-I C-terminal signal is much stronger than class-II + (paper AUC 0.91 vs 0.66); the pan models reflect this on these peptides.""" + pairs = list(CLASS_I_REF) + peptides = [p for p, _ in pairs] + c_flanks = [c for _, c in pairs] + s1 = [pp.cleavage.score + for pp in predictor_I.predict(peptides, c_flanks=c_flanks)] + s2 = [pp.endolysosomal_cleavage.score + for pp in predictor_II.predict(peptides, c_flanks=c_flanks)] + for a, b in zip(s1, s2): + assert a > b + + +@requires_netcleave +def test_predict_requires_c_flanks(predictor_I): + with pytest.raises(ValueError, match="c_flanks"): + predictor_I.predict(["SIINFEKL"]) + + +@requires_netcleave +def test_short_c_flank_yields_empty_result(predictor_I): + # 2-residue c_flank can't build the 4+3 site -> empty PeptideResult, + # but the result list stays aligned 1:1 with the input. + results = predictor_I.predict(["SIINFEKL"], c_flanks=["AB"]) + assert len(results) == 1 + assert results[0].preds == () + + +@requires_netcleave +def test_mixed_valid_and_short_flanks_stay_aligned(predictor_I): + results = predictor_I.predict( + ["SIINFEKL", "GILGFVFTL"], c_flanks=["AB", "QRS"]) + assert len(results) == 2 + assert results[0].preds == () # short flank + assert results[1].cleavage is not None # valid + + +@requires_netcleave +def test_predict_proteins(predictor_I): + protein = "MASIINFEKLDGHKQRLLWNGPMAVQRSTTT" # SIINFEKL@2, LLWNGPMAV@16 + results = predictor_I.predict_proteins({"p": protein}, peptide_lengths=[8, 9]) + assert "p" in results + scored = {pp.cleavage.peptide: pp.cleavage + for pp in results["p"] if pp.cleavage} + # peptide scored in real protein context matches the CLI reference + assert scored["SIINFEKL"].score == pytest.approx(0.842449, abs=1e-3) + assert scored["SIINFEKL"].offset == 2 + assert scored["SIINFEKL"].source_sequence_name == "p" + + +@requires_netcleave +def test_predict_dataframe_schema(predictor_II): + df = predictor_II.predict_dataframe(["SIINFEKL"], c_flanks=["DGH"]) + assert list(df.columns) == list(COLUMNS) + assert df["kind"].iloc[0] == Kind.endolysosomal_cleavage + assert df["c_flank"].iloc[0] == "DGH" + + +@requires_netcleave +def test_repeated_calls_consistent(predictor_I): + a = predictor_I.predict(["SIINFEKL"], c_flanks=["DGH"])[0].cleavage.score + b = predictor_I.predict(["SIINFEKL"], c_flanks=["DGH"])[0].cleavage.score + assert a == b + + +@requires_netcleave +def test_subclasses_set_class(predictor_I, predictor_II): + assert predictor_I.mhc_class == "I" + assert predictor_II.mhc_class == "II" + assert isinstance(predictor_I, NetCleave) + assert isinstance(predictor_II, NetCleave) From ec92f67fd0a0b3fec593c334d605a60d954cacc3 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 22:17:40 -0400 Subject: [PATCH 2/3] Fix NetCleave output-row mapping and temp-file collision (review) Addresses two bugs found in review (self-review + ultrareview): 1. predict_proteins crashed ("NetCleave returned N rows for M epitopes") on any protein containing a repeated peptide, and predict() failed when a peptide recurred within peptide+flank. NetCleave's pred_input 3 emits one row per regex occurrence of the epitope, so output rows are not 1:1 with input rows. Fix: score each peptide in its own surrogate protein (peptide + downstream flank) and de-multiplex output back to inputs by row id + expected cleavage site, instead of asserting positional 1:1. The surrogate yields identical scores (the site depends only on the C-terminal 4 + 3 downstream residues) and shrinks the predict_proteins CSV. 2. Concurrent NetCleave instances collided on a shared output/.csv (basename was pid + per-instance counter). Fix: use uuid4 in the basename. Also: document that an instance isn't thread-safe, and return [] from predict([]) instead of raising. Adds regression tests: repeated peptide in predict_proteins (both occurrences scored in their own downstream context, cross-checked against predict()), peptide recurring in flank, two-instance independence, and empty input. Claude-Session: https://claude.ai/code/session_01LZahFhBSCiehXTESCYQ7wG --- mhctools/netcleave.py | 130 ++++++++++++++++++++++++++++------------ tests/test_netcleave.py | 51 ++++++++++++++++ 2 files changed, 142 insertions(+), 39 deletions(-) diff --git a/mhctools/netcleave.py b/mhctools/netcleave.py index 7aa6dd7..ba00f2c 100644 --- a/mhctools/netcleave.py +++ b/mhctools/netcleave.py @@ -45,6 +45,8 @@ import subprocess import sys import tempfile +import uuid +from collections import defaultdict import pandas as pd @@ -121,6 +123,10 @@ class NetCleave(object): ----- NetCleave is distributed under GPL-v2; this wrapper only *runs* a user-provided installation and vendors none of it. + + A ``NetCleave`` instance is not safe to call from multiple threads + concurrently (each call shells out via a per-call temp file); use one + instance per thread. """ VALID_CLASSES = ("I", "II") @@ -188,26 +194,37 @@ def supported_kinds(self): # Subprocess # ------------------------------------------------------------------ - def _run_netcleave(self, epitopes, protein_seqs): - """Run NetCleave (pred_input 3) on parallel epitope/protein lists. + def _run_netcleave(self, requests): + """Run NetCleave (pred_input 3) on a list of scoring requests. + + Each request is ``(epitope, protein_seq, expected_site)``: the + peptide, the sequence to search it in, and the expected 7-residue + cleavage site (the peptide's C-terminal 4 residues + 3 downstream). + + NetCleave emits **one output row per regex occurrence** of the + epitope in ``protein_seq`` (and drops occurrences without 3 + downstream residues), so results are matched back to requests by + input-row id and cleavage site — never by row position. - Returns a list of scores (float or None) aligned to the inputs; - None means NetCleave could not build a cleavage site (e.g. missing - downstream context). + Returns a list of scores (float or None) aligned to *requests*; + None means no matching valid cleavage site was produced. """ + if not requests: + return [] self._call_counter += 1 - # Basename must contain no '.' before the extension — NetCleave + # Basename must be unique (a shared NetCleave output/ dir is used by + # every instance) and contain no '.' before the extension — NetCleave # derives the output filename via ``basename.split('.')[0]``. - basename = "mhctools_netcleave_%d_%d" % (os.getpid(), self._call_counter) + basename = "mhctools_netcleave_%d_%s" % (os.getpid(), uuid.uuid4().hex) tmp_dir = tempfile.mkdtemp(prefix="mhctools_netcleave_") input_csv = os.path.join(tmp_dir, basename + ".csv") output_csv = os.path.join( self.netcleave_dir, "output", basename + "_NetCleave.csv") pd.DataFrame({ - "epitope": [e.upper() for e in epitopes], - "protein_seq": [s.upper() for s in protein_seqs], - "protein_name": [str(i) for i in range(len(epitopes))], + "epitope": [r[0].upper() for r in requests], + "protein_seq": [r[1].upper() for r in requests], + "protein_name": [str(i) for i in range(len(requests))], }).to_csv(input_csv, index=False) try: @@ -224,7 +241,7 @@ def _run_netcleave(self, epitopes, protein_seqs): except subprocess.TimeoutExpired as e: raise RuntimeError( "NetCleave timed out after %d seconds on %d epitopes" - % (self.subprocess_timeout, len(epitopes))) from e + % (self.subprocess_timeout, len(requests))) from e except OSError as e: raise RuntimeError( "Could not run NetCleave with %r: %s" @@ -243,14 +260,23 @@ def _run_netcleave(self, epitopes, protein_seqs): % (output_csv, stderr_text)) out = pd.read_csv(output_csv) - if len(out) != len(epitopes): - raise RuntimeError( - "NetCleave returned %d rows for %d epitopes" - % (len(out), len(epitopes))) + # Map input-row id -> {cleavage_site (upper): score}. NetCleave + # carries our per-row protein_name through as ``uniprot_id`` and + # the 4+3 site as ``cleavage_site``; invalid sites read as NaN. + by_row = defaultdict(dict) + for rid, site, value in zip( + out["uniprot_id"], out["cleavage_site"], out["prediction"]): + try: + score = float(value) + except (TypeError, ValueError): + score = None + if score is not None and score != score: # NaN + score = None + by_row[str(rid)][str(site).upper()] = score + scores = [] - for value in out["prediction"].tolist(): - # NaN (no cleavage site) is not equal to itself. - scores.append(None if value != value else float(value)) + for i, (_epitope, _protein_seq, expected_site) in enumerate(requests): + scores.append(by_row.get(str(i), {}).get(expected_site.upper())) return scores finally: for path in (input_csv, output_csv): @@ -263,6 +289,14 @@ def _run_netcleave(self, epitopes, protein_seqs): except OSError: pass + @staticmethod + def _expected_site(peptide, c_flank): + """The 7-residue cleavage site NetCleave scores: peptide C-terminal + 4 residues + 3 downstream. Empty if there isn't enough context.""" + if len(peptide) < 4 or len(c_flank) < _C_FLANK_REQUIRED: + return "" + return (peptide[-4:] + c_flank[:_C_FLANK_REQUIRED]).upper() + # ------------------------------------------------------------------ # Public API # ------------------------------------------------------------------ @@ -290,6 +324,8 @@ def predict(self, peptides, n_flanks=None, c_flanks=None): """ peptide_list, n_flank_list, c_flank_list = _check_flank_inputs( peptides, n_flanks, c_flanks) + if not peptide_list: + return [] if c_flank_list is None: raise ValueError( "NetCleave.predict requires c_flanks (>= %d residues per " @@ -297,31 +333,28 @@ def predict(self, peptides, n_flanks=None, c_flanks=None): "peptides in a protein context, use predict_proteins()." % _C_FLANK_REQUIRED) - epitopes, proteins, meta = [], [], [] + requests, meta = [], [] for i, peptide in enumerate(peptide_list): n_flank = n_flank_list[i] if n_flank_list is not None else "" c_flank = c_flank_list[i] - if len(c_flank) < _C_FLANK_REQUIRED: + expected_site = self._expected_site(peptide, c_flank) + if not expected_site: logger.warning( "NetCleave: peptide %r has c_flank %r shorter than %d " "residues; cannot score its cleavage site", peptide, c_flank, _C_FLANK_REQUIRED) - meta.append(None) + meta.append((peptide, n_flank, c_flank, None)) continue - meta.append((peptide, n_flank, c_flank)) - epitopes.append(peptide) - proteins.append(n_flank + peptide + c_flank) + # Score each peptide in its own surrogate protein so its cleavage + # site is well-defined regardless of repeats elsewhere. + meta.append((peptide, n_flank, c_flank, len(requests))) + requests.append((peptide, peptide + c_flank, expected_site)) - scores = self._run_netcleave(epitopes, proteins) if epitopes else [] + scores = self._run_netcleave(requests) - results, score_idx = [], 0 - for entry in meta: - if entry is None: - results.append(PeptideResult(preds=())) - continue - peptide, n_flank, c_flank = entry - score = scores[score_idx] - score_idx += 1 + results = [] + for peptide, n_flank, c_flank, req_idx in meta: + score = None if req_idx is None else scores[req_idx] results.append(self._make_result( peptide, score, n_flank=n_flank, c_flank=c_flank)) return results @@ -351,15 +384,34 @@ def predict_proteins(self, sequence_dict, peptide_lengths=None, if isinstance(peptide_lengths, int): peptide_lengths = [peptide_lengths] + # Ensure at least 3 downstream residues are captured for the cleavage + # site, even if the caller asks for a shorter recorded flank. contexts = _peptide_contexts( - sequence_dict, peptide_lengths, flank_length) - epitopes = [c.peptide for c in contexts] - proteins = [sequence_dict[c.source_sequence_name] for c in contexts] - scores = self._run_netcleave(epitopes, proteins) if epitopes else [] + sequence_dict, peptide_lengths, flank_length, + n_flank_length=flank_length, + c_flank_length=max(flank_length, _C_FLANK_REQUIRED)) + + # Score each peptide in a surrogate protein (peptide + its downstream + # flank). The cleavage site depends only on those residues, so this + # matches the full-protein score while keeping each peptide unique in + # its row (NetCleave emits one row per regex occurrence). + requests, meta = [], [] + for context in contexts: + expected_site = self._expected_site(context.peptide, context.c_flank) + if not expected_site: + meta.append((context, None)) + continue + meta.append((context, len(requests))) + requests.append(( + context.peptide, + context.peptide + context.c_flank, + expected_site)) + + scores = self._run_netcleave(requests) - from collections import defaultdict results = defaultdict(list) - for context, score in zip(contexts, scores): + for context, req_idx in meta: + score = None if req_idx is None else scores[req_idx] results[context.source_sequence_name].append(self._make_result( context.peptide, score, n_flank=context.n_flank, c_flank=context.c_flank, diff --git a/tests/test_netcleave.py b/tests/test_netcleave.py index c595ff8..76cfda1 100644 --- a/tests/test_netcleave.py +++ b/tests/test_netcleave.py @@ -183,6 +183,51 @@ def test_predict_proteins(predictor_I): assert scored["SIINFEKL"].source_sequence_name == "p" +@requires_netcleave +def test_predict_proteins_repeated_peptide(predictor_I): + """Regression: a peptide occurring more than once in the protein must not + crash (NetCleave emits one row per occurrence). Each occurrence is scored + in its own downstream context.""" + # SIINFEKL appears at offset 2 (downstream DGH) and offset 16 (downstream QRS) + protein = "MASIINFEKLDGHKQRSIINFEKLQRSTTT" + results = predictor_I.predict_proteins({"p": protein}, peptide_lengths=[8]) + hits = [pp.cleavage for pp in results["p"] + if pp.cleavage and pp.cleavage.peptide == "SIINFEKL"] + assert len(hits) == 2 + by_offset = {h.offset: h for h in hits} + assert by_offset[2].c_flank == "DGH" + assert by_offset[16].c_flank == "QRS" + # Each occurrence matches predict() with the same downstream flank, and + # the differing downstream context gives differing scores. + assert by_offset[2].score == pytest.approx( + predictor_I.predict(["SIINFEKL"], c_flanks=["DGH"])[0].cleavage.score, + abs=1e-6) + assert by_offset[16].score == pytest.approx( + predictor_I.predict(["SIINFEKL"], c_flanks=["QRS"])[0].cleavage.score, + abs=1e-6) + assert by_offset[2].score != by_offset[16].score + + +@requires_netcleave +def test_predict_recurring_peptide_in_flank(predictor_I): + """Regression: a peptide that recurs within peptide+flank must not crash.""" + pp = predictor_I.predict(["AAAAAAAA"], c_flanks=["AAA"])[0] + assert pp.cleavage is not None + assert 0.0 <= pp.cleavage.score <= 1.0 + + +@requires_netcleave +def test_two_instances_do_not_interfere(predictor_I, predictor_II): + """Regression: distinct instances must not collide on the shared output + dir — each keeps its own class-correct score and Kind.""" + a = predictor_I.predict(["SIINFEKL"], c_flanks=["DGH"])[0].cleavage + b = predictor_II.predict(["SIINFEKL"], c_flanks=["DGH"])[0].endolysosomal_cleavage + assert a.kind == Kind.proteasome_cleavage + assert b.kind == Kind.endolysosomal_cleavage + assert a.score == pytest.approx(CLASS_I_REF[("SIINFEKL", "DGH")], abs=1e-3) + assert b.score == pytest.approx(CLASS_II_REF[("SIINFEKL", "DGH")], abs=1e-3) + + @requires_netcleave def test_predict_dataframe_schema(predictor_II): df = predictor_II.predict_dataframe(["SIINFEKL"], c_flanks=["DGH"]) @@ -198,6 +243,12 @@ def test_repeated_calls_consistent(predictor_I): assert a == b +@requires_netcleave +def test_predict_empty_returns_empty(predictor_I): + assert predictor_I.predict([], c_flanks=[]) == [] + assert predictor_I.predict([]) == [] # no c_flanks needed for empty input + + @requires_netcleave def test_subclasses_set_class(predictor_I, predictor_II): assert predictor_I.mhc_class == "I" From 954216ba49f85e5996e9f543b50c95cabc531fb4 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 22:44:31 -0400 Subject: [PATCH 3/3] Guard NetCleave row de-mux against epitope mismatch Belt-and-suspenders for the occurrence-mapping fix: when resolving each input row's score, verify NetCleave's output rows for that row id actually carry the requested epitope, so a future NetCleave output-contract change raises loudly instead of silently misaligning scores. Claude-Session: https://claude.ai/code/session_01LZahFhBSCiehXTESCYQ7wG --- mhctools/netcleave.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/mhctools/netcleave.py b/mhctools/netcleave.py index ba00f2c..f2b725c 100644 --- a/mhctools/netcleave.py +++ b/mhctools/netcleave.py @@ -264,8 +264,10 @@ def _run_netcleave(self, requests): # carries our per-row protein_name through as ``uniprot_id`` and # the 4+3 site as ``cleavage_site``; invalid sites read as NaN. by_row = defaultdict(dict) - for rid, site, value in zip( - out["uniprot_id"], out["cleavage_site"], out["prediction"]): + row_epitope = {} + for rid, epi, site, value in zip( + out["uniprot_id"], out["epitope"], + out["cleavage_site"], out["prediction"]): try: score = float(value) except (TypeError, ValueError): @@ -273,10 +275,19 @@ def _run_netcleave(self, requests): if score is not None and score != score: # NaN score = None by_row[str(rid)][str(site).upper()] = score + row_epitope.setdefault(str(rid), str(epi).upper()) scores = [] - for i, (_epitope, _protein_seq, expected_site) in enumerate(requests): - scores.append(by_row.get(str(i), {}).get(expected_site.upper())) + for i, (epitope, _protein_seq, expected_site) in enumerate(requests): + key = str(i) + # Guard against a NetCleave contract change silently + # misaligning rows: rows for this id must be our epitope. + seen = row_epitope.get(key) + if seen is not None and seen != epitope.upper(): + raise RuntimeError( + "NetCleave row %s epitope %r does not match the " + "requested peptide %r" % (key, seen, epitope)) + scores.append(by_row.get(key, {}).get(expected_site.upper())) return scores finally: for path in (input_csv, output_csv):