check for ring planarity

Check for and report non-planar rings in the modelling result. This check is
also incorporated in energy minimization, potentially increasing the number
of minimization steps.

CHANGELOG

@@ -5,6 +5,14 @@
 Changelog
 ================================================================================
 
+
+Release 3.x.x
+--------------------------------------------------------------------------------
+
+* Check for and report non-planar rings in the modelling result. This check is
+  also incorporated in energy minimization, potentially increasing the number
+  of minimization steps.
+
 Release 3.2.1
 --------------------------------------------------------------------------------
 

modelling/doc/model_checking.rst

@@ -35,6 +35,15 @@ Detecting Ring Punches
 
 .. autofunction:: FilterCandidatesWithSC
 
+
+Detecting Non-Planar Rings
+--------------------------------------------------------------------------------
+
+.. autofunction:: GetNonPlanarRings
+
+.. autofunction:: HasNonPlanarRings
+
+
 Model Checking With MolProbity
 --------------------------------------------------------------------------------
 

modelling/pymod/CMakeLists.txt

@@ -27,6 +27,7 @@ set(MODELLING_PYMOD
   _mhandle_helper.py
   _alignment_fiddling.py
   _raw_model.py
+  _planar_rings.py
 )
 
 pymod(NAME modelling

modelling/pymod/__init__.py

@@ -27,3 +27,4 @@ from ._fragger_handle import *
 from ._monte_carlo import *
 from ._mhandle_helper import *
 from ._raw_model import *
+from ._planar_rings import *

modelling/pymod/_pipeline.py

@@ -25,6 +25,7 @@ from ._modelling import *
 from ._reconstruct_sidechains import *
 from ._closegaps import *
 from ._ring_punches import *
+from ._planar_rings import *
 from ._mhandle_helper import *
 # external
 import ost
@@ -363,8 +364,9 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
             ost.PopVerbosityLevel()
         # checks above would remove bad atoms
         cur_energy = sim.GetEnergy()
-        if len(temp_ent_stereo_checked.Select("ele!=H").atoms) \
-           == len(temp_ent.Select("ele!=H").atoms):
+        n_checked = len(temp_ent_stereo_checked.Select("ele!=H").atoms)
+        n_tot = len(temp_ent.Select("ele!=H").atoms)
+        if n_checked == n_tot and HasNonPlanarRings(temp_ent) == False:
             ost.LogInfo("No more stereo-chemical problems "
                         "-> final energy: %g" % cur_energy)
             break
@@ -396,6 +398,7 @@ def CheckFinalModel(mhandle):
     - Remaining stereo-chemical problems after energy minimization. The affected
       residues will have the boolean property "stereo_chemical_problem_backbone"
       set to True, if the problem affects backbone atoms.
+    - Residues that contain non-planar rings.
     
     :param mhandle: Modelling handle for which to perform checks.
     :type mhandle:  :class:`ModellingHandle`
@@ -477,6 +480,13 @@ def CheckFinalModel(mhandle):
             ost.LogInfo(msg)
             AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res])
 
