From 23c7562a7ff6f38b6fd5103babeda656fd75c2eb Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Tue, 7 Apr 2026 12:07:10 -0400 Subject: [PATCH 1/3] Add unified prediction model (Pred, PeptidePreds, MultiSample) New data model with self-contained Pred objects, PeptidePreds grouping with best_affinity/best_presentation/best_stability helpers, and MultiSample for running across patients with different HLA genotypes. New API: - predictor.predict(peptides) -> list[PeptidePreds] - predictor.predict_proteins(seq_dict) -> {name: list[PeptidePreds]} - predictor.predict_dataframe() / predict_proteins_dataframe() - MultiSample(samples={name: alleles}, predictor_class=...) wrapper Header-driven parser (parse_netmhcpan_to_preds) returns Pred objects natively. NetMHCpan 4.1 emits both pMHC_affinity and pMHC_presentation per peptide-allele pair. Full backward compatibility: BindingPrediction, predict_peptides(), predict_subsequences() all still work. Conversion bridge via to_pred()/from_pred()/to_peptide_preds(). --- README.md | 246 +++++++++++++---- mhctools/__init__.py | 7 + mhctools/base_commandline_predictor.py | 89 ++++++- mhctools/base_predictor.py | 113 +++++++- mhctools/binding_prediction.py | 39 +++ mhctools/binding_prediction_collection.py | 21 ++ mhctools/netmhc_pan28.py | 3 +- mhctools/netmhc_pan3.py | 3 +- mhctools/netmhc_pan4.py | 3 +- mhctools/netmhc_pan41.py | 3 +- mhctools/parsing.py | 245 ++++++++++++----- mhctools/pred.py | 169 ++++++++++++ mhctools/sample.py | 94 +++++++ tests/test_known_class2_epitopes.py | 2 +- tests/test_mhc_formats.py | 107 ++++++++ tests/test_pred.py | 305 ++++++++++++++++++++++ 16 files changed, 1328 insertions(+), 121 deletions(-) create mode 100644 mhctools/pred.py create mode 100644 mhctools/sample.py create mode 100644 tests/test_pred.py diff --git a/README.md b/README.md index ac82c49..ee3f8f9 100644 --- a/README.md +++ b/README.md @@ -3,19 +3,188 @@ PyPI - - - # mhctools Python interface to running command-line and web-based MHC binding predictors. +## Installation + +```sh +pip install mhctools +``` + +For MHCflurry support, also run: + +```sh +mhcflurry-downloads fetch +``` + +## Quick start + +```python +from mhctools import NetMHCpan41 + +predictor = NetMHCpan41(alleles=["HLA-A*02:01", "HLA-B*07:02"]) + +# predict for specific peptides +results = predictor.predict(["SIINFEKL", "GILGFVFTL"]) + +# results is a list of PeptidePreds — one per peptide +for pp in results: + best = pp.best_affinity + if best: + print(f"{best.peptide} -> {best.allele} IC50={best.value:.1f}nM") +``` + +## Python API + +### Predicting peptides + +`predict()` takes a list of peptide sequences and returns a `list[PeptidePreds]`. +Each `PeptidePreds` contains `Pred` objects for every allele and measurement +kind the predictor supports. + +```python +from mhctools import NetMHCpan41 + +predictor = NetMHCpan41(alleles=["HLA-A*02:01", "HLA-B*07:02"]) +results = predictor.predict(["SIINFEKL", "GILGFVFTL"]) + +pp = results[0] +pp.best_affinity # Pred with highest affinity score +pp.best_affinity.allele # "HLA-A*02:01" +pp.best_affinity.value # IC50 in nM +pp.best_affinity.score # higher = better (~0-1) +pp.best_affinity.percentile_rank # lower = better (0-100) + +pp.best_affinity_by_rank # Pred with lowest percentile rank +pp.best_presentation # best EL/presentation score +pp.best_presentation_by_rank # best EL percentile rank +pp.best_stability # best pMHC stability (if available) +pp.best_stability_by_rank + +# filter by kind or allele +pp.filter(kind=Kind.pMHC_affinity) +pp.filter(allele="HLA-A*02:01") +``` + +NetMHCpan 4.1 automatically emits both `pMHC_affinity` and `pMHC_presentation` +predictions per peptide-allele pair. + +### Scanning proteins + +`predict_proteins()` takes a dictionary of protein sequences and returns +`{sequence_name: list[PeptidePreds]}`: + +```python +proteins = predictor.predict_proteins( + {"TP53": "MEEPQSDPSVEPPLSQETFS...", "KRAS": "MTEYKLVVVGAGGVGKS..."}, + peptide_lengths=[9, 10], +) + +for pp in proteins["TP53"]: + best = pp.best_affinity + if best and best.value < 500: + print(f" offset={best.offset} {best.peptide} IC50={best.value:.0f}") +``` + +### DataFrames + +Every level has a `_dataframe` variant that flattens to a pandas DataFrame +with consistent columns: + +```python +df = predictor.predict_dataframe(["SIINFEKL"], sample_name="pat001") +df = predictor.predict_proteins_dataframe({"TP53": "MEEPQ..."}, sample_name="pat001") +``` + +Columns: `sample_name`, `peptide`, `n_flank`, `c_flank`, +`source_sequence_name`, `offset`, `predictor_name`, `predictor_version`, +`allele`, `kind`, `score`, `value`, `percentile_rank`. + +### Multi-sample predictions + +`MultiSample` runs a predictor across multiple samples, each with its own +HLA genotype: + +```python +from mhctools import MultiSample, NetMHCpan41 + +ms = MultiSample( + samples={ + "pat001": ["HLA-A*02:01", "HLA-B*07:02"], + "pat002": ["HLA-A*01:01", "HLA-B*08:01"], + }, + predictor_class=NetMHCpan41, +) + +# {sample_name: list[PeptidePreds]} +results = ms.predict(["SIINFEKL", "GILGFVFTL"]) + +# {sample_name: {seq_name: list[PeptidePreds]}} +protein_results = ms.predict_proteins({"TP53": "MEEPQ..."}) + +# flat DataFrames with sample_name column +df = ms.predict_dataframe(["SIINFEKL"]) +df = ms.predict_proteins_dataframe({"TP53": "MEEPQ..."}) +``` + +### Measurement kinds + +The `Kind` enum describes what biological quantity a `Pred` measures: + +| Kind | Meaning | +|---|---| +| `pMHC_affinity` | Peptide-MHC binding affinity | +| `pMHC_presentation` | Likelihood of surface presentation (EL) | +| `pMHC_stability` | Peptide-MHC complex stability | +| `cellular_presentation` | Cross-allele presentation (e.g. MHCflurry) | +| `antigen_processing` | Combined processing score | +| `proteasome_cleavage` | Proteasomal cleavage score | +| `tap_transport` | TAP transport score | +| `erap_trimming` | ERAP trimming score | + +### The Pred object + +Every prediction is a frozen, self-contained `Pred` dataclass: + +```python +from mhctools import Pred, Kind + +pred = Pred( + kind=Kind.pMHC_affinity, + score=0.85, # ~0-1, higher = better + peptide="SIINFEKL", + allele="HLA-A*02:01", + value=120.5, # IC50 in nM + percentile_rank=0.8, + source_sequence_name="TP53", + offset=42, + predictor_name="netMHCpan", + predictor_version="4.1", +) +``` + +`score` is always higher-is-better. `value` is in native units (nM for +affinity, hours for stability). `percentile_rank` is always optional, +0-100, lower = stronger. + +## Supported predictors + +| Predictor | Kinds produced | Requires | +|---|---|---| +| `NetMHCpan` / `NetMHCpan41` | affinity + presentation | [NetMHCpan](http://www.cbs.dtu.dk/services/NetMHCpan/) | +| `NetMHCpan4` | affinity or presentation | NetMHCpan 4.0 | +| `NetMHCpan3` / `NetMHCpan28` | affinity | older NetMHCpan | +| `NetMHC` / `NetMHC3` / `NetMHC4` | affinity | [NetMHC](http://www.cbs.dtu.dk/services/NetMHC/) | +| `NetMHCIIpan` | affinity or presentation | [NetMHCIIpan](http://www.cbs.dtu.dk/services/NetMHCIIpan/) | +| `NetMHCcons` | affinity | [NetMHCcons](http://www.cbs.dtu.dk/services/NetMHCcons/) | +| `NetMHCstabpan` | stability | [NetMHCstabpan](http://www.cbs.dtu.dk/services/NetMHCstabpan/) | +| `MHCflurry` | affinity | `pip install mhcflurry` + `mhcflurry-downloads fetch` | +| `MixMHCpred` | presentation | [MixMHCpred](https://github.com/GfellerLab/MixMHCpred) | +| `RandomBindingPredictor` | affinity | (built-in) | +| `NetChop` | cleavage | [NetChop](http://www.cbs.dtu.dk/services/NetChop/) | + ## Commandline examples ### Prediction for user-supplied peptide sequences @@ -30,48 +199,27 @@ mhctools --sequence SIINFEKL SIINFEKLQ --mhc-predictor netmhc --mhc-alleles A020 mhctools --sequence AAAQQQSIINFEKL --extract-subsequences --mhc-peptide-lengths 8-10 --mhc-predictor mhcflurry --mhc-alleles A0201 ``` -## Python usage +## Legacy API + +The old `predict_peptides()` and `predict_subsequences()` methods still work +and return `BindingPredictionCollection` objects: ```python -from mhctools import NetMHCpan -# Run NetMHCpan for alleles HLA-A*01:01 and HLA-A*02:01 -predictor = NetMHCpan(alleles=["A*02:01", "hla-a0101"]) - -# scan the short proteins 1L2Y and 1L3Y for epitopes -protein_sequences = { - "1L2Y": "NLYIQWLKDGGPSSGRPPPS", - "1L3Y": "ECDTINCERYNGQVCGGPGRGLCFCGKCRCHPGFEGSACQA" -} - -binding_predictions = predictor.predict_subsequences(protein_sequences, peptide_lengths=[9]) - -# flatten binding predictions into a Pandas DataFrame -df = binding_predictions.to_dataframe() - -# epitope collection is sorted by percentile rank -# of binding predictions -for binding_prediction in binding_predictions: - if binding_prediction.affinity < 100: - print("Strong binder: %s" % (binding_prediction,)) +predictor = NetMHCpan(alleles=["A*02:01"]) +collection = predictor.predict_subsequences( + {"1L2Y": "NLYIQWLKDGGPSSGRPPPS"}, + peptide_lengths=[9], +) +df = collection.to_dataframe() + +for bp in collection: + if bp.affinity < 100: + print("Strong binder: %s" % bp) ``` -## API - -The following MHC binding predictors are available in `mhctools`: - -- `MHCflurry`: open source predictor installed by default with `mhctools`, requires the user run `mhcflurry-downloads fetch` first to download MHCflurry models -- `NetMHC3`: requires locally installed version of [NetMHC 3.x](http://www.cbs.dtu.dk/services/NetMHC-3.4/) -- `NetMHC4`: requires locally installed version of [NetMHC 4.x](http://www.cbs.dtu.dk/services/NetMHC/) -- `NetMHC`: a wrapper function to automatically use `NetMHC3` or `NetMHC4` depending on what's installed. -- `NetMHCpan`: requires locally installed version of [NetMHCpan](http://www.cbs.dtu.dk/services/NetMHCpan/) -- `NetMHCIIpan`: requires locally installed version of [NetMHCIIpan](http://www.cbs.dtu.dk/services/NetMHCIIpan/) -- `NetMHCcons`: requires locally installed version of [NetMHCcons](http://www.cbs.dtu.dk/services/NetMHCcons/) -- `IedbMhcClass1`: Uses IEDB's REST API for class I binding predictions. -- `IedbMhcClass2`: Uses IEDB's REST API for class II binding predictions. -- `RandomBindingPredictor`: Creates binding predictions with random IC50 and percentile rank values. - -Every binding predictor is constructed with an `alleles` argument specifying the HLA type for which to make predictions. Predictions are generated by calling the `predict` method with a dictionary mapping sequence IDs or names to amino acid sequences. +To convert legacy results to the new types: -Additionally there is a module for running the [NetChop](http://www.cbs.dtu.dk/services/NetChop) proteosomal cleavage predictor: - -- `NetChop`: requires locally installed version of [NetChop-3.1](http://www.cbs.dtu.dk/services/NetChop/) +```python +preds = collection.to_preds() # list of Pred +pp_list = collection.to_peptide_preds() # list of PeptidePreds +``` diff --git a/mhctools/__init__.py b/mhctools/__init__.py index 5469a39..03fbe4b 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -1,5 +1,7 @@ from .binding_prediction import BindingPrediction from .binding_prediction_collection import BindingPredictionCollection +from .pred import Pred, PeptidePreds, Kind, preds_from_rows +from .sample import MultiSample from .iedb import ( IedbNetMHCcons, IedbNetMHCpan, @@ -27,6 +29,11 @@ __version__ = "2.2.0" __all__ = [ + "Pred", + "PeptidePreds", + "Kind", + "preds_from_rows", + "MultiSample", "BindingPrediction", "BindingPredictionCollection", "IedbNetMHCcons", diff --git a/mhctools/base_commandline_predictor.py b/mhctools/base_commandline_predictor.py index 44c44f4..d4b1700 100644 --- a/mhctools/base_commandline_predictor.py +++ b/mhctools/base_commandline_predictor.py @@ -27,7 +27,9 @@ from .cleanup_context import CleanupFiles from .input_file_formats import create_input_peptides_files from .process_helpers import run_multiple_commands_redirect_stdout +from .binding_prediction import BindingPrediction from .binding_prediction_collection import BindingPredictionCollection +from .pred import Pred, PeptidePreds logger = logging.getLogger(__name__) @@ -53,7 +55,8 @@ def __init__( default_peptide_lengths=[9], group_peptides_by_length=False, min_peptide_length=8, - max_peptide_length=None,): + max_peptide_length=None, + parse_to_preds_fn=None,): """ Parameters ---------- @@ -146,6 +149,7 @@ def __init__( self.process_limit = process_limit self.parse_output_fn = parse_output_fn + self.parse_to_preds_fn = parse_to_preds_fn if isinstance(default_peptide_lengths, int): default_peptide_lengths = [default_peptide_lengths] @@ -284,6 +288,89 @@ def _run_commands_and_collect_predictions( logger.warning("No binding predictions from %s" % self.program_name) return BindingPredictionCollection(binding_predictions) + def _run_commands_and_collect_preds( + self, + commands, + input_filenames, + temp_dir_list, + sequence_key_mapping=None): + """Like _run_commands_and_collect_predictions but returns list[Pred].""" + if sequence_key_mapping is None: + sequence_key_mapping = defaultdict(lambda: "seq") + all_preds = [] + + filenames_to_delete = list(input_filenames) + [ + f.name for f in commands.keys()] + with CleanupFiles( + filenames=filenames_to_delete, + directories=temp_dir_list): + run_multiple_commands_redirect_stdout( + commands, + print_commands=True, + process_limit=self.process_limit) + for output_file, command in commands.items(): + output_file.close() + with open(output_file.name, 'r') as f: + file_contents = f.read() + all_preds.extend( + self.parse_to_preds_fn( + stdout=file_contents, + sequence_key_mapping=sequence_key_mapping, + predictor_name=self.program_name)) + + if len(all_preds) == 0: + logger.warning("No predictions from %s" % self.program_name) + + # Group by (peptide, offset, source) into PeptidePreds + groups = defaultdict(list) + for pred in all_preds: + key = (pred.peptide, pred.offset, pred.source_sequence_name) + groups[key].append(pred) + return [PeptidePreds(preds=tuple(preds)) for preds in groups.values()] + + def predict(self, peptides): + """ + Predict for a list of peptide sequences. + + Returns list of PeptidePreds. When a native parse_to_preds_fn is + available, parses directly to Pred objects. Otherwise falls back + to converting from BindingPrediction. + """ + if self.parse_to_preds_fn is None: + return super().predict(peptides) + + self._check_peptide_inputs(peptides) + input_filenames = create_input_peptides_files( + peptides, + max_peptides_per_file=self.max_peptides_per_file, + group_by_length=self.group_peptides_by_length) + commands = {} + dirs = [] + + for i, input_filename in enumerate(input_filenames): + for j, allele in enumerate(self.alleles): + if self.tempdir_flag: + temp_dirname = tempfile.mkdtemp( + prefix="tmp_%d_%d_%s" % (i, j, self.program_name), + suffix="XXXXXX") + dirs.append(temp_dirname) + else: + temp_dirname = None + output_file = tempfile.NamedTemporaryFile( + "w", + prefix="%s_output_length_%d_%d" % ( + self.program_name, i, j), + delete=False) + commands[output_file] = self._build_command( + input_filename=input_filename, + allele=allele, + peptide_mode=True, + temp_dirname=temp_dirname) + return self._run_commands_and_collect_preds( + commands=commands, + input_filenames=input_filenames, + temp_dir_list=dirs) + def predict_peptides(self, peptides): self._check_peptide_inputs(peptides) input_filenames = create_input_peptides_files( diff --git a/mhctools/base_predictor.py b/mhctools/base_predictor.py index 5a5b60b..1b686cc 100644 --- a/mhctools/base_predictor.py +++ b/mhctools/base_predictor.py @@ -18,6 +18,7 @@ from .unsupported_allele import UnsupportedAllele from .binding_prediction_collection import BindingPredictionCollection +from .pred import Pred, Kind, PeptidePreds logger = logging.getLogger(__name__) @@ -86,14 +87,117 @@ def __str__(self): self.alleles, self.default_peptide_lengths) + # --- new API --- + + def predict(self, peptides): + """ + Predict for a list of peptide sequences. + + Returns + ------- + list of PeptidePreds + """ + collection = self.predict_peptides(peptides) + return collection.to_peptide_preds(kind=self._default_pred_kind()) + + def predict_dataframe(self, peptides, sample_name=""): + """predict() flattened to a DataFrame.""" + import pandas as pd + dfs = [pp.to_dataframe(sample_name) for pp in self.predict(peptides)] + if not dfs: + from .pred import COLUMNS + return pd.DataFrame(columns=COLUMNS) + return pd.concat(dfs, ignore_index=True) + + def predict_proteins(self, sequence_dict, peptide_lengths=None): + """ + Scan protein sequences and predict for all subsequences. + + Parameters + ---------- + sequence_dict : dict or str + Mapping of sequence names to amino acid strings. + If a string, treated as {"seq": string}. + + peptide_lengths : list of int, optional + + Returns + ------- + dict mapping sequence_name -> list of PeptidePreds + """ + if isinstance(sequence_dict, str): + sequence_dict = {"seq": sequence_dict} + elif isinstance(sequence_dict, (list, tuple)): + sequence_dict = {seq: seq for seq in sequence_dict} + + peptide_lengths = self._check_peptide_lengths(peptide_lengths) + + peptide_set = set() + peptide_to_name_offset_pairs = defaultdict(list) + + for name, sequence in sequence_dict.items(): + for peptide_length in peptide_lengths: + for i in range(len(sequence) - peptide_length + 1): + peptide = sequence[i:i + peptide_length] + peptide_set.add(peptide) + peptide_to_name_offset_pairs[peptide].append((name, i)) + + peptide_list = sorted(peptide_set) + flat_preds = self.predict(peptide_list) + + # Expand: each PeptidePreds may map to multiple (name, offset) positions + results = defaultdict(list) + for pp in flat_preds: + if not pp.preds: + continue + peptide = pp.preds[0].peptide + for name, offset in peptide_to_name_offset_pairs.get(peptide, []): + relocated = PeptidePreds(preds=tuple( + Pred( + kind=p.kind, + score=p.score, + peptide=p.peptide, + allele=p.allele, + value=p.value, + percentile_rank=p.percentile_rank, + source_sequence_name=name, + offset=offset, + predictor_name=p.predictor_name, + predictor_version=p.predictor_version, + ) for p in pp.preds + )) + results[name].append(relocated) + return dict(results) + + def predict_proteins_dataframe(self, sequence_dict, peptide_lengths=None, sample_name=""): + """predict_proteins() flattened to a DataFrame.""" + import pandas as pd + dfs = [] + for name, pp_list in self.predict_proteins(sequence_dict, peptide_lengths).items(): + for pp in pp_list: + dfs.append(pp.to_dataframe(sample_name)) + if not dfs: + from .pred import COLUMNS + return pd.DataFrame(columns=COLUMNS) + return pd.concat(dfs, ignore_index=True) + + def _default_pred_kind(self): + """Override in subclasses to set the Kind for compat conversion.""" + return Kind.pMHC_affinity + + # --- deprecated API (still works) --- + def predict_peptides(self, peptides): """ - Given a list of peptide sequences, returns a BindingPredictionCollection + Deprecated: use predict() instead. + + Given a list of peptide sequences, returns a BindingPredictionCollection. """ raise NotImplementedError( "%s must implement predict_peptides" % (self.__class__.__name__,)) def predict_peptides_dataframe(self, peptides): + """Deprecated: use predict_dataframe() instead.""" return self.predict_peptides(peptides).to_dataframe() def _check_peptide_lengths(self, peptide_lengths=None): @@ -180,6 +284,8 @@ def predict_subsequences( sequence_dict, peptide_lengths=None): """ + Deprecated: use predict_proteins() instead. + Given a dictionary mapping sequence names to amino acid strings, and an optional list of peptide lengths, returns a BindingPredictionCollection. @@ -221,14 +327,11 @@ def predict_subsequences( alleles=self.alleles) return BindingPredictionCollection(results) - def predict(self, sequence_dict, peptide_lengths=None): - logger.warning("Deprecated method 'predict', use 'predict_subsequences") - return self.predict_subsequences(sequence_dict, peptide_lengths=None) - def predict_subsequences_dataframe( self, sequence_dict, peptide_lengths=None): + """Deprecated: use predict_proteins_dataframe() instead.""" return self.predict_subsequences( sequence_dict=sequence_dict, peptide_lengths=peptide_lengths).to_dataframe() diff --git a/mhctools/binding_prediction.py b/mhctools/binding_prediction.py index 930cc21..055431b 100644 --- a/mhctools/binding_prediction.py +++ b/mhctools/binding_prediction.py @@ -10,9 +10,13 @@ # See the License for the specific language governing permissions and # limitations under the License. +import warnings + import numpy as np from serializable import Serializable +from .pred import Pred, Kind + class BindingPrediction(Serializable): def __init__( self, @@ -158,3 +162,38 @@ def __hash__(self): def __lt__(self, other): return self.value < other.value + def to_pred(self, kind=Kind.pMHC_affinity): + """Convert to a Pred object. + + Parameters + ---------- + kind : Kind + What this prediction measures. Defaults to pMHC_affinity since + most BindingPrediction objects represent binding affinity. + """ + return Pred( + kind=kind, + score=self.score if self.score is not None else 0.0, + peptide=self.peptide, + allele=self.allele or "", + value=self.affinity, + percentile_rank=self.percentile_rank, + source_sequence_name=self.source_sequence_name, + offset=self.offset, + predictor_name=self.prediction_method_name, + ) + + @classmethod + def from_pred(cls, pred): + """Create a BindingPrediction from a Pred object.""" + return cls( + peptide=pred.peptide, + allele=pred.allele, + score=pred.score, + affinity=pred.value, + percentile_rank=pred.percentile_rank, + source_sequence_name=pred.source_sequence_name, + offset=pred.offset, + prediction_method_name=pred.predictor_name, + ) + diff --git a/mhctools/binding_prediction_collection.py b/mhctools/binding_prediction_collection.py index 18888fa..67d0a39 100644 --- a/mhctools/binding_prediction_collection.py +++ b/mhctools/binding_prediction_collection.py @@ -18,6 +18,7 @@ from sercol import Collection from .binding_prediction import BindingPrediction +from .pred import Kind, PeptidePreds class BindingPredictionCollection(Collection): def to_dataframe( @@ -29,3 +30,23 @@ def to_dataframe( return pd.DataFrame.from_records( [tuple([getattr(x, name) for name in columns]) for x in self], columns=columns) + + def to_preds(self, kind=Kind.pMHC_affinity): + """Convert all BindingPredictions to Pred objects. + + Returns a list of Pred (not grouped into PeptidePreds, since + BindingPredictionCollection has no peptide-position grouping). + """ + return [bp.to_pred(kind=kind) for bp in self] + + def to_peptide_preds(self, kind=Kind.pMHC_affinity): + """Convert to a list of PeptidePreds, grouped by (peptide, offset, source). + + Each PeptidePreds contains all alleles for one peptide position. + """ + from collections import defaultdict + groups = defaultdict(list) + for bp in self: + key = (bp.peptide, bp.offset, bp.source_sequence_name) + groups[key].append(bp.to_pred(kind=kind)) + return [PeptidePreds(preds=tuple(preds)) for preds in groups.values()] diff --git a/mhctools/netmhc_pan28.py b/mhctools/netmhc_pan28.py index 5ff9f31..e350d33 100644 --- a/mhctools/netmhc_pan28.py +++ b/mhctools/netmhc_pan28.py @@ -13,7 +13,7 @@ from __future__ import print_function, division, absolute_import from .base_commandline_predictor import BaseCommandlinePredictor -from .parsing import parse_netmhcpan28_stdout +from .parsing import parse_netmhcpan28_stdout, parse_netmhcpan_to_preds class NetMHCpan28(BaseCommandlinePredictor): def __init__( @@ -30,6 +30,7 @@ def __init__( default_peptide_lengths=default_peptide_lengths, group_peptides_by_length=True, parse_output_fn=parse_netmhcpan28_stdout, + parse_to_preds_fn=parse_netmhcpan_to_preds, supported_alleles_flag="-listMHC", input_file_flag="-f", length_flag="-l", diff --git a/mhctools/netmhc_pan3.py b/mhctools/netmhc_pan3.py index 8079405..f8dc527 100644 --- a/mhctools/netmhc_pan3.py +++ b/mhctools/netmhc_pan3.py @@ -12,7 +12,7 @@ from .base_commandline_predictor import BaseCommandlinePredictor -from .parsing import parse_netmhcpan3_stdout +from .parsing import parse_netmhcpan3_stdout, parse_netmhcpan_to_preds class NetMHCpan3(BaseCommandlinePredictor): def __init__( @@ -28,6 +28,7 @@ def __init__( alleles=alleles, default_peptide_lengths=default_peptide_lengths, parse_output_fn=parse_netmhcpan3_stdout, + parse_to_preds_fn=parse_netmhcpan_to_preds, supported_alleles_flag="-listMHC", input_file_flag="-f", length_flag="-l", diff --git a/mhctools/netmhc_pan4.py b/mhctools/netmhc_pan4.py index 0578a82..57cc4f5 100644 --- a/mhctools/netmhc_pan4.py +++ b/mhctools/netmhc_pan4.py @@ -13,7 +13,7 @@ from __future__ import print_function, division, absolute_import from .base_commandline_predictor import BaseCommandlinePredictor -from .parsing import parse_netmhcpan4_stdout, parse_netmhcpan41_stdout +from .parsing import parse_netmhcpan4_stdout, parse_netmhcpan41_stdout, parse_netmhcpan_to_preds from functools import partial class NetMHCpan4(BaseCommandlinePredictor): @@ -46,6 +46,7 @@ def __init__( alleles=alleles, default_peptide_lengths=default_peptide_lengths, parse_output_fn=partial(parse_netmhcpan4_stdout, mode=mode), + parse_to_preds_fn=parse_netmhcpan_to_preds, supported_alleles_flag="-listMHC", input_file_flag="-f", length_flag="-l", diff --git a/mhctools/netmhc_pan41.py b/mhctools/netmhc_pan41.py index a4b46f3..16e6355 100644 --- a/mhctools/netmhc_pan41.py +++ b/mhctools/netmhc_pan41.py @@ -11,7 +11,7 @@ # limitations under the License. from .base_commandline_predictor import BaseCommandlinePredictor -from .parsing import parse_netmhc41_stdout +from .parsing import parse_netmhc41_stdout, parse_netmhcpan_to_preds from functools import partial @@ -45,6 +45,7 @@ def __init__( alleles=alleles, default_peptide_lengths=default_peptide_lengths, parse_output_fn=partial(parse_netmhc41_stdout, mode=mode), + parse_to_preds_fn=parse_netmhcpan_to_preds, supported_alleles_flag="-listMHC", input_file_flag="-f", length_flag="-l", diff --git a/mhctools/parsing.py b/mhctools/parsing.py index fadd44b..ec586d3 100644 --- a/mhctools/parsing.py +++ b/mhctools/parsing.py @@ -20,6 +20,7 @@ from .allele_normalization import normalize_allele_name from .binding_prediction import BindingPrediction +from .pred import Pred, Kind logger = logging.getLogger(__name__) @@ -523,35 +524,54 @@ def _detect_netmhcpan_version_label(stdout, header_fields): return "NetMHCpan (unknown version)" -def parse_netmhcpan_stdout( +def _safe_float(fields, index): + """Extract a float from fields at index, or None if index is None.""" + if index is None: + return None + return float(fields[index]) + + +def _affinity_score(ic50, raw_score): + """Compute a higher-is-better score (~0-1) from IC50 or raw log score. + + If IC50 is valid, compute 1 - log(IC50)/log(50000). + Otherwise fall back to raw_score (already higher-is-better for NetMHC tools). + """ + if ic50 is not None and valid_affinity(ic50): + return 1.0 - (np.log(ic50) / np.log(50000)) + if raw_score is not None and np.isfinite(raw_score): + return raw_score + return 0.0 + + +def parse_netmhcpan_to_preds( stdout, - prediction_method_name="netmhcpan", - sequence_key_mapping=None, - mode=None): + predictor_name="netMHCpan", + predictor_version="", + sequence_key_mapping=None): """ - Auto-detecting parser for NetMHCpan output of any supported version - (2.8, 3.0, 4.0, 4.1). + Auto-detecting parser for NetMHCpan output. Returns list[Pred]. + + Reads the header line between dash separators to determine column + positions. Supports NetMHCpan 2.8, 3.0, 4.0, and 4.1. - Parses the header line between dash separators to determine column - positions, then extracts binding predictions accordingly. + For NetMHCpan 4.1 (which has both EL and BA columns), emits both + a pMHC_affinity Pred and a pMHC_presentation Pred per data row. Parameters ---------- stdout : str - Raw stdout from any version of NetMHCpan. - prediction_method_name : str + predictor_name : str - sequence_key_mapping : dict or None + predictor_version : str + If empty, auto-detected from stdout. - mode : str or None - One of "binding_affinity", "elution_score", or None. - Only relevant when the output contains both EL and BA columns - (NetMHCpan 4.1). Defaults to "binding_affinity". + sequence_key_mapping : dict or None Returns ------- - list of BindingPrediction + list of Pred """ check_stdout_error(stdout, "NetMHCpan") @@ -566,67 +586,170 @@ def parse_netmhcpan_stdout( version_label = _detect_netmhcpan_version_label(stdout, header_fields) logger.info("Detected %s output format", version_label) - # Position column (case varies across versions) + if not predictor_version: + match = re.search(r'# NetMHCpan version (\S+)', stdout) + predictor_version = match.group(1) if match else "" + + # --- locate columns from header --- + offset_index = field_index.get('Pos', field_index.get('pos')) if offset_index is None: - raise ValueError("No position column found in header: %s" % header_fields) + raise ValueError("No position column in header: %s" % header_fields) - # Allele column allele_index = field_index.get('HLA', field_index.get('MHC')) if allele_index is None: - raise ValueError("No allele column found in header: %s" % header_fields) + raise ValueError("No allele column in header: %s" % header_fields) - # Peptide column peptide_index = field_index.get('Peptide', field_index.get('peptide')) if peptide_index is None: - raise ValueError("No peptide column found in header: %s" % header_fields) + raise ValueError("No peptide column in header: %s" % header_fields) - # Identity/key column key_index = field_index.get('Identity') if key_index is None: - raise ValueError("No Identity column found in header: %s" % header_fields) + raise ValueError("No Identity column in header: %s" % header_fields) - # Versions with a Core column (3.0+) use 1-based offsets has_core = 'Core' in field_index - transforms = {} - if has_core: - transforms[offset_index] = lambda x: int(x) - 1 - - # Determine score, ic50, rank columns from what is present - if 'Score_EL' in field_index: - # NetMHCpan 4.1 format with separate EL and BA columns - effective_mode = mode or "binding_affinity" - if effective_mode == "binding_affinity" and 'Score_BA' in field_index: - score_index = field_index['Score_BA'] - rank_index = field_index.get('%Rank_BA') - ic50_index = field_index.get('Aff(nM)') + offset_is_one_based = has_core + + # Detect which column groups are available + has_el = 'Score_EL' in field_index + has_ba_separate = 'Score_BA' in field_index + has_log_score = '1-log50k(aff)' in field_index + has_plain_score = 'Score' in field_index + + # --- parse data rows --- + + preds = [] + for fields in split_stdout_lines(stdout): + # Strip optional trailing bind-level tokens (<=, WB, SB) + # These appear after the last numeric column and shift nothing + + offset = int(fields[offset_index]) + if offset_is_one_based: + offset -= 1 + + peptide = str(fields[peptide_index]) + allele = normalize_allele_name(str(fields[allele_index])) + + key = str(fields[key_index]) + if sequence_key_mapping: + key = sequence_key_mapping.get(key, key) + + shared = dict( + peptide=peptide, + allele=allele, + source_sequence_name=key, + offset=offset, + predictor_name=predictor_name, + predictor_version=predictor_version, + ) + + if has_el and has_ba_separate: + # NetMHCpan 4.1: emit BOTH affinity and presentation Preds + el_score = _safe_float(fields, field_index.get('Score_EL')) + el_rank = _safe_float(fields, field_index.get('%Rank_EL')) + ba_score = _safe_float(fields, field_index.get('Score_BA')) + ba_rank = _safe_float(fields, field_index.get('%Rank_BA')) + ic50 = _safe_float(fields, field_index.get('Aff(nM)')) + + preds.append(Pred( + kind=Kind.pMHC_affinity, + score=_affinity_score(ic50, ba_score), + value=ic50, + percentile_rank=ba_rank, + **shared)) + preds.append(Pred( + kind=Kind.pMHC_presentation, + score=el_score if el_score is not None else 0.0, + percentile_rank=el_rank, + **shared)) + + elif has_el and not has_ba_separate: + # NetMHCpan 4.1 EL-only (no -BA flag) or 4.0 EL mode + el_score = _safe_float(fields, field_index.get('Score_EL', field_index.get('Score'))) + el_rank = _safe_float(fields, field_index.get('%Rank_EL', field_index.get('%Rank'))) + preds.append(Pred( + kind=Kind.pMHC_presentation, + score=el_score if el_score is not None else 0.0, + percentile_rank=el_rank, + **shared)) + + elif has_log_score: + # NetMHCpan 2.8: 1-log50k(aff), Affinity(nM), %Rank + raw_score = _safe_float(fields, field_index['1-log50k(aff)']) + ic50 = _safe_float(fields, field_index.get('Affinity(nM)')) + rank = _safe_float(fields, field_index.get('%Rank')) + preds.append(Pred( + kind=Kind.pMHC_affinity, + score=_affinity_score(ic50, raw_score), + value=ic50, + percentile_rank=rank, + **shared)) + + elif has_plain_score: + # NetMHCpan 3.0 or 4.0 BA mode: Score, Aff(nM), %Rank + raw_score = _safe_float(fields, field_index['Score']) + ic50 = _safe_float(fields, field_index.get('Aff(nM)')) + rank = _safe_float(fields, field_index.get('%Rank')) + + if ic50 is not None: + # BA mode — has affinity + preds.append(Pred( + kind=Kind.pMHC_affinity, + score=_affinity_score(ic50, raw_score), + value=ic50, + percentile_rank=rank, + **shared)) + else: + # 4.0 EL mode — Score and %Rank only, no Aff(nM) + preds.append(Pred( + kind=Kind.pMHC_presentation, + score=raw_score if raw_score is not None else 0.0, + percentile_rank=rank, + **shared)) else: - score_index = field_index['Score_EL'] - rank_index = field_index.get('%Rank_EL') - ic50_index = None - elif '1-log50k(aff)' in field_index: - # NetMHCpan 2.8 format - score_index = field_index['1-log50k(aff)'] - rank_index = field_index.get('%Rank') - ic50_index = field_index.get('Affinity(nM)') - else: - # NetMHCpan 3.0 or 4.0 format - score_index = field_index.get('Score') - rank_index = field_index.get('%Rank') - ic50_index = field_index.get('Aff(nM)') + raise ValueError( + "Could not determine score columns from header: %s" % header_fields) - return parse_stdout( + return preds + + +def parse_netmhcpan_stdout( + stdout, + prediction_method_name="netmhcpan", + sequence_key_mapping=None, + mode=None): + """ + Auto-detecting parser for NetMHCpan output. Returns list[BindingPrediction] + for backward compatibility. + + For new code, use parse_netmhcpan_to_preds() instead. + """ + preds = parse_netmhcpan_to_preds( stdout=stdout, - prediction_method_name=prediction_method_name, - sequence_key_mapping=sequence_key_mapping, - key_index=key_index, - offset_index=offset_index, - peptide_index=peptide_index, - allele_index=allele_index, - score_index=score_index, - rank_index=rank_index, - ic50_index=ic50_index, - transforms=transforms) + predictor_name=prediction_method_name, + sequence_key_mapping=sequence_key_mapping) + + # For compat, pick one Pred per data row based on mode + if mode == "elution_score": + keep_kind = Kind.pMHC_presentation + else: + keep_kind = Kind.pMHC_affinity + + # Group by row key, prefer keep_kind + from collections import defaultdict + by_row = defaultdict(list) + for pred in preds: + row_key = (pred.peptide, pred.allele, pred.offset, pred.source_sequence_name) + by_row[row_key].append(pred) + + results = [] + for row_preds in by_row.values(): + preferred = [p for p in row_preds if p.kind == keep_kind] + chosen = preferred[0] if preferred else row_preds[0] + results.append(BindingPrediction.from_pred(chosen)) + + return results def parse_netmhccons_stdout( diff --git a/mhctools/pred.py b/mhctools/pred.py new file mode 100644 index 0000000..b688af5 --- /dev/null +++ b/mhctools/pred.py @@ -0,0 +1,169 @@ +# 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. + +from __future__ import annotations + +from dataclasses import dataclass, field +from enum import Enum +from typing import Optional + +import pandas as pd + + +class Kind(Enum): + """What biological quantity is being predicted.""" + # Per-allele + pMHC_affinity = "pMHC_affinity" + pMHC_presentation = "pMHC_presentation" + pMHC_stability = "pMHC_stability" + # Cell-level + cellular_presentation = "cellular_presentation" + # Processing pathway + antigen_processing = "antigen_processing" + proteasome_cleavage = "proteasome_cleavage" + tap_transport = "tap_transport" + erap_trimming = "erap_trimming" + + +COLUMNS = ( + "sample_name", + "peptide", + "n_flank", + "c_flank", + "source_sequence_name", + "offset", + "predictor_name", + "predictor_version", + "allele", + "kind", + "score", + "value", + "percentile_rank", +) + + +@dataclass(frozen=True) +class Pred: + """Single prediction from one model on one peptide. Self-contained.""" + kind: Kind + score: float + peptide: str = "" + allele: str = "" + n_flank: str = "" + c_flank: str = "" + value: Optional[float] = None + percentile_rank: Optional[float] = None + source_sequence_name: Optional[str] = None + offset: int = 0 + predictor_name: str = "" + predictor_version: str = "" + + def to_row(self, sample_name=""): + return { + "sample_name": sample_name, + "peptide": self.peptide, + "n_flank": self.n_flank, + "c_flank": self.c_flank, + "source_sequence_name": self.source_sequence_name, + "offset": self.offset, + "predictor_name": self.predictor_name, + "predictor_version": self.predictor_version, + "allele": self.allele, + "kind": self.kind.value, + "score": self.score, + "value": self.value, + "percentile_rank": self.percentile_rank, + } + + +@dataclass +class PeptidePreds: + """All Preds for one peptide from one predictor.""" + preds: tuple[Pred, ...] = () + + # --- best allele by score (higher = better) --- + + @property + def best_affinity(self) -> Optional[Pred]: + return self._best_by_score(Kind.pMHC_affinity) + + @property + def best_presentation(self) -> Optional[Pred]: + return self._best_by_score(Kind.pMHC_presentation) + + @property + def best_stability(self) -> Optional[Pred]: + return self._best_by_score(Kind.pMHC_stability) + + # --- best allele by rank (lower = better) --- + + @property + def best_affinity_by_rank(self) -> Optional[Pred]: + return self._best_by_rank(Kind.pMHC_affinity) + + @property + def best_presentation_by_rank(self) -> Optional[Pred]: + return self._best_by_rank(Kind.pMHC_presentation) + + @property + def best_stability_by_rank(self) -> Optional[Pred]: + return self._best_by_rank(Kind.pMHC_stability) + + # --- filtering --- + + def filter(self, kind=None, allele=None): + """Filter preds. None means don't filter on that field.""" + return [p for p in self.preds + if (kind is None or p.kind == kind) + and (allele is None or p.allele == allele)] + + # --- dataframe --- + + def to_dataframe(self, sample_name=""): + rows = [p.to_row(sample_name) for p in self.preds] + if not rows: + return pd.DataFrame(columns=COLUMNS) + return pd.DataFrame(rows, columns=COLUMNS) + + # --- internals --- + + def _best_by_score(self, kind) -> Optional[Pred]: + candidates = [p for p in self.preds if p.kind == kind and p.allele] + return max(candidates, key=lambda p: p.score) if candidates else None + + def _best_by_rank(self, kind) -> Optional[Pred]: + candidates = [p for p in self.preds + if p.kind == kind and p.allele + and p.percentile_rank is not None] + return min(candidates, key=lambda p: p.percentile_rank) if candidates else None + + +def preds_from_rows(rows, **shared): + """Build a PeptidePreds from dicts, with shared fields filled in. + + Example:: + + preds_from_rows( + [ + dict(kind=Kind.pMHC_affinity, allele="HLA-A*02:01", + score=0.85, value=120.5, percentile_rank=0.8), + dict(kind=Kind.pMHC_presentation, allele="HLA-A*02:01", + score=0.92, percentile_rank=0.3), + ], + peptide="SIINFEKL", + predictor_name="netMHCpan", + predictor_version="4.1", + ) + """ + return PeptidePreds(preds=tuple( + Pred(**{**shared, **row}) for row in rows + )) diff --git a/mhctools/sample.py b/mhctools/sample.py new file mode 100644 index 0000000..a2f00ca --- /dev/null +++ b/mhctools/sample.py @@ -0,0 +1,94 @@ +# 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 pandas as pd + +from .pred import PeptidePreds, COLUMNS + + +class MultiSample: + """ + Run a predictor across multiple samples, each with its own alleles. + + Parameters + ---------- + samples : dict + Mapping of sample_name -> list of allele strings. + Example: {"pat001": ["HLA-A*02:01", "HLA-B*07:02"], + "pat002": ["HLA-A*01:01", "HLA-B*08:01"]} + + predictor_class : class + A predictor class (e.g. NetMHCpan41) that accepts an ``alleles`` + keyword argument. + + **predictor_kwargs + Additional keyword arguments forwarded to the predictor constructor. + """ + + def __init__(self, samples, predictor_class, **predictor_kwargs): + self.samples = samples + self.predictor_class = predictor_class + self.predictor_kwargs = predictor_kwargs + + def _make_predictor(self, alleles): + return self.predictor_class(alleles=alleles, **self.predictor_kwargs) + + # --- peptide predictions --- + + def predict(self, peptides): + """ + Returns + ------- + dict mapping sample_name -> list of PeptidePreds + """ + results = {} + for sample_name, alleles in self.samples.items(): + predictor = self._make_predictor(alleles) + results[sample_name] = predictor.predict(peptides) + return results + + def predict_dataframe(self, peptides): + """predict() flattened to a DataFrame with sample_name column.""" + dfs = [] + for sample_name, pp_list in self.predict(peptides).items(): + 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) + + # --- protein scanning --- + + def predict_proteins(self, sequence_dict, peptide_lengths=None): + """ + Returns + ------- + dict mapping sample_name -> {sequence_name: list of PeptidePreds} + """ + results = {} + for sample_name, alleles in self.samples.items(): + predictor = self._make_predictor(alleles) + results[sample_name] = predictor.predict_proteins( + sequence_dict, peptide_lengths=peptide_lengths) + return results + + def predict_proteins_dataframe(self, sequence_dict, peptide_lengths=None): + """predict_proteins() flattened to a DataFrame with sample_name column.""" + dfs = [] + for sample_name, seq_dict in self.predict_proteins( + sequence_dict, peptide_lengths).items(): + for seq_name, pp_list in seq_dict.items(): + 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) diff --git a/tests/test_known_class2_epitopes.py b/tests/test_known_class2_epitopes.py index c04bcd8..5d54bc0 100644 --- a/tests/test_known_class2_epitopes.py +++ b/tests/test_known_class2_epitopes.py @@ -18,7 +18,7 @@ def expect_binder(mhc_model, peptide): - prediction = mhc_model.predict(peptide)[0] + prediction = mhc_model.predict_subsequences(peptide)[0] if prediction.value: assert prediction.value < 500, "Expected %s to have IC50 < 500nM, got %s" % ( peptide, prediction) diff --git a/tests/test_mhc_formats.py b/tests/test_mhc_formats.py index db0e832..5b125f3 100644 --- a/tests/test_mhc_formats.py +++ b/tests/test_mhc_formats.py @@ -14,9 +14,11 @@ parse_netmhcpan28_stdout, parse_netmhcpan3_stdout, parse_netmhcpan_stdout, + parse_netmhcpan_to_preds, parse_netmhc3_stdout, parse_netmhc4_stdout, ) +from mhctools.pred import Kind def test_netmhc3_stdout(): """ @@ -345,3 +347,108 @@ def test_auto_detect_netmhcpan41_el(): assert abs(entry.score - 0.0100620) < 1e-6 # Score_EL assert abs(entry.percentile_rank - 6.723) < 0.01 # %Rank_EL assert entry.value is None # ic50 not used in EL mode + + +# --- Tests for parse_netmhcpan_to_preds (new-style Pred output) --- + +def test_to_preds_netmhcpan28(): + """New-style parser returns Pred objects for NetMHCpan 2.8.""" + output = """ + ---------------------------------------------------x + pos HLA peptide Identity 1-log50k(aff) Affinity(nM) %Rank BindLevel + ---------------------------------------------------------------------------- + 0 HLA-A*02:03 QQQQQYFPE id0 0.024 38534.25 50.00 + 11 HLA-A*02:03 HIIIASSSL id0 0.515 189.74 4.00 <= WB + """ + preds = parse_netmhcpan_to_preds(output) + assert len(preds) == 2 + for p in preds: + assert p.kind == Kind.pMHC_affinity + assert p.allele == "HLA-A*02:03" + last = [p for p in preds if p.peptide == "HIIIASSSL"][0] + assert abs(last.value - 189.74) < 0.01 + assert last.percentile_rank == 4.00 + assert last.score > 0.4 # 1 - log(189.74)/log(50000) + + +def test_to_preds_netmhcpan41_emits_both_kinds(): + """NetMHCpan 4.1 emits both affinity and presentation Preds per row.""" + output = """ +# NetMHCpan version 4.1b + +# Make both EL and BA predictions + +--------------------------------------------------------------------------------------------------------------------------- + Pos MHC Peptide Core Of Gp Gl Ip Il Icore Identity Score_EL %Rank_EL Score_BA %Rank_BA Aff(nM) BindLevel +--------------------------------------------------------------------------------------------------------------------------- + 1 HLA-A*02:01 SIINFEKL SII-NFEKL 0 0 0 3 1 SIINFEKL PEPLIST 0.0100620 6.723 0.110414 20.171 15140.42 +--------------------------------------------------------------------------------------------------------------------------- +""" + preds = parse_netmhcpan_to_preds(output) + assert len(preds) == 2 # one affinity + one presentation + + affinity = [p for p in preds if p.kind == Kind.pMHC_affinity] + presentation = [p for p in preds if p.kind == Kind.pMHC_presentation] + assert len(affinity) == 1 + assert len(presentation) == 1 + + aff = affinity[0] + assert aff.allele == "HLA-A*02:01" + assert aff.peptide == "SIINFEKL" + assert abs(aff.value - 15140.42) < 0.1 # IC50 in nM + assert abs(aff.percentile_rank - 20.171) < 0.01 + assert aff.score > 0 # higher-is-better transform of IC50 + + pres = presentation[0] + assert abs(pres.score - 0.0100620) < 1e-6 + assert abs(pres.percentile_rank - 6.723) < 0.01 + assert pres.value is None # presentation has no native-unit value + + +def test_to_preds_netmhcpan4_el(): + """NetMHCpan 4.0 EL mode (no Aff column) returns presentation Preds.""" + output = """ +# NetMHCpan version 4.0 +----------------------------------------------------------------------------------- + Pos HLA Peptide Core Of Gp Gl Ip Il Icore Identity Score %Rank BindLevel +----------------------------------------------------------------------------------- + 1 HLA-A*02:01 SIINFEKL SIINF-EKL 0 0 0 5 1 SIINFEKL PEPLIST 0.3456780 5.1230 +----------------------------------------------------------------------------------- +""" + preds = parse_netmhcpan_to_preds(output) + assert len(preds) == 1 + p = preds[0] + assert p.kind == Kind.pMHC_presentation # no Aff(nM) → EL mode + assert abs(p.score - 0.3456780) < 1e-6 + + +def test_to_preds_auto_detects_version(): + """Predictor version auto-detected from stdout.""" + output = """ +# NetMHCpan version 4.1b +--------------------------------------------------------------------------------------------------------------------------- + Pos MHC Peptide Core Of Gp Gl Ip Il Icore Identity Score_EL %Rank_EL Score_BA %Rank_BA Aff(nM) BindLevel +--------------------------------------------------------------------------------------------------------------------------- + 1 HLA-A*02:01 SIINFEKL SII-NFEKL 0 0 0 3 1 SIINFEKL PEPLIST 0.0100620 6.723 0.110414 20.171 15140.42 +--------------------------------------------------------------------------------------------------------------------------- +""" + preds = parse_netmhcpan_to_preds(output) + assert preds[0].predictor_version == "4.1b" + + +def test_to_preds_self_contained(): + """Each Pred carries full context — peptide, allele, source, offset.""" + output = """ + ---------------------------------------------------x + pos HLA peptide Identity 1-log50k(aff) Affinity(nM) %Rank BindLevel + ---------------------------------------------------------------------------- + 5 HLA-A*02:03 YFPEITHII id0 0.231 4123.85 15.00 + """ + preds = parse_netmhcpan_to_preds(output, predictor_name="netMHCpan", predictor_version="2.8") + p = preds[0] + assert p.peptide == "YFPEITHII" + assert p.allele == "HLA-A*02:03" + assert p.source_sequence_name == "id0" + assert p.offset == 5 + assert p.predictor_name == "netMHCpan" + assert p.predictor_version == "2.8" diff --git a/tests/test_pred.py b/tests/test_pred.py new file mode 100644 index 0000000..bf1f2a4 --- /dev/null +++ b/tests/test_pred.py @@ -0,0 +1,305 @@ +# 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. + +from mhctools.pred import Pred, PeptidePreds, Kind, preds_from_rows, COLUMNS +from mhctools.sample import MultiSample +from mhctools.binding_prediction import BindingPrediction +from mhctools.binding_prediction_collection import BindingPredictionCollection +from mhctools.random_predictor import RandomBindingPredictor + + +# -- Pred -- + +def test_pred_basic(): + p = Pred( + kind=Kind.pMHC_affinity, + score=0.85, + peptide="SIINFEKL", + allele="HLA-A*02:01", + value=120.5, + percentile_rank=0.8, + ) + assert p.score == 0.85 + assert p.value == 120.5 + assert p.peptide == "SIINFEKL" + assert p.allele == "HLA-A*02:01" + assert p.kind == Kind.pMHC_affinity + + +def test_pred_frozen(): + p = Pred(kind=Kind.pMHC_affinity, score=0.5) + try: + p.score = 0.9 + assert False, "should be frozen" + except AttributeError: + pass + + +def test_pred_to_row(): + p = Pred( + kind=Kind.pMHC_affinity, + score=0.85, + peptide="SIINFEKL", + allele="HLA-A*02:01", + ) + row = p.to_row(sample_name="pat001") + assert row["sample_name"] == "pat001" + assert row["peptide"] == "SIINFEKL" + assert row["kind"] == "pMHC_affinity" + assert set(row.keys()) == set(COLUMNS) + + +def test_pred_defaults(): + p = Pred(kind=Kind.antigen_processing, score=0.7) + assert p.peptide == "" + assert p.allele == "" + assert p.n_flank == "" + assert p.c_flank == "" + assert p.value is None + assert p.percentile_rank is None + + +# -- PeptidePreds -- + +def _make_pred_set(): + return preds_from_rows( + [ + dict(kind=Kind.pMHC_affinity, allele="HLA-A*02:01", + score=0.85, value=120.5, percentile_rank=0.8), + dict(kind=Kind.pMHC_affinity, allele="HLA-B*07:02", + score=0.42, value=5000.0, percentile_rank=15.0), + dict(kind=Kind.pMHC_presentation, allele="HLA-A*02:01", + score=0.92, percentile_rank=0.3), + dict(kind=Kind.pMHC_presentation, allele="HLA-B*07:02", + score=0.15, percentile_rank=12.0), + dict(kind=Kind.antigen_processing, score=0.85), + ], + peptide="SIINFEKL", + predictor_name="mhcflurry", + predictor_version="2.1", + ) + + +def test_preds_from_rows(): + ps = _make_pred_set() + assert len(ps.preds) == 5 + for p in ps.preds: + assert p.peptide == "SIINFEKL" + assert p.predictor_name == "mhcflurry" + + +def test_best_affinity(): + ps = _make_pred_set() + best = ps.best_affinity + assert best is not None + assert best.allele == "HLA-A*02:01" + assert best.score == 0.85 + + +def test_best_affinity_by_rank(): + ps = _make_pred_set() + best = ps.best_affinity_by_rank + assert best is not None + assert best.allele == "HLA-A*02:01" + assert best.percentile_rank == 0.8 + + +def test_best_presentation(): + ps = _make_pred_set() + best = ps.best_presentation + assert best is not None + assert best.allele == "HLA-A*02:01" + assert best.score == 0.92 + + +def test_best_presentation_by_rank(): + ps = _make_pred_set() + best = ps.best_presentation_by_rank + assert best is not None + assert best.allele == "HLA-A*02:01" + assert best.percentile_rank == 0.3 + + +def test_best_stability_empty(): + ps = _make_pred_set() + assert ps.best_stability is None + assert ps.best_stability_by_rank is None + + +def test_filter(): + ps = _make_pred_set() + affinity_preds = ps.filter(kind=Kind.pMHC_affinity) + assert len(affinity_preds) == 2 + + a1_preds = ps.filter(allele="HLA-A*02:01") + assert len(a1_preds) == 2 + + processing = ps.filter(kind=Kind.antigen_processing) + assert len(processing) == 1 + assert processing[0].allele == "" + + +def test_to_dataframe(): + ps = _make_pred_set() + df = ps.to_dataframe(sample_name="pat001") + assert list(df.columns) == list(COLUMNS) + assert len(df) == 5 + assert (df["sample_name"] == "pat001").all() + assert (df["peptide"] == "SIINFEKL").all() + + +def test_empty_pred_set(): + ps = PeptidePreds() + assert ps.best_affinity is None + df = ps.to_dataframe() + assert list(df.columns) == list(COLUMNS) + assert len(df) == 0 + + +# -- predict() on a predictor -- + +def test_predict_returns_peptide_preds(): + predictor = RandomBindingPredictor( + alleles=["HLA-A*02:01", "HLA-B*07:02"], + default_peptide_lengths=[9]) + results = predictor.predict(["SIINFEKLL", "GILGFVFTL"]) + assert isinstance(results, list) + assert all(isinstance(pp, PeptidePreds) for pp in results) + assert len(results) == 2 + for pp in results: + assert pp.best_affinity is not None + + +def test_predict_dataframe(): + predictor = RandomBindingPredictor( + alleles=["HLA-A*02:01"], + default_peptide_lengths=[9]) + df = predictor.predict_dataframe(["SIINFEKLL"], sample_name="pat001") + assert list(df.columns) == list(COLUMNS) + assert (df["sample_name"] == "pat001").all() + + +def test_predict_proteins(): + predictor = RandomBindingPredictor( + alleles=["HLA-A*02:01"], + default_peptide_lengths=[9]) + result = predictor.predict_proteins({"TP53": "SIINFEKLLAA"}) + assert "TP53" in result + assert isinstance(result["TP53"], list) + assert all(isinstance(pp, PeptidePreds) for pp in result["TP53"]) + for pp in result["TP53"]: + for pred in pp.preds: + assert pred.source_sequence_name == "TP53" + + +def test_predict_proteins_dataframe(): + predictor = RandomBindingPredictor( + alleles=["HLA-A*02:01"], + default_peptide_lengths=[9]) + df = predictor.predict_proteins_dataframe( + {"TP53": "SIINFEKLLAA"}, sample_name="pat001") + assert list(df.columns) == list(COLUMNS) + assert (df["source_sequence_name"] == "TP53").all() + + +# -- MultiSample -- + +def test_multi_sample_predict(): + ms = MultiSample( + samples={ + "pat001": ["HLA-A*02:01"], + "pat002": ["HLA-B*07:02"], + }, + predictor_class=RandomBindingPredictor, + default_peptide_lengths=[9], + ) + results = ms.predict(["SIINFEKLL"]) + assert "pat001" in results + assert "pat002" in results + assert all(isinstance(pp, PeptidePreds) for pp in results["pat001"]) + + +def test_multi_sample_predict_dataframe(): + ms = MultiSample( + samples={ + "pat001": ["HLA-A*02:01"], + "pat002": ["HLA-B*07:02"], + }, + predictor_class=RandomBindingPredictor, + default_peptide_lengths=[9], + ) + df = ms.predict_dataframe(["SIINFEKLL"]) + assert set(df["sample_name"]) == {"pat001", "pat002"} + + +# -- Compat layer -- + +def test_binding_prediction_to_pred(): + bp = BindingPrediction( + peptide="SIINFEKL", + allele="HLA-A*02:01", + affinity=200.0, + percentile_rank=0.3, + source_sequence_name="seq", + offset=5, + prediction_method_name="netMHCpan", + ) + pred = bp.to_pred() + assert pred.kind == Kind.pMHC_affinity + assert pred.peptide == "SIINFEKL" + assert pred.allele == "HLA-A*02:01" + assert pred.value == 200.0 + assert pred.percentile_rank == 0.3 + assert pred.source_sequence_name == "seq" + assert pred.offset == 5 + assert pred.predictor_name == "netMHCpan" + + +def test_binding_prediction_from_pred(): + pred = Pred( + kind=Kind.pMHC_affinity, + score=0.85, + peptide="SIINFEKL", + allele="HLA-A*02:01", + value=120.5, + percentile_rank=0.8, + predictor_name="mhcflurry", + ) + bp = BindingPrediction.from_pred(pred) + assert bp.peptide == "SIINFEKL" + assert bp.allele == "HLA-A*02:01" + assert bp.affinity == 120.5 + assert bp.percentile_rank == 0.8 + assert bp.score == 0.85 + assert bp.prediction_method_name == "mhcflurry" + + +def test_collection_to_preds(): + bps = BindingPredictionCollection([ + BindingPrediction(peptide="SIINFEKL", allele="HLA-A*02:01", affinity=200.0), + BindingPrediction(peptide="SIINFEKL", allele="HLA-B*07:02", affinity=5000.0), + ]) + preds = bps.to_preds() + assert len(preds) == 2 + assert all(isinstance(p, Pred) for p in preds) + + +def test_collection_to_peptide_preds(): + bps = BindingPredictionCollection([ + BindingPrediction(peptide="SIINFEKL", allele="HLA-A*02:01", affinity=200.0, offset=0), + BindingPrediction(peptide="SIINFEKL", allele="HLA-B*07:02", affinity=5000.0, offset=0), + BindingPrediction(peptide="GILGFVFTL", allele="HLA-A*02:01", affinity=50.0, offset=10), + ]) + pp_list = bps.to_peptide_preds() + assert len(pp_list) == 2 # two distinct peptide positions + sizes = sorted(len(pp.preds) for pp in pp_list) + assert sizes == [1, 2] From b5f4ee500d7ff312a5d26f75dc47203642f90789 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Tue, 7 Apr 2026 12:18:52 -0400 Subject: [PATCH 2/3] Switch from pylint to ruff, fix all lint errors Replace pylint (crashing on modern type annotations) with ruff. Fix pre-existing issues: bare excepts, ambiguous variable names, unused imports, None comparisons, type() comparisons. --- lint.sh | 10 ++------- mhctools/__init__.py | 3 +++ mhctools/base_commandline_predictor.py | 5 ++--- mhctools/binding_prediction.py | 2 -- mhctools/cleanup_context.py | 8 +++---- mhctools/iedb.py | 8 +++---- mhctools/input_file_formats.py | 2 +- mhctools/netmhc_pan4.py | 2 +- mhctools/parsing.py | 12 +++++------ mhctools/pred.py | 2 +- mhctools/sample.py | 2 +- pyproject.toml | 2 +- tests/test_mhcflurry.py | 2 +- tests/test_netmhcii_pan.py | 2 +- tests/test_netmhcii_pan32.py | 30 +++++++++++++++----------- 15 files changed, 46 insertions(+), 46 deletions(-) diff --git a/lint.sh b/lint.sh index 4e59af7..12f1465 100755 --- a/lint.sh +++ b/lint.sh @@ -1,12 +1,6 @@ #!/bin/bash set -o errexit -# getting false positives due to this issue with pylint: -# https://bitbucket.org/logilab/pylint/issues/701/false-positives-with-not-an-iterable-and +python -m ruff check mhctools tests -find mhctools tests -name '*.py' \ - | xargs python -m pylint \ - --errors-only \ - --disable=unsubscriptable-object,not-an-iterable - -echo 'Passes pylint check' +echo 'Passes ruff check' diff --git a/mhctools/__init__.py b/mhctools/__init__.py index 03fbe4b..d092f91 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -62,6 +62,9 @@ "NetMHCIIpan4", "NetMHCIIpan4_BA", "NetMHCIIpan4_EL", + "NetMHCIIpan43", + "NetMHCIIpan43_BA", + "NetMHCIIpan43_EL", "NetMHCstabpan", "RandomBindingPredictor", "UnsupportedAllele", diff --git a/mhctools/base_commandline_predictor.py b/mhctools/base_commandline_predictor.py index d4b1700..3d6f0e4 100644 --- a/mhctools/base_commandline_predictor.py +++ b/mhctools/base_commandline_predictor.py @@ -27,9 +27,8 @@ from .cleanup_context import CleanupFiles from .input_file_formats import create_input_peptides_files from .process_helpers import run_multiple_commands_redirect_stdout -from .binding_prediction import BindingPrediction from .binding_prediction_collection import BindingPredictionCollection -from .pred import Pred, PeptidePreds +from .pred import PeptidePreds logger = logging.getLogger(__name__) @@ -166,7 +165,7 @@ def __init__( # it's present try: run_command([self.program_name]) - except: + except Exception: raise SystemError("Failed to run %s" % self.program_name) valid_alleles = None diff --git a/mhctools/binding_prediction.py b/mhctools/binding_prediction.py index 055431b..d1556b9 100644 --- a/mhctools/binding_prediction.py +++ b/mhctools/binding_prediction.py @@ -10,8 +10,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -import warnings - import numpy as np from serializable import Serializable diff --git a/mhctools/cleanup_context.py b/mhctools/cleanup_context.py index eeb9274..7adb160 100644 --- a/mhctools/cleanup_context.py +++ b/mhctools/cleanup_context.py @@ -49,24 +49,24 @@ def __exit__(self, type, value, traceback): logger.debug("Cleaning up %s", name) try: f.close() - except: + except Exception: pass try: os.remove(f.name) - except: + except Exception: pass for name in self.filenames: logger.debug("Cleaning up %s", name) try: os.remove(name) - except: + except Exception: pass for dirname in self.directories: logger.debug("Removing directory %s", dirname) try: shutil.rmtree(dirname) - except: + except Exception: pass diff --git a/mhctools/iedb.py b/mhctools/iedb.py index 223cf08..a950114 100644 --- a/mhctools/iedb.py +++ b/mhctools/iedb.py @@ -85,7 +85,7 @@ def _parse_iedb_response(response): df = pd.read_csv(io.BytesIO(response), delim_whitespace=True, header=0) # pylint doesn't realize that df is a DataFrame, so tell is - assert type(df) == pd.DataFrame + assert isinstance(df, pd.DataFrame) df = pd.DataFrame(df) if len(df) == 0: @@ -218,9 +218,9 @@ def predict_subsequences(self, sequence_dict, peptide_lengths=None): normalized_alleles = [] for key, amino_acid_sequence in sequence_dict.items(): - for l in peptide_lengths: - for i in range(len(amino_acid_sequence) - l + 1): - expected_peptides.add(amino_acid_sequence[i:i + l]) + for plen in peptide_lengths: + for i in range(len(amino_acid_sequence) - plen + 1): + expected_peptides.add(amino_acid_sequence[i:i + plen]) self._check_peptide_inputs(expected_peptides) for allele in self.alleles: # IEDB MHCII predictor expects DRA1 to be omitted. diff --git a/mhctools/input_file_formats.py b/mhctools/input_file_formats.py index d65085d..4b1984b 100644 --- a/mhctools/input_file_formats.py +++ b/mhctools/input_file_formats.py @@ -30,7 +30,7 @@ def create_input_peptides_files( """ if group_by_length: peptide_lengths = {len(p) for p in peptides} - peptide_groups = {l: [] for l in peptide_lengths} + peptide_groups = {plen: [] for plen in peptide_lengths} for p in peptides: peptide_groups[len(p)].append(p) else: diff --git a/mhctools/netmhc_pan4.py b/mhctools/netmhc_pan4.py index 57cc4f5..3954556 100644 --- a/mhctools/netmhc_pan4.py +++ b/mhctools/netmhc_pan4.py @@ -13,7 +13,7 @@ from __future__ import print_function, division, absolute_import from .base_commandline_predictor import BaseCommandlinePredictor -from .parsing import parse_netmhcpan4_stdout, parse_netmhcpan41_stdout, parse_netmhcpan_to_preds +from .parsing import parse_netmhcpan4_stdout, parse_netmhcpan_to_preds from functools import partial class NetMHCpan4(BaseCommandlinePredictor): diff --git a/mhctools/parsing.py b/mhctools/parsing.py index ec586d3..cc7ab8d 100644 --- a/mhctools/parsing.py +++ b/mhctools/parsing.py @@ -54,23 +54,23 @@ def split_stdout_lines(stdout): # all the NetMHC formats use lines full of dashes before any actual # binding results seen_dash = False - for l in stdout.split("\n"): - l = l.strip() + for line in stdout.split("\n"): + line = line.strip() # wait for a line like '----------' before trying to parse entries # have to include multiple dashes here since NetMHC 4.0 sometimes # gives negative positions in its "peptide" input mode - if l.startswith("---"): + if line.startswith("---"): seen_dash = True continue if not seen_dash: continue # ignore empty lines and comments - if not l or l.startswith("#"): + if not line or line.startswith("#"): continue # beginning of headers in NetMHC - if any(l.startswith(word) for word in NETMHC_TOKENS): + if any(line.startswith(word) for word in NETMHC_TOKENS): continue - yield l.split() + yield line.split() def clean_fields(fields, ignored_value_indices, transforms): diff --git a/mhctools/pred.py b/mhctools/pred.py index b688af5..438f353 100644 --- a/mhctools/pred.py +++ b/mhctools/pred.py @@ -12,7 +12,7 @@ from __future__ import annotations -from dataclasses import dataclass, field +from dataclasses import dataclass from enum import Enum from typing import Optional diff --git a/mhctools/sample.py b/mhctools/sample.py index a2f00ca..4217109 100644 --- a/mhctools/sample.py +++ b/mhctools/sample.py @@ -12,7 +12,7 @@ import pandas as pd -from .pred import PeptidePreds, COLUMNS +from .pred import COLUMNS class MultiSample: diff --git a/pyproject.toml b/pyproject.toml index 8383a4c..2737c63 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,7 +36,7 @@ dependencies = [ [project.optional-dependencies] dev = [ "build", - "pylint>=3.0", + "ruff", "pytest", "pytest-cov", "twine", diff --git a/tests/test_mhcflurry.py b/tests/test_mhcflurry.py index 91244df..098b7f2 100644 --- a/tests/test_mhcflurry.py +++ b/tests/test_mhcflurry.py @@ -10,7 +10,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -import sys + from .common import eq_ from numpy import testing diff --git a/tests/test_netmhcii_pan.py b/tests/test_netmhcii_pan.py index 0882c3c..60c1e2b 100644 --- a/tests/test_netmhcii_pan.py +++ b/tests/test_netmhcii_pan.py @@ -140,7 +140,7 @@ def test_netmhcii_pan(): # TODO: add more actual tests to this def check_netmhcii_pan(program_name, fail_if_no_such_program=True): try: - predictor = NetMHCIIpan( + NetMHCIIpan( alleles=[DEFAULT_ALLELE], program_name=program_name) except FileNotFoundError: if fail_if_no_such_program: diff --git a/tests/test_netmhcii_pan32.py b/tests/test_netmhcii_pan32.py index 577e41e..e935ca0 100644 --- a/tests/test_netmhcii_pan32.py +++ b/tests/test_netmhcii_pan32.py @@ -2,20 +2,26 @@ def test_netmhciipan43(): - predictor=NetMHCIIpan43(alleles=['DRB1_0101']) - predictions=predictor.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe() - assert predictions.shape==(2, 9) - assert all([pp==None for pp in predictions.affinity])==True ## no binding predictions + predictor = NetMHCIIpan43(alleles=['DRB1_0101']) + predictions = predictor.predict_subsequences( + ["AAGAARIIAVDINKD", "AAGIVESVGEGVTTV"]).to_dataframe() + assert predictions.shape == (2, 9) + assert all(pp is None for pp in predictions.affinity) + def test_netmhciipan43_ba(): - predictor_ba=NetMHCIIpan43_BA(alleles=['DRB1_0101',"HLA-DQA1*05:11-DQB1*03:02"]) - binding_predictions=predictor_ba.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe() - assert binding_predictions.shape==(4, 9) - assert all([pp==None for pp in binding_predictions.affinity])==False # output should preturn binding predictions + predictor_ba = NetMHCIIpan43_BA( + alleles=['DRB1_0101', "HLA-DQA1*05:11-DQB1*03:02"]) + binding_predictions = predictor_ba.predict_subsequences( + ["AAGAARIIAVDINKD", "AAGIVESVGEGVTTV"]).to_dataframe() + assert binding_predictions.shape == (4, 9) + assert not all(pp is None for pp in binding_predictions.affinity) def test_netmhciipan43_el(): - predictor_el=NetMHCIIpan43_EL(alleles=['DRB1_0101',"HLA-DQA1*05:11-DQB1*03:02"]) - EL_predictions=predictor_el.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe() - assert EL_predictions.shape==(4, 9) - assert all([pp==None for pp in EL_predictions.affinity])==True # no affinity predictions + predictor_el = NetMHCIIpan43_EL( + alleles=['DRB1_0101', "HLA-DQA1*05:11-DQB1*03:02"]) + el_predictions = predictor_el.predict_subsequences( + ["AAGAARIIAVDINKD", "AAGIVESVGEGVTTV"]).to_dataframe() + assert el_predictions.shape == (4, 9) + assert all(pp is None for pp in el_predictions.affinity) From 068cb4dbf14fab035ec2275b9cd1dcf39bbf1f07 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Tue, 7 Apr 2026 12:28:39 -0400 Subject: [PATCH 3/3] Add NetMHCpan42 class, make factory auto-detect any version Add NetMHCpan42/NetMHCpan42_BA/NetMHCpan42_EL classes (same output format as 4.1). Rewrite NetMHCpan factory to parse the version string and look up in a version map. Unknown versions >= 4.1 fall back to the latest known class with the header-driven auto-detecting parser instead of crashing. --- mhctools/__init__.py | 4 ++ mhctools/netmhc_pan.py | 67 ++++++++++++++++++++++------- mhctools/netmhc_pan42.py | 93 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 148 insertions(+), 16 deletions(-) create mode 100644 mhctools/netmhc_pan42.py diff --git a/mhctools/__init__.py b/mhctools/__init__.py index d092f91..aff4444 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -21,6 +21,7 @@ from .netmhc_pan3 import NetMHCpan3 from .netmhc_pan4 import NetMHCpan4, NetMHCpan4_BA, NetMHCpan4_EL from .netmhc_pan41 import NetMHCpan41, NetMHCpan41_BA, NetMHCpan41_EL +from .netmhc_pan42 import NetMHCpan42, NetMHCpan42_BA, NetMHCpan42_EL from .netmhcii_pan import NetMHCIIpan, NetMHCIIpan3, NetMHCIIpan4, NetMHCIIpan4_BA, NetMHCIIpan4_EL, NetMHCIIpan43, NetMHCIIpan43_BA, NetMHCIIpan43_EL from .random_predictor import RandomBindingPredictor from .netmhcstabpan import NetMHCstabpan @@ -57,6 +58,9 @@ "NetMHCpan4_EL", "NetMHCpan41_BA", "NetMHCpan41_EL", + "NetMHCpan42", + "NetMHCpan42_BA", + "NetMHCpan42_EL", "NetMHCIIpan", "NetMHCIIpan3", "NetMHCIIpan4", diff --git a/mhctools/netmhc_pan.py b/mhctools/netmhc_pan.py index 1076d93..8c5f834 100644 --- a/mhctools/netmhc_pan.py +++ b/mhctools/netmhc_pan.py @@ -11,6 +11,7 @@ # limitations under the License. import logging +import re from subprocess import check_output import os @@ -18,9 +19,32 @@ from .netmhc_pan3 import NetMHCpan3 from .netmhc_pan4 import NetMHCpan4 from .netmhc_pan41 import NetMHCpan41 +from .netmhc_pan42 import NetMHCpan42 logger = logging.getLogger(__name__) +# Maps (major, minor) version tuple to the predictor class. +# Checked in order; the last entry is the fallback for any 4.x not +# explicitly listed (future-proofing). +_VERSION_MAP = [ + ((2, 8), NetMHCpan28), + ((3, 0), NetMHCpan3), + ((4, 0), NetMHCpan4), + ((4, 1), NetMHCpan41), + ((4, 2), NetMHCpan42), +] + + +def _parse_version(version_str): + """Parse '4.2c' into (4, 2). Returns None on failure.""" + # strip trailing letter suffixes like 'b', 'c' + cleaned = re.sub(r'[a-zA-Z]+$', '', version_str) + parts = cleaned.split('.') + try: + return (int(parts[0]), int(parts[1])) + except (IndexError, ValueError): + return None + def NetMHCpan( alleles, @@ -29,17 +53,23 @@ def NetMHCpan( default_peptide_lengths=[9], extra_flags=[]): """ - This function wraps NetMHCpan28 and NetMHCpan3 to automatically detect which class - to use, with the help of the miraculous and strange '--version' netmhcpan argument. + Auto-detecting wrapper for any installed version of NetMHCpan. + + Runs ``netMHCpan --version`` to detect the installed version and returns + the appropriate predictor class. For unrecognized versions >= 4.1, + falls back to the latest known class (which uses the header-driven + auto-detecting parser). """ - # convert to str since Python3 returns a `bytes` object. - # The '_MHCTOOLS_VERSION_SNIFFING' here is meaningless, but it is necessary - # to call `netmhcpan --version` with some argument, otherwise it hangs. with open(os.devnull, 'w') as devnull: output = check_output([ program_name, "--version", "_MHCTOOLS_VERSION_SNIFFING"], stderr=devnull) output_str = output.decode("ascii", "ignore") + + match = re.search(r'# NetMHCpan version (\S+)', output_str) + version_str = match.group(1) if match else "" + version_tuple = _parse_version(version_str) if version_str else None + common_kwargs = { "alleles": alleles, "default_peptide_lengths": default_peptide_lengths, @@ -47,14 +77,19 @@ def NetMHCpan( "process_limit": process_limit, "extra_flags": extra_flags, } - if "NetMHCpan version 2.8" in output_str: - return NetMHCpan28(**common_kwargs) - elif "NetMHCpan version 3.0" in output_str: - return NetMHCpan3(**common_kwargs) - elif "NetMHCpan version 4.0" in output_str: - return NetMHCpan4(**common_kwargs) - elif "NetMHCpan version 4.1" in output_str: - return NetMHCpan41(**common_kwargs) - else: - raise RuntimeError( - "This software expects NetMHCpan version 2.8, 3.0, or 4.0, or 4.1") + + # Exact match + if version_tuple: + for (major_minor, cls) in _VERSION_MAP: + if version_tuple == major_minor: + logger.info("Detected NetMHCpan %s, using %s", version_str, cls.__name__) + return cls(**common_kwargs) + + # Fallback: use the latest known class (header-driven parser handles + # any output format with Score_EL / Score_BA columns) + fallback_cls = _VERSION_MAP[-1][1] + logger.warning( + "NetMHCpan version %s not explicitly supported, falling back to %s " + "(header-driven auto-detecting parser)", + version_str or "unknown", fallback_cls.__name__) + return fallback_cls(**common_kwargs) diff --git a/mhctools/netmhc_pan42.py b/mhctools/netmhc_pan42.py new file mode 100644 index 0000000..1425b7a --- /dev/null +++ b/mhctools/netmhc_pan42.py @@ -0,0 +1,93 @@ +# 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. + +from .base_commandline_predictor import BaseCommandlinePredictor +from .parsing import parse_netmhcpan41_stdout, parse_netmhcpan_to_preds +from functools import partial + + +class NetMHCpan42(BaseCommandlinePredictor): + def __init__( + self, + alleles, + default_peptide_lengths=[9], + program_name="netMHCpan", + process_limit=-1, + mode="binding_affinity", + extra_flags=[]): + """ + Wrapper for NetMHCpan 4.2. + + Output format is identical to 4.1 (Score_EL, %Rank_EL, Score_BA, + %Rank_BA, Aff(nM) columns). + + The mode argument should be one of "binding_affinity" (default) or + "elution_score". + """ + if mode == "binding_affinity": + flags = ["-BA"] + elif mode == "elution_score": + flags = [] + else: + raise ValueError("Unsupported mode", mode) + + BaseCommandlinePredictor.__init__( + self, + program_name=program_name, + alleles=alleles, + default_peptide_lengths=default_peptide_lengths, + parse_output_fn=partial(parse_netmhcpan41_stdout, mode=mode), + parse_to_preds_fn=parse_netmhcpan_to_preds, + supported_alleles_flag="-listMHC", + input_file_flag="-f", + length_flag="-l", + allele_flag="-a", + extra_flags=flags + extra_flags, + process_limit=process_limit) + + +class NetMHCpan42_EL(NetMHCpan42): + """NetMHCpan 4.2 in elution score mode.""" + def __init__( + self, + alleles, + default_peptide_lengths=[9], + program_name="netMHCpan", + process_limit=-1, + extra_flags=[]): + NetMHCpan42.__init__( + self, + alleles=alleles, + default_peptide_lengths=default_peptide_lengths, + program_name=program_name, + process_limit=process_limit, + mode="elution_score", + extra_flags=extra_flags) + + +class NetMHCpan42_BA(NetMHCpan42): + """NetMHCpan 4.2 in binding affinity mode.""" + def __init__( + self, + alleles, + default_peptide_lengths=[9], + program_name="netMHCpan", + process_limit=-1, + extra_flags=[]): + NetMHCpan42.__init__( + self, + alleles=alleles, + default_peptide_lengths=default_peptide_lengths, + program_name=program_name, + process_limit=process_limit, + mode="binding_affinity", + extra_flags=extra_flags)