diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 30f2b73..f9ad91e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -55,6 +55,7 @@ jobs: tests/test_processing_predictor.py \ tests/test_pepsickle.py \ tests/test_bigmhc.py \ + tests/test_unparseable_alleles.py \ tests/test_random.py integration-netmhc: diff --git a/mhctools/__init__.py b/mhctools/__init__.py index b30aa3d..91aaf39 100644 --- a/mhctools/__init__.py +++ b/mhctools/__init__.py @@ -78,7 +78,7 @@ def __getattr__(name): raise AttributeError( "module %r has no attribute %r" % (__name__, name)) -__version__ = "3.18.0" +__version__ = "3.19.0" __all__ = [ "Prediction", diff --git a/mhctools/allele_normalization.py b/mhctools/allele_normalization.py index 55d666e..397e764 100644 --- a/mhctools/allele_normalization.py +++ b/mhctools/allele_normalization.py @@ -266,3 +266,33 @@ def normalize_allele_name(raw_allele, omit_dra1=False, infer_class2_pair=True): normalized = "-".join(normalized_parts) _normalized_allele_cache[cache_key] = normalized return normalized + + +def normalize_allele_name_or_raw(name): + """ + Normalize an allele name via mhcgnomes, or fall back to a canonical "raw" + spelling for names it can't parse. + + Some alleles appear in a predictor's own supported list (e.g. netMHCpan's + ``-listMHC``) yet are rejected by mhcgnomes: exotic non-human alleles such + as ``H-2-Qa1`` or ``BoLA-amani.1``, and HLA low/null-expression variants + such as ``HLA-A30:14L``. For these we key on the predictor's own spelling + so that a requested allele and the same allele echoed back in the + predictor's output resolve to a single identity. + + The raw fallback strips surrounding whitespace and the ``*`` gene/allele + separator (which a predictor may add or omit — netMHCpan prints + ``HLA-A*30:14L`` for a requested ``HLA-A30:14L``), preserving case to match + the predictor's ``-listMHC`` / output spelling. + + Note: because parseable names are uppercased before normalization but the + raw fallback preserves case, un-normalizable exotic alleles are effectively + **case-sensitive** — request them exactly as the predictor's ``-listMHC`` + spells them (e.g. ``H-2-Qa1``, not ``h-2-qa1``), or validation against the + raw supported list won't match. + """ + stripped = str(name).strip() + try: + return normalize_allele_name(stripped.upper()) + except AlleleParseError: + return stripped.replace("*", "") diff --git a/mhctools/base_commandline_predictor.py b/mhctools/base_commandline_predictor.py index 1b5f8d0..410cef5 100644 --- a/mhctools/base_commandline_predictor.py +++ b/mhctools/base_commandline_predictor.py @@ -257,7 +257,10 @@ def __init__( valid_alleles=None, default_peptide_lengths=default_peptide_lengths, min_peptide_length=min_peptide_length, - max_peptide_length=max_peptide_length) + max_peptide_length=max_peptide_length, + # Non-human alleles the normalizer can't parse (e.g. H-2-Qa1) are + # kept verbatim and validated against the tool's raw -listMHC list. + keep_unparseable_alleles=True) self._resolve_supported_alleles() diff --git a/mhctools/base_predictor.py b/mhctools/base_predictor.py index 78c4889..a197097 100644 --- a/mhctools/base_predictor.py +++ b/mhctools/base_predictor.py @@ -15,7 +15,10 @@ from collections import defaultdict, namedtuple from typechecks import require_iterable_of -from .allele_normalization import normalize_allele_name +from .allele_normalization import ( + normalize_allele_name, + normalize_allele_name_or_raw, +) from .unsupported_allele import UnsupportedAllele from .binding_prediction_collection import BindingPredictionCollection @@ -115,7 +118,8 @@ def __init__( min_peptide_length=8, max_peptide_length=None, allow_X_in_peptides=False, - allow_lowercase_in_peptides=False): + allow_lowercase_in_peptides=False, + keep_unparseable_alleles=False): """ Parameters ---------- @@ -143,12 +147,19 @@ def __init__( allow_lowercase_in_peptides : bool Allow lowercase letters in peptide sequences + + keep_unparseable_alleles : bool + If True, keep allele names our normalizer can't parse verbatim + instead of raising (used by command-line predictors that validate + against the tool's own supported-allele list, so that non-human + alleles like H-2-Qa1 or BoLA-amani.1 can still be requested). """ # I find myself often constructing a predictor with just one allele # so as a convenience, allow user to not wrap that allele as a list if type(alleles) is str: alleles = alleles.split(',') - self.alleles = self._check_hla_alleles(alleles, valid_alleles) + self.alleles = self._check_hla_alleles( + alleles, valid_alleles, keep_unparseable=keep_unparseable_alleles) if type(default_peptide_lengths) is int: default_peptide_lengths = [default_peptide_lengths] @@ -515,20 +526,31 @@ def predict_subsequences_dataframe( @staticmethod def _check_hla_alleles( alleles, - valid_alleles=None): + valid_alleles=None, + keep_unparseable=False): """ Given a list of HLA alleles and an optional list of valid HLA alleles, return a set of alleles that we will pass into the MHC binding predictor. + + When keep_unparseable is True, allele names the normalizer can't + parse are kept verbatim (original case) rather than raising, so a + caller that validates against the tool's own supported list can still + accept non-human alleles like H-2-Qa1 or BoLA-amani.1. """ require_iterable_of(alleles, str, "HLA alleles") - # Don't run the MHC predictor twice for homozygous alleles, - # only run it for unique alleles - alleles = { - normalize_allele_name(allele.strip().upper()) - for allele in alleles - } + # Keep only unique alleles (don't run the predictor twice for a + # homozygous genotype). When keep_unparseable is set (command-line + # predictors, which validate against the tool's own -listMHC list), + # normalize_allele_name_or_raw retains names mhcgnomes can't parse + # using the same canonical fallback the output parser applies, so a + # requested allele and its echoed form share one identity. Otherwise + # normalize_allele_name raises on an unparseable name, as before. + if keep_unparseable: + alleles = {normalize_allele_name_or_raw(a) for a in alleles} + else: + alleles = {normalize_allele_name(a.strip().upper()) for a in alleles} if valid_alleles: # For some reason netMHCpan drops the '*' in names, so # 'HLA-A*03:01' becomes 'HLA-A03:01' diff --git a/mhctools/parsing.py b/mhctools/parsing.py index a41833d..9502bd4 100644 --- a/mhctools/parsing.py +++ b/mhctools/parsing.py @@ -15,7 +15,7 @@ import numpy as np -from .allele_normalization import normalize_allele_name +from .allele_normalization import normalize_allele_name_or_raw from .binding_prediction import BindingPrediction from .pred import Prediction, Kind @@ -223,7 +223,7 @@ def parse_stdout( source_sequence_name=original_key, offset=offset, peptide=peptide, - allele=normalize_allele_name(allele), + allele=normalize_allele_name_or_raw(allele), score=score, affinity=ic50, percentile_rank=rank, @@ -644,7 +644,7 @@ def parse_netmhcpan_to_preds( offset -= 1 peptide = str(fields[peptide_index]) - allele = normalize_allele_name(str(fields[allele_index])) + allele = normalize_allele_name_or_raw(str(fields[allele_index])) key = str(fields[key_index]) if sequence_key_mapping: diff --git a/tests/test_netmhc_pan.py b/tests/test_netmhc_pan.py index 92f7a21..675c3f4 100644 --- a/tests/test_netmhc_pan.py +++ b/tests/test_netmhc_pan.py @@ -94,6 +94,32 @@ def test_netmhc_pan_multiple_alleles(): "Expected both alleles, got %s" % (observed_alleles,) +def test_netmhc_pan_exotic_unnormalizable_alleles(): + """netMHCpan lists non-human alleles (H-2-Qa1, BoLA-amani.1) and HLA + low/null-expression variants (HLA-A30:14L) that mhcgnomes can't normalize + (issue #220). They must round-trip: request them, run, and get predictions + back keyed by the requested identity. netMHCpan echoes HLA-A30:14L back + with a '*' (HLA-A*30:14L), so both spellings must resolve to one identity. + """ + cases = [ + ("H-2-Qa1", "H-2-Qa1"), + ("BoLA-amani.1", "BoLA-amani.1"), + ("HLA-A30:14L", "HLA-A30:14L"), # requested without '*' + ("HLA-A*30:14L", "HLA-A30:14L"), # requested with '*', same identity + ] + peptides = ["SIINFEKLL", "GILGFVFTL"] + for requested, expected_identity in cases: + predictor = NetMHCpan(alleles=[requested]) + binding_predictions = predictor.predict_peptides(peptides) + assert len(binding_predictions) == len(peptides), \ + "Expected %d predictions for %s, got %s" % ( + len(peptides), requested, binding_predictions) + observed = {bp.allele for bp in binding_predictions} + assert observed == {expected_identity}, \ + "Expected identity %r for requested %r, got %s" % ( + expected_identity, requested, observed) + + def test_netmhc_pan_batched_matches_per_allele(): """Batching alleles into one `-a A,B,C` invocation must produce exactly the same scores as running one process per allele. Uses the real binary's diff --git a/tests/test_unparseable_alleles.py b/tests/test_unparseable_alleles.py new file mode 100644 index 0000000..6a5ac1e --- /dev/null +++ b/tests/test_unparseable_alleles.py @@ -0,0 +1,131 @@ +# 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. + +""" +Unit tests for supporting alleles that netMHCpan lists but mhcgnomes can't +parse (issue #220): exotic non-human alleles (H-2-Qa1, BoLA-amani.1) and HLA +low/null-expression variants (HLA-A30:14L). + +The key invariant is that a requested allele and the same allele echoed back +in a predictor's output resolve to a single identity, so `_check_results` +matches them. Both go through `normalize_allele_name_or_raw`; when parsing +fails it strips whitespace and the '*' gene/allele separator (which netMHCpan +adds or omits) so e.g. requested "HLA-A30:14L" and echoed "HLA-A*30:14L" agree. + +No predictor binary is required. +""" + +import pytest + +from mhctools.allele_normalization import ( + normalize_allele_name_or_raw, + AlleleParseError, +) +from mhctools.base_predictor import BasePredictor +from mhctools.parsing import parse_netmhc4_stdout + + +# --------------------------------------------------------------------------- +# normalize_allele_name_or_raw: parse when possible, canonical raw fallback +# --------------------------------------------------------------------------- + +def test_parseable_allele_is_normalized(): + # Case-insensitive, and the canonical form keeps the '*'. + assert normalize_allele_name_or_raw("hla-a*02:01") == "HLA-A*02:01" + assert normalize_allele_name_or_raw("HLA-A02:01") == "HLA-A*02:01" + + +def test_unparseable_non_human_kept_verbatim(): + # These netMHCpan names don't parse at all; keep the tool's own spelling + # (including original case, which the output echoes). + for name in ["H-2-Qa1", "H-2-Qa2", "BoLA-amani.1", "BoLA-JSP.1", + "BoLA-T2a", "BoLA-gb1.7"]: + assert normalize_allele_name_or_raw(name) == name + + +def test_whitespace_is_stripped(): + assert normalize_allele_name_or_raw(" H-2-Qa1 \n") == "H-2-Qa1" + + +def test_star_stripped_from_unparseable_fallback(): + # netMHCpan prints "HLA-A*30:14L" for a requested "HLA-A30:14L"; both + # spellings must collapse to one identity so results match. + bare = normalize_allele_name_or_raw("HLA-A30:14L") + starred = normalize_allele_name_or_raw("HLA-A*30:14L") + assert bare == starred == "HLA-A30:14L" + + +def test_request_and_echo_share_identity(): + # The invariant _check_results relies on: whatever the user typed and + # whatever netMHCpan echoes back reduce to the same string. + for requested, echoed in [ + ("HLA-A30:14L", "HLA-A*30:14L"), # star added by netMHCpan + ("H-2-Qa1", "H-2-Qa1"), # verbatim + ("BoLA-amani.1", "BoLA-amani.1")]: + assert (normalize_allele_name_or_raw(requested) == + normalize_allele_name_or_raw(echoed)) + + +# --------------------------------------------------------------------------- +# _check_hla_alleles: keep_unparseable gate +# --------------------------------------------------------------------------- + +def test_check_hla_alleles_raises_by_default(): + # Default behavior (used by IEDB / in-process predictors) must still + # reject names it can't normalize. + with pytest.raises(AlleleParseError): + BasePredictor._check_hla_alleles(["H-2-Qa1"]) + + +def test_check_hla_alleles_keeps_unparseable_when_requested(): + result = set(BasePredictor._check_hla_alleles( + ["H-2-Qa1", "BoLA-amani.1"], keep_unparseable=True)) + assert result == {"H-2-Qa1", "BoLA-amani.1"} + + +def test_check_hla_alleles_still_normalizes_parseable_when_keeping(): + # Mixing parseable and unparseable: parseable ones are still canonicalized. + result = set(BasePredictor._check_hla_alleles( + ["hla-a*02:01", "H-2-Qa1"], keep_unparseable=True)) + assert result == {"HLA-A*02:01", "H-2-Qa1"} + + +def test_check_hla_alleles_dedupes_star_variants_of_unparseable(): + # Two spellings of the same un-normalizable allele collapse to one entry. + result = BasePredictor._check_hla_alleles( + ["HLA-A30:14L", "HLA-A*30:14L"], keep_unparseable=True) + assert result == ["HLA-A30:14L"] + + +# --------------------------------------------------------------------------- +# Output parser falls back to the raw name instead of raising (binary-free). +# --------------------------------------------------------------------------- + +def test_parser_keeps_unparseable_output_allele(): + # A netMHCpan/netMHC-style table whose allele column is an exotic name + # mhcgnomes can't parse. Before the fallback this raised in the parser; + # now the raw name is carried through verbatim. + stdout = """ +# NetMHC version 4.0 +----------------------------------------------------------------------------------- + pos HLA peptide Core Offset I_pos I_len D_pos D_len iCore Identity 1-log50k(aff) Affinity(nM) %Rank BindLevel +----------------------------------------------------------------------------------- + 0 H-2-Qa1 SIINFEKLL SIINFEKLL 0 0 0 0 0 SIINFEKLL SEQ_A 0.349 1147.39 4.50 +----------------------------------------------------------------------------------- + +Protein PEPLIST. Allele H-2-Qa1. Number of high binders 0. Number of weak binders 0. Number of peptides 1 +----------------------------------------------------------------------------------- +""" + preds = parse_netmhc4_stdout(stdout) + assert len(preds) == 1 + assert preds[0].allele == "H-2-Qa1" # raw name, not raised on + assert preds[0].peptide == "SIINFEKLL"