(protein) termini patching

.gitignore

@@ -24,3 +24,5 @@ tests/test_data/gen_seq/output/*.itp
 doc/build
 doc/source/api
 doc/source/.doctrees
+
+/polyply/data/martini3007/*

bin/polyply

@@ -77,6 +77,20 @@ def main(): # pylint: disable=too-many-locals,too-many-statements
     dna_group = parser_gen_itp.add_argument_group('DNA specifc options')
     dna_group.add_argument('-dsdna', dest='dsdna', action='store_true',
                            help='complement single sequence to dsDNA sequence')
+    modifications_group = parser_gen_itp.add_argument_group('Modifications group')
+    modifications_group.add_argument('-mods', dest='mods', action='append',
+                                     default=[], type=lambda s: s.split(":"),
+                                     help=('Add a modification to a residue. The format is '
+                                           '<resname><resid>:modification_name, '
+                                           'e.g. ASP1:N-ter')
+                                     )
+    modifications_group.add_argument('-protter', dest='protter',
+                                     action='store_true',
+                                     help=('Apply N-ter and C-ter modifications to a protein. '
+                                           'Alias for "-mods RES1:N-ter -mods RESN:C-ter", '
+                                           'where RES1 is the name and resid of the first '
+                                           'residue in the protein and RESN is the name and resid of the last')
+                                    )
 
     parser_gen_itp.set_defaults(func=gen_itp)
 

polyply/data/martini3/aminoacids.ff

@@ -803,59 +803,6 @@ CA  BB
 [ edges ]
 CA BB
 
-[ link ]
-resname $protein_resnames
-[ atoms ]
-CA { }
-BB { }
-+BB { }
-[ exclusions ]
-#meta {"comment": "CA-BB"}
-CA +BB
-[ edges ]
-BB +BB
-
-;; Protein termini. These links should be applied last.
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "Q5", "charge": 1}}
-[ non-edges ]
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "Q5", "charge": -1}}
-[ non-edges ]
-BB +BB
-
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "P6", "charge": 0}}
-[ non-edges ]
-BB +BB
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ molmeta ]
-neutral_termini true
-[ atoms ]
-BB {"replace": {"atype": "P5", "charge": 0}}
-[ non-edges ]
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ molmeta ]
-neutral_termini true
-[ atoms ]
-BB {"replace": {"atype": "P6", "charge": 0}}
-[ non-edges ]
-BB +BB
-
 ;; Cystein bridge
 [ link ]
 resname "CYS"

polyply/data/martini3/modifications.ff

@@ -0,0 +1,37 @@
+[ citations ]
+Martini3
+
+[ modification ]
+C-ter
+[ atoms ]
+BB {"replace": {"atype": "Q5", "charge": -1.0}}
+
+[ modification ]
+N-ter
+[ atoms ]
+BB {"replace": {"atype": "Q5", "charge": 1.0}}
+
+[ modification ]
+zwitter
+[ atoms ]
+BB {"replace": {"atype": "Q5"}}
+
+[ modification ]
+COOH-ter
+[ atoms ]
+BB {"replace": {"atype": "P6", "charge": 0.0}}
+
+[ modification ]
+NH2-ter
+[ atoms ]
+BB {"replace": {"atype": "P6", "charge": 0.0}}
+
+[ modification ]
+CCAP-ter
+[ atoms ]
+BB {"replace": {"atype": "P1", "charge": 0.0}}
+
+[ modification ]
+NCAP-ter
+[ atoms ]
+BB {"replace": {"atype": "P2", "charge": 0.0}}

polyply/src/apply_links.py

@@ -230,16 +230,16 @@ def match_link_and_residue_atoms(meta_molecule, link, link_to_resid):
         # relative resid has been asserted before so we can
         # exclude it here
         ignore = ['order', 'charge_group', 'replace', 'resid']
-        matchs = list(find_atoms(block, ignore=ignore, **attrs))
+        matches = list(find_atoms(block, ignore=ignore, **attrs))
 
-        if len(matchs) == 1:
-            link_to_mol[node] = matchs[0]
-        elif len(matchs) == 0:
-            msg = "Found no matchs for node {} in resiue {}. Cannot apply link."
+        if len(matches) == 1:
+            link_to_mol[node] = matches[0]
+        elif len(matches) == 0:
+            msg = "Found no matches for node {} in resiue {}. Cannot apply link."
             raise MatchError(msg.format(node, resid))
         else:
             msg = "Found {} matches for node {} in resiue {}. Cannot apply link."
-            raise MatchError(msg.format(len(matchs), node, resid))
+            raise MatchError(msg.format(len(matches), node, resid))
 
     return link_to_mol
 
@@ -475,8 +475,8 @@ def run_molecule(self, meta_molecule):
                                              res_link,
                                              node_match=_res_match,
                                              edge_match=_linktype_match)
-            raw_matchs = GM.subgraph_isomorphisms_iter()
-            for match in raw_matchs:
+            raw_matches = GM.subgraph_isomorphisms_iter()
+            for match in raw_matches:
                 nodes = match.keys()
                 resids =[meta_molecule.nodes[node]["resid"] for node in nodes]
                 orders = [ res_link.nodes[match[node]]["order"] for node in nodes]

polyply/src/apply_modifications.py

@@ -1,4 +1,4 @@
-# Copyright 2020 University of Groningen
+# Copyright 2024 University of Groningen
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -12,5 +12,96 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-class ApplyModifications():
-      pass
+from vermouth.log_helpers import StyleAdapter, get_logger
+from vermouth.processors.annotate_mut_mod import parse_residue_spec
+LOGGER = StyleAdapter(get_logger(__name__))
+from .processor import Processor
+import networkx as nx
+
+def _get_protein_termini(meta_molecule):
+    """
+    make a resspec for a protein with correct N-ter and C-ter
+    """
+    max_resid = meta_molecule.max_resid
+    last_node = max_resid - 1
+    last_resname = meta_molecule.nodes[last_node]['resname']
+
+    protein_termini = [({'resid': 1, 'resname': meta_molecule.nodes[0]['resname']}, 'N-ter'),
+                       ({'resid': max_resid, 'resname': last_resname}, 'C-ter')
+                       ]
+
+    return protein_termini
+
+def apply_mod(meta_molecule, modifications):
+    """
+
+    Apply a terminal modification.
+    Note this assumes that the modification is additive, ie. no atoms are
+    removed
+
+    Parameters
+    ----------
+    meta_molecule: :class:`polyply.src.meta_molecule.MetaMolecule`
+    modifications: list
+        list of (resspec, modification) pairs to apply
+
+    Returns
+    ----------
+    meta_molecule
+    """
+
+    molecule = meta_molecule.molecule
+
+    if not molecule.force_field.modifications:
+        LOGGER.warning('No modifications present in forcefield, none will be applied')
+        return meta_molecule
+
+    for target, desired_mod in modifications:
+
+        target_resid = target['resid']
+
+        mod_atoms = {}
+        for mod_atom in molecule.force_field.modifications[desired_mod].atoms:
+            if 'replace' in mod_atom:
+                mod_atoms[mod_atom['atomname']] = mod_atom['replace']
+            else:
+                mod_atoms[mod_atom['atomname']] = {}
+
+        anum_dict = {}
+        for atomid, node in enumerate(molecule.nodes):
+            if molecule.nodes[node]['resid'] == target_resid:
+                # add the modification from the modifications dict
+                aname = molecule.nodes[node]['atomname']
+                if aname in mod_atoms.keys():
+                    anum_dict[aname] = atomid
+                    for key, value in mod_atoms[aname].items():
+                        molecule.nodes[node][key] = value
+
+        mod_interactions = molecule.force_field.modifications[desired_mod].interactions
+        for i in mod_interactions:
+            for j in molecule.force_field.modifications[desired_mod].interactions[i]:
+                molecule.add_interaction(i,
+                                         (anum_dict[j.atoms[0]],
+                                          anum_dict[j.atoms[1]]),
+                                         j.parameters,
+                                         meta=j.meta)
+
+    return meta_molecule
+
+class ApplyModifications(Processor):
+    def __init__(self, modifications=None, protter=False):
+        if not modifications:
+            modifications = []
+        self.modifications = []
+        for resspec, val in modifications:
+            self.modifications.append((parse_residue_spec(resspec), val))
+        self.protter = protter
+
+    def run_molecule(self, meta_molecule):
+
+        if self.protter:
+            self.modifications = _get_protein_termini(meta_molecule)
+
+        apply_mod(meta_molecule, self.modifications)
+
+        return meta_molecule

polyply/src/gen_itp.py

@@ -34,6 +34,7 @@
 from polyply.src.graph_utils import find_missing_edges
 from .load_library import load_ff_library
 from .gen_dna import complement_dsDNA
+from .apply_modifications import ApplyModifications
 
 LOGGER = StyleAdapter(get_logger(__name__))
 
@@ -60,7 +61,10 @@ def split_seq_string(sequence):
         monomers.append(Monomer(resname=resname, n_blocks=n_blocks))
     return monomers
 
-def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[], lib=None, seq=None, seq_file=None, dsdna=False):
+def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[],
+               lib=None, seq=None, seq_file=None,
+               dsdna=False, mods=[], protter=False):
+
     """
     Top level function for running the polyply parameter generation.
     Parameters seq and seq_file are mutually exclusive. Set the other
@@ -107,6 +111,8 @@ def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[], lib=None,
     LOGGER.info("applying links between residues",  type="step")
     meta_molecule = ApplyLinks().run_molecule(meta_molecule)
 
+    meta_molecule = ApplyModifications(modifications=mods, protter=protter).run_molecule(meta_molecule)
+
     # Raise warning if molecule is disconnected
     msg = "Missing a link between residue {idxA} {resA} and residue {idxB} {resB}."
     for missing in find_missing_edges(meta_molecule, meta_molecule.molecule):