Add ERAMER ERAP1-trimming predictor (fills the erap_trimming gap)

.github/workflows/tests.yml

@@ -57,6 +57,7 @@ jobs:
             tests/test_mixmhc2pred2.py \
             tests/test_deeptap.py \
             tests/test_calis.py \
+            tests/test_eramer.py \
             tests/test_processing_predictor.py \
             tests/test_pepsickle.py \
             tests/test_bigmhc.py \

README.md

@@ -169,7 +169,7 @@ The canonical prediction kind strings are defined in `mhctools.pred.Kind`.
 | `proteasome_cleavage` | Proteasomal (MHC-I, cytosolic) C-terminal cleavage score |
 | `endolysosomal_cleavage` | Endolysosomal (MHC-II, cathepsin) C-terminal cleavage score |
 | `tap_transport` | TAP transport / binding score |
-| `erap_trimming` | ERAP trimming score (reserved, not yet used) |
+| `erap_trimming` | ERAP1 N-terminal trimming score |
 
 Predictors also expose `kind_support()` so downstream code can tell what MHC
 context is meaningful for each emitted kind:
@@ -210,6 +210,7 @@ Examples:
 | `NetCleave_I` | `proteasome_cleavage` | `none` | `I` |
 | `NetCleave_II` | `endolysosomal_cleavage` | `none` | `II` |
 | `DeepTAP` | `tap_transport` | `none` | `none` |
+| `ERAMER` | `erap_trimming` | `none` | `I` |
 | `NetTCR` | `pMHC_TCR_binding` | `none` | `I` |
 | `Tulip` | `pMHC_TCR_binding` | `single_allele` | `I` |
 | `BigMHC_IM` | `immunogenicity` | `single_allele` | `I` |
@@ -390,6 +391,36 @@ results[1].tap_transport.score             # 0-1, higher = stronger TAP binding
 > exists for any tool (true of the whole TAP field). Treat the score as a useful
 > pathway signal for prioritization, not a validated oracle.
 
+### ERAP1 trimming
+
+| Predictor | Kinds produced | Requires |
+|---|---|---|
+| `ERAMER` | ERAP1 trimming (`erap_trimming`) | [ERAMER](https://github.com/aalokaily/ERAMER) clone with `PWM.xlsx` (set `ERAMER_HOME`) + `openpyxl` |
+
+ERAP1 trims the N-termini of 9–16mer precursor peptides in the ER down to the
+8–10mers MHC-I presents — the step between TAP transport and MHC loading, and
+otherwise the last empty stage in the pathway. `ERAMER` scores a precursor by
+averaging a per-length position-weight-matrix specificity over each residue
+trimmed off as it is cut toward a target epitope length (allele-independent, one
+`erap_trimming` prediction per peptide; `score` roughly −1…1, higher = more
+likely trimmed).
+
+ERAMER is **GPLv3** and its PWM ships in a GPL-licensed `PWM.xlsx`, so mhctools
+vendors neither: this is a clean-room Python-3 reimplementation of the
+(Python-2.7) tool's trimming-cascade average that loads the PWM from a
+user-provided ERAMER checkout at runtime. Point at the clone with `ERAMER_HOME`.
+
+```python
+from mhctools import ERAMER
+
+predictor = ERAMER(epitope_length=8)       # resolves ERAMER_HOME / ~/ERAMER
+results = predictor.predict(["GGGGGVVVVVVAAAEE"])   # a 9-16mer precursor
+results[0].erap_trimming.score
+```
+
+> ⚠️ ERAMER's evaluation is self-reported and ERAP1 trimming is an intrinsically
+> noisy signal; treat the score as a pathway prior, not a validated oracle.
+
 ### Immunogenicity
 
 | Predictor | Kinds produced | Requires |

mhctools/__init__.py

@@ -28,6 +28,7 @@
 from .prime import PRIME
 from .deeptap import DeepTAP
 from .calis import Calis
+from .eramer import ERAMER
 from .processing_predictor import (
     ProcessingPredictor,
     SCORING_MODES,
@@ -84,7 +85,7 @@ def __getattr__(name):
     raise AttributeError(
         "module %r has no attribute %r" % (__name__, name))
 
-__version__ = "3.26.0"
+__version__ = "3.27.0"
 
 __all__ = [
     "Prediction",
@@ -115,6 +116,7 @@ def __getattr__(name):
     "PRIME",
     "DeepTAP",
     "Calis",
+    "ERAMER",
     "MHCflurry",
     "MHCflurry_Affinity",
     "ProcessingPredictor",

mhctools/annotate.py

@@ -72,6 +72,7 @@
     "stability": (Kind.pMHC_stability, "value"),
     "immunogenicity": (Kind.immunogenicity, "score"),
     "tap_transport": (Kind.tap_transport, "score"),
+    "erap_trimming": (Kind.erap_trimming, "score"),
 }
 
 

mhctools/cli/args.py

@@ -66,6 +66,7 @@
     PRIME,
     DeepTAP,
     Calis,
+    ERAMER,
 )
 
 
@@ -173,6 +174,7 @@ def __hash__(self):
     "prime": PRIME,
     "deeptap": DeepTAP,
     "calis": Calis,
+    "eramer": ERAMER,
 }
 
 

