Select Git revision
_pipeline.py
_pipeline.py 8.76 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.
'''
import ost
from ost import mol,conop
from ost.mol import mm
from promod3 import loop
from promod3 import sidechain
from . import _modelling
from . import _closegaps
def BuildSidechains(mhandle):
'''Build sidechains for model.
This is esentially a wrapper for :func:`promod3.sidechain.Reconstruct`.
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
'''
# reconstruct the sidechains
sidechain.Reconstruct(mhandle.model, keep_sidechains=True)
def MinimizeModelEnergy(mhandle, max_iterations=3, max_iter_sd=30,
max_iter_lbfgs=20):
'''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 iterations of a steepest descend method
- run up to *max_iter_lbfgs* minimization iterations 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.
: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`
'''
# 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.platform = mm.Platform.Reference
settings.nonbonded_method = mm.NonbondedMethod.CutoffNonPeriodic
settings.keep_ff_specific_naming = False
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 BuildFromRawModel(mhandle, fragment_db=None, structure_db=None,
torsion_sampler=None, merge_distance=4):
'''Build a model starting with a raw model (see :func:`BuildRawModel`).
This function will:
#. close small deletions and remove terminal gaps
(see :func:`CloseSmallDeletions` and :func:`RemoveTerminalGaps`)
#. iteratively: merge gaps of distance i and fill loops by database
(i from 0 to *merge_distance*-1)
(see :func:`MergeGapsByDistance` and :func:`FillLoopsByDatabase`)
#. close remaining gaps by Monte Carlo
(see :func:`FillLoopsByMonteCarlo`)
#. build sidechains
(see :func:`BuildSidechains`)
#. minimize energy for the final model, which is then returned
(see :func:`MinimizeModelEnergy`)
If 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:`~promod3.modelling.ModellingHandle`
:param fragment_db: A fragment database coupled to the *structure_db*.
Loads the default one shipped with |project| if omitted.
:type fragment_db: :class:`~promod3.loop.FragDB`
:param structure_db: The structural database. Loads the default one shipped
with |project| if omitted.
:type structure_db: :class:`~promod3.loop.StructureDB`
:param torsion_sampler: A sampler for torsion angles. Loads the default one
shipped with |project| if omitted.
:type torsion_sampler: :class:`~promod3.loop.TorsionSampler`
:param merge_distance: Maximal distance up to which 2 gaps are merged. The
gaps need to be immediate neighbours, distance is
counted by residues, including stems. So **A-A-A**
have a distance of 0 and therefore merging is only
applied until *merge_distance* - 1.
:type merge_distance: :class:`int`
: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
if fragment_db == None:
fragment_db = loop.LoadFragDB()
if structure_db == None:
structure_db = loop.LoadStructureDB()
if torsion_sampler == None:
torsion_sampler = loop.LoadTorsionSamplerCoil()
scorer = _closegaps.SetupBackboneScorer(mhandle)
# step 1a: close small deletions in the raw model
ost.LogInfo("Trying to close small deletions (no. of gaps: %d)." % \
len(mhandle.gaps))
_closegaps.CloseSmallDeletions(mhandle, scorer)
# step 1b: remove terminal gaps since we cannot deal with them below
num_gaps_removed = _modelling.RemoveTerminalGaps(mhandle)
ost.LogInfo("Removed %d terminal gaps." % num_gaps_removed)
# step 2: simple handling of further gaps: merge, then fill from db
for distance in range(merge_distance):
# step 2a: Combine gaps living close to each other
ost.LogInfo("Trying to merge %d gap(s) with distance %d." % \
(len(mhandle.gaps), distance))
_closegaps.MergeGapsByDistance(mhandle, distance)
# step 2b: fit loops into the model
ost.LogInfo("Trying to fill %d gap(s) by database." % len(mhandle.gaps))
_closegaps.FillLoopsByDatabase(mhandle, scorer, fragment_db,
structure_db, torsion_sampler)
# check if we can abort
ost.LogInfo("%d gap(s) left after database search." % len(mhandle.gaps))
if len(mhandle.gaps) == 0:
break
# step 3: close remaining gaps by Monte Carlo
if len(mhandle.gaps) > 0:
ost.LogInfo("Trying to fill %d gap(s) by Monte Carlo." % len(mhandle.gaps))
_closegaps.FillLoopsByMonteCarlo(mhandle, scorer, torsion_sampler)
# check if we succeeded
if len(mhandle.gaps) > 0:
ost.LogWarning("Failed to close %d gap(s). Returning incomplete model!" % \
len(mhandle.gaps))
# step 4: build sidechains
ost.LogInfo("Rebuilding sidechains.")
BuildSidechains(mhandle)
# step 5: minimize energy of final model using molecular mechanics
ost.LogInfo("Minimize energy.")
MinimizeModelEnergy(mhandle)
# done
return mhandle.model
# these methods will be exported into module
__all__ = ('BuildFromRawModel', 'BuildSidechains', 'MinimizeModelEnergy',)