diff --git a/modelling/doc/pipeline.rst b/modelling/doc/pipeline.rst index f5fa202666feef0d0c7bc1482cdc15f7238288a6..fb5b239c057bc2f6beda172c78a373f25a93588a 100644 --- a/modelling/doc/pipeline.rst +++ b/modelling/doc/pipeline.rst @@ -146,6 +146,23 @@ Build Raw Modelling Handle :type: :class:`~promod3.modelling.SidechainReconstructor` + .. attribute:: fragger_handles + + Optional attribute which is set in :meth:`SetFraggerHandles`. Use + :meth:`hasattr` to check if it's available. If it's set, it is used in + :meth:`BuildFromRawModel`. + + :type: :class:`list` of :class:`FraggerHandle` + + .. attribute:: modelling_issues + + Optional attribute which is set in :meth:`AddModellingIssue`. Use + :meth:`hasattr` to check if it's available. If it's set, it can be used to + check issues which occurred in :meth:`BuildFromRawModel` (see + :meth:`MinimizeModelEnergy` and :meth:`CheckFinalModel` for details). + + :type: :class:`list` of :class:`ModellingIssue` + .. method:: Copy() Generates a deep copy. Everything will be copied over to the returned @@ -422,3 +439,7 @@ Modelling Steps .. autofunction:: CheckFinalModel +.. autoclass:: ModellingIssue + :members: + +.. autofunction:: AddModellingIssue diff --git a/modelling/pymod/_mhandle_helper.py b/modelling/pymod/_mhandle_helper.py index 0e7fad147f5f749f47abd93e76c981228e228ecc..c7f1048a9b35fe7149b2c5fe54cdcf631857b102 100644 --- a/modelling/pymod/_mhandle_helper.py +++ b/modelling/pymod/_mhandle_helper.py @@ -13,9 +13,75 @@ # See the License for the specific language governing permissions and # limitations under the License. +class ModellingIssue: + """ + + :param text: Sets :attr:`text`. + :param severity: Sets :attr:`severity`. + :param residue_list: Sets :attr:`residue_list`. + + .. attribute:: text + + Description of issue. + + :type: :class:`str` + + .. attribute:: severity + + Severity of issue. + + :type: :class:`Severity` + + .. attribute:: residue_list + + List of residues affected by issue (or empty list if global issue). + + :type: :class:`list` of :class:`~ost.mol.ResidueHandle` / + :class:`~ost.mol.ResidueView` + """ + + class Severity: + """Enumerates severities.""" + #: Minor issues like remaining stereo-chemistry problems after MM. + MINOR = 0 + #: Major issues like MM-failures, incomplete models and ring punches. + MAJOR = 10 + + def __init__(self, text, severity, residue_list=[]): + self.text = text + self.severity = severity + self.residue_list = residue_list + + def is_major(self): + """ + :return: True if this is a major issue. + :rtype: :class:`bool` + """ + return self.severity == ModellingIssue.Severity.MAJOR + +def AddModellingIssue(mhandle, text, severity, residue_list=[]): + """Adds a new :class:`ModellingIssue` to + :attr:`~ModellingHandle.modelling_issues` in *mhandle*. + + If *mhandle* doesn't contain the :attr:`~ModellingHandle.modelling_issues` + attribute yet, it is added. + + :param mhandle: Will have the issue added to. + :type mhandle: :class:`ModellingHandle` + :param text: Sets :attr:`ModellingIssue.text`. + :param severity: Sets :attr:`ModellingIssue.severity`. + :param residue_list: Sets :attr:`ModellingIssue.residue_list`. + """ + modelling_issue = ModellingIssue(text, severity, residue_list) + if hasattr(mhandle, "modelling_issues"): + mhandle.modelling_issues.append(modelling_issue) + else: + mhandle.modelling_issues = [modelling_issue] + + def SetFraggerHandles(mhandle, fragger_handles): - """ Sets the :attr:`fragger_handles` in **mhandle** while ensuring - consistency with the :attr:`ModellingHandle.seqres` + """Sets :attr:`~ModellingHandle.fragger_handles` in *mhandle* while + ensuring consistency with the :attr:`ModellingHandle.seqres`. :param mhandle: Will have the fragger handles attached afterwards :param fragger_handles: The fragger handles to attach @@ -23,8 +89,8 @@ def SetFraggerHandles(mhandle, fragger_handles): :type mhandle: :class:`ModellingHandle` :type fragger_handles: :class:`list` of :class:`FraggerHandle` - :raises: :class:`ValueError` when the given **fragger_handles** are not - consistent with seqres in **mhandle** + :raises: :class:`ValueError` when the given *fragger_handles* are not + consistent with seqres in *mhandle* """ if len(mhandle.seqres) != len(fragger_handles): @@ -34,3 +100,6 @@ def SetFraggerHandles(mhandle, fragger_handles): raise RuntimeError("Sequence in FraggerHandle must match sequence "+ "in SEQRES!") mhandle.fragger_handles = fragger_handles + +# these classes/methods will be exported into module +__all__ = ('ModellingIssue', 'AddModellingIssue', 'SetFraggerHandles') diff --git a/modelling/pymod/_pipeline.py b/modelling/pymod/_pipeline.py index ec1e88f70ff43aa7ebf9d881ff4359133d0c4f9e..d8202b4266d5ae843331c02cda0070c10117f42f 100644 --- a/modelling/pymod/_pipeline.py +++ b/modelling/pymod/_pipeline.py @@ -25,6 +25,7 @@ from _modelling import * from _reconstruct_sidechains import * from _closegaps import * from _ring_punches import * +from _mhandle_helper import * # external import ost from ost import mol, conop @@ -272,6 +273,11 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20, env. variable ``PM3_OPENMM_CPU_THREADS`` to the number of desired threads. If the variable is not set, 1 thread will be used by default. + If the starting model is so bad that the energy is NaN or Inf from the start + (happens if atoms are on top of each other or almost), the energy + minimization is aborted. This issue is logged and added as a major issue to + :attr:`~ModellingHandle.modelling_issues` of *mhandle*. + :param mhandle: Modelling handle on which to apply change. :type mhandle: :class:`ModellingHandle` @@ -312,12 +318,16 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20, # check for certain failure -> we get NaN/Inf if atoms are on top each other cur_energy = sim.GetEnergy() if math.isnan(cur_energy): - ost.LogError("OpenMM could not minimize energy as atoms are on top of " - "each other!") + msg = "OpenMM could not minimize energy as atoms are on top of "\ + "each other!" + ost.LogError(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MAJOR) return if math.isinf(cur_energy): - ost.LogError("OpenMM could not minimize energy as atoms are almost " - "on top of each other!") + msg = "OpenMM could not minimize energy as atoms are almost "\ + "on top of each other!" + ost.LogError(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MAJOR) return # settings to check for stereochemical problems @@ -365,6 +375,21 @@ def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20, def CheckFinalModel(mhandle): '''Performs samity checks on final models and reports problems. + + Issues are logged and tracked in :attr:`~ModellingHandle.modelling_issues` + of *mhandle*. Major issues: + + - Chains with less than 3 residues (usually due to bad templates). + - Incomplete models (i.e. some gaps couldn't be closed). One issue is + created per unclosed gap and the stems of the gap are added to the issue. + - Complete models with sequence mismatches (should never happen). + - Residues with rings which have been punched by another bond. + + Minor issues: + + - 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. :param mhandle: Modelling handle for which to perform checks. :type mhandle: :class:`ModellingHandle` @@ -374,14 +399,19 @@ def CheckFinalModel(mhandle): # report incomplete models for chain in mhandle.model.chains: if chain.residue_count < 3: - ost.LogWarning("Chain %s of returned model contains only %d "\ - "residues! This typically indicates that the "\ - "template is mostly a Calpha trace or contains "\ - "too many D-peptides."\ - % (chain.name, chain.residue_count)) + msg = "Chain %s of returned model contains only %d residues! "\ + "This typically indicates that the template is mostly "\ + "a Calpha trace or contains too many D-peptides."\ + % (chain.name, chain.residue_count) + ost.LogWarning(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MAJOR) if len(mhandle.gaps) > 0: ost.LogWarning("Failed to close %d gap(s). Returning incomplete model!"\ % len(mhandle.gaps)) + for gap in mhandle.gaps: + AddModellingIssue(mhandle, "Failed to close gap (%s)." % (str(gap)), + ModellingIssue.Severity.MAJOR, + [gap.before, gap.after]) else: # check sequences for chain, seq in zip(mhandle.model.chains, mhandle.seqres): @@ -390,15 +420,18 @@ def CheckFinalModel(mhandle): expected_seq = seq[a-1:b] actual_seq = ''.join([r.one_letter_code for r in chain.residues]) if expected_seq != actual_seq: - ost.LogWarning("Sequence mismatch in chain %s!"\ - " Expected '%s'. Got '%s'" \ - % (chain.name, expected_seq, actual_seq)) + msg = "Sequence mismatch in chain %s! Expected '%s'. Got '%s'" \ + % (chain.name, expected_seq, actual_seq) + ost.LogWarning(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MAJOR) # report ring punchings rings = GetRings(mhandle.model) ring_punches = GetRingPunches(rings, mhandle.model) for res in ring_punches: - ost.LogWarning("Ring of " + str(res) + " has been punched!") + msg = "Ring of " + str(res) + " has been punched!" + ost.LogWarning(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MAJOR, [res]) # report stereo-chemical problems ost.PushVerbosityLevel(0) @@ -429,12 +462,14 @@ def CheckFinalModel(mhandle): for res in model_src.residues: if res.HasProp("stereo_chemical_problem_backbone") and\ res.GetBoolProp("stereo_chemical_problem_backbone"): - ost.LogInfo("Stereo-chemical problem in backbone " + \ - "of residue " + str(res)) + msg = "Stereo-chemical problem in backbone of residue " + str(res) + ost.LogInfo(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res]) elif res.HasProp("stereo_chemical_problem") and\ res.GetBoolProp("stereo_chemical_problem"): - ost.LogInfo("Stereo-chemical problem in sidechain " + \ - "of residue " + str(res)) + msg = "Stereo-chemical problem in sidechain of residue " + str(res) + ost.LogInfo(msg) + AddModellingIssue(mhandle, msg, ModellingIssue.Severity.MINOR, [res]) def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list(), model_termini=False): diff --git a/modelling/tests/test_pipeline.py b/modelling/tests/test_pipeline.py index 65a43b325782640bd7720cc197bdb80402c217ea..ac4e4bac43404c911e311fdffdaf554f1870a4a7 100644 --- a/modelling/tests/test_pipeline.py +++ b/modelling/tests/test_pipeline.py @@ -235,15 +235,27 @@ class PipelineTests(unittest.TestCase): self.assertEqual(log.messages['WARNING'][-1], "Failed to close %d gap(s). " \ "Returning incomplete model!" % exp_gaps) + self.assertGreaterEqual(len(mhandle.modelling_issues), exp_gaps) + for i in range(exp_gaps): + issue = mhandle.modelling_issues[i] + self.assertTrue(issue.text.startswith("Failed to close")) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 2) # fake remove gaps to get sequence mismatch warnings mhandle.gaps = modelling.StructuralGapList() + mhandle.modelling_issues = list() log.messages['WARNING'] = list() modelling.CheckFinalModel(mhandle) - # 1 warning per chain expected + # 1 warning/issue per chain expected exps = "Sequence mismatch in chain" self.assertGreaterEqual(len(log.messages['WARNING']), num_chains) + self.assertGreaterEqual(len(mhandle.modelling_issues), num_chains) for i in range(num_chains): self.assertTrue(log.messages['WARNING'][i].startswith(exps)) + issue = mhandle.modelling_issues[i] + self.assertTrue(issue.text.startswith(exps)) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 0) ####################################################################### @@ -305,7 +317,7 @@ class PipelineTests(unittest.TestCase): modelling.MinimizeModelEnergy(mhandle) self.compare(mhandle.model, 'data/1crn_final.pdb') - def testMinimizeModelEnergy(self): + def testMinimizeModelEnergyFails(self): '''Check energy minimization fails.''' # setup log log = _FetchLog() @@ -316,17 +328,27 @@ class PipelineTests(unittest.TestCase): log.messages['ERROR'] = list() modelling.MinimizeModelEnergy(mhandle) self.assertEqual(len(log.messages['ERROR']), 1) - self.assertEqual(log.messages['ERROR'][0], - "OpenMM could not minimize energy as atoms are on " - "top of each other!") + self.assertEqual(len(mhandle.modelling_issues), 1) + exp_msg = "OpenMM could not minimize energy as atoms are on "\ + "top of each other!" + self.assertEqual(log.messages['ERROR'][0], exp_msg) + issue = mhandle.modelling_issues[0] + self.assertEqual(issue.text, exp_msg) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 0) # catch atoms almost on top mhandle = self.getMockModel(io.LoadPDB('data/gly_almost_on_top.pdb')) log.messages['ERROR'] = list() modelling.MinimizeModelEnergy(mhandle) self.assertEqual(len(log.messages['ERROR']), 1) - self.assertEqual(log.messages['ERROR'][0], - "OpenMM could not minimize energy as atoms are " - "almost on top of each other!") + self.assertEqual(len(mhandle.modelling_issues), 1) + exp_msg = "OpenMM could not minimize energy as atoms are "\ + "almost on top of each other!" + self.assertEqual(log.messages['ERROR'][0], exp_msg) + issue = mhandle.modelling_issues[0] + self.assertEqual(issue.text, exp_msg) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 0) def testBuildFromRawModel(self): '''Check that BuildFromRawModel works.''' @@ -387,8 +409,14 @@ class PipelineTests(unittest.TestCase): log.messages['WARNING'] = list() modelling.CheckFinalModel(mhandle) self.assertEqual(len(log.messages['WARNING']), 1) - self.assertEqual(log.messages['WARNING'][0], - "Ring of A.TYR321 has been punched!") + self.assertGreaterEqual(len(mhandle.modelling_issues), 1) + exp_msg = "Ring of A.TYR321 has been punched!" + self.assertEqual(log.messages['WARNING'][0], exp_msg) + issue = mhandle.modelling_issues[0] + self.assertEqual(issue.text, exp_msg) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 1) + self.assertEqual(str(issue.residue_list[0]), "A.TYR321") def testCheckFinalModelClash(self): '''Check that we report models with stereo-chemical problems.''' @@ -402,12 +430,21 @@ class PipelineTests(unittest.TestCase): log.messages['INFO'] = list() modelling.CheckFinalModel(mhandle) self.assertEqual(len(log.messages['INFO']), 2) - self.assertEqual(log.messages['INFO'][0], - "Stereo-chemical problem in backbone of residue "\ - "A.GLY79") - self.assertEqual(log.messages['INFO'][1], - "Stereo-chemical problem in sidechain of residue "\ - "A.HIS98") + self.assertEqual(len(mhandle.modelling_issues), 2) + exp_msg = "Stereo-chemical problem in backbone of residue A.GLY79" + self.assertEqual(log.messages['INFO'][0], exp_msg) + issue = mhandle.modelling_issues[0] + self.assertEqual(issue.text, exp_msg) + self.assertFalse(issue.is_major()) + self.assertEqual(len(issue.residue_list), 1) + self.assertEqual(str(issue.residue_list[0]), "A.GLY79") + exp_msg = "Stereo-chemical problem in sidechain of residue A.HIS98" + self.assertEqual(log.messages['INFO'][1], exp_msg) + issue = mhandle.modelling_issues[1] + self.assertEqual(issue.text, exp_msg) + self.assertFalse(issue.is_major()) + self.assertEqual(len(issue.residue_list), 1) + self.assertEqual(str(issue.residue_list[0]), "A.HIS98") def testMolProbity(self): '''Check if binding to molprobity works.''' @@ -503,11 +540,15 @@ class PipelineTests(unittest.TestCase): log.messages['WARNING'] = list() modelling.CheckFinalModel(mhandle) self.assertEqual(len(log.messages['WARNING']), 1) - self.assertEqual(log.messages['WARNING'][-1], - "Chain A of returned model contains only 2 "\ - "residues! This typically indicates that the "\ - "template is mostly a Calpha trace or contains "\ - "too many D-peptides.") + self.assertEqual(len(mhandle.modelling_issues), 1) + exp_msg = "Chain A of returned model contains only 2 residues! "\ + "This typically indicates that the template is mostly a "\ + "Calpha trace or contains too many D-peptides." + self.assertEqual(log.messages['WARNING'][-1], exp_msg) + issue = mhandle.modelling_issues[0] + self.assertEqual(issue.text, exp_msg) + self.assertTrue(issue.is_major()) + self.assertEqual(len(issue.residue_list), 0) if __name__ == "__main__": from ost import testutils