mhctools/eramer.py

@@ -0,0 +1,281 @@
+# 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.
+
+"""ERAMER — ERAP1 trimming prediction, filling mhctools' `erap_trimming` gap.
+
+ERAP1 (ER aminopeptidase 1) trims the N-termini of 9-16mer precursor peptides in
+the ER down to the 8-10mers that MHC-I presents — the processing step between
+TAP transport and MHC loading, and previously the one empty ``Kind`` in
+mhctools. ERAMER (Al-okaily 2024) models ERAP1 specificity as a per-length
+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=``).
+
+Upstream: https://github.com/aalokaily/ERAMER
+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.
+"""
+
+import os
+from os.path import isdir, isfile, join
+
+import pandas as pd
+
+from .pred import COLUMNS, Kind, PeptideResult, Prediction
+
+# ERAP1 processes precursors in this length range; ERAMER ships one PWM sheet
+# per precursor length in [9, 16].
+ERAMER_MIN_PRECURSOR_LENGTH = 9
+ERAMER_MAX_PRECURSOR_LENGTH = 16
+_PWM_LENGTHS = tuple(range(
+    ERAMER_MIN_PRECURSOR_LENGTH, ERAMER_MAX_PRECURSOR_LENGTH + 1))
+_VALID_AMINO_ACIDS = frozenset("ACDEFGHIKLMNPQRSTVWY")
+
+
+_MISSING_PWM_HELP = (
+    "Set ERAMER_HOME to a clone of https://github.com/aalokaily/ERAMER "
+    "(or pass eramer_home= / pwm_path=). mhctools does not vendor the "
+    "GPL-licensed PWM.")
+
+
+def _find_pwm_path(eramer_home=None, pwm_path=None):
+    """Resolve the path to ERAMER's ``PWM.xlsx``.
+
+    An explicit *pwm_path* or *eramer_home* is honored exactly (and must exist —
+    no silent fallback). Otherwise the resolver tries ``$ERAMER_PWM``, then
+    ``$ERAMER_HOME/PWM.xlsx``, then ``~/ERAMER/PWM.xlsx``.
+    """
+    if pwm_path:
+        if isfile(pwm_path):
+            return pwm_path
+        raise FileNotFoundError(
+            "ERAMER PWM not found at pwm_path=%r." % pwm_path)
+    if eramer_home:
+        candidate = join(eramer_home, "PWM.xlsx")
+        if isfile(candidate):
+            return candidate
+        raise FileNotFoundError(
+            "PWM.xlsx not found in eramer_home=%r." % eramer_home)
+
+    candidates = []
+    if os.environ.get("ERAMER_PWM"):
+        candidates.append(os.environ["ERAMER_PWM"])
+    if os.environ.get("ERAMER_HOME"):
+        candidates.append(join(os.environ["ERAMER_HOME"], "PWM.xlsx"))
+    home = join(os.path.expanduser("~"), "ERAMER")
+    if isdir(home):
+        candidates.append(join(home, "PWM.xlsx"))
+    for candidate in candidates:
+        if isfile(candidate):
+            return candidate
+    raise FileNotFoundError("ERAMER's PWM.xlsx not found. " + _MISSING_PWM_HELP)
+
+
+def load_pwm(pwm_path):
+    """Load ERAMER's PWM into ``{length: {amino_acid: [weight_per_position]}}``.
+
+    Reads the ``Length 9`` .. ``Length 16`` sheets: each has 20 amino-acid rows
+    (letter in column B, per-position weights in the following columns).
+    """
+    try:
+        import openpyxl
+    except ImportError as e:
+        raise ImportError(
+            "Reading ERAMER's PWM.xlsx requires openpyxl (`pip install "
+            "openpyxl`).") from e
+
+    workbook = openpyxl.load_workbook(pwm_path, data_only=True, read_only=True)
+    weights_by_length = {}
+    for length in _PWM_LENGTHS:
+        sheet_name = "Length %d" % length
+        if sheet_name not in workbook.sheetnames:
+            raise ValueError(
+                "ERAMER PWM.xlsx is missing sheet %r" % sheet_name)
+        worksheet = workbook[sheet_name]
+        aa_weights = {}
+        # 20 amino-acid rows: openpyxl rows 2..21, letter in col 2, the L
+        # per-position weights in cols 3..(2 + length).
+        for row in range(2, 22):
+            aa = worksheet.cell(row=row, column=2).value
+            if not aa:
+                continue
+            aa = str(aa).strip().upper()
+            aa_weights[aa] = [
+                float(worksheet.cell(row=row, column=3 + pos).value)
+                for pos in range(length)]
+        missing = _VALID_AMINO_ACIDS - set(aa_weights)
+        if missing:
+            raise ValueError(
+                "ERAMER PWM sheet %r missing amino acid(s): %s"
+                % (sheet_name, "".join(sorted(missing))))
+        weights_by_length[length] = aa_weights
+    workbook.close()
+    return weights_by_length
+
+
+def _specificity(peptide, aa_weights):
+    """Mean PWM specificity of a peptide (one trimming intermediate)."""
+    length = len(peptide)
+    total = sum(aa_weights[peptide[i]][i] for i in range(length))
+    # Rounded per-residue-average, matching upstream ERAMER's compute_specificty.
+    return round(total / length, 6)
+
+
+def eramer_score(precursor, epitope_length, weights_by_length):
+    """ERAMER trimming score: mean specificity over the trimming cascade.
+
+    As the precursor is trimmed one residue at a time from length ``L`` down to
+    ``epitope_length + 1``, each intermediate C-terminal ``l``-mer is scored by
+    its PWM specificity; the returned score is their mean (higher = more likely
+    trimmed by ERAP1).
+
+    Returns ``None`` if the precursor is already at/below ``epitope_length``
+    (nothing to trim).
+    """
+    length = len(precursor)
+    scores = [
+        _specificity(precursor[-sub_len:], weights_by_length[sub_len])
+        for sub_len in range(length, epitope_length, -1)]
+    if not scores:
+        return None
+    return sum(scores) / len(scores)
+
+
+class ERAMER:
+    """ERAMER ERAP1-trimming predictor (clean-room reimplementation).
+
+    Allele-independent: ``predict()`` returns one ``Kind.erap_trimming``
+    prediction per peptide (the peptide is treated as an ERAP1 *precursor*),
+    with an empty ``allele``.
+
+    Parameters
+    ----------
+    epitope_length : int
+        Target epitope length the precursor is trimmed toward; the cascade runs
+        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``.
+    pwm_path : str, optional
+        Direct path to ``PWM.xlsx`` (overrides *eramer_home*).
+    """
+
+    def __init__(self, epitope_length=8, eramer_home=None, pwm_path=None):
+        if not (8 <= epitope_length <= ERAMER_MAX_PRECURSOR_LENGTH - 1):
+            raise ValueError(
+                "epitope_length must be in 8-%d (the shortest trimmed "
+                "intermediate must still be a >=9mer with a PWM); got %d"
+                % (ERAMER_MAX_PRECURSOR_LENGTH - 1, epitope_length))
+        self.epitope_length = epitope_length
+        self.pwm_path = _find_pwm_path(eramer_home, pwm_path)
+        self._weights_by_length = None
+
+    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)
+        return self._weights_by_length
+
+    def _check_peptides(self, peptides):
+        for peptide in peptides:
+            n = len(peptide)
+            if not (ERAMER_MIN_PRECURSOR_LENGTH <= n
+                    <= ERAMER_MAX_PRECURSOR_LENGTH):
+                raise ValueError(
+                    "ERAMER precursors must be %d-%d residues; got %r (length "
+                    "%d)" % (ERAMER_MIN_PRECURSOR_LENGTH,
+                             ERAMER_MAX_PRECURSOR_LENGTH, peptide, n))
+            if n <= self.epitope_length:
+                raise ValueError(
+                    "Precursor %r (length %d) is not longer than the target "
+                    "epitope_length %d, so there is nothing to trim"
+                    % (peptide, n, self.epitope_length))
+            invalid = set(peptide) - _VALID_AMINO_ACIDS
+            if invalid:
+                raise ValueError(
+                    "Peptide %r contains non-standard amino acid(s): %s"
+                    % (peptide, "".join(sorted(invalid))))
+
+    def predict(self, peptides):
+        """Predict ERAP1 trimming for a list of precursor peptides.
+
+        Returns
+        -------
+        list of PeptideResult
+            One entry per input peptide, each holding a single
+            ``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]
+        self._check_peptides(peptide_list)
+        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(
+                kind=Kind.erap_trimming,
+                score=float(score),
+                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)