Merge pull request #593 from csbrasnett/idp-go

Idp go

bin/martinize2

@@ -170,7 +170,7 @@ def pdb_to_universal(
     return canonicalized
 
 
-def martinize(system, mappings, to_ff, delete_unknown=False):
+def martinize(system, mappings, to_ff, delete_unknown=False, idrs=False, disordered_regions=None):
     """
     Convert a system from one force field to an other at lower resolution.
     """
@@ -185,6 +185,9 @@ def martinize(system, mappings, to_ff, delete_unknown=False):
     ).run_system(system)
     LOGGER.info("Averaging the coordinates.", type="step")
     vermouth.DoAverageBead(ignore_missing_graphs=True).run_system(system)
+    if idrs:
+        LOGGER.info("Annotating IDRs.", type="step")
+        vermouth.AnnotateIDRs(id_regions=[(int(start), int(stop)) for start, stop in disordered_regions]).run_system(system)
     LOGGER.info("Applying the links.", type="step")
     vermouth.DoLinks().run_system(system)
     LOGGER.info("Placing the charge dummies.", type="step")
@@ -612,6 +615,13 @@ def entry():
               "format: <start_resid_1>:<end_resid_1> <start_resid_2>:<end_resid_2>..."
              ),
     )
+    idr_tuning = water_group.add_argument(
+        "-idr-tune",
+        dest="idr_tune",
+        action="store_true",
+        default=False,
+        help=("Tune the idr regions with specific bonded potentials."),
+    )
 
 
     prot_group = parser.add_argument_group("Protein description")
@@ -922,6 +932,13 @@ def entry():
         )
     vermouth.SetMoleculeMeta(scfix=args.scfix).run_system(system)
 
+    vermouth.SetMoleculeMeta(idr=args.idr_tune).run_system(system)
+    if args.idr_tune:
+        if not target_ff.has_feature("idr"):
+            LOGGER.warning('Improved IDR potentials are not implemented '
+                           'for this force field',
+                           type="missing-feature")
+
     ss_sequence = list(
         itertools.chain(
             *(
@@ -943,6 +960,8 @@ def entry():
         mappings=known_mappings,
         to_ff=known_force_fields[args.to_ff],
         delete_unknown=True,
+        idrs=args.idr_tune,
+        disordered_regions=args.water_idrs
     )
 
     # Apply position restraints if required.
@@ -1051,6 +1070,7 @@ def entry():
             vermouth.rcsu.go_vs_includes.VirtualSiteCreator().run_system(system)
             itp_paths = {"atomtypes": "virtual_sites_atomtypes.itp",
                          "nonbond_params": "virtual_sites_nonbond_params.itp"}
+
         # now we add a bias by defining specific virtual-site water interactions
         vermouth.processors.ComputeWaterBias(args.water_bias,
                                              { s:float(eps) for s, eps in args.water_bias_eps},

vermouth/data/force_fields/martini3001/aminoacids.ff

@@ -962,6 +962,73 @@ scfix true
 [ angles ]
 BB +BB +SC1 10 100 15 {"group": "SC-BB-BB and BB-BB-SC scFix", "comment": "BB-BB-SC"}
 
+;;; IDR specific links
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ angles ]
+BB +BB +SC1 10 85 10 {"group": "idp-fix"}
+
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ angles ]
+-SC1 -BB BB 10 85 10 {"group": "idp-fix"}
+
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ dihedrals ]
+BB +BB ++BB +++BB 9 -120 -1 1 {"group": "idp-fix", "comment": "BB-BB-BB-BB-v1", "version":1}
+BB +BB ++BB +++BB 9 -120 -1 2 {"group": "idp-fix", "comment": "BB-BB-BB-BB-v2", "version":2}
+
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ dihedrals ]
+SC1 BB +BB +SC1 9 -130 -1.5 1 {"group": "idp-fix", "version": 1, "comment": "SC1-BB-BB-SC1-v1"}
+SC1 BB +BB +SC1 9  100 -1.5 2 {"group": "idp-fix", "version": 2, "comment": "SC1-BB-BB-SC1-v2"}
+
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ dihedrals ]
+-SC1 -BB BB +BB 1 115 -4.5 1 {"group": "idp-fix", "comment": "SC1-BB-BB(GLY)-BB"}
+[ non-edges ]
+BB SC1
+
+[ link ]
+resname $protein_resnames
+cgidr true
+[ features ]
+idr
+[ molmeta ]
+idr true
+[ dihedrals ]
+-BB BB +BB +SC1 1 0 -2.0 1 {"group": "idp-fix", "comment": "BB-BB(GLY)-BB-SC1"}
+[ non-edges ]
+BB SC1
+
 [ link ]
 resname $protein_resnames
 [ angles ]

vermouth/processors/__init__.py

@@ -44,3 +44,4 @@
 from .merge_all_molecules import MergeAllMolecules
 from .annotate_mut_mod import AnnotateMutMod
 from .water_bias import ComputeWaterBias
+from .annotate_idrs import AnnotateIDRs