+    # report non-planar rings
+    non_planar_rings = GetNonPlanarRings(mhandle.model)
+    for res in non_planar_rings:
+        msg = "Ring of " + str(res) + " is non-planar!"
+        ost.LogWarning(msg)
+        AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res])
+
 def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list(),
                       model_termini=False):
     '''Build a model starting with a raw model (see :func:`BuildRawModel`).

modelling/pymod/_planar_rings.py

+# Copyright (c) 2013-2020, SIB - Swiss Institute of Bioinformatics and
+#                          Biozentrum - University of Basel
+# 
+# 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.
+
+'''Helper functions to deal with non-planar rings'''
+import ost
+from ost import geom
+
+from promod3 import core
+
+def _GetMaxDist(res):
+    """Get max distance of any atom from ring plane
+
+    None is returned if res.one_letter_code is not in  ['Y', 'F', 'W', 'H']
+    or when not all expected ring atoms are present
+
+    :param res:  Residue for which you want to get the max dist to ring plane
+    :type res: :class:`~ost.mol.ResidueHandle` or :class:`~ost.mol.ResidueView`
+    :returns: Max distance of any ring atom to optimal ring plane or None if
+              residue contains no ring or not all expected ring atoms are
+              present.
+    """
+    prof = core.StaticRuntimeProfiler.StartScoped("planar_rings::GetMaxDist")
+    if res.one_letter_code in "YF":
+        atom_names = ["CG", "CD1", "CD2", "CE1", "CE2", "CZ"]
+    elif res.one_letter_code == 'W':
+        atom_names = ["CG", "CD1", "NE1", "CD2", "CE2", "CE3", "CZ2", "CZ3",
+                      "CH2"]
+    elif res.one_letter_code == 'H':
+        atom_names = ["CG", "CD2", "ND1", "CE1", "NE2"]
+    else:
+        return None
+    positions = geom.Vec3List()
+    for atom_name in atom_names:
+        a = res.FindAtom(atom_name)
+        if a.IsValid():
+            positions.append(a.GetPos())
+        else:
+            return None
+    origin=positions.center
+    normal=positions.principal_axes.GetRow(0)
+    plane = geom.Plane(origin, normal)
+    return max([abs(geom.Distance(plane, pos)) for pos in positions])
+
+def GetNonPlanarRings(ent, max_dist_thresh = 0.1):
+    '''Get list of residues with rings that non-planar in the given structure.
+
+    Only residues with res.one_letter_code in ['Y', 'F', 'W', 'H'] that contain
+    all expected ring atoms are considered.
+
+    Planarity is defined by the maximum distance of any ring atom to the optimal
+    plane describing the ring. This plane is constructed by defining a point on
+    the plane, here the geometric center of the ring atoms and a normal. We
+    construct an orthogonal basis [e1, e2, e3]. e1 points in direction with
+    lowest variance of ring atoms and e3 in direction with highest variance.
+    e1 is thus the plane normal. Internally this is done using a singular value
+    decomposition algorithm.
+
+    To put the default threshold of 0.1 in perspective: if you calculate these
+    distances on the same non-redundant set of experimental structures as
+    |project| used to derive its rotamer libraries, 99.9 % of the residues are
+    within 0.065 (HIS), 0.075 (TRP), 0.057 (TYR), 0.060 (PHE).
+
+    :param ent: Structure for which to detect non-planar rings.
+    :type ent:  :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
+    :param max_dist_thresh: A residue that contains a ring is considered
+                            non-planar if the max distance of any ring-atom to
+                            the optimal ring plane is above this threshold.
+    :type max_dist_thresh: :class:`float`
+
+    :return: :class:`list` of residues (:class:`~ost.mol.ResidueHandle`/
+             :class:`~ost.mol.ResidueView`) which have a non-planar ring.
+    '''
+    prof_name = 'ring_punches::GetNonPlanarRings'
+    prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
+    dists = [_GetMaxDist(res) for res in ent.residues]
+    mdt = max_dist_thresh
+    return [r for d, r in zip(dists, ent.residues) if d is not None and d > mdt]
+
+def HasNonPlanarRings(ent, max_dist_thresh = 0.1):
+    '''Check if any ring is non-planar in the given structure.
+
+    Calls :func:`GetNonPlanarRings` with given parameters and returns if any
+    residue is considered non-planar.
+
+    :param ent: Structure for which to detect non-planar rings
+    :type ent:  :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
+    :param max_dist_thresh: A residue that contains a ring is considered
+                            non-planar if the max distance of any ring-atom to
+                            the optimal ring plane is above this threshold.
+    :type max_dist_thresh: :class:`float`
+
+
+    :return: True, if any ring is non-planar
+    :rtype:  :class:`bool`
+    '''
+    prof_name = 'ring_punches::HasNonPlanarRings'
+    prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
+    return len(GetNonPlanarRings(ent, max_dist_thresh)) > 0
+
+# these methods will be exported into module
+__all__ = ('GetNonPlanarRings', 'HasNonPlanarRings')