Merge pull request #168 from fkgruber/master

added support to NetMHCIIpan 4.3

Author: iskandr

mhctools/__init__.py

@@ -19,7 +19,7 @@
 from .netmhc_pan3 import NetMHCpan3
 from .netmhc_pan4 import NetMHCpan4, NetMHCpan4_BA, NetMHCpan4_EL
 from .netmhc_pan41 import NetMHCpan41, NetMHCpan41_BA, NetMHCpan41_EL
-from .netmhcii_pan import NetMHCIIpan, NetMHCIIpan3, NetMHCIIpan4, NetMHCIIpan4_BA, NetMHCIIpan4_EL
+from .netmhcii_pan import NetMHCIIpan, NetMHCIIpan3, NetMHCIIpan4, NetMHCIIpan4_BA, NetMHCIIpan4_EL, NetMHCIIpan43, NetMHCIIpan43_BA, NetMHCIIpan43_EL
 from .random_predictor import RandomBindingPredictor
 from .netmhcstabpan import NetMHCstabpan
 from .unsupported_allele import UnsupportedAllele
@@ -48,7 +48,7 @@
     "NetMHCpan41",
     "NetMHCpan4_BA",
     "NetMHCpan4_EL",
-    "NetMHCpan41_BA", 
+    "NetMHCpan41_BA",
     "NetMHCpan41_EL",
     "NetMHCIIpan",
     "NetMHCIIpan3",

mhctools/base_commandline_predictor.py

