Add auto-detecting parser for NetMHCpan output

mhctools/__init__.py

@@ -24,7 +24,7 @@
 from .netmhcstabpan import NetMHCstabpan
 from .unsupported_allele import UnsupportedAllele
 
-__version__ = "2.1.0"
+__version__ = "2.2.0"
 
 __all__ = [
     "BindingPrediction",

mhctools/parsing.py

@@ -12,12 +12,17 @@
 
 from __future__ import print_function, division, absolute_import
 
+import logging
+import re
+
 import numpy as np
 
 from .allele_normalization import normalize_allele_name
 
 from .binding_prediction import BindingPrediction
 
+logger = logging.getLogger(__name__)
+
 
 NETMHC_TOKENS = {
     "pos",
@@ -475,6 +480,155 @@ def parse_netmhcpan41_stdout(
         transforms=transforms)
 
 
+def _find_header_fields(stdout):
+    """
+    Find the header line in NetMHCpan stdout and return its whitespace-split
+    fields. The header is the first non-dash, non-empty line that follows a
+    line of dashes and contains known column names.
+    """
+    prev_was_dash = False
+    for line in stdout.split("\n"):
+        stripped = line.strip()
+        if stripped.startswith("---"):
+            prev_was_dash = True
+            continue
+        if prev_was_dash and stripped:
+            fields = stripped.split()
+            if any(f.lower() in ('pos', 'hla', 'mhc') for f in fields):
+                return fields
+            prev_was_dash = False
+    return None
+
+
+def _detect_netmhcpan_version_label(stdout, header_fields):
+    """
+    Best-effort detection of NetMHCpan version for logging.
+    Checks for an explicit version string first, then infers from header
+    column names.
+    """
+    match = re.search(r'# NetMHCpan version (\S+)', stdout)
+    if match:
+        return "NetMHCpan %s" % match.group(1)
+
+    field_set = set(header_fields)
+    if 'Score_EL' in field_set:
+        return "NetMHCpan 4.1"
+    elif '1-log50k(aff)' in field_set and 'Core' not in field_set:
+        return "NetMHCpan 2.8"
+    elif 'Core' in field_set and 'Aff(nM)' in field_set:
+        return "NetMHCpan 3.x/4.x (binding affinity mode)"
+    elif 'Core' in field_set:
+        return "NetMHCpan 4.x (elution score mode)"
+    else:
+        return "NetMHCpan (unknown version)"
+
+
+def parse_netmhcpan_stdout(
+        stdout,
+        prediction_method_name="netmhcpan",
+        sequence_key_mapping=None,
+        mode=None):
+    """
+    Auto-detecting parser for NetMHCpan output of any supported version
+    (2.8, 3.0, 4.0, 4.1).
+
+    Parses the header line between dash separators to determine column
+    positions, then extracts binding predictions accordingly.
+
+    Parameters
+    ----------
+    stdout : str
+        Raw stdout from any version of NetMHCpan.
+
+    prediction_method_name : str
+
+    sequence_key_mapping : dict or None
+
+    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".
+
+    Returns
+    -------
+    list of BindingPrediction
+    """
+    check_stdout_error(stdout, "NetMHCpan")
+
+    header_fields = _find_header_fields(stdout)
+    if header_fields is None:
+        raise ValueError(
+            "Could not find header line in NetMHCpan output. "
+            "Expected a header row between lines of dashes.")
+
+    field_index = {name: i for i, name in enumerate(header_fields)}
+
+    version_label = _detect_netmhcpan_version_label(stdout, header_fields)
+    logger.info("Detected %s output format", version_label)
+
+    # Position column (case varies across versions)
+    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)
+
+    # 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)
+
+    # 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)
+
+    # 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)
+
+    # 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)')
+        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)')
+
+    return parse_stdout(
+        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)
+
+
 def parse_netmhccons_stdout(
         stdout,
         prediction_method_name="netmhccons",

pyproject.toml

@@ -40,6 +40,7 @@ dev = [
   "pytest",
   "pytest-cov",
   "twine",
+  "wheel",
 ]
 
 [project.scripts]

tests/test_mhc_formats.py

@@ -13,6 +13,7 @@
 from mhctools.parsing import (
   parse_netmhcpan28_stdout,
   parse_netmhcpan3_stdout,
+  parse_netmhcpan_stdout,
   parse_netmhc3_stdout,
   parse_netmhc4_stdout,
 )
@@ -185,3 +186,162 @@ def test_mhcpan3_stdout():
             # expect the epitopes to be sorted in increasing IC50
             assert entry.value == 18234.7, entry
             assert entry.percentile_rank == 10.00, entry
+
+
+def test_auto_detect_netmhcpan28():
+    """Auto-detecting parser should produce identical results to parse_netmhcpan28_stdout."""
+    netmhcpan28_output = """
+    # Affinity Threshold for Strong binding peptides  50.000',
+    # Affinity Threshold for Weak binding peptides 500.000',
+    # Rank Threshold for Strong binding peptides   0.500',
+    # Rank Threshold for Weak binding peptides   2.000',
+    ---------------------------------------------------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
+      1  HLA-A*02:03    QQQQYFPEI   id0         0.278      2461.53   15.00
+     11  HLA-A*02:03    HIIIASSSL   id0         0.515       189.74    4.00 <= WB
+    """
+    results = parse_netmhcpan_stdout(netmhcpan28_output)
+    assert len(results) == 3
+    for entry in results:
+        assert entry.allele == 'HLA-A*02:03'
+    last = [e for e in results if e.peptide == "HIIIASSSL"][0]
+    assert last.value == 189.74
+    assert last.percentile_rank == 4.00
+    assert last.offset == 11
+
+
+def test_auto_detect_netmhcpan3():
+    """Auto-detecting parser should produce identical results to parse_netmhcpan3_stdout."""
+    netmhcpan3_output = """
+    # Rank Threshold for Strong binding peptides   0.500
+    # Rank Threshold for Weak binding peptides   2.000
+    -----------------------------------------------------------------------------------
+      Pos          HLA         Peptide       Core Of Gp Gl Ip Il        Icore        Identity   Score Aff(nM)   %Rank  BindLevel
+    -----------------------------------------------------------------------------------
+        1  HLA-B*18:01        QQQQQYFP  QQQQQYFP-  0  0  0  8  1     QQQQQYFP             id0 0.06456 24866.4   17.00
+        9  HLA-B*18:01        EITHIIAS  -EITHIIAS  0  0  0  0  1     EITHIIAS             id0 0.09323 18234.7   10.00
+    """
+    results = parse_netmhcpan_stdout(netmhcpan3_output)
+    assert len(results) == 2
+    for entry in results:
+        assert entry.allele == 'HLA-B*18:01'
+    first = results[0]
+    assert first.peptide == "QQQQQYFP"
+    assert first.offset == 0  # 1-based converted to 0-based
+    assert first.percentile_rank == 17.00
+
+    second = results[1]
+    assert second.value == 18234.7
+    assert second.percentile_rank == 10.00
+    assert second.offset == 8  # 9 -> 8 (1-based to 0-based)
+
+
+def test_auto_detect_netmhcpan4_ba():
+    """Auto-detecting parser for NetMHCpan 4.0 binding affinity mode."""
+    netmhcpan4_ba_output = """
+# NetMHCpan version 4.0
+
+# Input is in PEPTIDE format
+
+# Make binding affinity predictions
+
+HLA-A02:01 : Distance to training data  0.000 (using nearest neighbor HLA-A02:01)
+
+# Rank Threshold for Strong binding peptides   0.500
+# Rank Threshold for Weak binding peptides   2.000
+-----------------------------------------------------------------------------------
+  Pos          HLA         Peptide       Core Of Gp Gl Ip Il        Icore        Identity     Score Aff(nM)   %Rank  BindLevel
+-----------------------------------------------------------------------------------
+    1  HLA-A*02:01        SIINFEKL  SIINF-EKL  0  0  0  5  1     SIINFEKL         PEPLIST 0.1141340 14543.1 18.9860
+-----------------------------------------------------------------------------------
+"""
+    results = parse_netmhcpan_stdout(netmhcpan4_ba_output)
+    assert len(results) == 1
+    entry = results[0]
+    assert entry.allele == 'HLA-A*02:01'
+    assert entry.peptide == "SIINFEKL"
+    assert entry.offset == 0  # 1-based to 0-based
+    assert abs(entry.score - 0.1141340) < 1e-6
+    assert abs(entry.value - 14543.1) < 0.1
+    assert abs(entry.percentile_rank - 18.9860) < 0.001
+
+
+def test_auto_detect_netmhcpan4_el():
+    """Auto-detecting parser for NetMHCpan 4.0 elution score mode (no Aff column)."""
+    netmhcpan4_el_output = """
+# NetMHCpan version 4.0
+
+# Input is in PEPTIDE format
+
+HLA-A02:01 : Distance to training data  0.000 (using nearest neighbor HLA-A02:01)
+
+# Rank Threshold for Strong binding peptides   0.500
+# Rank Threshold for Weak binding peptides   2.000
+-----------------------------------------------------------------------------------
+  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
+-----------------------------------------------------------------------------------
+"""
+    results = parse_netmhcpan_stdout(netmhcpan4_el_output)
+    assert len(results) == 1
+    entry = results[0]
+    assert entry.allele == 'HLA-A*02:01'
+    assert entry.peptide == "SIINFEKL"
+    assert entry.offset == 0
+    assert abs(entry.score - 0.3456780) < 1e-6
+    assert entry.value is None  # no Aff(nM) column
+    assert abs(entry.percentile_rank - 5.1230) < 0.001
+
+
+def test_auto_detect_netmhcpan41_ba():
+    """Auto-detecting parser for NetMHCpan 4.1 in binding affinity mode."""
+    netmhcpan41_output = """
+# NetMHCpan version 4.1b
+
+# Input is in PEPTIDE format
+
+# Make both EL and BA predictions
+
+HLA-A02:01 : Distance to training data  0.000 (using nearest neighbor HLA-A02:01)
+
+# Rank Threshold for Strong binding peptides   0.500
+# Rank Threshold for Weak binding peptides   2.000
+---------------------------------------------------------------------------------------------------------------------------
+ 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
+---------------------------------------------------------------------------------------------------------------------------
+"""
+    results = parse_netmhcpan_stdout(netmhcpan41_output, mode="binding_affinity")
+    assert len(results) == 1
+    entry = results[0]
+    assert entry.allele == 'HLA-A*02:01'
+    assert entry.peptide == "SIINFEKL"
+    assert entry.offset == 0
+    assert abs(entry.score - 0.110414) < 1e-5  # Score_BA
+    assert abs(entry.percentile_rank - 20.171) < 0.01  # %Rank_BA
+    assert abs(entry.value - 15140.42) < 0.1  # Aff(nM)
+
+
+def test_auto_detect_netmhcpan41_el():
+    """Auto-detecting parser for NetMHCpan 4.1 in elution score mode."""
+    netmhcpan41_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
+---------------------------------------------------------------------------------------------------------------------------
+"""
+    results = parse_netmhcpan_stdout(netmhcpan41_output, mode="elution_score")
+    assert len(results) == 1
+    entry = results[0]
+    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