move all operations, that should lead to a gapless model into one function

The idea is that this whole procedure can then be applied to only
a part of the model by providing chain_idx and resnum ranges.

doc/tests/scripts/modelling_steps.py

@@ -13,32 +13,20 @@ aln = io.LoadAlignment('data/1crn.fasta')
 aln.AttachView(1, tpl.CreateFullView())
 mhandle = modelling.BuildRawModel(aln)
 
-# perform loop modelling to close all gaps
-modelling.SetupDefaultBackboneScorer(mhandle)
+# we're not modelling termini
 modelling.RemoveTerminalGaps(mhandle)
-modelling.CloseSmallDeletions(mhandle)
-for distance in range(merge_distance):
-    modelling.MergeGapsByDistance(mhandle, distance)
-    modelling.FillLoopsByDatabase(mhandle, fragment_db,
-                                  structure_db, torsion_sampler,
-                                  min_loops_required=-1,
-                                  max_res_extension=6)
-# if above fails, try DB-fill with less restrictions
-modelling.FillLoopsByDatabase(mhandle, fragment_db,
-                              structure_db, torsion_sampler,
-                              min_loops_required=-1)
-modelling.FillLoopsByDatabase(mhandle, fragment_db,
-                              structure_db, torsion_sampler)
-# if above fails on some gaps, use Monte Carlo
-modelling.FillLoopsByMonteCarlo(mhandle, torsion_sampler)
-# as a last resort, try to close large deletions
-modelling.CloseLargeDeletions(mhandle, structure_db)
+
+# perform loop modelling to close all gaps
+modelling.CloseGaps(mhandle)
+
 # build sidechains
 modelling.BuildSidechains(mhandle, merge_distance,
                           fragment_db, structure_db,
                           torsion_sampler)
+
 # minimize energy of final model using molecular mechanics
 modelling.MinimizeModelEnergy(mhandle)
+
 # check final model and report issues
 modelling.CheckFinalModel(mhandle)
 

modelling/doc/pipeline.rst

@@ -9,7 +9,9 @@ A protein homology modelling pipeline has the following main steps:
 - Perform loop modelling to close (or remove) all gaps (see functions
   :func:`CloseSmallDeletions`, :func:`RemoveTerminalGaps`,
   :func:`MergeGapsByDistance`, :func:`FillLoopsByDatabase`,
-  :func:`FillLoopsByMonteCarlo`, :func:`CloseLargeDeletions`)
+  :func:`FillLoopsByMonteCarlo`, :func:`CloseLargeDeletions` or
+  :func:`CloseGaps` that calls all these functions using predefined
+  heuristics)
 - Build sidechains (see :func:`BuildSidechains` function)
 - Minimize energy of final model using molecular mechanics
   (see :func:`MinimizeModelEnergy` function)
@@ -23,13 +25,6 @@ purposes:
 
 .. literalinclude:: ../../../tests/doc/scripts/modelling_steps.py
 
-In the default pipeline above, we call :func:`FillLoopsByDatabase` multiple
-times. First, we try to close "easy" gaps which require few extensions (we wish
-to limit the damage we do on the template) and for which we have plenty of loop
-candidates. If some gaps cannot be closed like this, we try less restrictive
-options. This approach is helpful if neighboring gaps are close together and the
-one closer to the C-terminus is easier to close. Several variants of the
-pipeline were evaluated on 1752 target-template-pairs and this one worked best.
 
 Build Raw Modelling Handle
 --------------------------------------------------------------------------------
@@ -207,6 +202,8 @@ Modelling Steps
   :return: Number of gaps which were removed.
   :rtype:  :class:`int`
 
+.. autofunction:: CloseGaps
+
 .. autofunction:: CloseSmallDeletions
 
 .. autofunction:: MergeGapsByDistance

modelling/pymod/_closegaps.py

