Select Git revision
_pipeline.py 13.09 KiB
'''High-level functionality for modelling module to build pipelines. Added in
the __init__.py file. To be used directly by passing a ModellingHandle instance
as argument.
'''
# internal
from promod3 import loop
from promod3 import sidechain
from _modelling import *
from _closegaps import *
from _ring_punches import *
# external
import ost
from ost import mol
from ost.mol import mm
import os
def BuildSidechains(mhandle, merge_distance=4, scorer=None, fragment_db=None,
structure_db=None, torsion_sampler=None):
'''Build sidechains for model.
This is a wrapper for :func:`promod3.sidechain.Reconstruct`, followed
by a check for ring punches. If ring punches are found it introduces gaps
for the residues with punched rings and tries to fill them with
:func:`FillLoopsByDatabase` with *ring_punch_detection=2*.
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
:param merge_distance: Used as parameter for :func:`MergeGapsByDistance`
if ring punches are found.
:type merge_distance: :class:`int`
:param scorer: Used as parameter for :func:`FillLoopsByDatabase`
if ring punches are found. A default one is created
if None.
:type scorer: :class:`~promod3.loop.BackboneLoopScorer`
:param fragment_db: Used as parameter for :func:`FillLoopsByDatabase`
if ring punches are found. A default one is loaded
if None.
:type fragment_db: :class:`~promod3.loop.FragDB`
:param structure_db: Used as parameter for :func:`FillLoopsByDatabase`
if ring punches are found. A default one is loaded
if None.
:type structure_db: :class:`~promod3.loop.StructureDB`
:param torsion_sampler: Used as parameter for :func:`FillLoopsByDatabase`
if ring punches are found. A default one is loaded
if None.
:type torsion_sampler: :class:`~promod3.loop.TorsionSampler`
'''
ost.LogInfo("Rebuilding sidechains.")
sidechain.Reconstruct(mhandle.model, keep_sidechains=True)
# check for ring punches
rings = GetRings(mhandle.model)
ring_punches = GetRingPunches(rings, mhandle.model)
# try to fix them
if len(ring_punches) > 0:
ost.LogInfo("Trying to fix %d ring punch(es)." % len(ring_punches))
# backup old gaps
old_gaps = [g.Copy() for g in mhandle.gaps]
# new gaps for mhandle
mhandle.gaps = StructuralGapList()
for res in ring_punches:
mygap = StructuralGap(res.prev, res.next, res.one_letter_code)
mhandle.gaps.append(mygap)
# load stuff if needed
if scorer is None:
scorer = SetupBackboneScorer(mhandle)
if fragment_db is None:
fragment_db = loop.LoadFragDB()
if structure_db is None:
structure_db = loop.LoadStructureDB()
if torsion_sampler is None:
torsion_sampler = loop.LoadTorsionSamplerCoil()
# fix it
MergeGapsByDistance(mhandle, merge_distance)
FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
torsion_sampler, ring_punch_detection=2)
# re-build sidechains
sidechain.Reconstruct(mhandle.model, keep_sidechains=True)
# restore gaps
mhandle.gaps = StructuralGapList()
for g in old_gaps:
mhandle.gaps.append(g)
def MinimizeModelEnergy(mhandle, max_iterations=12, max_iter_sd=20,
max_iter_lbfgs=10):
'''Minimize energy of final model using molecular mechanics.
Uses :mod:`ost.mol.mm` to perform energy minimization.
It will iteratively (at most *max_iterations* times):
- run up to *max_iter_sd* minimization iter. of a steepest descend method
- run up to *max_iter_lbfgs* minimization iter. of a Limited-memory
Broyden-Fletcher-Goldfarb-Shanno method
- abort if no stereochemical problems found
The idea is that we don't want to minimize "too much". So, we iteratively
minimize until there are no stereochemical problems and not more.
To speed things up, this can run on multiple CPU threads by setting the
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.
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
:param max_iterations: Max. number of iterations for SD+LBFGS
:type max_iterations: :class:`int`
:param max_iter_sd: Max. number of iterations within SD method
:type max_iter_sd: :class:`int`
:param max_iter_lbfgs: Max. number of iterations within LBFGS method
:type max_iter_lbfgs: :class:`int`
'''
ost.LogInfo("Minimize energy.")
# ignore LogInfo in stereochemical problems if output up to info done
ignore_stereo_log = (ost.GetVerbosityLevel() == 3)
# setup mm simulation
settings = mm.Settings()
settings.integrator = mm.LangevinIntegrator(310, 1, 0.002)
settings.init_temperature = 0
settings.forcefield = mm.LoadCHARMMForcefield()
settings.nonbonded_method = mm.NonbondedMethod.CutoffNonPeriodic
settings.keep_ff_specific_naming = False
# switch this to "mm.Platform.Reference" for debugging
settings.platform = mm.Platform.CPU
if os.getenv('PM3_OPENMM_CPU_THREADS') is None:
settings.cpu_properties["CpuThreads"] = "1"
else:
settings.cpu_properties["CpuThreads"] = os.getenv('PM3_OPENMM_CPU_THREADS')
sim = mm.Simulation(mhandle.model,settings)
# settings to check for stereochemical problems
clashing_distances = mol.alg.DefaultClashingDistances()
bond_stereo_chemical_param = mol.alg.DefaultBondStereoChemicalParams()
angle_stereo_chemical_param = mol.alg.DefaultAngleStereoChemicalParams()
for i in range(max_iterations):
# update atoms
ost.LogInfo("Perform energy minimization (iteration %d)" % (i+1))
sim.ApplySD(tolerance = 1.0, max_iterations = max_iter_sd)
sim.ApplyLBFGS(tolerance = 1.0, max_iterations = max_iter_lbfgs)
sim.UpdatePositions()
# check for stereochemical problems
if ignore_stereo_log:
ost.PushVerbosityLevel(2)
temp_ent = sim.GetEntity()
temp_ent = temp_ent.Select("aname!=OXT")
temp_ent_clash_filtered = mol.alg.FilterClashes(\
temp_ent, clashing_distances)[0]
# note: 10,10 parameters below are hard coded bond-/angle-tolerances
temp_ent_stereo_checked = mol.alg.CheckStereoChemistry(\
temp_ent_clash_filtered,
bond_stereo_chemical_param,
angle_stereo_chemical_param,
10, 10)[0]
if ignore_stereo_log:
ost.PopVerbosityLevel()
# checks above would remove bad atoms
if len(temp_ent_stereo_checked.Select("ele!=H").atoms) \
== len(temp_ent.Select("ele!=H").atoms):
break
else:
ost.LogInfo("Stereo-chemical problems found, need more energy"
+ " minimization (iteration %d)" % (i+1))
# update model
simulation_ent = sim.GetEntity()
mhandle.model = mol.CreateEntityFromView(\
simulation_ent.Select("peptide=true and ele!=H"),
True)
def CheckFinalModel(mhandle):
'''Performs samity checks on final models and reports problems.
:param mhandle: Modelling handle for which to perform checks.
:type mhandle: :class:`ModellingHandle`
'''
# report incomplete models
if len(mhandle.gaps) > 0:
ost.LogWarning("Failed to close %d gap(s). Returning incomplete model!" % \
len(mhandle.gaps))
else:
# check sequences
for chain, seq in zip(mhandle.model.chains, mhandle.seqres):
a = chain.residues[0].GetNumber().GetNum()
b = chain.residues[-1].GetNumber().GetNum()
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))
# 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!")
# report stereo-chemical problems
ost.PushVerbosityLevel(0)
clashing_distances = mol.alg.DefaultClashingDistances()
bond_stereo_chemical_param = mol.alg.DefaultBondStereoChemicalParams()
angle_stereo_chemical_param = mol.alg.DefaultAngleStereoChemicalParams()
# extract problems
model_src = mhandle.model.Select("aname!=OXT")
clash_info = mol.alg.FilterClashes(model_src, clashing_distances)[1]
# note: 10,10 parameters below are hard coded bond-/angle-tolerances
stereo_info = mol.alg.CheckStereoChemistry(model_src,
bond_stereo_chemical_param,
angle_stereo_chemical_param,
10, 10)[1]
ost.PopVerbosityLevel()
# set bool props in model-residues
atoms = [e.GetFirstAtom() for e in clash_info.GetClashList()]\
+ [e.GetFirstAtom() for e in stereo_info.GetBondViolationList()]\
+ [e.GetSecondAtom() for e in stereo_info.GetAngleViolationList()]
for atomui in atoms:
res = model_src.FindResidue(atomui.GetChainName(), atomui.GetResNum())
res.SetBoolProp("stereo_chemical_problem", True)
if atomui.GetAtomName() in ["CA", "N", "O", "C"]:
res.SetBoolProp("stereo_chemical_problem_backbone", True)
# report bad residues
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))
elif res.HasProp("stereo_chemical_problem") and\
res.GetBoolProp("stereo_chemical_problem"):
ost.LogInfo("Stereo-chemical problem in sidechain " + \
"of residue " + str(res))
def BuildFromRawModel(mhandle):
'''Build a model starting with a raw model (see :func:`BuildRawModel`).
This function implements a recommended pipeline to generate complete models
from a raw model. The steps are shown in detail in the code example
:ref:`above <modelling_steps_example>`. If you wish to use your own
pipeline, you can use that code as a starting point for your own custom
modelling pipeline. For reproducibility, we recommend that you keep copies
of custom pipelines.
If the function fails to close all gaps, it will produce a warning and
return an incomplete model.
:param mhandle: The prepared template coordinates loaded with the input
alignment.
:type mhandle: :class:`ModellingHandle`
:return: Delivers the model as an |ost_s| entity.
:rtype: :class:`Entity <ost.mol.EntityHandle>`
'''
ost.LogInfo("Starting modelling based on a raw model.")
# a bit of setup
fragment_db = loop.LoadFragDB()
structure_db = loop.LoadStructureDB()
torsion_sampler = loop.LoadTorsionSamplerCoil()
merge_distance = 4
scorer = SetupBackboneScorer(mhandle)
# remove terminal gaps and close small deletions
RemoveTerminalGaps(mhandle)
CloseSmallDeletions(mhandle, scorer)
# iteratively merge gaps of distance i and fill loops by database
for distance in range(merge_distance):
MergeGapsByDistance(mhandle, distance)
FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
torsion_sampler, min_loops_required=-1,
max_res_extension=6)
FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
torsion_sampler, min_loops_required=-1)
FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
torsion_sampler)
# close remaining gaps by Monte Carlo
FillLoopsByMonteCarlo(mhandle, scorer, torsion_sampler)
CloseLargeDeletions(mhandle, scorer, structure_db)
# build sidechains
BuildSidechains(mhandle, merge_distance, scorer, fragment_db,
structure_db, torsion_sampler)
# minimize energy of final model using molecular mechanics
MinimizeModelEnergy(mhandle)
# sanity checks
CheckFinalModel(mhandle)
# done
return mhandle.model
# these methods will be exported into module
__all__ = ('BuildFromRawModel', 'BuildSidechains', 'MinimizeModelEnergy',
'CheckFinalModel')