SCHWED-843: allow to call Phenix.MolProbity from PM3

modelling/doc/index.rst

@@ -115,6 +115,12 @@ Modelling Pipeline
 
 .. autofunction:: CheckFinalModel
 
+.. autofunction:: RunMolProbity
+
+.. autofunction:: RunMolProbityEntity
+
+.. autofunction:: ReportMolProbityScores
+
 
 Closing Gaps
 --------------------------------------------------------------------------------

modelling/pymod/CMakeLists.txt

@@ -9,6 +9,7 @@ export_extension_scheme.cc
 set(MODELLING_PYMOD
   __init__.py
   _closegaps.py
+  _molprobity.py
   _pipeline.py
   _ring_punches.py
 )

modelling/pymod/__init__.py

 """Initialise the modelling module."""
 from _modelling import *
 from _closegaps import *
+from _molprobity import *
 from _pipeline import *
 from _ring_punches import *

modelling/pymod/_molprobity.py

+'''Binding to external MolProbity tool to get scores.'''
+
+import ost
+from ost import settings
+import subprocess, tempfile, os
+
+def RunMolProbity(target_pdb, molprobity_bin=None):
+    '''Run ``MolProbity`` from ``Phenix`` on a given PDB file.
+
+    MolProbity score computation: (taken from molprobity source code)
+
+    .. code-block:: python
+
+        clashscore = result["Clashscore"]
+        rota_out = result["Rotamer outliers"]
+        rama_iffy = 100. - result["Ramachandran favored"]
+        mpscore = (( 0.426 * math.log(1 + clashscore) ) +
+                 ( 0.33 * math.log(1 + max(0, rota_out - 1)) ) +
+                 ( 0.25 * math.log(1 + max(0, rama_iffy - 2)) )) + 0.5
+    
+    :param target_pdb: Path to PDB file on which to do analysis.
+    :type target_pdb:  :class:`str`
+    
+    :param molprobity_bin: Path to ``phenix.molprobity`` executable. If None, it
+                           searches for it in the ``PATH`` or (if set) in the
+                           env. variable ``MOLPROBITY_EXECUTABLE``.
+                           The function was tested with ``Phenix 1.9-1692`` and
+                           with ``MolProbity 4.2`` which also includes it.
+    :type molprobity_bin:  :class:`str`
+
+    :return: Dictionary with scores produced by MolProbity. Entries:
+
+             - "Ramachandran outliers" (percentage [0,100] as :class:`float`)
+             - "Ramachandran favored" (percentage [0,100] as :class:`float`)
+             - "Rotamer outliers" (percentage [0,100] as :class:`float`)
+             - "C-beta deviations" (:class:`integer`)
+             - "Clashscore" (:class:`float`)
+             - "MolProbity score" (:class:`float`)
+             - "RMS(bonds)" (:class:`float`)
+             - "RMS(angles)" (:class:`float`)
+
+    :rtype: :class:`dict`
+
+    :raises: :class:`~ost.settings.FileNotFound` if the "phenix.molprobity"
+             executable is not found.
+    '''
+    # HELPER
+    def GetStringValue(my_line):
+        """Parse line formatted as ' bla  = X ...' and return 'X' as string."""
+        return my_line.split('=')[1].strip().split(' ')[0]
+
+    # locate molprobity
+    molprobity_exec = settings.Locate("phenix.molprobity",
+                                      explicit_file_name=molprobity_bin,
+                                      env_name='MOLPROBITY_EXECUTABLE')
+    # run molprobity avoiding any file output
+    cmd = [molprobity_exec, target_pdb, "output.quiet=True", \
+           "output.coot=False", "output.probe_dots=False"]
+    job = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+    sout, _ = job.communicate()
+    # parse
+    result = dict()
+    first_found = False
+    lines = sout.splitlines()
+    for line in lines:
+        my_line = line.strip().lower()
+        if my_line.startswith("ramachandran outliers"):
+            first_found = True
+            result["Ramachandran outliers"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("favored"):
+            result["Ramachandran favored"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("rotamer outliers"):
+            result["Rotamer outliers"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("clashscore"):
+            result["Clashscore"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("c-beta deviations"):
+            result["C-beta deviations"] = int(GetStringValue(line))
+        elif first_found and my_line.startswith("molprobity score"):
+            result["MolProbity score"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("rms(bonds)"):
+            result["RMS(bonds)"] = float(GetStringValue(line))
+        elif first_found and my_line.startswith("rms(angles)"):
+            result["RMS(angles)"] = float(GetStringValue(line))
+    return result
+
+
+def RunMolProbityEntity(ost_ent, molprobity_bin=None):
+    '''Run molprobity from phenix on given OST entity.
+    
+    See :func:`RunMolProbity` for details.
+
+    :param ost_ent: OST entity on which to do analysis.
+    :type ost_ent:  :class:`Entity <ost.mol.EntityHandle>`
+    '''
+    # store as PDB in a tmp. file
+    _, tmpfile = tempfile.mkstemp(suffix='.pdb')
+    ost.PushVerbosityLevel(0)
+    ost.io.SavePDB(ost_ent, tmpfile)
+    ost.PopVerbosityLevel()
+    # get result
+    result = RunMolProbity(tmpfile, molprobity_bin)
+    # clean up
+    os.remove(tmpfile)
+    return result
+
+def ReportMolProbityScores(scores):
+    '''Print MolProbity score and its components to LogInfo.
+    
+    :param scores: MolProbity scores as generated by :func:`RunMolProbity`.
+    :type scores:  :class:`dict`
+    '''
+    ost.LogInfo("MolProbity score: "  + str(scores["MolProbity score"]))
+    ost.LogInfo("- Clashscore          : " + str(scores["Clashscore"]))
+    ost.LogInfo("- Rotamer outliers    : " + str(scores["Rotamer outliers"]) + " %")
+    ost.LogInfo("- Ramachandran favored: " + str(scores["Ramachandran favored"]) + " %")
+
+__all__ = ('RunMolProbity', 'RunMolProbityEntity', 'ReportMolProbityScores')

modelling/pymod/_pipeline.py

@@ -113,6 +113,9 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
 
     :param max_iter_lbfgs: Max. number of iterations within LBFGS method
     :type max_iter_lbfgs:  :class:`int`
+
+    :return: The model including all oxygens as used in the minimizer.
+    :rtype:  :class:`Entity <ost.mol.EntityHandle>`
     '''
     ost.LogInfo("Minimize energy.")
     # ignore LogInfo in stereochemical problems if output up to info done
@@ -178,6 +181,8 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
     mhandle.model = mol.CreateEntityFromView(\
                             simulation_ent.Select("peptide=true and ele!=H"),
                             True)
+    # return model with oxygens
+    return mol.CreateEntityFromView(simulation_ent.Select("peptide=true"), True)
 
 def CheckFinalModel(mhandle):
     '''Performs samity checks on final models and reports problems.

modelling/tests/test_pipeline.py

@@ -5,7 +5,7 @@ We go from a fake cut in crambin to reconstruct that same protein.
 import unittest
 import math
 import ost
-from ost import io, mol, seq
+from ost import io, mol, seq, settings
 from promod3 import loop, modelling
 
 ################################################################
@@ -303,6 +303,27 @@ class PipelineTests(unittest.TestCase):
         mhandle = self.getRawModelOligo("data/5d52-1_cut")
         self.checkFinalModel(mhandle, exp_gaps=2, num_chains=2)
 
+    def testMolProbity(self):
+        '''Check if binding to molprobity works.'''
+        # do we have it?
+        try:
+            settings.Locate("phenix.molprobity",
+                            env_name='MOLPROBITY_EXECUTABLE')
+        except settings.FileNotFound:
+            print 'phenix.molprobity is missing. Please install it and make '+\
+                  'it available in your PATH to execute test testMolProbity.'
+            return
+        # load 1crn and evaluate it
+        prot = io.LoadPDB('1crn', remote=True)
+        results = modelling.RunMolProbityEntity(prot)
+        # check output
+        self.assertAlmostEqual(results['Ramachandran favored'], 97.73, places=2)
+        self.assertAlmostEqual(results['MolProbity score'], 0.56, places=2)
+        self.assertEqual(results['C-beta deviations'], 0)
+        self.assertAlmostEqual(results['Clashscore'], 0.0, places=2)
+        self.assertAlmostEqual(results['Ramachandran outliers'], 0.0, places=2)
+        self.assertAlmostEqual(results['Rotamer outliers'], 0.0, places=2)
+
 if __name__ == "__main__":
     from ost import testutils
     testutils.RunTests()