@@ -393,7 +393,8 @@ def CloseSmallDeletions(mhandle, max_extension=9, clash_thresh=1.0,
         if not success:
             current_gap_index += 1
 
-def MergeGapsByDistance(mhandle, distance):
+def MergeGapsByDistance(mhandle, distance, chain_idx = None, 
+                        resnum_range = None):
     '''Merge 2 neighbouring gaps by deleting residues in-between.
 
     Check if two neighbouring gaps are at max. *distance* residues apart from
@@ -411,6 +412,14 @@ def MergeGapsByDistance(mhandle, distance):
     :param distance: The max. no. of residues between two gaps up to which
                      merge happens.
     :type distance: :class:`int`
+
+    :param chain_idx: If not None, only gaps from chain with given index get
+                      processed
+    :type chain_idx:  :class:`int`
+
+    :param resnum_range: If not None, two gaps only get merged if they're
+                         both in this resnum range.
+    :type resnum_range: :class:`tuple` containing two :class:`int`
     '''
 
     prof_name = 'closegaps::MergeGapsByDistance'
@@ -438,6 +447,10 @@ def MergeGapsByDistance(mhandle, distance):
             # check that we are on the same chain
             if current_gap.GetChain() != next_gap.GetChain():
                 continue
+            # check for range (if given)
+            if not (_InRange(current_gap, chain_idx, resnum_range) and\
+                    _InRange(next_gap, chain_idx, resnum_range)):
+                continue
             # no merging of gaps at the end AND the start :)
             if current_gap.IsNTerminal() and next_gap.IsCTerminal():
                 continue
@@ -1013,7 +1026,8 @@ def ModelTermini(mhandle, torsion_sampler, fragger_handles=None,
 
 def CloseLargeDeletions(mhandle, structure_db, linker_length=8,
                         num_fragments=500, use_scoring_extender=True,
-                        use_full_extender=True):
+                        use_full_extender=True, chain_idx=None,
+                        resnum_range=None):
 
     '''Try to close large deletions.
 
@@ -1046,6 +1060,14 @@ def CloseLargeDeletions(mhandle, structure_db, linker_length=8,
                               :class:`GapExtender`.
                               See :func:`CloseSmallDeletions`.
     :type use_full_extender:  :class:`bool`
+
+    :param chain_idx: If not None, only gaps from chain with given index get
+                      processed
+    :type chain_idx:  :class:`int`
+
+    :param resnum_range: If not None, only gaps within this resnum range get
+                         processed
+    :type resnum_range: :class:`tuple` containing two :class:`int`
     '''
 
     prof_name = 'closegaps::CloseLargeDeletions'
@@ -1067,6 +1089,12 @@ def CloseLargeDeletions(mhandle, structure_db, linker_length=8,
 
         # keep copy of original gap
         gap_orig = mhandle.gaps[gap_idx].Copy()
+
+        # check whether we are in the desired range
+        if not _InRange(gap_orig, chain_idx, resnum_range):
+            gap_idx += 1
+            continue
+
         # get chain for gap
         actual_chain_idx = gap_orig.GetChainIndex()
         actual_chain = mhandle.model.chains[actual_chain_idx]

modelling/pymod/_pipeline.py

@@ -162,6 +162,99 @@ def _GetSimEntity(sim):
     return ent
 ###############################################################################
 
+def CloseGaps(mhandle, merge_distance=4, fragment_db=None, structure_db=None, 
+              torsion_sampler=None, chain_idx = None, resnum_range = None):
+    '''
+    Tries to close all gaps in a model, except termini. It will go through
+    following steps:
+
+    - Try to close small deletions by relaxing them 
+      (see :func:`CloseSmallDeletions`)
+    - Iteratively merge gaps up to a distance **merge_distance**
+      (see :func:`MergeGapsByDistance`) and try to fill them with a database 
+      approach (see :func:`FillLoopsByDatabase`)
+    - Try to fill remaining gaps using a Monte Carlo approach
+      (see :func:`FillLoopsByMonteCarlo`)
+    - Large deletions get closed using a last resort approach
+      (see :func:`CloseLargeDeletions`)
+
+    :param mhandle: Modelling handle on which to apply change.
+    :type mhandle:  :class:`ModellingHandle`
+
+    :param merge_distance:  Max. merge distance when performing the database 
+                            approach
+    :type merge_distance:   :class:`int`
+    :param fragment_db:     Database for searching fragments in database 
+                            approach, must be consistent with provided
+                            **structure_db**. A default is loaded if None.
+    :type fragment_db:      :class:`~promod3.loop.FragDB`
+    :param structure_db:    Structure db from which the **fragment_db** gets
+                            it's structural information. A default is loaded 
+                            if None.
+    :type structure_db:     :class:`~promod3.loop.StructureDB`
+    :param torsion_sampler: Used as parameter for :func:`FillLoopsByDatabase`
+                            and :func:`FillLoopsByMonteCarlo` A default one is 
+                            loaded if None.
+    :param chain_idx: If not None, only gaps from chain with given index get
+                      processed
+    :type chain_idx:  :class:`int`
+
+    :param resnum_range: If not None, only gaps within this resnum range get
+                         processed.
+    :type resnum_range: :class:`tuple` containing two :class:`int`  
+
+    '''
+
+
+    # load stuff if needed
+    if not IsBackboneScorerSet(mhandle) or not IsBackboneScorerEnvSet(mhandle):
+        SetupDefaultBackboneScorer(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()
+
+    #try to close small deletions by relaxing them
+    CloseSmallDeletions(mhandle, chain_idx = chain_idx, 
+                        resnum_range = resnum_range)
+
+    # iteratively merge gaps of distance i and fill loops by database
+    for distance in range(merge_distance):
+        MergeGapsByDistance(mhandle, distance, chain_idx = chain_idx,
+                            resnum_range = resnum_range)
+        FillLoopsByDatabase(mhandle, fragment_db, structure_db,
+                            torsion_sampler, min_loops_required=-1,
+                            max_res_extension=6, chain_idx = chain_idx,
+                            resnum_range = resnum_range)
+        
+    # if above fails, try DB-fill with less restrictions
+    FillLoopsByDatabase(mhandle, fragment_db, structure_db,
+                        torsion_sampler, min_loops_required=-1,
+                        chain_idx = chain_idx, resnum_range = resnum_range)
+    FillLoopsByDatabase(mhandle, fragment_db, structure_db,
+                        torsion_sampler, chain_idx = chain_idx,
+                        resnum_range = resnum_range)
+
+    # close remaining gaps by Monte Carlo
+    FillLoopsByMonteCarlo(mhandle, torsion_sampler, chain_idx = chain_idx,
+                          resnum_range = resnum_range)
+
+    # last resort approach to close large deletions
+    CloseLargeDeletions(mhandle, structure_db, chain_idx = chain_idx, 
+                        resnum_range = resnum_range)
+
+    #In the function above, we call :func:`FillLoopsByDatabase` multiple
+    #times. First, we try to close "easy" gaps which require few extensions 
+    #(we wish to limit the damage we do on the template) and for which we have 
+    #plenty of loop candidates. If some gaps cannot be closed like this, we try 
+    #less restrictive options. This approach is helpful if neighboring gaps are 
+    #close together and the one closer to the C-terminus is easier to close. 
+    #Several variants were evaluated on 1752 target-template-pairs and this one 
+    #worked best.
+
+
 def BuildSidechains(mhandle, merge_distance=4, fragment_db=None,
                     structure_db=None, torsion_sampler=None):
     '''Build sidechains for model.
@@ -467,24 +560,13 @@ def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     else:
         SetupDefaultBackboneScorer(mhandle)
 
-    # remove terminal gaps and close small deletions
+    # remove terminal gaps
     RemoveTerminalGaps(mhandle)
-    CloseSmallDeletions(mhandle)
-
-    # iteratively merge gaps of distance i and fill loops by database
-    for distance in range(merge_distance):
-        MergeGapsByDistance(mhandle, distance)
-        FillLoopsByDatabase(mhandle, fragment_db, structure_db,
-                            torsion_sampler, min_loops_required=-1,
-                            max_res_extension=6)
-    FillLoopsByDatabase(mhandle, fragment_db, structure_db,
-                        torsion_sampler, min_loops_required=-1)
-    FillLoopsByDatabase(mhandle, fragment_db, structure_db,
-                        torsion_sampler)
 
-    # close remaining gaps by Monte Carlo
-    FillLoopsByMonteCarlo(mhandle, torsion_sampler)
-    CloseLargeDeletions(mhandle, structure_db)
+    #close gaps
+    CloseGaps(mhandle, merge_distance = merge_distance, 
+              fragment_db = fragment_db, structure_db = structure_db, 
+              torsion_sampler = torsion_sampler)
 
     # build sidechains
     BuildSidechains(mhandle, merge_distance, fragment_db,
@@ -501,5 +583,5 @@ def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     return mhandle.model
 
 # these methods will be exported into module
-__all__ = ('BuildFromRawModel', 'BuildSidechains', 'MinimizeModelEnergy',
-           'CheckFinalModel')
+__all__ = ('BuildFromRawModel', 'CloseGaps','BuildSidechains', 
+           'MinimizeModelEnergy', 'CheckFinalModel')