diff --git a/mhctools/__init__.py b/mhctools/__init__.py index c04c8a7..a9f0147 100644 --- a/mhctools/__init__.py +++ b/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", diff --git a/mhctools/netmhc_pan.py b/mhctools/netmhc_pan.py index 8c5f834..e82a673 100644 --- a/mhctools/netmhc_pan.py +++ b/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 "" diff --git a/mhctools/parsing.py b/mhctools/parsing.py index d6942ea..e812c6f 100644 --- a/mhctools/parsing.py +++ b/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,6 +980,8 @@ 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 @@ -968,6 +989,8 @@ def parse_netmhciipan43_stdout( 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. @@ -975,6 +998,17 @@ def parse_netmhciipan43_stdout( 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) diff --git a/tests/test_mhc_formats.py b/tests/test_mhc_formats.py index 5b125f3..8384734 100644 --- a/tests/test_mhc_formats.py +++ b/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