vermouth/processors/annotate_idrs.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# 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.
+# You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# 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.
+"""
+Provides processors that can add and remove IDR specific bonds
+"""
+
+from .processor import Processor
+from ..rcsu.go_utils import _in_resid_region
+from ..log_helpers import StyleAdapter, get_logger
+LOGGER = StyleAdapter(get_logger(__name__))
+
+
+def annotate_disorder(molecule, id_regions, annotation="cgidr"):
+    """
+    Annotate the disordered regions of the molecule
+
+    molecule: :class:`vermouth.molecule.Molecule`
+            the molecule
+    idr_regions: list
+        list of tuples defining the resids of the idrs in the molecule
+    annotation: str
+        name of the annotation in the node
+    """
+
+    for key, node in molecule.nodes.items():
+        _old_resid = node['_old_resid']
+        if _in_resid_region(_old_resid, id_regions):
+            molecule.nodes[key][annotation] = True
+
+class AnnotateIDRs(Processor):
+    """
+    Processor to annotate intrinsically disordered regions of a molecule.
+
+    This processor is designed primarily for the work described in the reference
+    M3_GO, but is generally applicable for such circumstances where extra
+    addition/removals are necessary.
+
+    """
+
+    def __init__(self, id_regions=None):
+        """
+        Parameters
+        ----------
+        id_regions:
+            regions defining the IDRs
+        """
+        self.id_regions = id_regions
+
+    def run_molecule(self, molecule):
+        """
+        Assign disordered regions for a single molecule
+        """
+
+        annotate_disorder(molecule, self.id_regions)
+
+        return molecule
+
+    def run_system(self, system):
+        """
+        Assign the water bias of the Go model to file. Biasing
+        is always molecule specific i.e. no two different
+        vermouth molecules can have the same bias.
+
+        Parameters
+        ----------
+        system: :class:`vermouth.system.System`
+        """
+        if not self.id_regions:
+            return system
+        LOGGER.info("Annotating disordered regions", type="step")
+        super().run_system(system)

vermouth/processors/water_bias.py

@@ -16,6 +16,7 @@
 from ..graph_utils import make_residue_graph
 from ..rcsu.go_utils import get_go_type_from_attributes, _get_bead_size, _in_resid_region
 from ..gmx.topology import NonbondParam
+import numpy as np
 
 class ComputeWaterBias(Processor):
     """
@@ -54,8 +55,8 @@ def __init__(self,
         water_bias: dict[str, float]
             a dict of secondary structure codes and
             epsilon value for the water bias in kJ/mol
-        idr_regions:
-            regions defining the IDRs
+        idr_regions: list
+            list of tuples of residue regions defining the IDRs
         prefix: str
             prefix of the Go virtual-site atomtypes
         system: vermouth.system.System
@@ -118,6 +119,41 @@ def assign_residue_water_bias(self, molecule, res_graph):
                                       meta={"comment": ["water bias", sec_struc]})
             self.system.gmx_topology_params["nonbond_params"].append(water_bias)
 
+    def remove_cross_nb_interactions(self, molecule, res_graph):
+        """
+        Remove Go bonds between folded and disordered regions of a molecule
+
+        Parameters
+        ----------
+        molecule: :class:`vermouth.molecule.Molecule`
+            the molecule
+        res_graph: :class:`vermouth.molecule.Molecule`
+            the residue graph of the molecule
+        """
+        #list of all the Go pairs in the molecule
+        all_go_pairs = np.array([list(i.atoms) for i in self.system.gmx_topology_params["nonbond_params"] if 'W' not in list(i.atoms)])
+        # list to record which items we don't want. cross = go potential between folded and disordered domain.
+        all_cross_pairs = []
+
+        for res_node in res_graph.nodes:
+            resid = res_graph.nodes[res_node]['resid']
+            _old_resid = res_graph.nodes[res_node]['_old_resid']
+            chain = res_graph.nodes[res_node]['chain']
+
+            if _in_resid_region(_old_resid, self.idr_regions):
+                vs_go_node = next(get_go_type_from_attributes(res_graph.nodes[res_node]['graph'],
+                                                              resid=resid,
+                                                              chain=chain,
+                                                              prefix=molecule.meta.get('moltype')))
+                all_cross_pairs.append(np.where(all_go_pairs == vs_go_node)[0]) #just need the first one
+
+        # make sure we only have one entry in case a site has more than one interaction
+        all_cross_pairs = np.unique([x for xs in all_cross_pairs for x in xs])
+        # delete the folded-disordered Go interactions from the list going backwards.
+        # otherwise list order gets messed up.
+        for i in reversed(all_cross_pairs):
+            del self.system.gmx_topology_params["nonbond_params"][i]
+
     def run_molecule(self, molecule):
         """
         Assign water bias for a single molecule
@@ -134,6 +170,7 @@ def run_molecule(self, molecule):
             res_graph = make_residue_graph(molecule)
 
         self.assign_residue_water_bias(molecule, res_graph)
+        self.remove_cross_nb_interactions(molecule, res_graph)
 
         return molecule
 

vermouth/tests/data/integration_tests/tier-1/hst5/README

@@ -6,3 +6,4 @@
 
 # -ff martini 3001
 # BB-W bias for idr
+# add extra IDR potentials

vermouth/tests/data/integration_tests/tier-1/hst5/martinize2/command

@@ -7,3 +7,4 @@ martinize2
 -water-bias
 -water-bias-eps idr:0.5
 -id-regions 1:24
+-idr-tune