From 34a865e5fe5ff949c999e42221f3fa4b13089bc8 Mon Sep 17 00:00:00 2001
From: Gerardo Tauriello <gerardo.tauriello@unibas.ch>
Date: Fri, 11 Dec 2015 18:50:22 +0100
Subject: [PATCH] added fill gaps by Monte Carlo (with docs and tests)
---
core/tests/test_helper.py | 3 +-
loop/doc/loop_candidate.rst | 201 ++++++------
modelling/doc/index.rst | 2 +
modelling/pymod/_closegaps.py | 489 ++++++++++++++++++++++-------
modelling/src/model.cc | 2 +-
modelling/tests/test_close_gaps.py | 212 +++++++++++--
modelling/tests/test_modelling.py | 10 +
pipeline/pymod/run_engine.py | 20 +-
8 files changed, 672 insertions(+), 267 deletions(-)
diff --git a/core/tests/test_helper.py b/core/tests/test_helper.py
index bf570fb3..9005334a 100644
--- a/core/tests/test_helper.py
+++ b/core/tests/test_helper.py
@@ -9,7 +9,8 @@ class _FetchLog(ost.LogSink):
self.messages = dict()
def LogMessage(self, message, severity):
- levels = ['ERROR', 'WARNING', 'INFO', 'VERBOSE', 'DEBUG', 'TRACE']
+ levels = ['ERROR', 'WARNING', 'SCRIPT', 'INFO', 'VERBOSE', 'DEBUG',
+ 'TRACE']
level = levels[severity]
if not level in self.messages.keys():
self.messages[level] = list()
diff --git a/loop/doc/loop_candidate.rst b/loop/doc/loop_candidate.rst
index 8ddf9789..e6e889f1 100644
--- a/loop/doc/loop_candidate.rst
+++ b/loop/doc/loop_candidate.rst
@@ -94,29 +94,38 @@ Loop Candidates
.. class:: LoopCandidates(n_stem,c_stem,seq)
- The *LoopCandidates* is the basic object used for loop modelling. It contains the positions of the two anchor residues, the *n_stem* and *c_stem* residues, and the sequence of the residues to be modelled in between these anchor residues.
- It also contains a list of :class:`LoopCandidate`, which are possible conformations of the backbone between the stem residues.
+ The *LoopCandidates* is the basic object used for loop modelling.
+ It contains the positions of the two anchor residues, the *n_stem* and *c_stem*
+ residues, and the sequence of the residues to be modelled in between these
+ anchor residues.
+ It also contains a list of :class:`LoopCandidate`, which are possible
+ conformations of the backbone between the stem residues.
:param n_stem: The residue at the N-terminal end of the loop
:param c_stem: The residue at the C-terminal end of the loop
- :param seq: The sequence of residues to be added between the *n_stem* and *c_stem*
+ :param seq: The sequence of residues to be added between the
+ *n_stem* and *c_stem*
:type n_stem: :class:`ost.mol.ResidueHandle`
:type c_stem: :class:`ost.mol.ResidueHandle`
:type sequence: :class:`str`
- .. method:: FillFromDatabase(n_stem, c_stem, seq, frag_db, structural_db, [extended_search=False])
+ .. staticmethod:: FillFromDatabase(n_stem, c_stem, seq, frag_db, \
+ structural_db, extended_search=False)
Searches for loop candidates matching the length (number of residues in *seq*) and geometry (of *n_stem* and *c_stem*) of the loop to be modelled in a fragment database.
:param n_stem: The residue at the N-terminal end of the loop
:param c_stem: The residue at the C-terminal end of the loop
- :param seq: The sequence of residues to be added between the *n_stem* and *c_stem*
+ :param seq: The sequence of residues to be added between the
+ *n_stem* and *c_stem*
:param frag_db: The fragment database
:param structural_db: The according structural database
- :param extended_search: Wether search should be extended to fragments matching the geometry of the *n_stem* and *c_stem* a bit less precisely.
+ :param extended_search: Wether search should be extended to fragments
+ matching the geometry of the *n_stem* and *c_stem*
+ a bit less precisely.
:type n_stem: :class:`ost.mol.ResidueHandle`
:type c_stem: :class:`ost.mol.ResidueHandle`
@@ -128,134 +137,108 @@ Loop Candidates
:returns: A list of loop candidates
:rtype: :class:`LoopCandidates`
- .. method:: FillFromMonteCarloSampler(n_stem, c_stem, initial_bb, seq, num_loops, steps, sampler, closer, scorer, cooler)
+ .. staticmethod::
+ FillFromMonteCarloSampler(n_stem, c_stem, seq, num_loops, \
+ steps, sampler, closer, scorer, cooler)
+ FillFromMonteCarloSampler(n_stem, c_stem, initial_bb, seq, num_loops, \
+ steps, sampler, closer, scorer, cooler)
Uses Monte Carlo simulated annealing to sample the loop to be modelled.
- :param n_stem: The residue at the N-terminal end of the loop
- :param c_stem: The residue at the C-terminal end of the loop
- :param initial_bb: Initial configuration used for the sampling
- :param seq: The sequence of residues to be added between the *n_stem* and *c_stem*
- :param num_loops: Number of loop candidates to return
- :param steps: Number of MC steps to perform for each loop candidate generated
- :param sampler: Used to generate a new configuration at each MC step
- :param closer: Used to close the loop (make it match the *n_stem* and *c_stem* geometry) after each MC step
- :param scorer: Used to score the generated configurations at each MC step
- :param cooler: Controls the temperature profile of the simulated annealing
-
- :type n_stem: :class:`ost.mol.ResidueHandle`
- :type c_stem: :class:`ost.mol.ResidueHandle`
- :type seq: :class:`str`
- :type initial_bb: :class:`BackboneList`
- :type num_loops: :class:`int`
- :type steps: :class:`int`
- :type sampler: :ref:`mc-sampler-object`
- :type closer: :ref:`mc-closer-object`
- :type scorer: :ref:`mc-scorer-object`
- :type cooler: :ref:`mc-cooler-object`
-
- :returns: A list of loop candidates
- :rtype: :class:`LoopCandidates`
-
- .. method:: FillFromMonteCarloSampler(n_stem, c_stem, seq, num_loops, steps, sampler, closer, scorer, cooler)
-
- Uses Monte Carlo simulated annealing to sample the loop to be modelled.
-
- :param n_stem: The residue at the N-terminal end of the loop
- :param c_stem: The residue at the C-terminal end of the loop
- :param seq: The sequence of residues to be added between the *n_stem* and *c_stem*
- :param num_loops: Number of loop candidates to return
- :param steps: Number of MC steps to perform for each loop candidate generated
- :param sampler: Used to generate a new configuration at each MC step
- :param closer: Used to close the loop (make it match the *n_stem* and *c_stem* geometry) after each MC step
- :param scorer: Used to score the generated configurations at each MC step
- :param cooler: Controls the temperature profile of the simulated annealing
-
- :type n_stem: :class:`ost.mol.ResidueHandle`
- :type c_stem: :class:`ost.mol.ResidueHandle`
- :type seq: :class:`str`
- :type num_loops: :class:`int`
- :type steps: :class:`int`
- :type sampler: :ref:`mc-sampler-object`
- :type closer: :ref:`mc-closer-object`
- :type scorer: :ref:`mc-scorer-object`
- :type cooler: :ref:`mc-cooler-object`
-
- :returns: A list of loop candidates
- :rtype: :class:`LoopCandidates`
+ :param n_stem: The residue at the N-terminal end of the loop
+ :param c_stem: The residue at the C-terminal end of the loop
+ :param initial_bb: Initial configuration used for the sampling (optional)
+ :param seq: The sequence of residues to be added between the
+ *n_stem* and *c_stem*
+ :param num_loops: Number of loop candidates to return
+ :param steps: Number of MC steps to perform for each loop candidate
+ generated
+ :param sampler: Used to generate a new configuration at each MC step
+ :param closer: Used to close the loop (make it match the *n_stem*
+ and *c_stem* geometry) after each MC step
+ :param scorer: Used to score the generated configurations at each
+ MC step
+ :param cooler: Controls the temperature profile of the simulated annealing
+
+ :type n_stem: :class:`ost.mol.ResidueHandle`
+ :type c_stem: :class:`ost.mol.ResidueHandle`
+ :type initial_bb: :class:`BackboneList`
+ :type seq: :class:`str`
+ :type num_loops: :class:`int`
+ :type steps: :class:`int`
+ :type sampler: :ref:`mc-sampler-object`
+ :type closer: :ref:`mc-closer-object`
+ :type scorer: :ref:`mc-scorer-object`
+ :type cooler: :ref:`mc-cooler-object`
+
+ :returns: A list of loop candidates
+ :rtype: :class:`LoopCandidates`
+
+ :raises: A :exc:`RuntimeError`, if *initial_bb* is not given and the Monte
+ Carlo sampler failed to initialize (i.e. stems are too far apart)
+ :raises: A :exc:`RuntimeError`, if *initial_bb* is given and it is
+ inconsistent with *seq*
.. method:: ClusterCandidates()
- Clusters the loop candidates according to their pairwise CA-RMSD using a FLAME (fuzzy clustering by local approximation fo memberships) clustering algorithm.
-
- :returns: A list of :class:`LoopCandidates`. Each element in the list corresponds to the candidates in one cluster.
-
- .. method:: ApplyCCD(torsion_sampler, [max_iterations=1000, rmsd_cutoff=0.1, keep_non_converged=false, random_seed=0, step_size=1.0])
-
- Closes all :class:`LoopCandidate` in :class:`LoopCandidates` (i.e. modifies the :class:`LoopCandidate` so that its stem residues match those of the :class:`LoopCandidates`, which are the stem residues of the loop being modelled), using the CCD algorithm. CCD (cyclic coordinate descent) is an iterative minimization algorithm. At each step of the closing the *torsion_sampler* is used to calculate the probability of the proposed move and accepted or not depending on a metropolis criterium.
-
- :param torsion_sampler: The torsion sampler
- :param max_iterations: Maximum number of iteration
- :param rmsd_cutoff: Cutoff in stem residue RMSD used to determine convergence
- :param keep_non_converged: Wether to keep loop candidates for which the closing did not converged
+ Clusters the loop candidates according to their pairwise CA-RMSD using a
+ FLAME (fuzzy clustering by local approximation fo memberships) clustering
+ algorithm.
+
+ :returns: A list of :class:`LoopCandidates`. Each element in the list
+ corresponds to the candidates in one cluster.
+
+ .. method:: ApplyCCD(max_iterations=1000, rmsd_cutoff=0.1, \
+ keep_non_converged=false, random_seed=0, step_size=1.0)
+ ApplyCCD(torsion_sampler, max_iterations=1000, rmsd_cutoff=0.1, \
+ keep_non_converged=false, random_seed=0, step_size=1.0)
+ ApplyCCD(torsion_samplers, max_iterations=1000, rmsd_cutoff=0.1, \
+ keep_non_converged=false, random_seed=0, step_size=1.0)
+
+ Closes all :class:`LoopCandidate` in :class:`LoopCandidates` using the CCD
+ algorithm (i.e. modifies the :class:`LoopCandidate` so that its stem residues
+ match those of the :class:`LoopCandidates`, which are the stem residues of the
+ loop being modelled).
+ CCD (cyclic coordinate descent, see :class:`~promod3.loop.CCD`) is an
+ iterative minimization algorithm.
+
+ If *torsion_sampler* or *torsion_samplers* is given, it is used at each step
+ of the closing to calculate the probability of the proposed move, which is
+ then accepted (or not) depending on a metropolis criterium.
+
+ :param torsion_sampler: the torsion sampler
+ :param torsion_samplers: a list containing one torsion sampler for each
+ residue in the loop.
+ :param max_iterations: maximum number of iteration
+ :param rmsd_cutoff: cutoff in stem residue RMSD used to determine
+ convergence
+ :param keep_non_converged: whether to keep loop candidates for which the
+ closing did not converge. If this is False, it
+ erases loop candidates which didn't converge!
:param random_seed: seed for random number generator
:param step_size: size of the steps taken during loop closing.
:type torsion_sampler: :class:`TorsionSampler`
+ :type torsion_samplers: :class:`list` of :class:`TorsionSampler`
:type max_iterations: :class:`int`
:type rmsd_cutoff: :class:`float`
:type keep_non_converged: :class:`bool`
:type random_seed: :class:`int`
- :type step_size: :class:`float`
-
- .. method:: ApplyCCD(torsion_samplers, [max_iterations=1000, rmsd_cutoff=0.1, keep_non_converged=false, random_seed=0, step_size=1.0])
-
- Closes all :class:`LoopCandidate` in :class:`LoopCandidates` (i.e. modifies the :class:`LoopCandidate` so that its stem residues match those of the :class:`LoopCandidates`, which are the stem residues of the loop being modelled) using the CCD algorithm. CCD (cyclic coordinate descent) is an iterative minimization algorithm. At each step of the closing the *torsion_sampler* is used to calculate the probability of the proposed move and accepted or not depending on a metropolis criterium.
-
- :param torsion_sampler: A list containing one torsion sampler for each residue in the loop.
- :param max_iterations: Maximum number of iteration
- :param rmsd_cutoff: Cutoff in stem residue RMSD used to determine convergence
- :param keep_non_converged: Wether to keep loop candidates for which the closing did not converged
- :param random_seed: seed for random number generator
- :param step_size: size of the steps taken during loop closing.
-
- :type torsion_samplers: :class:`TorsionSampler`
- :type max_iterations: :class:`int`
- :type rmsd_cutoff: :class:`float`
- :type keep_non_converged: :class:`bool`
- :type random_seed: :class:`int`
- :type step_size: :class:`float`
-
-
- .. method:: ApplyCCD([max_iterations=1000, rmsd_cutoff=0.1, keep_non_converged=false, random_seed=0, step_size=1.0])
-
- Closes all :class:`LoopCandidate` in :class:`LoopCandidates` (i.e. modifies the :class:`LoopCandidate` so that its stem residues match those of the :class:`LoopCandidates`, which are the stem residues of the loop being modelled), using the CCD algorithm. CCD (cyclic coordinate descent) is an iterative minimization algorithm.
-
- :param max_iterations: Maximum number of iteration
- :param rmsd_cutoff: Cutoff in stem residue RMSD used to determine convergence
- :param keep_non_converged: Wether to keep loop candidates for which the closing did not converged
- :param random_seed: seed for random number generator
- :param step_size: size of the steps taken during loop closing.
-
- :type max_iterations: :class:`int`
- :type rmsd_cutoff: :class:`float`
- :type keep_non_converged: :class:`bool`
- :type random_seed: :class:`int`
- :type step_size: :class:`float`
+ :type step_size: :class:`float`
.. method:: ApplyKIC(pivot_one,pivot_two,pivot_three)
Closes all :class:`LoopCandidate` in :class:`LoopCandidates` (i.e. modifies the :class:`LoopCandidate` so that its stem residues match those of the :class:`LoopCandidates`, which are the stem residues of the loop being modelled), using the KIC algorithm. This algorithm finds analytical solutions for closing the loop by modifying the torsion angles of just three pivot residues. Due to the underlying mathematical formalism in KIC, up to 16 solutions can be found for every candidate. This leads to an increase in number of loops.
- :param pivot_one: Maximum number of iteration
- :param pivot_two: Cutoff in stem residue RMSD used
+ :param pivot_one: Maximum number of iteration
+ :param pivot_two: Cutoff in stem residue RMSD used
:param pivot_three: Wether to keep loop candidates for
:type max_iterations: :class:`int`
:type rmsd_cutoff: :class:`float`
:type keep_non_converged: :class:`bool`
:type random_seed: :class:`int`
- :type step_size: :class:`float`
+ :type step_size: :class:`float`
.. method:: ToEntity(index)
diff --git a/modelling/doc/index.rst b/modelling/doc/index.rst
index 390a66d7..eceed77d 100644
--- a/modelling/doc/index.rst
+++ b/modelling/doc/index.rst
@@ -93,6 +93,8 @@ Modelling API
.. autofunction:: FillLoopsByDatabase
+.. autofunction:: FillLoopsByMonteCarlo
+
.. function:: ClearGaps(mhandle, gap)
Removes all gaps from mhandle which are fully enclosed by given gap.
diff --git a/modelling/pymod/_closegaps.py b/modelling/pymod/_closegaps.py
index f32d39f9..b1c48dbc 100644
--- a/modelling/pymod/_closegaps.py
+++ b/modelling/pymod/_closegaps.py
@@ -1,5 +1,5 @@
-'''Functions added as high-level functionality to modelling module in the
-__init__.py file. To be used directly by passing a ModellingHandle instance
+'''Functions added as high-level functionality to modelling module in the
+__init__.py file. To be used directly by passing a ModellingHandle instance
as argument.
'''
@@ -7,6 +7,7 @@ import ost
#pylint: disable=no-name-in-module
from . import _modelling as modelling
from promod3 import loop
+import sys
def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
e_thresh=200):
@@ -21,33 +22,33 @@ def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
Closed gaps are removed from :attr:`mhandle.gaps`.
.. testcode:: closesmalldel
- :hide:
+ :hide:
- from promod3 import modelling
- from promod3 import loop
+ from promod3 import modelling
+ from promod3 import loop
- tpl = ost.io.LoadPDB('../tests/modelling/data/gly.pdb')
- aln = ost.seq.CreateAlignment(ost.seq.CreateSequence('trg', 'GGG-GGG'),
- ost.seq.CreateSequence('tpl', 'GGGAGGG'))
- aln.AttachView(1, tpl.CreateFullView())
- mhandle = modelling.BuildRawModel(aln)
- assert len(mhandle.gaps) == 1
- scorer = loop.SetupBackboneScorer(mhandle)
- modelling.CloseSmallDeletions(mhandle, scorer)
- assert len(mhandle.gaps) == 0
+ tpl = ost.io.LoadPDB('../tests/modelling/data/gly.pdb')
+ aln = ost.seq.CreateAlignment(ost.seq.CreateSequence('trg', 'GGG-GGG'),
+ ost.seq.CreateSequence('tpl', 'GGGAGGG'))
+ aln.AttachView(1, tpl.CreateFullView())
+ mhandle = modelling.BuildRawModel(aln)
+ assert len(mhandle.gaps) == 1
+ scorer = loop.SetupBackboneScorer(mhandle)
+ modelling.CloseSmallDeletions(mhandle, scorer)
+ assert len(mhandle.gaps) == 0
.. doctest:: closesmalldel
- import ost
- from promod3 import modelling
- from promod3 import loop
+ import ost
+ from promod3 import modelling
+ from promod3 import loop
- tpl = ost.io.LoadPDB('gly.pdb')
- aln = ost.io.LoadAlignment('seq.fasta')
- aln.AttachView(1, tpl.CreateFullView())
- mhandle = modelling.BuildRawModel(aln)
- scorer = loop.SetupBackboneScorer(mhandle)
- modelling.CloseSmallDeletions(mhandle, scorer)
+ tpl = ost.io.LoadPDB('gly.pdb')
+ aln = ost.io.LoadAlignment('seq.fasta')
+ aln.AttachView(1, tpl.CreateFullView())
+ mhandle = modelling.BuildRawModel(aln)
+ scorer = loop.SetupBackboneScorer(mhandle)
+ modelling.CloseSmallDeletions(mhandle, scorer)
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
@@ -68,8 +69,7 @@ def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
:param e_thresh: Potential energy should be lower than this.
:type e_thresh: :class:`float`
- :return: If gaps are deleted, the scorer needs to be updated and is
- returned.
+ :return: If gaps are deleted, the `scorer` is updated and returned.
:rtype: :class:`~promod3.loop.BackboneLoopScorer`
'''
current_gap_index = 0
@@ -84,7 +84,7 @@ def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
# A deletion is a gap of size 0 in the template, this means that to
# transform the template sequence into the target sequence, aa's vanish,
- # so the target sequence has gap caracters, the template not.
+ # so the target sequence has gap characters, the template not.
# If we are not looking at a deletion, do nothing.
if len(mhandle.gaps[current_gap_index].seq) == 0:
# work on a copy of the gap, if not closed in the end, no harm done
@@ -106,9 +106,7 @@ def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
current_resnum = current_gap.before.GetNumber()
after_resnum = current_gap.after.GetNumber()
found_residues = True
- while True:
- if current_resnum > after_resnum:
- break
+ while current_resnum <= after_resnum:
res = current_chain.FindResidue(current_resnum)
if not res.IsValid():
found_residues = False
@@ -133,7 +131,7 @@ def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0,
if score < clash_thresh and \
potential_e < e_thresh and \
scorer.TransOmegaTorsions(bb_list):
- ost.LogVerbose("Closed: %s by relaxing %s" % \
+ ost.LogInfo("Closed: %s by relaxing %s" % \
(mhandle.gaps[current_gap_index], current_gap))
chain = current_gap.before.GetChain()
bb_list.InsertInto(chain, current_gap.before.GetNumber().GetNum(),
@@ -162,39 +160,37 @@ def MergeGapsByDistance(mhandle, distance):
removed. Stem residues count to the gap, so **A-A-A** has a distance of 0.
.. testcode:: mergegapsbydist
- :hide:
+ :hide:
- from promod3 import modelling
+ from promod3 import modelling
- tpl = ost.io.LoadPDB('../tests/modelling/data/1mcg.pdb')
- aln = ost.seq.CreateAlignment(ost.seq.CreateSequence('trg',
- 'DDFAGDTKNLGHN'),
- ost.seq.CreateSequence('tpl',
- 'NN----A----LF'))
- aln.AttachView(1, tpl.CreateFullView())
- mhandle = modelling.BuildRawModel(aln)
- assert len(mhandle.gaps) == 2
- modelling.MergeGapsByDistance(mhandle, 0)
- assert len(mhandle.gaps) == 1
+ tpl = ost.io.LoadPDB('../tests/modelling/data/1mcg.pdb')
+ aln = ost.seq.CreateAlignment(ost.seq.CreateSequence('trg',
+ 'DDFAGDTKNLGHN'),
+ ost.seq.CreateSequence('tpl',
+ 'NN----A----LF'))
+ aln.AttachView(1, tpl.CreateFullView())
+ mhandle = modelling.BuildRawModel(aln)
+ assert len(mhandle.gaps) == 2
+ modelling.MergeGapsByDistance(mhandle, 0)
+ assert len(mhandle.gaps) == 1
.. doctest:: mergegapsbydist
- import ost
- from promod3 import modelling
+ import ost
+ from promod3 import modelling
- tpl = ost.io.LoadPDB('1mcg.pdb')
- aln = ost.io.LoadAlignment('1mcg_aln.fasta')
- aln.AttachView(1, tpl.CreateFullView())
- mhandle = modelling.BuildRawModel(aln)
- modelling.MergeGapsByDistance(mhandle, 0)
+ tpl = ost.io.LoadPDB('1mcg.pdb')
+ aln = ost.io.LoadAlignment('1mcg_aln.fasta')
+ aln.AttachView(1, tpl.CreateFullView())
+ mhandle = modelling.BuildRawModel(aln)
+ modelling.MergeGapsByDistance(mhandle, 0)
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
:param distance: The max. no. of residues between two gaps up to which
merge happens.
:type distance: :class:`int`
-
- :return: Nothing.
'''
# IMPORTANT: the assumption is that ModellingHandle stores gaps
# sequentially
@@ -224,7 +220,7 @@ def MergeGapsByDistance(mhandle, distance):
if dist <= distance:
# gaps are close enough, combine! combine!
modelling.MergeGaps(mhandle, i)
- ost.LogVerbose("Merged gap %s and %s into %s" % \
+ ost.LogInfo("Merged gap %s and %s into %s" % \
(current_gap, next_gap, mhandle.gaps[i]))
try_again = True
break
@@ -237,42 +233,44 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
Usually this will extend the gaps a bit to match candidates from the
database. Do not expect a gap being filled in between its actual stem
residues.
+ This function cannot fill gaps at C- or N-terminal.
.. testcode:: fillloopsbydb
- :hide:
+ :hide:
- from promod3 import modelling
- from promod3 import loop
+ from promod3 import modelling
+ from promod3 import loop
- tpl = ost.io.LoadPDB('../tests/modelling/data/2dbs.pdb')
- aln = ost.seq.CreateAlignment(
- ost.seq.CreateSequence('trg', 'TLNGFTVPAGNTLV--LNPDKGATVTMA'),
- ost.seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSNVHIRVE'))
+ tpl = ost.io.LoadPDB('../tests/modelling/data/2dbs.pdb')
+ aln = ost.seq.CreateAlignment(
+ ost.seq.CreateSequence('trg', 'TLNGFTVPAGNTLV--LNPDKGATVTMA'),
+ ost.seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSNVHIRVE'))
- aln.AttachView(1, tpl.CreateFullView())
- mhandle = modelling.BuildRawModel(aln)
- assert len(mhandle.gaps) == 1
- scorer = loop.SetupBackboneScorer(mhandle)
- modelling.FillLoopsByDatabase(mhandle, scorer, loop.LoadFragDB(),
- loop.LoadStructureDB(),
- loop.LoadTorsionSamplerCoil())
- assert len(mhandle.gaps) == 0
+ aln.AttachView(1, tpl.CreateFullView())
+ mhandle = modelling.BuildRawModel(aln)
+ assert len(mhandle.gaps) == 1
+ scorer = loop.SetupBackboneScorer(mhandle)
+ modelling.FillLoopsByDatabase(mhandle, scorer, loop.LoadFragDB(),
+ loop.LoadStructureDB(),
+ loop.LoadTorsionSamplerCoil())
+ assert len(mhandle.gaps) == 0
.. doctest:: fillloopsbydb
- from promod3 import modelling
- from promod3 import loop
-
- tpl = io.LoadPDB('2dbs.pdb')
- aln = ost.io.LoadAlignment('2dbs.fasta')
- aln.AttachView(1, tpl.CreateFullView())
-
- mhandle = modelling.BuildRawModel(aln)
- scorer = loop.SetupBackboneScorer(mhandle)
- scorer = modelling.FillLoopsByDatabase(mhandle, scorer, loop.LoadFragDB(),
- loop.LoadStructureDB(),
- loop.LoadTorsionSamplerCoil())
+ from promod3 import modelling
+ from promod3 import loop
+
+ tpl = io.LoadPDB('2dbs.pdb')
+ aln = ost.io.LoadAlignment('2dbs.fasta')
+ aln.AttachView(1, tpl.CreateFullView())
+
+ mhandle = modelling.BuildRawModel(aln)
+ scorer = loop.SetupBackboneScorer(mhandle)
+ scorer = modelling.FillLoopsByDatabase(mhandle, scorer,
+ loop.LoadFragDB(),
+ loop.LoadStructureDB(),
+ loop.LoadTorsionSamplerCoil())
:param mhandle: Modelling handle on which to apply change.
:type mhandle: :class:`ModellingHandle`
@@ -298,8 +296,7 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
in the database.
:type max_db_loop_len: :class:`int`
- :return: If gaps are deleted, the scorer needs to be updated and is
- returned.
+ :return: If gaps are deleted, the scorer is updated returned.
:rtype: :class:`~promod3.loop.BackboneLoopScorer`
'''
#pylint: disable=too-many-branches,too-many-statements
@@ -334,6 +331,7 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
score += ex*coef
return score
+ # main code
svm_scaler = _Scaler()
svm_predictor = _Predictor()
# point to the current gap
@@ -343,77 +341,92 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
# an updated list. gap_idx is only increased when necessary, e.g. current
# gap could not be removed.
while gap_idx < len(mhandle.gaps):
- # count how many loops have been found/ searched so far, needed as a
- # stop criterion at some point to start investigating found loops.
- found_loops = 0
+
+ # keep copy of original gap
+ gap_orig = mhandle.gaps[gap_idx].Copy()
+ # get chain for gap
+ actual_chain_idx = gap_orig.GetChainIndex()
+ actual_chain = mhandle.model.chains[actual_chain_idx]
# terminal gaps are no loops, database needs 2 stems to work on
- if mhandle.gaps[gap_idx].IsNTerminal() or mhandle.gaps[gap_idx].IsCTerminal():
+ if gap_orig.IsTerminal():
gap_idx += 1
continue
- # get info on current gap
- actual_gap = mhandle.gaps[gap_idx].Copy()
- actual_length = actual_gap.length
- actual_extender = modelling.GapExtender(actual_gap)
- min_before_resnum = 100000000000
- max_after_resnum = -100000000000
+ ##################################
+ # find loop candidates
+ ##################################
+ # min/max stem res. numbers
+ min_before_resnum = sys.maxint
+ max_after_resnum = -sys.maxint-1
+ # list of LoopCandidates to fill gap
actual_candidates = list()
+ # list of StructuralGap corresponding to actual_candidates
actual_extended_gaps = list()
- actual_chain_idx = actual_gap.GetChainIndex()
- actual_chain = mhandle.model.chains[actual_chain_idx]
-
- # Find loop candidates.
- while actual_length <= max_db_loop_len:
- # extend the loop
+ # number of loops found (stop if >= max_loops_to_search)
+ found_loops = 0
+ # currently extended gap and gap-extender
+ actual_gap = gap_orig.Copy()
+ actual_extender = modelling.GapExtender(actual_gap)
+ last_gap_length = actual_gap.length
+ # iteratively extend actual_gap until we have enough loops
+ while last_gap_length <= max_db_loop_len:
+ # extend the gap
if not actual_extender.Extend():
break
+ # ensure both stems are valid
if not actual_gap.before.IsValid() or \
not actual_gap.after.IsValid():
break
-
- # check if we found 'enough' loops plus room for the newly fitted
- # stem residues
- if found_loops >= max_loops_to_search and \
- actual_gap.length > actual_length:
+ # check if we have enough and ensure that we add all candidates with same length
+ if found_loops >= max_loops_to_search and actual_gap.length > last_gap_length:
break
+
# get candidates for the current loop
candidates = loop.LoopCandidates.FillFromDatabase(
actual_gap.before, actual_gap.after, actual_gap.full_seq,
fragment_db, structure_db, False)
- # update gap info
- actual_length = actual_gap.length
+ # TODO: WTF
+ last_gap_length = actual_gap.length
min_before_resnum = min(min_before_resnum,
actual_gap.before.GetNumber().GetNum())
max_after_resnum = max(max_after_resnum,
actual_gap.after.GetNumber().GetNum())
+ # skip gaps with no loop candidates
if len(candidates) == 0:
continue
+ # try to close loops
try:
#pylint: disable=broad-except
+ candidates.ApplyCCD(torsion_sampler)
+ except Exception:
# ccd requires the residues before and after the stems
# to be valid, it could be that we extend too much...
# just neglect if that happens, further extension will
# fail anyway and causes the extension process to stop/
- candidates.ApplyCCD(torsion_sampler)
- except Exception:
continue
+ # skip if no loop was successfully closed
if len(candidates) == 0:
continue
+ # update candidate lists
actual_candidates.append(candidates)
actual_extended_gaps.append(actual_gap.Copy())
found_loops += len(candidates)
+ # check result
if len(actual_candidates) == 0:
+ ost.LogInfo("Failed at loop insterion (%s)" % str(gap_orig))
gap_idx += 1
continue
+ ##################################
# all loop candidates will be scored along the full max. extension ever
# reached in the search before, so we build an overall frame, where we
# insert the loops
+ ##################################
frame_seq = mhandle.seqres[actual_chain_idx]\
[min_before_resnum-1:max_after_resnum]
frame_backbone_list = loop.BackboneList(frame_seq)
@@ -433,29 +446,31 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
overall_after_residue,
frame_seq)
+ ##################################
# prepare loop candidates for scoring
+ ##################################
back_mapper = list()
for i, loop_candidates in enumerate(actual_candidates):
start_index = loop_candidates.GetNStem().GetNumber().GetNum()\
- min_before_resnum
for j, loop_candidate in enumerate(loop_candidates):
actual_frame_backbone_list = frame_backbone_list
- actual_frame_backbone_list.ReplaceFragment(\
- loop_candidate.bb_list,
- start_index,
- False)
+ actual_frame_backbone_list.ReplaceFragment(
+ loop_candidate.bb_list, start_index, False)
final_loop_candidates.Add(actual_frame_backbone_list)
back_mapper.append((i, j))
final_loop_candidates.AttachScorer(scorer)
final_loop_candidates.CalculateClashScores(actual_chain_idx)
- final_loop_candidates.CalculateCBetaScores()
+ final_loop_candidates.CalculateCBetaScores(actual_chain_idx)
final_loop_candidates.CalculateTorsionScores(actual_chain_idx)
final_loop_candidates.CalculateCBPackingScores(actual_chain_idx)
final_loop_candidates.CalculateHBondScores(actual_chain_idx)
+ ##################################
# compare scores
- min_score = 10000000000
+ ##################################
+ min_score = sys.maxint
optimal_candidate = -1
for i in range(len(final_loop_candidates)):
features = [final_loop_candidates[i].torsion_score,
@@ -469,7 +484,9 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
optimal_candidate = i
min_score = score
+ ##################################
# finally resolve loop
+ ##################################
if optimal_candidate != -1:
idx_a = back_mapper[optimal_candidate][0]
idx_b = back_mapper[optimal_candidate][1]
@@ -479,21 +496,261 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
bb_list,
actual_candidates[idx_a].GetNStem().GetNumber().GetNum(),
actual_chain_idx)
- ost.LogVerbose("Resolved %s by filling %s" % \
- (str(actual_gap), str(actual_extended_gaps[idx_a])))
+ ost.LogInfo("Resolved %s by filling %s" % \
+ (str(gap_orig), str(actual_extended_gaps[idx_a])))
next_gap = modelling.ClearGaps(mhandle, actual_extended_gaps[idx_a])
if next_gap == -1:
break
else:
gap_idx = next_gap
else:
- ost.LogVerbose("Failed at loop insterion (%s)" % str(actual_gap))
+ ost.LogInfo("Failed at loop insterion (%s)" % str(gap_orig))
+ gap_idx += 1
+
+ return scorer
+
+def FillLoopsByMonteCarlo(mhandle, scorer, torsion_sampler, max_loops_to_search=6,
+ max_extension=30, mc_num_loops=2, mc_steps=5000):
+ '''Try to fill up loops with Monte Carlo sampling.
+
+ This is meant as a "last-resort" approach when it is not possible to fill
+ the loops from the database with :func:`FillLoopsByDatabase`.
+
+ .. doctest:: fillloopsbymc
+
+ from promod3 import modelling
+ from promod3 import loop
+
+ tpl = io.LoadPDB('2dbs.pdb')
+ aln = ost.io.LoadAlignment('2dbs.fasta')
+ aln.AttachView(1, tpl.CreateFullView())
+
+ mhandle = modelling.BuildRawModel(aln)
+ scorer = loop.SetupBackboneScorer(mhandle)
+ scorer = modelling.FillLoopsByMonteCarlo(mhandle, scorer,
+ loop.LoadTorsionSamplerCoil())
+
+ :param mhandle: Modelling handle on which to apply change.
+ :type mhandle: :class:`ModellingHandle`
+
+ :param scorer: A scorer dedicated to this raw model.
+ :type scorer: :class:`~promod3.loop.BackboneLoopScorer`
+
+ :param torsion_sampler: A sampler for torsion angles.
+ :type torsion_sampler: :class:`~promod3.loop.TorsionSampler`
+
+ :param max_loops_to_search: Maximal number of loops to consider.
+ :type max_loops_to_search: :class:`int`
+
+ :param max_extension: Maximal number of gap extension steps to perform.
+ :type max_extension: :class:`int`
+
+ :param mc_num_loops: Number of loop candidates to consider for each gap
+ (see :meth:`promod3.loop.FillFromMonteCarloSampler`)
+ :param mc_steps: Number of MC steps to perform for each loop candidate
+ (see :meth:`promod3.loop.FillFromMonteCarloSampler`)
+
+ :return: If gaps are deleted, the scorer is updated returned.
+ :rtype: :class:`~promod3.loop.BackboneLoopScorer`
+ '''
+
+ # point to the current gap
+ gap_idx = 0
+
+ # Iterate all gaps. Since the number of gaps may change, always compare to
+ # an updated list. gap_idx is only increased when necessary, e.g. current
+ # gap could not be removed.
+ while gap_idx < len(mhandle.gaps):
+
+ # keep copy of original gap
+ gap_orig = mhandle.gaps[gap_idx].Copy()
+ # get chain for gap
+ actual_chain_idx = gap_orig.GetChainIndex()
+ actual_chain = mhandle.model.chains[actual_chain_idx]
+
+ # ignore terminal gaps
+ if gap_orig.IsTerminal():
gap_idx += 1
+ continue
+
+ ##################################
+ # find loop candidates
+ ##################################
+ # min/max stem res. numbers
+ min_before_resnum = sys.maxint
+ max_after_resnum = -sys.maxint-1
+ # list of LoopCandidates to fill gap
+ actual_candidates = list()
+ # list of StructuralGap corresponding to actual_candidates
+ actual_extended_gaps = list()
+ # number of loops found (stop if >= max_loops_to_search)
+ found_loops = 0
+ # currently extended gap and gap-extender
+ actual_gap = gap_orig.Copy()
+ actual_extender = modelling.GapExtender(actual_gap)
+
+ # make sure, that the loop seq has at least length 3
+ while (len(actual_gap.seq) < 3):
+ if not actual_extender.Extend():
+ break
+
+ # iteratively extend actual_gap
+ for _ in range(max_extension):
+ # extend the gap
+ if not actual_extender.Extend():
+ break
+ # ensure both stems are valid
+ if not actual_gap.before.IsValid() or \
+ not actual_gap.after.IsValid():
+ break
+
+ # setup sampler, closer, scorer and cooler for MC
+ try:
+ mc_sampler = loop.PhiPsiSampler(actual_gap.full_seq,
+ torsion_sampler)
+
+ mc_closer = loop.DirtyCCDCloser(actual_gap.before,
+ actual_gap.after)
+
+ weights = dict()
+ weights["cbeta"]=5.0
+ weights["cb_packing"]=1.0
+ weights["clash"]=0.05
+ mc_scorer = loop.LinearScorer(scorer, actual_gap.before.GetNumber().GetNum(),
+ actual_chain_idx, weights)
+
+ mc_cooler = loop.ExponentialCooler(200,100,0.9)
+ except:
+ # could throw an exception upon construction when the neighbouring
+ # residues of the stem residues are invalid
+ ost.LogError('Failed to set up MC components')
+ # TODO - ignore for now...??
+ continue
+
+ # try to get candidates for the current loop
+ try:
+ ost.LogVerbose("firing MC for " + str(actual_gap))
+ candidates = loop.LoopCandidates.FillFromMonteCarloSampler(
+ actual_gap.before, actual_gap.after, actual_gap.full_seq,
+ mc_num_loops, mc_steps, mc_sampler, mc_closer,
+ mc_scorer, mc_cooler)
+ except RuntimeError:
+ # monte carlo cannot be initialized when the stems are too far apart
+ # => we need further loop extension
+ ost.LogVerbose("failed in setting up " + str(actual_gap) +
+ "stems might be too far away...")
+ continue
+ # skip if nothing found
+ if len(candidates) == 0:
+ continue
+
+ # update candidate lists
+ min_before_resnum = min(min_before_resnum,
+ actual_gap.before.GetNumber().GetNum())
+ max_after_resnum = max(max_after_resnum,
+ actual_gap.after.GetNumber().GetNum())
+ actual_candidates.append(candidates)
+ actual_extended_gaps.append(actual_gap.Copy())
+ found_loops += len(candidates)
+
+ # abort if we found enough loops
+ if found_loops >= max_loops_to_search:
+ break
+
+ # check result
+ if len(actual_candidates) == 0:
+ ost.LogInfo("Failed at loop insterion (%s)" % str(gap_orig))
+ gap_idx += 1
+ continue
+
+ ##################################
+ # all loop candidates will be scored along the full max. extension ever
+ # reached in the search before, so we build an overall frame, where we
+ # insert the loops (NOTE: identical to DB)
+ ##################################
+ frame_seq = mhandle.seqres[actual_chain_idx]\
+ [min_before_resnum-1:max_after_resnum]
+ frame_backbone_list = loop.BackboneList(frame_seq)
+ actual_res_num = ost.mol.ResNum(min_before_resnum)
+ for i in range(len(frame_seq)):
+ actual_res = actual_chain.FindResidue(actual_res_num)
+ if actual_res.IsValid():
+ frame_backbone = loop.Backbone(actual_res, frame_seq[i])
+ frame_backbone_list[i] = frame_backbone
+ actual_res_num += 1
+
+ overall_before_residue = actual_chain.FindResidue(
+ ost.mol.ResNum(min_before_resnum))
+ overall_after_residue = actual_chain.FindResidue(
+ ost.mol.ResNum(max_after_resnum))
+ final_loop_candidates = loop.LoopCandidates(overall_before_residue,
+ overall_after_residue,
+ frame_seq)
+
+ ##################################
+ # prepare loop candidates for scoring (NOTE: identical to DB)
+ ##################################
+ back_mapper = list()
+ for i, loop_candidates in enumerate(actual_candidates):
+ start_index = loop_candidates.GetNStem().GetNumber().GetNum()\
+ - min_before_resnum
+ for j, loop_candidate in enumerate(loop_candidates):
+ actual_frame_backbone_list = frame_backbone_list
+ actual_frame_backbone_list.ReplaceFragment(
+ loop_candidate.bb_list, start_index, False)
+ final_loop_candidates.Add(actual_frame_backbone_list)
+ back_mapper.append((i,j))
+
+ final_loop_candidates.AttachScorer(scorer)
+ final_loop_candidates.CalculateClashScores(actual_chain_idx)
+ final_loop_candidates.CalculateCBetaScores(actual_chain_idx)
+ final_loop_candidates.CalculateTorsionScores(actual_chain_idx)
+ final_loop_candidates.CalculateCBPackingScores(actual_chain_idx)
+ final_loop_candidates.CalculateHBondScores(actual_chain_idx)
+
+ ##################################
+ # compare scores
+ ##################################
+ min_score = sys.maxint
+ optimal_candidate = -1
+ for i in range(len(final_loop_candidates)):
+ score = 0.0
+ score += 5.0*final_loop_candidates[i].cbeta_score
+ score += 1.0*final_loop_candidates[i].packing_score
+ score += 0.2*final_loop_candidates[i].clash_score
+ score += 0.004*final_loop_candidates[i].hbond_score
+ if score < min_score:
+ optimal_candidate = i
+ min_score = score
+
+ ##################################
+ # finally resolve loop (NOTE: identical to DB)
+ ##################################
+ if optimal_candidate != -1:
+ idx_a = back_mapper[optimal_candidate][0]
+ idx_b = back_mapper[optimal_candidate][1]
+ actual_candidates[idx_a].InsertInto(mhandle.model, idx_b)
+ bb_list = actual_candidates[idx_a][idx_b].bb_list
+ scorer.SetEnvironment(\
+ bb_list,
+ actual_candidates[idx_a].GetNStem().GetNumber().GetNum(),
+ actual_chain_idx)
+ ost.LogInfo("Resolved %s by filling %s" % \
+ (str(gap_orig), str(actual_extended_gaps[idx_a])))
+ next_gap = modelling.ClearGaps(mhandle, actual_extended_gaps[idx_a])
+ if next_gap == -1:
+ break
+ else:
+ gap_idx = next_gap
+ else:
+ ost.LogInfo("Failed at loop insterion (%s)" % str(gap_orig))
+ gap_idx += 1
+
return scorer
# these methods will be exported into module
__all__ = ('CloseSmallDeletions', 'MergeGapsByDistance',
- 'FillLoopsByDatabase',)
+ 'FillLoopsByDatabase', 'FillLoopsByMonteCarlo')
# LocalWords: modeling stereochemically param idx init
# LocalWords: py ost pylint modelling promod CloseSmallDeletions testcode
diff --git a/modelling/src/model.cc b/modelling/src/model.cc
index 12379bca..15d91575 100644
--- a/modelling/src/model.cc
+++ b/modelling/src/model.cc
@@ -655,7 +655,7 @@ void BuildRawChain(const ost::seq::AlignmentHandle& aln,
// that's an insertion
gap_list.push_back(StructuralGap(before, dst_res, gap_seq));
gap_seq="";
- } else if (last_index+1!=i) {
+ } else if (last_index > -1 && last_index+1 != i) {
// that's a deletion
gap_list.push_back(StructuralGap(before, dst_res, ""));
//form_peptide_bond=true;
diff --git a/modelling/tests/test_close_gaps.py b/modelling/tests/test_close_gaps.py
index d833591e..c30e2565 100644
--- a/modelling/tests/test_close_gaps.py
+++ b/modelling/tests/test_close_gaps.py
@@ -11,7 +11,8 @@ class _FetchLog(ost.LogSink):
self.messages = dict()
def LogMessage(self, message, severity):
- levels = ['ERROR', 'WARNING', 'DEBUG', 'INFO', 'VERBOSE', 'TRACE']
+ levels = ['ERROR', 'WARNING', 'SCRIPT', 'INFO', 'VERBOSE', 'DEBUG',
+ 'TRACE']
level = levels[severity]
if not level in self.messages.keys():
self.messages[level] = list()
@@ -23,7 +24,14 @@ class CloseGapsTests(unittest.TestCase):
ost.PushLogSink(self.log)
ost.PushVerbosityLevel(4)
- def testClosesmalldel(self):
+ @classmethod
+ def setUpClass(cls):
+ # load dbs etc here for all tests
+ cls.frag_db = loop.LoadFragDB()
+ cls.structure_db = loop.LoadStructureDB()
+ cls.torsion_sampler = loop.LoadTorsionSamplerCoil()
+
+ def testCloseSmallDel(self):
# check that very small gaps are closed
# create a raw model to work with (which actually is a bit dangerous:
# the PDB file has nothing to do with the alignment...)
@@ -36,14 +44,13 @@ class CloseGapsTests(unittest.TestCase):
self.assertEqual(len(rmodel.gaps), 1)
# obtain the scorer
scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
modelling.CloseSmallDeletions(rmodel, scorer)
self.assertEqual(len(rmodel.gaps), 0)
- self.assertEqual(self.log.messages['VERBOSE'],
- ['Assigning MOL_IDs',
- 'selected 1 chain(s), 20 residue(s), 80 atom(s) 79 '+
- 'bond(s)',
- 'Closed: A.GLY3-()-A.GLY4 by relaxing '+
- 'A.GLY3-(GG)-A.GLY6'])
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Closed: A.GLY3-()-A.GLY4 by relaxing A.GLY3-(GG)-A.GLY6')
def testMergeGapsByDistance(self):
# test that merging two close gaps works
@@ -55,16 +62,15 @@ class CloseGapsTests(unittest.TestCase):
self.assertEqual(len(rmodel.gaps), 2)
self.assertEqual(str(rmodel.gaps[0]), 'A.ASP2-(FAGD)-A.THR7')
self.assertEqual(str(rmodel.gaps[1]), 'A.THR7-(KNLG)-A.HIS12')
- modelling.MergeGapsByDistance(rmodel,0)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.MergeGapsByDistance(rmodel, 0)
self.assertEqual(len(rmodel.gaps), 1)
self.assertEqual(str(rmodel.gaps[0]), 'A.ASP2-(FAGDTKNLG)-A.HIS12')
- self.assertEqual(self.log.messages['VERBOSE'],
- ['Assigning MOL_IDs',
- 'selected 1 chain(s), 5 residue(s), 40 atom(s) 40 '+
- 'bond(s)',
- 'Merged gap A.ASP2-(FAGD)-A.THR7 and '+
- 'A.THR7-(KNLG)-A.HIS12 into '+
- 'A.ASP2-(FAGDTKNLG)-A.HIS12'])
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Merged gap A.ASP2-(FAGD)-A.THR7 and A.THR7-(KNLG)-A.HIS12'+
+ ' into A.ASP2-(FAGDTKNLG)-A.HIS12')
def testMergeGapsByDistanceBothTerminals(self):
# test that we do not delete the whole thing for gaps at terminals
@@ -74,7 +80,7 @@ class CloseGapsTests(unittest.TestCase):
aln.AttachView(1, tpl.CreateFullView())
rmodel = modelling.BuildRawModel(aln)
self.assertEqual(len(rmodel.gaps), 2)
- modelling.MergeGapsByDistance(rmodel,4)
+ modelling.MergeGapsByDistance(rmodel, 4)
self.assertEqual(len(rmodel.gaps), 2)
def testMergeGapsByDistanceOneTerminal(self):
@@ -85,11 +91,36 @@ class CloseGapsTests(unittest.TestCase):
aln.AttachView(1, tpl.CreateFullView())
rmodel = modelling.BuildRawModel(aln)
self.assertEqual(len(rmodel.gaps), 2)
- modelling.MergeGapsByDistance(rmodel,2)
+ modelling.MergeGapsByDistance(rmodel, 2)
self.assertEqual(len(rmodel.gaps), 1)
def testFillGapsByDatabase(self):
- # get rid of gaps by db
+ # get rid of 2 gaps by db
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLV--LNPDKGATVTMAAAA'),
+ seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSN---VHIRVE'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 2)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByDatabase(rmodel, scorer, self.frag_db,
+ self.structure_db,
+ self.torsion_sampler)
+ self.assertEqual(len(rmodel.gaps), 0)
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 2)
+ self.assertEqual(self.log.messages['INFO'][-2],
+ 'Resolved A.VAL14-()-A.LEU15 by filling '+
+ 'A.ASN11-(TLVLNP)-A.ASP18')
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.GLY20-(ATV)-A.THR24 by filling '+
+ 'A.GLY20-(ATVTMAA)-A.ALA28')
+
+ def testFillGapsByDatabaseDeletion(self):
+ # get rid of deletion by db
tpl = io.LoadPDB('data/2dbs.pdb')
aln = seq.CreateAlignment(
seq.CreateSequence('trg', 'TLNGFTVPAGNTLV--LNPDKGATVTMA'),
@@ -99,18 +130,139 @@ class CloseGapsTests(unittest.TestCase):
rmodel = modelling.BuildRawModel(aln)
self.assertEqual(len(rmodel.gaps), 1)
scorer = loop.SetupBackboneScorer(rmodel)
- modelling.FillLoopsByDatabase(rmodel, scorer, loop.LoadFragDB(),
- loop.LoadStructureDB(),
- loop.LoadTorsionSamplerCoil())
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByDatabase(rmodel, scorer, self.frag_db,
+ self.structure_db,
+ self.torsion_sampler)
self.assertEqual(len(rmodel.gaps), 0)
- self.assertEqual(self.log.messages['VERBOSE'],
- ['Assigning MOL_IDs',
- 'selected 1 chain(s), 28 residue(s), 218 atom(s) '+
- '222 bond(s)',
- 'Resolved A.PHE5-(TVPAGNTLV)-A.LEU15 by filling '+
- 'A.GLY10-(NTLVLNPD)-A.LYS19'])
-
-# fill db: terminal loops are not altered
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.VAL14-()-A.LEU15 by filling '+
+ 'A.GLY10-(NTLVLNPD)-A.LYS19')
+
+ def testFillGapsByDatabaseInsertion(self):
+ # get rid of insertion by db
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVAALNPDKGGAGATVTMA'),
+ seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSN---VHIRVE'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 1)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByDatabase(rmodel, scorer, self.frag_db,
+ self.structure_db,
+ self.torsion_sampler)
+ self.assertEqual(len(rmodel.gaps), 0)
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.GLY22-(GAG)-A.ALA26 by filling '+
+ 'A.GLY22-(GAGATVT)-A.MET30')
+
+ def testFillGapsByDatabaseTerminal(self):
+ # terminal loops are not altered
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVAALNPDKGGAGATVTMA'),
+ seq.CreateSequence('tpl', '-NGGTLLIPNGTYHFLGIQMKSNVHIRVE--'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 2)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ modelling.FillLoopsByDatabase(rmodel, scorer, self.frag_db,
+ self.structure_db,
+ self.torsion_sampler)
+ self.assertEqual(len(rmodel.gaps), 2)
+
+ def testFillGapsByMonteCarlo(self):
+ # get rid of 2 gaps by db
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVA-LNPDKGATVTMAA'),
+ seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSN-VHIRVE'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 2)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByMonteCarlo(rmodel, scorer,
+ self.torsion_sampler,
+ mc_steps=100)
+ self.assertEqual(len(rmodel.gaps), 0)
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 2)
+ self.assertEqual(self.log.messages['INFO'][-2],
+ 'Resolved A.ALA15-()-A.LEU16 by filling '+
+ 'A.LEU13-(VALNP)-A.ASP19')
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.GLY21-(A)-A.THR23 by filling '+
+ 'A.PRO18-(DKGA)-A.THR23')
+
+ def testFillGapsByMonteCarloDeletion(self):
+ # get rid of deletion by db
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVA-LNPDKGATVTMA'),
+ seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSNVHIRVE'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 1)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByMonteCarlo(rmodel, scorer,
+ self.torsion_sampler,
+ mc_steps=100)
+ self.assertEqual(len(rmodel.gaps), 0)
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.ALA15-()-A.LEU16 by filling '+
+ 'A.ASN11-(TLVAL)-A.ASN17')
+
+ def testFillGapsByMonteCarloInsertion(self):
+ # get rid of insertion by db
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVAALNPDKGAGTVTMA'),
+ seq.CreateSequence('tpl', 'NGGTLLIPNGTYHFLGIQMKSN-VHIRVE'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 1)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ nlogs = len(self.log.messages['INFO'])
+ modelling.FillLoopsByMonteCarlo(rmodel, scorer,
+ self.torsion_sampler,
+ mc_steps=100)
+ self.assertEqual(len(rmodel.gaps), 0)
+ # check logs
+ self.assertEqual(len(self.log.messages['INFO'])-nlogs, 1)
+ self.assertEqual(self.log.messages['INFO'][-1],
+ 'Resolved A.GLY22-(A)-A.GLY24 by filling '+
+ 'A.PRO19-(DKGA)-A.GLY24')
+
+ def testFillGapsByMonteCarloTerminal(self):
+ # terminal loops are not altered
+ tpl = io.LoadPDB('data/2dbs.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', 'TLNGFTVPAGNTLVAALNPDKGGAGATVTMA'),
+ seq.CreateSequence('tpl', '-NGGTLLIPNGTYHFLGIQMKSNVHIRVE--'))
+
+ aln.AttachView(1, tpl.CreateFullView())
+ rmodel = modelling.BuildRawModel(aln)
+ self.assertEqual(len(rmodel.gaps), 2)
+ scorer = loop.SetupBackboneScorer(rmodel)
+ modelling.FillLoopsByMonteCarlo(rmodel, scorer,
+ self.torsion_sampler,
+ mc_steps=100)
+ self.assertEqual(len(rmodel.gaps), 2)
if __name__ == "__main__":
from ost import testutils
diff --git a/modelling/tests/test_modelling.py b/modelling/tests/test_modelling.py
index 69c79bf9..2e6eff48 100644
--- a/modelling/tests/test_modelling.py
+++ b/modelling/tests/test_modelling.py
@@ -68,6 +68,16 @@ class ModellingTests(unittest.TestCase):
self.assertEqual(result.gaps[1].after, mol.ResidueHandle())
self.assertEqual(result.gaps[1].seq, 'G')
+ def testTerTrg(self):
+ # test if raw model ignores terminal gaps in target sequence
+ tpl = io.LoadPDB('data/gly.pdb')
+ aln = seq.CreateAlignment(
+ seq.CreateSequence('trg', '--GG--'),
+ seq.CreateSequence('tpl', 'GGGGGG'))
+ aln.AttachView(1, tpl.CreateFullView())
+ result = modelling.BuildRawModel(aln)
+ self.assertEqual(len(result.gaps), 0)
+
def testModified(self):
# test if we correctly strip off modifications
tpl = io.LoadPDB('data/sep.pdb')
diff --git a/pipeline/pymod/run_engine.py b/pipeline/pymod/run_engine.py
index 0acb61b2..791ea643 100644
--- a/pipeline/pymod/run_engine.py
+++ b/pipeline/pymod/run_engine.py
@@ -39,7 +39,7 @@ def BuildFromRawModel(mhandle, structure_db=None, fragment_db=None,
:return: Delivers the model as an |ost_s| entity.
:rtype: :class:`Entity <ost.mol.EntityHandle>`
'''
- ost.LogVerbose("Starting modelling based on a raw model")
+ ost.LogInfo("Starting modelling based on a raw model")
# a bit of setup
if structure_db == None:
@@ -52,22 +52,22 @@ def BuildFromRawModel(mhandle, structure_db=None, fragment_db=None,
scorer = loop.SetupBackboneScorer(mhandle)
# step 1: close small deletions in the raw model
- ost.LogVerbose("Trying to close small deletions (no. of gaps: %d)." % \
+ ost.LogInfo("Trying to close small deletions (no. of gaps: %d)." % \
len(mhandle.gaps))
- scorer = modelling.CloseSmallDeletions(mhandle, scorer)
+ modelling.CloseSmallDeletions(mhandle, scorer)
# 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.LogVerbose("Trying to merge gaps (%d) with distance %d." % \
+ ost.LogInfo("Trying to merge gaps (%d) with distance %d." % \
(len(mhandle.gaps), distance))
modelling.MergeGapsByDistance(mhandle, distance)
- ost.LogVerbose("%d gap(s) left after merging.\n" % len(mhandle.gaps)+
+ ost.LogInfo("%d gap(s) left after merging.\n" % len(mhandle.gaps)+
"Trying to fill loops by database")
# step 2b: fit loops into the model
- scorer = modelling.FillLoopsByDatabase(mhandle, scorer, fragment_db,
- structure_db, torsion_sampler,
- max_loops_to_search=30)
- ost.LogVerbose("%d gap(s) left after database search." % \
+ modelling.FillLoopsByDatabase(mhandle, scorer, fragment_db,
+ structure_db, torsion_sampler,
+ max_loops_to_search=30)
+ ost.LogInfo("%d gap(s) left after database search." % \
len(mhandle.gaps))
if len(mhandle.gaps) == 0:
break
@@ -77,6 +77,6 @@ def BuildFromRawModel(mhandle, structure_db=None, fragment_db=None,
__all__ = ('BuildProMod3Model')
# LocalWords: BuildProMod ProMod sticked ost promod BuildFromRawModel param
-# LocalWords: rtype bool LoadStructureDB LogVerbose LoadFragDB
+# LocalWords: rtype bool LoadStructureDB LogInfo LoadFragDB
# LocalWords: LoadTorsionSamplerCoil SetupBackboneScorer
# LocalWords: CloseSmallDeletions
--
GitLab