From 3a487caf3a3803eda9e305f6490c27fd285b5230 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Fri, 4 Sep 2026 09:10:43 -0400 Subject: [PATCH] Add complete MixMHCpred 3.0 support --- README.md | 67 ++++ mhctools/__init__.py | 12 +- mhctools/mixmhcpred.py | 801 ++++++++++++++++++++++++++++++++++----- pyproject.toml | 9 +- tests/test_mixmhcpred.py | 332 +++++++++++++++- 5 files changed, 1106 insertions(+), 115 deletions(-) diff --git a/README.md b/README.md index a1b0ad6..891f39a 100644 --- a/README.md +++ b/README.md @@ -302,6 +302,73 @@ affinity, hours for stability). `percentile_rank` is always optional, | `IedbNetMHCpan` / `IedbSMM` / `IedbNetMHCIIpan` | affinity | IEDB web API | | `RandomBindingPredictor` | affinity | (built-in) | +`MixMHCpred` 3.0 predicts **class-I presentation** for peptides of length +8-14. Version 3.0 adds pan-allele inference, MHC-I sequence alignment and +sequence-driven prediction, and optional binding-motif/peptide-length plots. +mhctools exposes all per-allele scores and percentile ranks through the +canonical prediction API. `predict_detailed` also retains MixMHCpred's raw +`Score_bestAllele`, `BestAllele`, and `%Rank_bestAllele` columns plus each +allele's closest training allele, sequence distance, and pan-allele status. + +MixMHCpred 3.0 is licensed for academic, non-commercial research and prohibits +redistribution without written permission, so its approximately 200 MB of +code, models, and reference data are not included in mhctools. Review the +[upstream license and installation guide](https://github.com/GfellerLab/MixMHCpred) +before downloading the official tagged release: + +```sh +git clone --branch v3.0 --depth 1 \ + https://github.com/GfellerLab/MixMHCpred.git +chmod +x MixMHCpred/MixMHCpred +export MIXMHCPRED_PATH="$PWD/MixMHCpred" +pip install "mhctools[mixmhcpred]" +``` + +The `mixmhcpred` extra installs the upstream Python dependencies. Sequence +alignment additionally needs the `mafft` executable. The upstream +`install_packages` script is another way to install both sets of dependencies. + +```python +from mhctools import MixMHCpred + +predictor = MixMHCpred( + alleles=["HLA-A*02:01", "HLA-A*01:02"], # A*01:02 uses v3 pan inference +) + +# Canonical mhctools output: one pMHC_presentation Prediction per allele. +results = predictor.predict(["SIINFEKL"]) +results[0].presentation.score + +# Complete native output and v3 quality/provenance metadata. +detailed = predictor.predict_detailed(["SIINFEKL"]) +detailed.table[["Score_bestAllele", "BestAllele", "%Rank_bestAllele"]] +detailed.allele_info[1].closest_training_allele +detailed.allele_info[1].distance +detailed.allele_info[1].pan_allele + +# Retain Binding_predictions.txt, PWM/PLD files and images, and the HTML view. +motifs = predictor.predict_detailed( + ["SIINFEKL"], output_dir="mixmhcpred-output", output_motifs=True) +motifs.artifacts.files + +# Align novel MHC-I sequences, then optionally predict and render their motifs. +sequence_result = predictor.predict_allele_sequences( + "unaligned-mhc-i.fasta", + peptides=["SIINFEKL"], + output_dir="mixmhcpred-sequence-output", + output_motifs=True, +) +sequence_result.aligned_sequences +sequence_result.table +sequence_result.allele_info[0].closest_database_allele +sequence_result.artifacts.files +``` + +Both artifact APIs require a new output path: the wrapper refuses an existing +path because MixMHCpred itself deletes and recreates its output directory. +`exclude_peptides_with_cysteine=True` is implemented by mhctools before the +external call, including under v3.0 where the legacy `-c` option was removed. + `MixMHC2pred` is a pan-allele **class-II** presentation predictor and a strong complement to `NetMHCIIpan` (independently co-best in the Frontiers in Immunology 2024 class-II benchmark). It emits one `pMHC_presentation` diff --git a/mhctools/__init__.py b/mhctools/__init__.py index 2bb7785..cb06dfb 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -24,7 +24,12 @@ IedbSMM_PMBEC, IedbNetMHCIIpan, ) -from .mixmhcpred import MixMHCpred +from .mixmhcpred import ( + MixMHCpred, + MixMHCpredAlleleInfo, + MixMHCpredArtifacts, + MixMHCpredResult, +) from .mixmhc2pred import MixMHC2pred from .prime import PRIME from .deeptap import DeepTAP @@ -88,7 +93,7 @@ def __getattr__(name): raise AttributeError( "module %r has no attribute %r" % (__name__, name)) -__version__ = "3.31.7" +__version__ = "3.31.8" __all__ = [ "Prediction", @@ -116,6 +121,9 @@ def __getattr__(name): "IedbSMM_PMBEC", "IedbNetMHCIIpan", "MixMHCpred", + "MixMHCpredAlleleInfo", + "MixMHCpredArtifacts", + "MixMHCpredResult", "MixMHC2pred", "PRIME", "DeepTAP", diff --git a/mhctools/mixmhcpred.py b/mhctools/mixmhcpred.py index 3aad510..8650add 100644 --- a/mhctools/mixmhcpred.py +++ b/mhctools/mixmhcpred.py @@ -10,129 +10,738 @@ # See the License for the specific language governing permissions and # limitations under the License. -import pandas as pd -from tempfile import mkdtemp, NamedTemporaryFile -from os.path import join, exists -from os import remove +"""Wrapper for MixMHCpred 3.0 class-I ligand presentation predictions. + +MixMHCpred 3.0 is a pan-allele predictor which can also align user-provided +MHC-I sequences and render binding-motif artifacts. Its academic license does +not permit redistribution, so mhctools shells out to a user-provided official +checkout and does not vendor MixMHCpred code, models, or reference data. + +The canonical :meth:`MixMHCpred.predict` API returns one +``Kind.pMHC_presentation`` prediction per peptide and allele. The +:meth:`MixMHCpred.predict_detailed` API additionally preserves MixMHCpred's +complete table, including its independently-computed best score, best allele, +best rank, closest training allele, distance score, and pan-allele status. +Sequence alignment, sequence-driven prediction, and motif/HTML generation are +available through :meth:`MixMHCpred.predict_allele_sequences`. + +Upstream: https://github.com/GfellerLab/MixMHCpred +""" -from .allele_normalization import normalize_allele_name +from __future__ import annotations + +import os +import re +import shutil +from dataclasses import dataclass +from pathlib import Path +from subprocess import STDOUT, CalledProcessError, check_output +from tempfile import TemporaryDirectory + +import pandas as pd -from .base_predictor import BasePredictor +from .allele_normalization import normalize_allele_name_or_raw +from .base_predictor import BasePredictor, _check_flank_inputs from .binding_prediction import BindingPrediction from .binding_prediction_collection import BindingPredictionCollection +from .pred import Kind, PeptideResult, Prediction from .process_helpers import run_command -from .cleanup_context import CleanupFiles + +_VERSION_RE = re.compile(r"MixMHCpred\s*(?:\(v)?([0-9]+(?:\.[0-9]+)+)") +_CLOSEST_RE = re.compile(r"(.+?)\s*\(([-+0-9.eE]+)\)\s*$") +_SEQUENCE_QUALITY_RE = re.compile( + r"

(.*?)

\s*" + r"

closest Allele = (.*?), distance = (.*?)

\s*" + r"

closest Allele from database = (.*?), distance = (.*?)

", + re.IGNORECASE | re.DOTALL, +) + + +@dataclass(frozen=True) +class MixMHCpredAlleleInfo: + """Quality/provenance information for one MixMHCpred output allele.""" + + allele: str + native_allele: str + closest_training_allele: str = "" + distance: float | None = None + pan_allele: bool = False + closest_database_allele: str = "" + database_distance: float | None = None + + +@dataclass(frozen=True) +class MixMHCpredArtifacts: + """Persistent files produced by a MixMHCpred 3.0 artifact run.""" + + output_dir: str + files: tuple + binding_predictions: str = "" + alignment: str = "" + overview_html: str = "" + + +@dataclass +class MixMHCpredResult: + """Complete parsed output from one MixMHCpred invocation. + + ``table`` is the exact upstream tabular output, retaining + ``Score_bestAllele``, ``BestAllele``, ``%Rank_bestAllele``, and every + allele-specific score/rank column. ``predictions`` maps those per-allele + columns into mhctools' canonical API. + """ + + table: pd.DataFrame + allele_info: tuple + version: str = "" + comments: tuple = () + stdout: str = "" + aligned_sequences: tuple = () + artifacts: MixMHCpredArtifacts | None = None + + @property + def predictions(self): + """Return the detailed table as canonical ``PeptideResult`` objects.""" + if self.table.empty: + return [] + + peptide_index = self.table.columns.get_loc("Peptide") + allele_columns = tuple( + ( + info, + self.table.columns.get_loc(f"Score_{info.native_allele}"), + self.table.columns.get_loc(f"%Rank_{info.native_allele}"), + ) + for info in self.allele_info + ) + results = [] + for row in self.table.itertuples(index=False, name=None): + peptide = str(row[peptide_index]) + preds = [] + for info, score_index, rank_index in allele_columns: + preds.append(Prediction( + kind=Kind.pMHC_presentation, + score=float(row[score_index]), + peptide=peptide, + allele=info.allele, + percentile_rank=float(row[rank_index]), + predictor_name="mixmhcpred", + predictor_version=self.version, + )) + results.append(PeptideResult(preds=tuple(preds))) + return results + + +def resolve_mixmhcpred_path(program_name=None): + """Resolve a MixMHCpred executable. + + Resolution order is an explicit ``program_name``, ``MIXMHCPRED_PATH``, + then ``MixMHCpred`` on ``PATH``. A path may name either the executable or + the root of an official checkout. + """ + candidate = program_name or os.environ.get("MIXMHCPRED_PATH") + if not candidate: + candidate = shutil.which("MixMHCpred") + if not candidate: + raise FileNotFoundError( + "MixMHCpred was not found. Download the official v3.0 release " + "from https://github.com/GfellerLab/MixMHCpred, accept its " + "academic/non-commercial license, and set MIXMHCPRED_PATH to " + "the checkout or executable.") + + path = Path(candidate).expanduser() + if path.is_dir(): + path = path / "MixMHCpred" + if not path.exists() and os.sep not in str(candidate): + resolved = shutil.which(str(candidate)) + if resolved: + path = Path(resolved) + if not path.is_file(): + raise FileNotFoundError( + f"MixMHCpred executable does not exist: {path}") + if " " in str(path.resolve()): + raise ValueError( + "MixMHCpred 3.0 does not support spaces in its installation path: " + f"{path}") + return str(path.resolve()) + + +def mixmhcpred_version(program_name): + """Return the version reported by ``MixMHCpred --help``.""" + program = resolve_mixmhcpred_path(program_name) + try: + output = check_output( + [program, "--help"], stderr=STDOUT, text=True) + except (CalledProcessError, OSError) as e: + raise RuntimeError( + f"Could not query MixMHCpred version from {program}: {e}") + match = _VERSION_RE.search(output) + return match.group(1) if match else "" + + +def _native_alleles_from_table(table): + return tuple( + column[len("Score_"):] + for column in table.columns + if column.startswith("Score_") and column != "Score_bestAllele" + ) + + +def _normalize_output_allele(allele): + try: + return normalize_allele_name_or_raw(allele) + except (TypeError, ValueError): + return allele + + +def _parse_comments(filename): + comments = [] + with open(filename) as f: + for line in f: + if line.startswith("#"): + comments.append(line.rstrip("\n")) + return tuple(comments) + + +def _comment_value(comments, labels): + for comment in comments: + text = comment.lstrip("#").strip() + for label in labels: + if text.lower().startswith(label.lower()): + return text[len(label):].strip() + return "" + + +def _sequence_quality_from_html(path): + if not path or not Path(path).is_file(): + return {} + text = Path(path).read_text() + result = {} + for match in _SEQUENCE_QUALITY_RE.finditer(text): + allele, closest, distance, database, database_distance = match.groups() + result[allele] = ( + closest, + float(distance), + database, + float(database_distance), + ) + return result + + +def parse_mixmhcpred_output( + filename, + alleles=None, + artifacts=None, + stdout="", + aligned_sequences=(), + sequence_mode=False): + """Parse the complete output of MixMHCpred 2.x or 3.0. + + The returned :class:`MixMHCpredResult` preserves the raw table and header + metadata. If ``alleles`` are supplied, their caller-facing spellings are + retained while native names are taken from the output columns. Set + ``sequence_mode=True`` for predictions generated from user-provided MHC-I + sequences so their pan-allele provenance is retained. + """ + comments = _parse_comments(filename) + table = pd.read_csv(filename, comment="#", sep="\t") + required = { + "Peptide", "Score_bestAllele", "BestAllele", "%Rank_bestAllele", + } + missing = required - set(table.columns) + if missing: + raise ValueError( + "Unexpected MixMHCpred output; missing columns " + f"{sorted(missing)} (columns: {list(table.columns)})") + + native_alleles = _native_alleles_from_table(table) + for native in native_alleles: + rank_column = f"%Rank_{native}" + if rank_column not in table.columns: + raise ValueError( + f"MixMHCpred output missing column {rank_column!r} " + f"(columns: {list(table.columns)})") + + if alleles is None: + public_alleles = tuple( + _normalize_output_allele(name) for name in native_alleles) + else: + public_alleles = tuple(alleles) + if len(public_alleles) != len(native_alleles): + raise ValueError( + f"MixMHCpred returned {len(native_alleles)} allele(s), " + f"expected {len(public_alleles)}: {native_alleles}") + + version = "" + for comment in comments: + match = _VERSION_RE.search(comment) + if match: + version = match.group(1) + break + + allele_comment = _comment_value( + comments, ("Alleles:", "Predictions for Alleles :")) + pan_native = set() + if " - predicted motif for " in allele_comment: + _, pan_text = allele_comment.split(" - predicted motif for ", 1) + pan_native = {x.strip() for x in pan_text.split(",") if x.strip()} + + closest_text = _comment_value( + comments, ("Closest Allele (distance score):",)) + closest = [] + if closest_text: + for item in closest_text.split(" -- "): + match = _CLOSEST_RE.match(item) + if not match: + raise ValueError( + "Could not parse MixMHCpred closest-allele metadata: " + f"{item!r}") + closest.append((match.group(1), float(match.group(2)))) + if len(closest) != len(native_alleles): + raise ValueError( + f"MixMHCpred returned {len(closest)} closest-allele entries " + f"for {len(native_alleles)} alleles") + + sequence_quality = _sequence_quality_from_html( + artifacts.overview_html if artifacts else "") + allele_info = [] + for i, (public, native) in enumerate(zip(public_alleles, native_alleles)): + closest_name = "" + distance = None + database_name = "" + database_distance = None + if closest: + closest_name, distance = closest[i] + elif native in sequence_quality: + (closest_name, distance, + database_name, database_distance) = sequence_quality[native] + allele_info.append(MixMHCpredAlleleInfo( + allele=public, + native_allele=native, + closest_training_allele=_normalize_output_allele(closest_name), + distance=distance, + pan_allele=sequence_mode or native in pan_native, + closest_database_allele=_normalize_output_allele(database_name), + database_distance=database_distance, + )) + + return MixMHCpredResult( + table=table, + allele_info=tuple(allele_info), + version=version, + comments=comments, + stdout=stdout, + aligned_sequences=tuple(aligned_sequences), + artifacts=artifacts, + ) + + +def parse_mixmhcpred_results(filename): + """Parse MixMHCpred output into legacy ``BindingPrediction`` objects.""" + result = parse_mixmhcpred_output(filename) + return [ + BindingPrediction( + peptide=str(peptide), + allele=_normalize_output_allele(allele), + score=float(score), + percentile_rank=float(percentile_rank), + prediction_method_name="mixmhcpred", + ) + for peptide, allele, score, percentile_rank in zip( + result.table["Peptide"], + result.table["BestAllele"], + result.table["Score_bestAllele"], + result.table["%Rank_bestAllele"], + ) + ] + + +def _validate_prediction_rows(result, peptides): + observed = result.table["Peptide"].astype(str).tolist() + if observed != list(peptides): + raise RuntimeError( + "MixMHCpred output peptides do not match its inputs. " + f"Expected {list(peptides)}, got {observed}") + return result + + +def _sequence_allele_info(names, artifacts=None): + quality = _sequence_quality_from_html( + artifacts.overview_html if artifacts else "") + result = [] + for name in names: + closest, distance, database, database_distance = quality.get( + name, ("", None, "", None)) + result.append(MixMHCpredAlleleInfo( + allele=name, + native_allele=name, + closest_training_allele=_normalize_output_allele(closest), + distance=distance, + pan_allele=True, + closest_database_allele=_normalize_output_allele(database), + database_distance=database_distance, + )) + return tuple(result) + + +def _read_fasta(path): + records = [] + name = None + sequence = [] + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + if line.startswith(">"): + if name is not None: + records.append((name, "".join(sequence))) + name = line[1:] + sequence = [] + elif name is None: + raise ValueError("Sequence data appears before the first FASTA header") + else: + sequence.append(line) + if name is not None: + records.append((name, "".join(sequence))) + return tuple(records) + + +def _write_fasta(path, records): + with open(path, "w") as f: + f.write("".join(f">{name}\n{sequence}\n" for name, sequence in records)) + + +def _collect_artifacts(output_dir): + root = Path(output_dir).resolve() + files = tuple(sorted( + str(path.resolve()) for path in root.rglob("*") if path.is_file())) + + def existing(name): + path = root / name + return str(path) if path.is_file() else "" + + return MixMHCpredArtifacts( + output_dir=str(root), + files=files, + binding_predictions=existing("Binding_predictions.txt"), + alignment=existing("final_alignment.fasta"), + overview_html=existing("Data_overview.html"), + ) + class MixMHCpred(BasePredictor): + """MixMHCpred class-I ligand-presentation predictor. + + Parameters + ---------- + alleles : sequence of str, optional + MHC-I alleles. Version 3.0 accepts its training alleles plus MHC-I + alleles present in its sequence database for pan-allele inference. + default_peptide_lengths : sequence of int + program_name : str, optional + Executable or official checkout directory. If omitted, resolve + ``MIXMHCPRED_PATH`` and then ``MixMHCpred`` on ``PATH``. + exclude_peptides_with_cysteine : bool + Exclude cysteine-containing peptides in mhctools before invocation. + This is version-independent; MixMHCpred 3.0 removed the old ``-c`` + option. + """ + + mhc_class = "I" + def __init__( self, - alleles, - default_peptide_lengths=[9], - program_name="MixMHCpred", + alleles=None, + default_peptide_lengths=(9,), + program_name=None, exclude_peptides_with_cysteine=False): - """ - Wrapper for MixMHCpred - - Parameters - ---------- - alleles : list of str - - default_peptide_lengths : list of int - - program_name : str - - exclude_peptides_with_cysteine : bool - If True then drop peptides which contain 'C' from predictions, - default is True. - """ + if isinstance(alleles, str): + alleles = alleles.split(",") + normalized = [] + for allele in alleles or (): + value = normalize_allele_name_or_raw(allele) + if value not in normalized: + normalized.append(value) BasePredictor.__init__( self, - alleles=alleles, + alleles=normalized, default_peptide_lengths=default_peptide_lengths, min_peptide_length=8, max_peptide_length=14, allow_X_in_peptides=False, - allow_lowercase_in_peptides=False) + allow_lowercase_in_peptides=False, + keep_unparseable_alleles=True, + ) + # BasePredictor historically de-duplicates through a set. Restore the + # caller's stable order so native output columns have deterministic + # identities. + self.alleles = normalized self.program_name = program_name self.exclude_peptides_with_cysteine = exclude_peptides_with_cysteine + self._version = None + + @property + def executable(self): + return resolve_mixmhcpred_path(self.program_name) + + @property + def version(self): + if self._version is None: + self._version = mixmhcpred_version(self.executable) + return self._version + + def _require_v3(self, feature): + version = self.version + if not version or int(version.split(".", 1)[0]) < 3: + raise RuntimeError( + f"{feature} requires MixMHCpred 3.0 or newer; detected " + f"{version or 'an unknown version'} at {self.executable}") + + @staticmethod + def _validate_output_dir(output_dir): + path = Path(output_dir).expanduser().resolve() + if path.exists(): + raise FileExistsError( + "Refusing to let MixMHCpred replace existing output path: " + f"{path}") + if " " in str(path): + raise ValueError( + "MixMHCpred 3.0 does not support spaces in output paths: " + f"{path}") + return str(path) + + def _run_command(self, args, stdout_path): + try: + with open(stdout_path, "w") as stdout_file: + run_command( + args, + suppress_stderr=False, + redirect_stdout_file=stdout_file, + ) + except (CalledProcessError, OSError) as e: + stdout = ( + Path(stdout_path).read_text() + if Path(stdout_path).exists() else "") + raise RuntimeError( + f"MixMHCpred failed: {e}\n{stdout.strip()}") from e + return Path(stdout_path).read_text() + + def predict_detailed(self, peptides, output_dir=None, output_motifs=False): + """Run allele-based prediction and retain every MixMHCpred output. + + Set ``output_motifs=True`` and provide a new ``output_dir`` to retain + v3.0's motif matrices/images, peptide-length distributions, HTML, and + binding-prediction table. The method refuses an existing output path + because upstream deletes its output directory. + """ + peptide_list, _, _ = _check_flank_inputs(peptides) + self._check_peptide_inputs(peptide_list) if self.exclude_peptides_with_cysteine: - self.extra_commandline_args = ["-c"] + peptide_list = [p for p in peptide_list if "C" not in p] + if not peptide_list or not self.alleles: + return MixMHCpredResult( + table=pd.DataFrame(columns=( + "Peptide", "Score_bestAllele", "BestAllele", + "%Rank_bestAllele")), + allele_info=(), + ) + if output_motifs and output_dir is None: + raise ValueError("output_dir is required when output_motifs=True") + if output_motifs: + self._require_v3("Motif generation") + + persistent_dir = ( + self._validate_output_dir(output_dir) if output_motifs else None) + with TemporaryDirectory(prefix="mhctools-mixmhcpred-") as temp_dir: + input_path = Path(temp_dir) / "peptides.txt" + stdout_path = Path(temp_dir) / "stdout.txt" + input_path.write_text("\n".join(peptide_list) + "\n") + output_path = ( + Path(persistent_dir) if output_motifs + else Path(temp_dir) / "predictions.txt") + args = [ + self.executable, + "-i", str(input_path), + "-o", str(output_path), + "-a", ",".join(self.alleles), + ] + if output_motifs: + args.extend(("-m", "1")) + stdout = self._run_command(args, stdout_path) + prediction_path = ( + output_path / "Binding_predictions.txt" + if output_motifs else output_path) + if not prediction_path.is_file(): + raise RuntimeError( + "MixMHCpred exited without creating predictions at " + f"{prediction_path}. stdout: {stdout.strip()}") + artifacts = ( + _collect_artifacts(persistent_dir) if output_motifs else None) + result = parse_mixmhcpred_output( + prediction_path, + alleles=self.alleles, + artifacts=artifacts, + stdout=stdout, + ) + if not result.version: + result.version = self.version + return _validate_prediction_rows(result, peptide_list) + + def predict(self, peptides, n_flanks=None, c_flanks=None): + """Predict class-I presentation, preserving one result per input.""" + peptide_list, _, _ = _check_flank_inputs( + peptides, n_flanks, c_flanks) + self._check_peptide_inputs(peptide_list) + if not peptide_list: + return [] + if not self.alleles: + return [PeptideResult() for _ in peptide_list] + + if self.exclude_peptides_with_cysteine: + allowed = [p for p in peptide_list if "C" not in p] else: - self.extra_commandline_args = [] + allowed = peptide_list + predicted = iter(self.predict_detailed(allowed).predictions) + results = [] + for peptide in peptide_list: + if self.exclude_peptides_with_cysteine and "C" in peptide: + results.append(PeptideResult()) + else: + results.append(next(predicted)) + return results def predict_peptides(self, peptides): - """ + """Legacy ``BindingPredictionCollection`` prediction API.""" + return BindingPredictionCollection([ + BindingPrediction.from_pred(pred) + for peptide_result in self.predict(peptides) + for pred in peptide_result.preds + ]) + + def predict_allele_sequences( + self, + allele_sequences, + peptides=None, + output_dir=None, + output_motifs=False): + """Align MHC-I sequences and optionally predict peptides/motifs. Parameters ---------- - peptides : list of str + allele_sequences : mapping or str/path-like + Mapping of allele identifiers to unaligned MHC-I sequences, or a + FASTA filename. + peptides : sequence of str, optional + Peptides to score against sequence-derived motifs. Omit for + alignment and/or motif generation only. + output_dir : str/path-like, optional + New directory in which to retain all upstream artifacts. Required + for ``output_motifs=True``. If omitted for alignment/prediction, + files are parsed from a temporary directory and only structured + Python output is retained. + output_motifs : bool + Generate PWM tables/images, peptide-length distributions, and the + HTML overview. + """ + if hasattr(allele_sequences, "items"): + records = tuple( + (str(name), str(sequence)) + for name, sequence in allele_sequences.items()) + else: + records = _read_fasta(allele_sequences) + if not records: + raise ValueError("At least one MHC-I allele sequence is required") + names = [name for name, _ in records] + if any(not name for name in names) or len(set(names)) != len(names): + raise ValueError("MHC-I FASTA identifiers must be non-empty and unique") + for name, sequence in records: + if not sequence or not sequence.isalpha() or not sequence.isupper(): + raise ValueError( + f"Invalid MHC-I sequence for {name!r}: expected uppercase " + "amino acids") + peptide_list = None + if peptides is not None: + peptide_list, _, _ = _check_flank_inputs(peptides) + self._check_peptide_inputs(peptide_list) + if self.exclude_peptides_with_cysteine: + peptide_list = [p for p in peptide_list if "C" not in p] + if not peptide_list: + peptide_list = None + if output_motifs and output_dir is None: + raise ValueError("output_dir is required when output_motifs=True") + self._require_v3("Allele-sequence alignment and prediction") - Returns - ------- - list of BindingPrediction - """ - self._check_peptide_inputs(peptides) - results = [] - for allele in self.alleles: - - temp_dir = mkdtemp(prefix="mhctools", suffix="mixmhcpred") - input_file_path = join(temp_dir, "mixmhcpred_inputs.txt") - output_file_path = join(temp_dir, "mixmhcpred_outputs.txt") - - with open(input_file_path, "w") as f: - for i, p in enumerate(peptides): - f.write(p) - if i < len(peptides) - 1: - f.write("\n") - with CleanupFiles( - filenames=[input_file_path, output_file_path], - directories=[temp_dir]): - with NamedTemporaryFile(prefix="MixMHCpred_stdout", mode="w", delete=False) as stdout_file: - stdout_file_name = stdout_file.name - run_command([ - self.program_name, - "-i", input_file_path, - "-o", output_file_path, - "-a", normalize_allele_name(allele)] + self.extra_commandline_args, - suppress_stderr=False, - redirect_stdout_file=stdout_file) - if exists(output_file_path): - results.extend(parse_mixmhcpred_results(output_file_path)) - else: - with open(stdout_file_name, "r") as f: - stdout = f.read().strip() - raise ValueError( - "MixMHCpred failed on allele '%s' with stdout '%s'" % (allele, stdout)) - remove(stdout_file_name) - return BindingPredictionCollection(results) + persistent_dir = ( + self._validate_output_dir(output_dir) + if output_dir is not None else None) + with TemporaryDirectory( + prefix="mhctools-mixmhcpred-sequences-") as temp_dir: + sequence_path = Path(temp_dir) / "alleles.fasta" + peptide_path = Path(temp_dir) / "peptides.txt" + stdout_path = Path(temp_dir) / "stdout.txt" + _write_fasta(sequence_path, records) + if peptide_list is not None: + peptide_path.write_text("\n".join(peptide_list) + "\n") + run_output_dir = Path(persistent_dir or (Path(temp_dir) / "output")) + args = [ + self.executable, + "-s", str(sequence_path), + "-o", str(run_output_dir), + ] + if peptide_list is not None: + args.extend(("-i", str(peptide_path), "-p", "1")) + if output_motifs: + args.extend(("-m", "1")) + stdout = self._run_command(args, stdout_path) -def parse_mixmhcpred_results(filename): - """ - Parses output files of MixMHCpred that are expected to look like: + alignment_path = run_output_dir / "final_alignment.fasta" + if not alignment_path.is_file() or alignment_path.stat().st_size == 0: + raise RuntimeError( + "MixMHCpred exited without a non-empty sequence alignment. " + "Ensure MAFFT is installed and runnable. stdout: " + f"{stdout.strip()}") + aligned = _read_fasta(alignment_path) + artifacts = ( + _collect_artifacts(persistent_dir) if persistent_dir else None) - Peptide Score_bestAllele BestAllele %Rank_bestAllele Score_A0201 %Rank_A0201 - MLDDFSAGA 0.182093 A0201 0.3 0.182093 0.3 - SPEGEETII -0.655341 A0201 51.0 -0.655341 51.0 - ILDRIITNA 0.203906 A0201 0.3 0.203906 0.3 + if peptide_list is None: + return MixMHCpredResult( + table=pd.DataFrame(columns=( + "Peptide", "Score_bestAllele", "BestAllele", + "%Rank_bestAllele")), + allele_info=_sequence_allele_info(names, artifacts), + version=self.version, + stdout=stdout, + aligned_sequences=aligned, + artifacts=artifacts, + ) - Parameters - ---------- - filename : str - - Returns list of BindingPrediction - """ - df = pd.read_csv(filename, comment="#", sep="\t") - binding_predictions = [] - for peptide, allele, score, pr in zip( - df["Peptide"], - df["BestAllele"], - df["Score_bestAllele"], - df["%Rank_bestAllele"]): - binding_predictions.append(BindingPrediction( - peptide=peptide, - allele=normalize_allele_name(allele), - score=score, - percentile_rank=pr, - prediction_method_name="mixmhcpred")) - return binding_predictions + prediction_path = run_output_dir / "Binding_predictions.txt" + if not prediction_path.is_file(): + raise RuntimeError( + "MixMHCpred exited without sequence-driven predictions at " + f"{prediction_path}. stdout: {stdout.strip()}") + result = parse_mixmhcpred_output( + prediction_path, + alleles=names, + artifacts=artifacts, + stdout=stdout, + aligned_sequences=aligned, + sequence_mode=True, + ) + if not result.version: + result.version = self.version + return _validate_prediction_rows(result, peptide_list) + + def _default_pred_kind(self): + return Kind.pMHC_presentation + + def kind_support(self): + return { + Kind.pMHC_presentation: { + "mhc_dependence": "single_allele", + "mhc_class": "I", + }, + } diff --git a/pyproject.toml b/pyproject.toml index 7dfeb16..f898fad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,12 @@ dependencies = [ pepsickle = [ "pepsickle", ] +mixmhcpred = [ + # Runtime packages used by the separately licensed upstream v3 checkout. + "scipy", + "logomaker", + "matplotlib", +] nettcr = [ # A TFLite runtime for NetTCR's bundled models. ai-edge-litert is the # lightweight option; tensorflow also works. The NetTCR-2.2 repository @@ -44,7 +50,8 @@ nettcr = [ ] dev = [ "build", - "ruff", + # Ruff 0.16 changed its default rule set; keep lint reproducible until #268. + "ruff==0.15.21", "pytest", "pytest-cov", "twine", diff --git a/tests/test_mixmhcpred.py b/tests/test_mixmhcpred.py index 277b2df..d1e28cc 100644 --- a/tests/test_mixmhcpred.py +++ b/tests/test_mixmhcpred.py @@ -10,22 +10,322 @@ # See the License for the specific language governing permissions and # limitations under the License. -from tempfile import NamedTemporaryFile -from mhctools.mixmhcpred import parse_mixmhcpred_results -from .common import eq_ +import os +import shutil +import sys +from pathlib import Path -example_output = """Peptide\tScore_bestAllele\tBestAllele\t%Rank_bestAllele\tScore_A0201\t%Rank_A0201 +import pytest + +from mhctools import MixMHCpred +from mhctools.mixmhcpred import ( + mixmhcpred_version, + parse_mixmhcpred_output, + parse_mixmhcpred_results, + resolve_mixmhcpred_path, +) +from mhctools.pred import COLUMNS, Kind + +V3_OUTPUT = """#################### +# Output from MixMHCpred (v3.0) +# Alleles: A0201, A0102 - predicted motif for A0102 +# Closest Allele (distance score): A0201 (0.0) -- A0101 (0.037) +# Input file: peptides.txt +#################### +Peptide\tScore_bestAllele\tBestAllele\t%Rank_bestAllele\tScore_A0201\t%Rank_A0201\tScore_A0102\t%Rank_A0102 +SIINFEKL\t-1.119737\tA0201\t2.698623\t-1.119737\t2.698623\t-3.902923\t41.899867 +MLDDFSAGA\t0.182093\tA0201\t0.3\t0.182093\t0.3\t-2.0\t20.0 +""" + +LEGACY_OUTPUT = """Peptide\tScore_bestAllele\tBestAllele\t%Rank_bestAllele\tScore_A0201\t%Rank_A0201 MLDDFSAGA\t0.182093\tA0201\t0.3\t0.182093\t0.3 SPEGEETII\t-0.655341\tA0201\t51.0\t-0.655341\t51.0 -ILDRIITNA\t0.203906\tA0201\t0.3\t0.203906\t0.3""" - -def test_parse_mixmhcpred_results(): - with NamedTemporaryFile(mode="r+") as f: - f.write(example_output) - f.flush() - binding_results = parse_mixmhcpred_results(f.name) - - eq_(len(binding_results), 3) - eq_(binding_results[0].peptide, "MLDDFSAGA") - eq_(binding_results[1].peptide, "SPEGEETII") - eq_(binding_results[2].peptide, "ILDRIITNA") +ILDRIITNA\t0.203906\tA0201\t0.3\t0.203906\t0.3 +""" + +SINGLE_OUTPUT = """# Output from MixMHCpred (v3.0) +Peptide\tScore_bestAllele\tBestAllele\t%Rank_bestAllele\tScore_A0201\t%Rank_A0201 +SIINFEKL\t-1.119737\tA0201\t2.698623\t-1.119737\t2.698623 +""" + +SEQUENCE_OUTPUT = """#################### +# Output from MixMHCpred (v3.0) +# Predictions for Alleles : novel-one, novel-two +#################### +Peptide\tScore_bestAllele\tBestAllele\t%Rank_bestAllele\tScore_novel-one\t%Rank_novel-one\tScore_novel-two\t%Rank_novel-two +SIINFEKL\t-0.25\tnovel-two\t0.23\t-3.5\t32.6\t-0.25\t0.23 +""" + + +def _write(path, text): + Path(path).write_text(text) + return path + + +def test_parse_v3_preserves_all_outputs(tmp_path): + path = _write(tmp_path / "v3.txt", V3_OUTPUT) + result = parse_mixmhcpred_output( + path, alleles=["HLA-A*02:01", "HLA-A*01:02"]) + + assert result.version == "3.0" + assert list(result.table.columns) == [ + "Peptide", "Score_bestAllele", "BestAllele", "%Rank_bestAllele", + "Score_A0201", "%Rank_A0201", "Score_A0102", "%Rank_A0102", + ] + assert result.table.loc[0, "Score_bestAllele"] == pytest.approx(-1.119737) + assert result.table.loc[0, "BestAllele"] == "A0201" + assert result.table.loc[0, "%Rank_bestAllele"] == pytest.approx(2.698623) + + a0201, a0102 = result.allele_info + assert a0201.allele == "HLA-A*02:01" + assert a0201.native_allele == "A0201" + assert a0201.closest_training_allele == "HLA-A*02:01" + assert a0201.distance == 0.0 + assert not a0201.pan_allele + assert a0102.closest_training_allele == "HLA-A*01:01" + assert a0102.distance == pytest.approx(0.037) + assert a0102.pan_allele + + predictions = result.predictions + assert len(predictions) == 2 + assert len(predictions[0].preds) == 2 + assert predictions[0].presentation.allele == "HLA-A*02:01" + assert all(p.kind == Kind.pMHC_presentation for p in predictions[0].preds) + assert all(p.predictor_version == "3.0" for p in predictions[0].preds) + assert predictions[0].filter(allele="HLA-A*01:02")[0].score == pytest.approx( + -3.902923) + + +def test_legacy_parser_remains_compatible(tmp_path): + path = _write(tmp_path / "legacy.txt", LEGACY_OUTPUT) + predictions = parse_mixmhcpred_results(path) + assert len(predictions) == 3 + assert [p.peptide for p in predictions] == [ + "MLDDFSAGA", "SPEGEETII", "ILDRIITNA"] + assert all(p.allele == "HLA-A*02:01" for p in predictions) + + v3_path = _write(tmp_path / "v3.txt", V3_OUTPUT) + v3_predictions = parse_mixmhcpred_results(v3_path) + assert len(v3_predictions) == 2 + assert all(p.allele == "HLA-A*02:01" for p in v3_predictions) + + +def test_parser_rejects_missing_per_allele_rank(tmp_path): + bad = V3_OUTPUT.replace("\t%Rank_A0102", "") + path = _write(tmp_path / "bad.txt", bad) + with pytest.raises(ValueError, match="%Rank_A0102"): + parse_mixmhcpred_output(path) + + +def test_resolve_path_from_environment_or_checkout(tmp_path, monkeypatch): + executable = _write(tmp_path / "MixMHCpred", "#!/bin/sh\n") + monkeypatch.setenv("MIXMHCPRED_PATH", str(tmp_path)) + assert resolve_mixmhcpred_path() == str(Path(executable).resolve()) + assert resolve_mixmhcpred_path(tmp_path) == str(Path(executable).resolve()) + + +def test_resolve_path_error_has_install_guidance(monkeypatch): + monkeypatch.delenv("MIXMHCPRED_PATH", raising=False) + monkeypatch.setattr(shutil, "which", lambda _: None) + with pytest.raises(FileNotFoundError, match="official v3.0 release"): + resolve_mixmhcpred_path() + + +def test_version_detection(tmp_path): + executable = _write( + tmp_path / "MixMHCpred", "#!/bin/sh\necho MixMHCpred3.0\n") + Path(executable).chmod(0o755) + assert mixmhcpred_version(executable) == "3.0" + + +def test_predict_batches_alleles_and_filters_cysteine(tmp_path, monkeypatch): + predictor = MixMHCpred( + alleles=["HLA-A*02:01", "HLA-A*01:02"], + program_name=sys.executable, + exclude_peptides_with_cysteine=True, + ) + calls = [] + + def fake_run(args, stdout_path): + calls.append(args) + input_path = Path(args[args.index("-i") + 1]) + assert input_path.read_text() == "SIINFEKL\n" + output_path = Path(args[args.index("-o") + 1]) + _write(output_path, V3_OUTPUT.splitlines()[0] + "\n" + "\n".join( + line for line in V3_OUTPUT.splitlines()[1:] + if not line.startswith("MLDDFSAGA")) + "\n") + _write(stdout_path, "DONE\n") + return "DONE\n" + + monkeypatch.setattr(predictor, "_run_command", fake_run) + results = predictor.predict(["SIINFEKL", "CILGFVFTL"]) + + assert len(calls) == 1 + assert calls[0][calls[0].index("-a") + 1] == ( + "HLA-A*02:01,HLA-A*01:02") + assert "-c" not in calls[0] + assert len(results) == 2 + assert len(results[0].preds) == 2 + assert results[1].preds == () + assert predictor.supported_kinds == (Kind.pMHC_presentation,) + assert list(results[0].to_dataframe().columns) == list(COLUMNS) + assert list(predictor.predict_dataframe([]).columns) == list(COLUMNS) + + +def test_predict_rejects_silent_or_stale_output(monkeypatch): + predictor = MixMHCpred( + alleles=["HLA-A*02:01"], program_name=sys.executable) + + def fake_run(args, stdout_path): + _write(Path(args[args.index("-o") + 1]), LEGACY_OUTPUT) + _write(stdout_path, "DONE\n") + return "DONE\n" + + monkeypatch.setattr(predictor, "_run_command", fake_run) + with pytest.raises(RuntimeError, match="do not match its inputs"): + predictor.predict(["SIINFEKL"]) + + +def test_motif_artifacts_are_returned_and_existing_path_is_safe( + tmp_path, monkeypatch): + predictor = MixMHCpred( + alleles=["HLA-A*02:01"], program_name=sys.executable) + predictor._version = "3.0" + output_dir = tmp_path / "motifs" + + def fake_run(args, stdout_path): + destination = Path(args[args.index("-o") + 1]) + destination.mkdir() + (destination / "Motifs").mkdir() + _write(destination / "Binding_predictions.txt", SINGLE_OUTPUT) + _write(destination / "Data_overview.html", "") + _write(destination / "Motifs" / "A0201_PWM9.png", "fake") + _write(stdout_path, "DONE\n") + return "DONE\n" + + monkeypatch.setattr(predictor, "_run_command", fake_run) + result = predictor.predict_detailed( + ["SIINFEKL"], output_dir=output_dir, output_motifs=True) + assert result.artifacts.output_dir == str(output_dir.resolve()) + assert result.artifacts.binding_predictions.endswith( + "Binding_predictions.txt") + assert result.artifacts.overview_html.endswith("Data_overview.html") + assert len(result.artifacts.files) == 3 + with pytest.raises(FileExistsError, match="Refusing"): + predictor.predict_detailed( + ["SIINFEKL"], output_dir=output_dir, output_motifs=True) + + +def test_sequence_prediction_captures_alignment_quality_and_artifacts( + tmp_path, monkeypatch): + predictor = MixMHCpred(program_name=sys.executable) + predictor._version = "3.0" + output_dir = tmp_path / "sequence-output" + sequences = { + "novel-one": "ACDEFGHIKLMNPQRSTVWY", + "novel-two": "YWVTSRQPNMLKIHGFEDCA", + } + + def fake_run(args, stdout_path): + destination = Path(args[args.index("-o") + 1]) + destination.mkdir() + (destination / "Motifs").mkdir() + sequence_path = Path(args[args.index("-s") + 1]) + shutil.copyfile(sequence_path, destination / "final_alignment.fasta") + _write(destination / "Binding_predictions.txt", SEQUENCE_OUTPUT) + _write( + destination / "Data_overview.html", + "

novel-one

" + "

closest Allele = A0201, distance = 0.1

" + "

closest Allele from database = A0202, distance = 0.01

" + "

novel-two

" + "

closest Allele = B0702, distance = 0.2

" + "

closest Allele from database = B0703, distance = 0.02

", + ) + _write(destination / "Motifs" / "PLD_sequences.txt", "fake") + _write(stdout_path, "DONE\n") + return "DONE\n" + + monkeypatch.setattr(predictor, "_run_command", fake_run) + result = predictor.predict_allele_sequences( + sequences, + peptides=["SIINFEKL"], + output_dir=output_dir, + output_motifs=True, + ) + + assert result.aligned_sequences == tuple(sequences.items()) + assert result.artifacts.alignment.endswith("final_alignment.fasta") + assert result.artifacts.overview_html.endswith("Data_overview.html") + assert result.predictions[0].filter(allele="novel-two")[0].score == -0.25 + first, second = result.allele_info + assert first.pan_allele and second.pan_allele + assert first.closest_training_allele == "HLA-A*02:01" + assert first.closest_database_allele == "HLA-A*02:02" + assert first.distance == 0.1 + assert first.database_distance == 0.01 + + +def test_alignment_only_retains_sequence_provenance(monkeypatch): + predictor = MixMHCpred(program_name=sys.executable) + predictor._version = "3.0" + + def fake_run(args, stdout_path): + destination = Path(args[args.index("-o") + 1]) + destination.mkdir() + shutil.copyfile( + Path(args[args.index("-s") + 1]), + destination / "final_alignment.fasta", + ) + _write(stdout_path, "DONE\n") + return "DONE\n" + + monkeypatch.setattr(predictor, "_run_command", fake_run) + result = predictor.predict_allele_sequences( + {"novel-one": "ACDEFGHIKLMNPQRSTVWY"}) + + assert result.table.empty + assert len(result.allele_info) == 1 + assert result.allele_info[0].allele == "novel-one" + assert result.allele_info[0].native_allele == "novel-one" + assert result.allele_info[0].pan_allele + + +MIXMHCPRED_V3 = os.environ.get("MIXMHCPRED_V3_PATH") +requires_mixmhcpred_v3 = pytest.mark.skipif( + not MIXMHCPRED_V3, + reason="set MIXMHCPRED_V3_PATH to an official MixMHCpred v3.0 checkout", +) + + +@requires_mixmhcpred_v3 +def test_mixmhcpred_v3_end_to_end_known_and_pan_alleles(): + predictor = MixMHCpred( + alleles=["HLA-A*02:01", "HLA-A*01:02"], + program_name=MIXMHCPRED_V3, + ) + result = predictor.predict_detailed(["SIINFEKL"]) + assert result.version == "3.0" + assert result.predictions[0].filter(allele="HLA-A*02:01")[0].score == ( + pytest.approx(-1.119737, abs=1e-6)) + assert result.predictions[0].filter(allele="HLA-A*01:02")[0].score == ( + pytest.approx(-3.902923, abs=1e-6)) + assert [info.pan_allele for info in result.allele_info] == [False, True] + assert result.allele_info[1].closest_training_allele == "HLA-A*01:01" + + +@requires_mixmhcpred_v3 +def test_mixmhcpred_v3_end_to_end_sequence_prediction(): + executable = Path(resolve_mixmhcpred_path(MIXMHCPRED_V3)) + upstream_fixture = executable.parent / "input" / "To_align_sequences.fasta" + if not upstream_fixture.is_file() or shutil.which("mafft") is None: + pytest.skip("official sequence fixture and MAFFT are required") + + result = MixMHCpred( + program_name=MIXMHCPRED_V3).predict_allele_sequences( + upstream_fixture, peptides=["SIINFEKL"]) + assert [name for name, _ in result.aligned_sequences] == [ + "HLA-A*01:07", "Mafa-B*008:02", "Mamu-B*030:02:01:08"] + assert len(result.predictions) == 1 + assert len(result.predictions[0].preds) == 3 + assert result.table.loc[0, "BestAllele"] == "Mamu-B*030:02:01:08"