Surface MHCflurry's antigen-processing (cleavage) score (#235)

README.md

@@ -204,6 +204,7 @@ Examples:
 | `MHCflurry` | `pMHC_affinity` | `single_allele` | `I` |
 | `MHCflurry` haplotype mode | `pMHC_presentation` | `haplotype` | `I` |
 | `MHCflurry` per-allele panel mode | `pMHC_presentation` | `single_allele` | `I` |
+| `MHCflurry` | `antigen_processing` | `none` | `none` |
 | `Pepsickle` | `proteasome_cleavage` | `none` | `none` |
 | `NetCleave_I` | `proteasome_cleavage` | `none` | `I` |
 | `NetCleave_II` | `endolysosomal_cleavage` | `none` | `II` |
@@ -285,7 +286,7 @@ affinity, hours for stability). `percentile_rank` is always optional,
 | `NetMHCIIpan` / `NetMHCIIpan43` | affinity or presentation | [NetMHCIIpan](https://services.healthtech.dtu.dk/services/NetMHCIIpan-4.3/) |
 | `NetMHCcons` | affinity | [NetMHCcons](https://services.healthtech.dtu.dk/services/NetMHCcons-1.1/) |
 | `NetMHCstabpan` | stability | [NetMHCstabpan](https://services.healthtech.dtu.dk/services/NetMHCstabpan-1.0/) |
-| `MHCflurry` | affinity + presentation | `pip install mhcflurry` + `mhcflurry-downloads fetch` |
+| `MHCflurry` | affinity + presentation + processing | `pip install mhcflurry` + `mhcflurry-downloads fetch` |
 | `MHCflurry_Affinity` | affinity | `pip install mhcflurry` + `mhcflurry-downloads fetch` |
 | `BigMHC` | presentation or immunogenicity | [BigMHC](https://github.com/KarchinLab/bigmhc) clone (set `BIGMHC_DIR`) |
 | `MixMHCpred` | presentation | [MixMHCpred](https://github.com/GfellerLab/MixMHCpred) |

mhctools/__init__.py

@@ -80,7 +80,7 @@ def __getattr__(name):
     raise AttributeError(
         "module %r has no attribute %r" % (__name__, name))
 
-__version__ = "3.21.0"
+__version__ = "3.22.0"
 
 __all__ = [
     "Prediction",

mhctools/annotate.py

@@ -68,6 +68,7 @@
     "rank": (None, "percentile_rank"),
     "value": (Kind.pMHC_affinity, "value"),
     "presentation": (Kind.pMHC_presentation, "score"),
+    "processing": (Kind.antigen_processing, "score"),
     "stability": (Kind.pMHC_stability, "value"),
     "immunogenicity": (Kind.immunogenicity, "score"),
 }

mhctools/mhcflurry.py

@@ -102,8 +102,11 @@ class MHCflurry(BasePredictor):
     Produces per-allele ``pMHC_affinity`` predictions. For presentation,
     ``presentation_allele_mode`` controls whether mhctools treats the allele
     set as one class-I haplotype or as a panel of independent one-allele
-    samples. The legacy ``predict_peptides`` method returns BindingPrediction
-    objects based on affinity values for backward compat.
+    samples. It also surfaces MHCflurry's ``antigen_processing`` (cleavage)
+    score — computed by the presentation model from the peptide + flanks and
+    allele-independent — as one allele-less prediction per peptide. The legacy
+    ``predict_peptides`` method returns BindingPrediction objects based on
+    affinity values for backward compat.
 
     See https://github.com/openvax/mhcflurry
     """
@@ -317,6 +320,11 @@ def predict(self, peptides, n_flanks=None, c_flanks=None):
                     expected_presentation_rows))
 
         pres_by_peptide_index = {i: [] for i in range(len(peptide_list))}
+        # MHCflurry's presentation predict() also returns a processing_score
+        # (its antigen-processing / cleavage head). It depends only on the
+        # peptide + flanks, not the allele, so it's identical across the
+        # per-allele rows of a peptide; we keep the first seen per peptide.
+        processing_by_peptide_index = {}
         seen_presentation_keys = set()
         for row_position, row in enumerate(pres_df.itertuples(index=False)):
             row_index = int(getattr(row, "peptide_num", row_position))
@@ -339,6 +347,8 @@ def predict(self, peptides, n_flanks=None, c_flanks=None):
                 row.presentation_percentile,
                 allele,
             ))
+            processing_by_peptide_index.setdefault(
+                row_index, getattr(row, "processing_score", None))
 
         groups = [list() for _ in peptide_list]
         for row_index, row in zip(batch_indices, aff_df.itertuples(index=False)):
@@ -391,6 +401,21 @@ def predict(self, peptides, n_flanks=None, c_flanks=None):
                     predictor_name="mhcflurry",
                 ))
 
+            # Surface MHCflurry's antigen-processing (cleavage) score, which
+            # its presentation predictor already computes from the peptide +
+            # flanks. Allele-independent, so emit once per peptide, allele-less.
+            processing_score = processing_by_peptide_index.get(row_index)
+            if processing_score is not None:
+                groups[row_index].append(Prediction(
+                    kind=Kind.antigen_processing,
+                    score=processing_score,
+                    peptide=pep,
+                    allele="",
+                    n_flank=n_flank,
+                    c_flank=c_flank,
+                    predictor_name="mhcflurry",
+                ))
+
         return [PeptideResult(preds=tuple(preds)) for preds in groups]
 
     def predict_with_flanks(self, peptides, n_flanks, c_flanks):
@@ -416,6 +441,10 @@ def kind_support(self):
                 "mhc_dependence": presentation_dependence,
                 "mhc_class": "I",
             },
+            Kind.antigen_processing: {
+                "mhc_dependence": "none",
+                "mhc_class": "none",
+            },
         }
 
 

mhctools/pred.py

@@ -272,6 +272,11 @@ def immunogenicity(self) -> Optional[Prediction]:
         """Best immunogenicity prediction, or None."""
         return self.best_by_score(Kind.immunogenicity)
 
+    @property
+    def processing(self) -> Optional[Prediction]:
+        """Best antigen-processing prediction, or None."""
+        return self.best_by_score(Kind.antigen_processing)
+
     @property
     def cleavage(self) -> Optional[Prediction]:
         """Best proteasomal cleavage prediction, or None."""

tests/test_annotate_table.py

@@ -345,6 +345,21 @@ def test_peptide_matched_ignoring_surrounding_whitespace():
 
 # --- allele-free predictor path (by-peptide lookup) -------------------------
 
+def test_processing_field_token_pins_antigen_processing_kind():
+    # the 'processing' field selects Kind.antigen_processing (higher-better)
+    spec = AnnotationSpec(_factory, "col", field="processing")
+    assert spec.kind == Kind.antigen_processing
+    assert spec.prediction_field == "score"
+    df = pd.DataFrame({"peptide": ["SIINFEKL", "GILGFVFTL", "WWWWWWWWW"]})
+    out = annotate_table(
+        df,
+        [AnnotationSpec(lambda alleles: _ProcessingFixturePredictor(),
+                        "proc", field="processing")])
+    assert out.iloc[0]["proc"] == 0.80
+    assert out.iloc[1]["proc"] == 0.30
+    assert math.isnan(out.iloc[2]["proc"])
+
+
 def test_allele_free_predictor_uses_by_peptide_lookup():
     df = pd.DataFrame({"peptide": ["SIINFEKL", "GILGFVFTL", "WWWWWWWWW"]})
     out = annotate_table(

tests/test_mhcflurry.py

@@ -165,14 +165,97 @@ def test_mhcflurry_multiple_alleles():
     eq_(1, len(results), "Expected one PeptideResult")
     r = results[0]
 
-    # Two per-allele affinity predictions plus one haplotype-level
-    # presentation prediction with best_allele attribution.
-    eq_(3, len(r.preds), "Expected 3 predictions")
+    # Two per-allele affinity predictions, one haplotype-level presentation
+    # prediction with best_allele attribution, and one allele-independent
+    # antigen-processing prediction.
+    eq_(4, len(r.preds), "Expected 4 predictions")
     eq_(2, len(r.filter(kind=Kind.pMHC_affinity)))
     eq_(1, len(r.filter(kind=Kind.pMHC_presentation)))
+    eq_(1, len(r.filter(kind=Kind.antigen_processing)))
 
-    # Both alleles should be present through affinity predictions.
+    # Both alleles should be present through affinity predictions
+    # (the processing prediction is allele-less, so it doesn't add one).
     assert r.alleles == set(alleles)
 
-    # Both kinds should be present
-    assert r.kinds == {Kind.pMHC_affinity, Kind.pMHC_presentation}
+    # All three kinds should be present
+    assert r.kinds == {
+        Kind.pMHC_affinity, Kind.pMHC_presentation, Kind.antigen_processing}
+
+
+def test_mhcflurry_processing_score():
+    """MHCflurry surfaces its antigen-processing (cleavage) score.
+
+    The score is allele-independent (one per peptide, no allele) and matches
+    the ``processing_score`` MHCflurry's presentation predictor computes.
+    """
+    from mhcflurry import Class1PresentationPredictor
+
+    alleles = ["HLA-A*02:01", "HLA-B*07:02"]
+    peptides = ["SIINFEKL", "GILGFVFTL"]
+    n_flanks = ["AAA", "CCC"]
+    c_flanks = ["KKK", "DDD"]
+
+    predictor = MHCflurry(alleles=alleles)
+    results = predictor.predict(peptides, n_flanks=n_flanks, c_flanks=c_flanks)
+
+    # Ground truth straight from MHCflurry.
+    raw = Class1PresentationPredictor.load().predict(
+        peptides=peptides,
+        alleles={a: [a] for a in alleles},
+        n_flanks=n_flanks,
+        c_flanks=c_flanks,
+        include_affinity_percentile=False,
+        verbose=0)
+    expected = dict(zip(raw["peptide"], raw["processing_score"]))
+
+    for r in results:
+        processing = r.filter(kind=Kind.antigen_processing)
+        eq_(1, len(processing),
+            "Expected exactly one processing prediction per peptide")
+        pred = processing[0]
+        assert pred.allele == "", "Processing score is allele-independent"
+        assert pred.predictor_name == "mhcflurry"
+        # flanks are carried through as provenance
+        idx = peptides.index(r.peptide)
+        assert pred.n_flank == n_flanks[idx]
+        assert pred.c_flank == c_flanks[idx]
+        testing.assert_allclose(pred.score, expected[r.peptide], rtol=1e-5)
+        # the convenience accessor resolves to the same prediction
+        assert r.processing is not None
+        assert r.processing.kind == Kind.antigen_processing
+
+    # antigen_processing is advertised as an allele-independent supported kind
+    support = predictor.kind_support()[Kind.antigen_processing]
+    eq_("none", support["mhc_dependence"])
+    eq_("none", support["mhc_class"])
+
+
+def test_mhcflurry_processing_score_without_flanks():
+    """The processing score is still emitted when no flanks are provided."""
+    predictor = MHCflurry(alleles=["HLA-A*02:01"])
+    r = predictor.predict(["SIINFEKL"])[0]
+    processing = r.filter(kind=Kind.antigen_processing)
+    eq_(1, len(processing), "Expected one processing prediction")
+    pred = processing[0]
+    assert pred.allele == ""
+    assert pred.n_flank == "" and pred.c_flank == ""
+    assert pred.score is not None
+
+
+def test_mhcflurry_legacy_predict_peptides_unchanged():
+    """Surfacing processing on predict() must not change the legacy path.
+
+    The CLI (--output-csv) and predict_peptides go through the affinity-only
+    BindingPredictionCollection path, which should stay one affinity record
+    per (peptide, allele) with no antigen_processing rows.
+    """
+    alleles = ["HLA-A*02:01", "HLA-B*07:02"]
+    peptides = ["SIINFEKL", "GILGFVFTL"]
+    predictor = MHCflurry(alleles=alleles)
+    collection = predictor.predict_peptides(peptides)
+
+    eq_(len(peptides) * len(alleles), len(collection),
+        "Legacy path should emit one affinity record per (peptide, allele)")
+    for bp in collection:
+        assert bp.affinity is not None
+        assert bp.prediction_method_name == "mhcflurry"