(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,13 @@ 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')
+ )
 
  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.
@@ -11,6 +11,116 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+import vermouth.molecule
+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
 
-class ApplyModifications():
- pass
+protein_resnames = "GLY|ALA|CYS|VAL|LEU|ILE|MET|PRO|HYP|ASN|GLN|ASP|ASP0|GLU|GLU0|THR|SER|LYS|LYS0|ARG|ARG0|HIS|HISH|PHE|TYR|TRP"
+
+def _patch_protein_termini(meta_molecule, ter_mods=['N-ter', 'C-ter']):
+ """
+ make a resspec for a protein with correct terminal modification
+ """
+ protein_termini = [({'resid': 1, 'resname': meta_molecule.nodes[0]['resname']}, ter_mods[0])]
+ max_resid = meta_molecule.max_resid
+ last_node = max_resid - 1
+ last_resname = meta_molecule.nodes[last_node]['resname']
+ if len(ter_mods) > 1:
+ last_mod = ({'resid': max_resid, 'resname': last_resname}, ter_mods[1])
+ protein_termini.append(last_mod)
+ else:
+ # if only one mod in ter_mods, apply the mod to both start and end residue
+ LOGGER.info("Only one terminal modification specified. "
+ f"Will apply {ter_mods[0]} to both {meta_molecule.nodes[0]['resname']}1 and {last_resname}{max_resid}")
+ protein_termini.append(({'resid': max_resid, 'resname': last_resname}, ter_mods[0]))
+
+ 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']] = {}
+
+ target_residue = meta_molecule.nodes[target_resid - 1]
+ # takes care to skip all residues that come from an itp file
+ if not target_residue.get('from_itp', 'False'):
+ LOGGER.warning("meta_molecule has come from itp. Will not attempt to modify.")
+ continue
+ # checks that the resname is a protein resname as defined above
+ if not vermouth.molecule.attributes_match(target_residue,
+ {'resname': vermouth.molecule.Choice(protein_resnames.split("|"))}):
+ LOGGER.warning("The resname of your target residue is not recognised a protein resname. "
+ "Will not attempt to modify.")
+ continue
+
+ anum_dict = {}
+ # this gives you the correct node indices
+ for node in target_residue['graph'].nodes:
+ aname = molecule.nodes[node]['atomname']
+ if aname in mod_atoms.keys():
+ anum_dict[aname] = node + 1
+ 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]]-1,
+ anum_dict[j.atoms[1]]-1),
+ j.parameters,
+ meta=j.meta)
+
+ return meta_molecule
+
+class ApplyModifications(Processor):
+ """
+ This processor takes a class:`polyply.src.MetaMolecule` and
+ based on modifications defined in the `force-field` attribute of the
+ MetaMolecule applies them when appropriate.
+
+ """
+ def __init__(self, meta_molecule, modifications=[]):
+ self.target_mods = []
+ for resspec, val in modifications:
+ self.target_mods.append((parse_residue_spec(resspec), val))
+ if len(self.target_mods) == 0:
+ self.target_mods = _patch_protein_termini(meta_molecule)
+
+ def run_molecule(self, meta_molecule):
+
+ apply_mod(meta_molecule, self.target_mods)
+ 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,9 @@ 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,
+ meta_molecule=meta_molecule).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):

polyply/tests/example_fixtures.py

@@ -47,6 +47,13 @@ def example_meta_molecule():
  3 8 12
  """
  force_field = vermouth.forcefield.ForceField("test")
+
+ # force_field.modifications.update({'N-ter':{'PTM_atom': False,
+ # 'replace': {'atype': 'Q5', 'charge': 1.0},
+ # 'order': 0, 'atomname': 'BB'}}
+ # )
+
+
  block_A = Block(force_field=force_field)
  block_A.add_nodes_from([(0 , {'resid': 1, 'atomname': 'BB', 'atype': 'A', 'charge': 0.0, 'other': 'A', 'resname': 'A'}),
  (1 , {'resid': 1, 'atomname': 'BB1', 'atype': 'B', 'charge': 0.0, 'other': 'A', 'resname': 'A'}),