@@ -301,7 +301,8 @@ def predict_peptides(self, peptides):
                         prefix="tmp_%d_%d_%s" % (
                             i,
                             j,
-                            self.program_name))
+                            self.program_name),
+                    suffix="XXXXXX")
                     logger.debug(
                         "Created temporary directory %s for allele %s",
                         temp_dirname,

mhctools/netmhcii_pan.py

@@ -18,7 +18,7 @@
 from mhcnames import parse_classi_or_classii_allele_name
 
 from .base_commandline_predictor import BaseCommandlinePredictor
-from .parsing import parse_netmhciipan_stdout, parse_netmhciipan4_stdout
+from .parsing import parse_netmhciipan_stdout, parse_netmhciipan4_stdout, parse_netmhciipan43_stdout
 
 logger = logging.getLogger(__name__)
 
@@ -216,3 +216,73 @@ def NetMHCIIpan(
         return NetMHCIIpan3(**kwargs)
     else:
         raise ValueError("This software expects NetMHCIIpan version 3.x or 4.0")
+
+
+class NetMHCIIpan43(NetMHCIIpanBase):
+    """
+    Wrapper for NetMHCIIpan 4.3, using a different parser.
+    """
+    def __init__(
+            self,
+            alleles,
+            default_peptide_lengths=[15, 16, 17, 18, 19, 20],
+            program_name="netMHCIIpan",
+            process_limit=-1,
+            mode="elution_score",
+            extra_flags=[]):
+
+        if mode not in ['binding_affinity', 'elution_score']:
+            raise ValueError("Unsupported mode", mode)
+
+        # Always include binding affinity data (-BA flag), though the main score and %rank will
+        # still be EL-based. This gives us access to the BA-based score and %rank columns.
+
+        NetMHCIIpanBase.__init__(
+            self,
+            alleles=alleles,
+            program_name=program_name,
+            process_limit=process_limit,
+            parse_output_fn=partial(parse_netmhciipan43_stdout, mode=mode),
+            default_peptide_lengths=default_peptide_lengths,
+            extra_flags=['-BA'] + extra_flags)
+
+class NetMHCIIpan43_EL(NetMHCIIpan43):
+    """
+    Wrapper for NetMHCIIpan43 when the preferred mode is elution score
+    """
+    def __init__(
+            self,
+            alleles,
+            default_peptide_lengths=[15, 16, 17, 18, 19, 20],
+            program_name="netMHCIIpan",
+            process_limit=-1,
+            extra_flags=[]):
+        NetMHCIIpan43.__init__(
+            self,
+            alleles=alleles,
+            default_peptide_lengths=default_peptide_lengths,
+            program_name=program_name,
+            process_limit=process_limit,
+            mode="elution_score",
+            extra_flags=extra_flags)
+
+
+class NetMHCIIpan43_BA(NetMHCIIpan43):
+    """
+    Wrapper for NetMHCIIpan43 when the preferred mode is binding affinity
+    """
+    def __init__(
+            self,
+            alleles,
+            default_peptide_lengths=[15, 16, 17, 18, 19, 20],
+            program_name="netMHCIIpan",
+            process_limit=-1,
+            extra_flags=[]):
+        NetMHCIIpan43.__init__(
+            self,
+            alleles=alleles,
+            default_peptide_lengths=default_peptide_lengths,
+            program_name=program_name,
+            process_limit=process_limit,
+            mode="binding_affinity",
+            extra_flags=extra_flags)

mhctools/parsing.py

@@ -565,13 +565,13 @@ def parse_netmhciipan4_stdout(
     --------------------------------------------------------------------------------------------------------------------------------------------
      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      APAPSWPLSSSVPSQ    3   PSWPLSSSV     0.333            test      0.001268    71.87       NA      0.346949       1171.30    50.15       
-       3     DRB1_0101      PAPSWPLSSSVPSQK    4   WPLSSSVPS     0.713            test      0.002836    54.45       NA      0.412004        579.40    36.66       
-       4     DRB1_0101      APSWPLSSSVPSQKT    3   WPLSSSVPS     0.773            test      0.003677    49.14       NA      0.448939        388.53    29.75       
-       5     DRB1_0101      PSWPLSSSVPSQKTY    2   WPLSSSVPS     0.407            test      0.001602    66.79       NA      0.470979        306.09    25.98       
-       6     DRB1_0101      SWPLSSSVPSQKTYQ    3   LSSSVPSQK     0.633            test      0.001671    65.82       NA      0.476222        289.21    25.07       
-       7     DRB1_0101      WPLSSSVPSQKTYQG    3   SSSVPSQKT     0.553            test      0.001697    65.45       NA      0.447217        395.83    30.05    
+       1     DRB1_0101      PAPAPSWPLSSSVPS    4   PSWPLSSSV     0.327            test      0.000857    79.79       NA      0.327674       1442.91    54.35
+       2     DRB1_0101      APAPSWPLSSSVPSQ    3   PSWPLSSSV     0.333            test      0.001268    71.87       NA      0.346949       1171.30    50.15
+       3     DRB1_0101      PAPSWPLSSSVPSQK    4   WPLSSSVPS     0.713            test      0.002836    54.45       NA      0.412004        579.40    36.66
+       4     DRB1_0101      APSWPLSSSVPSQKT    3   WPLSSSVPS     0.773            test      0.003677    49.14       NA      0.448939        388.53    29.75
+       5     DRB1_0101      PSWPLSSSVPSQKTY    2   WPLSSSVPS     0.407            test      0.001602    66.79       NA      0.470979        306.09    25.98
+       6     DRB1_0101      SWPLSSSVPSQKTYQ    3   LSSSVPSQK     0.633            test      0.001671    65.82       NA      0.476222        289.21    25.07
+       7     DRB1_0101      WPLSSSVPSQKTYQG    3   SSSVPSQKT     0.553            test      0.001697    65.45       NA      0.447217        395.83    30.05
 
 
     """
@@ -643,3 +643,53 @@ def parse_netmhcstabpan(
         rank_index=6,
         transforms=transforms)
 
+def parse_netmhciipan43_stdout(
+        stdout,
+        prediction_method_name="netmhciipan",
+        sequence_key_mapping=None,
+        mode="elution_score"):
+    """
+    # Threshold for Strong binding peptides (%Rank)  1.00%
+    # Threshold for Weak binding peptides (%Rank)    5.00%
+
+    # DRB1_0101 : Distance to training data  0.000 (using nearest neighbor DRB1_0101)
+
+    # Allele: DRB1_0101
+    --------------------------------------------------------------------------------------------------------------------------------------------
+     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        AAAGAEAGKATTE    1   AAGAEAGKA     0.740        0        Sequence      0.000143    72.70     0.000      0.086339    95.13      19645.62
+       2         DRB1_0101      AALAAAAGVPPADKY    2   LAAAAGVPP     0.950        0        Sequence      0.030661    15.86     0.300      0.616438     7.87         63.44
+       3         DRB1_0101         EKPGNRNPYENL    3   GNRNPYENL     0.680        0        Sequence      0.000000   100.00     0.670      0.030746    98.27      35850.53
+       4         DRB1_0101      STWLLKPGAGIMIFD    2   WLLKPGAGI     0.420        0        Sequence      0.027506    16.58     0.000      0.708645     2.85         23.39
+       5         DRB1_0101      KSVPLEMLLINLTTI    4   LEMLLINLT     0.980        0        Sequence      0.007346    26.97     0.560      0.662093     4.95         38.71
+       6         DRB1_0101      KTKEDLFGKKNLIPS    5   LFGKKNLIP     0.560        0        Sequence      0.000144    72.62     0.670      0.298919    55.60       1969.51
+       7         DRB1_0101      KIYHKCDNACIGSIR    2   YHKCDNACI     0.830        0        Sequence      0.003393    34.40     0.940      0.407849    34.13        606.04
+       8         DRB1_0101      PCLFMRTVSHVILHG    3   FMRTVSHVI     0.960        0        Sequence      0.047565    13.05     0.345      0.664932     4.80         37.54
+       9         DRB1_0101      GAATVAAGAATTAAG    4   VAAGAATTA     0.920        0        Sequence      0.164981     6.81     0.000      0.514261    17.98        191.63
+  """
+    check_stdout_error(stdout, "NetMHCIIpan")
+
+    if mode not in ["elution_score", "binding_affinity"]:
+        raise ValueError("Mode is %s but must be one of: elution_score, binding affinity" % mode)
+
+    # 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,
+    }
+
+    # 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(
+        stdout=stdout,
+        prediction_method_name=prediction_method_name,
+        sequence_key_mapping=sequence_key_mapping,
+        key_index=7,
+        offset_index=0,
+        peptide_index=2,
+        allele_index=1,
+        ic50_index=13 if mode == "binding_affinity" else None,
+        rank_index=9 if mode == "elution_score" else 12,
+        score_index=8 if mode == "elution_score" else 11,
+        transforms=transforms)

tests/test_netmhcii_pan32.py

@@ -0,0 +1,21 @@
+from mhctools import NetMHCIIpan43, NetMHCIIpan43_EL, NetMHCIIpan43_BA
+
+
+def test_netmhciipan43():
+    predictor=NetMHCIIpan43(alleles=['DRB1_0101'])
+    predictions=predictor.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe()
+    assert predictions.shape==(2, 9)
+    assert all([pp==None for pp in predictions.affinity])==True ## no binding predictions
+
+def test_netmhciipan43_ba():
+    predictor_ba=NetMHCIIpan43_BA(alleles=['DRB1_0101',"HLA-DQA1*05:11-DQB1*03:02"])
+    binding_predictions=predictor_ba.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe()
+    assert binding_predictions.shape==(4, 9)
+    assert all([pp==None for pp in binding_predictions.affinity])==False # output should preturn binding predictions
+
+
+def test_netmhciipan43_el():
+    predictor_el=NetMHCIIpan43_EL(alleles=['DRB1_0101',"HLA-DQA1*05:11-DQB1*03:02"])
+    EL_predictions=predictor_el.predict_subsequences(["AAGAARIIAVDINKD","AAGIVESVGEGVTTV"]).to_dataframe()
+    assert EL_predictions.shape==(4, 9)
+    assert all([pp==None for pp in EL_predictions.affinity])==True # no affinity predictions