Fix NetMHCpan version sniffing and NetMHCIIpan parsing errors

mhctools/__init__.py

@@ -40,7 +40,7 @@
 from .netmhcstabpan import NetMHCstabpan
 from .unsupported_allele import UnsupportedAllele
 
-__version__ = "3.0.1"
+__version__ = "3.0.2"
 
 __all__ = [
     "Pred",

mhctools/netmhc_pan.py

@@ -11,9 +11,9 @@
 # limitations under the License.
 
 import logging
-import re
-from subprocess import check_output
 import os
+import re
+import subprocess
 
 from .netmhc_pan28 import NetMHCpan28
 from .netmhc_pan3 import NetMHCpan3
@@ -60,11 +60,14 @@ def NetMHCpan(
     falls back to the latest known class (which uses the header-driven
     auto-detecting parser).
     """
-    with open(os.devnull, 'w') as devnull:
-        output = check_output([
-            program_name, "--version", "_MHCTOOLS_VERSION_SNIFFING"],
-            stderr=devnull)
-    output_str = output.decode("ascii", "ignore")
+    # Pass /dev/null as the input file so netMHCpan doesn't try to
+    # open a nonexistent file (it hangs without any argument).
+    # Use run() instead of check_output() because netMHCpan may exit
+    # non-zero even though the version string is still in stdout.
+    result = subprocess.run(
+        [program_name, "--version", os.devnull],
+        capture_output=True)
+    output_str = result.stdout.decode("ascii", "ignore")
 
     match = re.search(r'# NetMHCpan version (\S+)', output_str)
     version_str = match.group(1) if match else ""

mhctools/parsing.py

@@ -185,17 +185,17 @@ def parse_stdout(
         if score_index is None:
             score = None
         else:
-            score = float(fields[score_index])
+            score = _try_float(fields[score_index])
 
         if rank_index is None:
             rank = None
         else:
-            rank = float(fields[rank_index])
+            rank = _try_float(fields[rank_index])
 
         if ic50_index is None:
             ic50 = None
         else:
-            ic50 = float(fields[ic50_index])
+            ic50 = _try_float(fields[ic50_index])
 
         key = str(fields[key_index])
         if sequence_key_mapping:
@@ -522,6 +522,14 @@ def _detect_netmhcpan_version_label(stdout, header_fields):
         return "NetMHCpan (unknown version)"
 
 
+def _try_float(value):
+    """Convert *value* to float, returning NaN for non-numeric strings like 'NA'."""
+    try:
+        return float(value)
+    except (ValueError, TypeError):
+        return float('nan')
+
+
 def _safe_float(fields, index):
     """Extract a float from fields at index, or None if index is None."""
     if index is None:
@@ -619,8 +627,9 @@ def parse_netmhcpan_to_preds(
 
     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
+        # Optional trailing bind-level tokens (<=, WB, SB) may appear
+        # after the last numeric column; they are harmless because the
+        # header-driven indices never reach them.
 
         offset = int(fields[offset_index])
         if offset_is_one_based:
@@ -861,6 +870,15 @@ def parse_netmhciipan4_stdout(
         0: lambda x: int(x) - 1,
     }
 
+    # Strip optional BindLevel indicators (index 13 for v4.0 layout).
+    ignored = {
+        "<=SB": 13,
+        "<=WB": 13,
+        "<=": 13,
+        "SB": 14,
+        "WB": 14,
+    }
+
     # we're running NetMHCIIpan 4 with -BA every time so both EL and BA are available, but only
     # return one of them depending on the input mode
     return parse_stdout(
@@ -874,6 +892,7 @@ def parse_netmhciipan4_stdout(
         ic50_index=11 if mode == "binding_affinity" else None,
         rank_index=8 if mode == "elution_score" else 12,
         score_index=7 if mode == "elution_score" else 10,
+        ignored_value_indices=ignored,
         transforms=transforms)
 
 def parse_netmhcstabpan(
@@ -961,20 +980,35 @@ def parse_netmhciipan43_stdout(
         ba_score_index = 11
         ba_rank_index = 12
         ba_ic50_index = 13
+        # BindLevel (e.g. "<=SB", "<=WB") is optional at index 14
+        bind_level_index = 14
     else:
         key_index = 6
         el_score_index = 7
         el_rank_index = 8
         ba_score_index = 10
         ba_rank_index = 12
         ba_ic50_index = 11
+        # BindLevel is optional at index 13
+        bind_level_index = 13
 
     # the offset specified in "pos" (at index 0) is 1-based instead of 0-based. we adjust it to be
     # 0-based, as in all the other netmhc predictors supported by this library.
     transforms = {
         0: lambda x: int(x) - 1,
     }
 
+    # Strip optional BindLevel indicators that only appear on some rows.
+    # These are at the end of the line and may be one token ("<=SB") or
+    # two ("<=", "SB") depending on the netMHCIIpan version.
+    ignored = {
+        "<=SB": bind_level_index,
+        "<=WB": bind_level_index,
+        "<=": bind_level_index,
+        "SB": bind_level_index + 1,
+        "WB": bind_level_index + 1,
+    }
+
     # we're running NetMHCIIpan 4.3 with -BA every time so both EL and BA are available, but only
     # return one of them depending on the input mode
     return parse_stdout(
@@ -988,4 +1022,5 @@ def parse_netmhciipan43_stdout(
         ic50_index=ba_ic50_index if mode == "binding_affinity" else None,
         rank_index=el_rank_index if mode == "elution_score" else ba_rank_index,
         score_index=el_score_index if mode == "elution_score" else ba_score_index,
+        ignored_value_indices=ignored,
         transforms=transforms)

tests/test_mhc_formats.py

@@ -17,6 +17,8 @@
   parse_netmhcpan_to_preds,
   parse_netmhc3_stdout,
   parse_netmhc4_stdout,
+  parse_netmhciipan43_stdout,
+  parse_netmhciipan4_stdout,
 )
 from mhctools.pred import Kind
 
@@ -452,3 +454,76 @@ def test_to_preds_self_contained():
     assert p.offset == 5
     assert p.predictor_name == "netMHCpan"
     assert p.predictor_version == "2.8"
+
+
+# ---------- NetMHCIIpan 4.3 ----------
+
+def test_netmhciipan43_el_with_bind_level():
+    """Rows with <=SB BindLevel should parse without error (GH-169)."""
+    output = """
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     Pos               MHC              Peptide   Of        Core  Core_Rel Inverted        Identity      Score_EL %Rank_EL  Exp_Bind      Score_BA %Rank_BA  Affinity(nM)  BindLevel
+    --------------------------------------------------------------------------------------------------------------------------------------------
+       1         DRB1_0101      GAATVAAGAATTAAG    4   VAAGAATTA     0.920        0        Sequence      0.164981     6.81     0.000      0.514261    17.98        191.63
+       2         DRB1_0101      KSVPLEMLLINLTTI    4   LEMLLINLT     0.980        0        Sequence      0.807346     0.50     0.560      0.662093     4.95         38.71 <=SB
+    """
+    results = parse_netmhciipan43_stdout(output, mode="elution_score")
+    assert len(results) == 2
+    assert results[0].peptide == "GAATVAAGAATTAAG"
+    assert abs(results[0].percentile_rank - 6.81) < 0.01
+    assert results[1].peptide == "KSVPLEMLLINLTTI"
+    assert abs(results[1].percentile_rank - 0.50) < 0.01
+
+
+def test_netmhciipan43_ba_with_bind_level():
+    """BA mode should parse rows with <=WB BindLevel (GH-169)."""
+    output = """
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     Pos               MHC              Peptide   Of        Core  Core_Rel Inverted        Identity      Score_EL %Rank_EL  Exp_Bind      Score_BA %Rank_BA  Affinity(nM)  BindLevel
+    --------------------------------------------------------------------------------------------------------------------------------------------
+       1         DRB1_0101      GAATVAAGAATTAAG    4   VAAGAATTA     0.920        0        Sequence      0.164981     6.81     0.000      0.514261    17.98        191.63 <=WB
+    """
+    results = parse_netmhciipan43_stdout(output, mode="binding_affinity")
+    assert len(results) == 1
+    assert abs(results[0].affinity - 191.63) < 0.01
+    assert abs(results[0].percentile_rank - 17.98) < 0.01
+
+
+def test_netmhciipan43_exp_bind_na():
+    """Exp_Bind = NA should not crash float conversion (GH-169)."""
+    output = """
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     Pos               MHC              Peptide   Of        Core  Core_Rel Inverted        Identity      Score_EL %Rank_EL  Exp_Bind      Score_BA %Rank_BA  Affinity(nM)  BindLevel
+    --------------------------------------------------------------------------------------------------------------------------------------------
+       1         DRB1_0101      PAPAPSWPLSSSVPS    4   PSWPLSSSV     0.327        0        Sequence      0.000857    79.79       NA      0.327674    54.35       1442.91
+    """
+    results = parse_netmhciipan43_stdout(output, mode="elution_score")
+    assert len(results) == 1
+    assert abs(results[0].score - 0.000857) < 0.0001
+
+
+def test_netmhciipan43_two_token_bind_level():
+    """'<= SB' (with space) is two tokens — both should be stripped (GH-169)."""
+    output = """
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     Pos               MHC              Peptide   Of        Core  Core_Rel Inverted        Identity      Score_EL %Rank_EL  Exp_Bind      Score_BA %Rank_BA  Affinity(nM)  BindLevel
+    --------------------------------------------------------------------------------------------------------------------------------------------
+       1         DRB1_0101      KSVPLEMLLINLTTI    4   LEMLLINLT     0.980        0        Sequence      0.807346     0.50     0.560      0.662093     4.95         38.71 <= SB
+    """
+    results = parse_netmhciipan43_stdout(output, mode="binding_affinity")
+    assert len(results) == 1
+    assert abs(results[0].affinity - 38.71) < 0.01
+
+
+def test_netmhciipan4_with_bind_level():
+    """v4.0 layout with BindLevel should also parse (GH-169)."""
+    output = """
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     Pos           MHC              Peptide   Of        Core  Core_Rel        Identity      Score_EL %Rank_EL Exp_Bind      Score_BA  Affinity(nM) %Rank_BA  BindLevel
+    --------------------------------------------------------------------------------------------------------------------------------------------
+       1     DRB1_0101      PAPAPSWPLSSSVPS    4   PSWPLSSSV     0.327            test      0.000857    79.79       NA      0.327674       1442.91    54.35
+       2     DRB1_0101      GAATVAAGAATTAAG    4   VAAGAATTA     0.920            test      0.800000     0.50     0.000      0.514261        191.63    17.98 <=SB
+    """
+    results = parse_netmhciipan4_stdout(output, mode="elution_score")
+    assert len(results) == 2
+    assert abs(results[1].score - 0.800000) < 0.0001