(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,104 @@
 # 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
+ """
+
+ resids = nx.get_node_attributes(meta_molecule.molecule, 'resid')
+ last_resid = max([val for val in resids.values()])
+
+ resnames = nx.get_node_attributes(meta_molecule.molecule, 'resname')
+ last_resname = [val for val in resnames.values()][-1]
+
+ protein_termini = [({'resid': 1, 'resname': resnames[0]}, 'N-ter'),
+ ({'resid': last_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
+ """
+
+ if not modifications:
+ return 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(meta_molecule.molecule.nodes):
+ if meta_molecule.molecule.nodes[node]['resid'] == target_resid:
+ resname = meta_molecule.molecule.nodes[node]['resname']
+ meta_molecule.molecule.nodes[node]['resname'] = resname + 'M'
+ # add the modification from the modifications dict
+ aname = meta_molecule.molecule.nodes[node]['atomname']
+ if aname in mod_atoms.keys():
+ anum_dict[aname] = atomid
+ for key, value in mod_atoms[aname].items():
+ meta_molecule.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):