From 98809cc4a3400276ba7fb682bf92fd85c6ccc462 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Tue, 15 Nov 2016 16:17:20 +0100
Subject: [PATCH] give alternative ModellingHandles as optional argument in
 BuildFromRawModel

---
 modelling/pymod/_closegaps.py |  95 +++++++++++++++++++++++++++-
 modelling/pymod/_mtm.py       |  75 ++++++++++++----------
 modelling/pymod/_pipeline.py  | 113 ++++++----------------------------
 3 files changed, 154 insertions(+), 129 deletions(-)

diff --git a/modelling/pymod/_closegaps.py b/modelling/pymod/_closegaps.py
index 2fd43444..651c2bac 100644
--- a/modelling/pymod/_closegaps.py
+++ b/modelling/pymod/_closegaps.py
@@ -918,6 +918,98 @@ def FillLoopsByMonteCarlo(mhandle, torsion_sampler, max_loops_to_search=6,
             gap_idx = new_idx
 
 
+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 ModelTermini(mhandle, torsion_sampler, fragger_handles=None,
                  mc_num_loops=20, mc_steps=5000):
     '''Try to model termini with Monte Carlo sampling.
@@ -1242,4 +1334,5 @@ def CloseLargeDeletions(mhandle, structure_db, linker_length=8,
 
 # these methods will be exported into module
 __all__ = ('CloseSmallDeletions', 'MergeGapsByDistance', 'FillLoopsByDatabase',
-           'FillLoopsByMonteCarlo', 'ModelTermini', 'CloseLargeDeletions')
+           'FillLoopsByMonteCarlo', 'CloseGaps', 'ModelTermini', 
+           'CloseLargeDeletions')
diff --git a/modelling/pymod/_mtm.py b/modelling/pymod/_mtm.py
index 6d287b60..deee37f0 100644
--- a/modelling/pymod/_mtm.py
+++ b/modelling/pymod/_mtm.py
@@ -1,10 +1,9 @@
+from promod3 import loop
+from _modelling import *
+from _pipeline import *
 import sys
 import numpy as np
 from ost import seq, geom, mol
-from promod3 import loop
-import _modelling as modelling
-from _pipeline import CloseGaps
-from promod3.loop import FraggerHandle
 
 def _CheckAlternativeRawmodels(seed_rawmodel, alternative_rawmodels):
 
@@ -231,13 +230,13 @@ def _GetRigidBlocksSingleChain(rawmodel_one, rawmodel_two,
         return (list(), list(), list())
   
     if invert_transforms:
-        rigid_block_result = modelling.RigidBlocks(ca_positions_two_raw, 
-                                                   ca_positions_one_raw, 
-                                                   window_size, 20, 2.0, 0.9)
+        rigid_block_result = RigidBlocks(ca_positions_two_raw, 
+                                         ca_positions_one_raw, 
+                                         window_size, 20, 2.0, 0.9)
     else:
-        rigid_block_result = modelling.RigidBlocks(ca_positions_one_raw, 
-                                                   ca_positions_two_raw, 
-                                                   window_size, 20, 2.0, 0.9)
+        rigid_block_result = RigidBlocks(ca_positions_one_raw, 
+                                         ca_positions_two_raw, 
+                                         window_size, 20, 2.0, 0.9)
   
     rigid_blocks = rigid_block_result[0]
     rigid_block_transformations = rigid_block_result[1]
@@ -286,14 +285,14 @@ class ExtensionResult:
 
 def _ApplyExtensionResult(seed_rawmodel, seed_ch_idx, extension_result):
 
-    modelling.MergeMHandle(extension_result.arm, 
-                           seed_rawmodel, 
-                           extension_result.arm_ch_idx, 
-                           seed_ch_idx, 
-                           extension_result.start_rnum, 
-                           extension_result.end_rnum, 
-                           extension_result.transform,
-                           reset_scoring_env = True)
+    MergeMHandle(extension_result.arm, 
+                 seed_rawmodel, 
+                 extension_result.arm_ch_idx, 
+                 seed_ch_idx, 
+                 extension_result.start_rnum, 
+                 extension_result.end_rnum, 
+                 extension_result.transform,
+                 reset_scoring_env = True)
     # do linker if present
     if extension_result.linker_bb_list != None:
         linker = extension_result.linker_bb_list
@@ -337,7 +336,7 @@ def _NonOverlapLinkerSampler(seed_rawmodel, arm, seed_ch_idx, arm_ch_idx,
     # setup everything still being None
     if termini_type == "nter":
         anchor = seed_chain.FindResidue(mol.ResNum(linker_range[1]))
-        closer = modelling.NTerminalCloser(anchor)
+        closer = NTerminalCloser(anchor)
         arm_residues = list()
         for num in range(arm_range[0], arm_range[1] + 1):
             arm_res = arm_chain.FindResidue(mol.ResNum(num))
@@ -351,7 +350,7 @@ def _NonOverlapLinkerSampler(seed_rawmodel, arm, seed_ch_idx, arm_ch_idx,
   
     if termini_type == "cter":
         anchor = seed_chain.FindResidue(mol.ResNum(linker_range[0]))
-        closer = modelling.CTerminalCloser(anchor)
+        closer = CTerminalCloser(anchor)
         arm_residues = list()
         for num in range(arm_range[0], arm_range[1] + 1):
             arm_res = arm_chain.FindResidue(mol.ResNum(num))
@@ -365,19 +364,19 @@ def _NonOverlapLinkerSampler(seed_rawmodel, arm, seed_ch_idx, arm_ch_idx,
         full_bb_list.ReplaceFragment(arm_bb_list, arm_replace_idx, True)
 
     # do the sampling
-    num_fragments = 1000
+    num_fragments = 10000
     psipred_prediction = None
     profile = None
     if len(seed_rawmodel.profiles) > 0:
         profile = seed_rawmodel.profiles[seed_ch_idx]
     if len(seed_rawmodel.psipred_predictions) > 0:
         psipred_prediction = seed_rawmodel.psipred_predictions[seed_ch_idx]
-    fragger_handle = FraggerHandle(seed_rawmodel.seqres[seed_ch_idx],
-                                   profile = profile,
-                                   psipred_pred = psipred_prediction,
-                                   fragment_length = len(linker_seq),
-                                   fragments_per_position = num_fragments,
-                                   structure_db = structure_db)
+    fragger_handle = loop.FraggerHandle(seed_rawmodel.seqres[seed_ch_idx],
+                                        profile = profile,
+                                        psipred_pred = psipred_prediction,
+                                        fragment_length = len(linker_seq),
+                                        fragments_per_position = num_fragments,
+                                        structure_db = structure_db)
 
     fragger = fragger_handle.Get(linker_range[0] - 1)
 
@@ -537,7 +536,7 @@ def _FindNonOverlappingExtension(seed_rawmodel, arm,
     # Before we can extract a BackboneList from the arm we have to close
     # all gaps in the arm of the current chain in the region of interest
     copied_arm = arm.Copy()
-    modelling.SetupDefaultBackboneScorer(copied_arm)
+    SetupDefaultBackboneScorer(copied_arm)
     CloseGaps(copied_arm, chain_idx = arm_ch_idx, resnum_range = arm_range)
 
     transform = None
@@ -647,7 +646,7 @@ def _FindOverlappingExtension(seed_rawmodel, arm,
     #we do not operate on the arm itself, but rather on a copy of it
     copied_arm = arm.Copy()
     arm_chain = arm.model.chains[arm_ch_idx]
-    modelling.SetupDefaultBackboneScorer(copied_arm)
+    SetupDefaultBackboneScorer(copied_arm)
   
     if termini_type == "nter":
         CloseGaps(copied_arm, chain_idx = arm_ch_idx,
@@ -712,7 +711,7 @@ def _FindOverlappingExtension(seed_rawmodel, arm,
 
 
 def ModelExtensions(seed_rawmodel, alternative_rawmodels,
-                    scoring_weights, extension_strategy = "overlapping",
+                    scoring_weights = None, extension_strategy = "overlapping",
                     max_linker_extension = 0, structure_db = None):
     '''
     Tries to extend every chain of the **seed_rawmodel** towards nter and cter 
@@ -768,10 +767,20 @@ def ModelExtensions(seed_rawmodel, alternative_rawmodels,
         raise ValueError("Invalid extension_strategy")
 
     # make sure the seed has a scorer
-    if not (modelling.IsBackboneScorerSet(seed_rawmodel) and\
-            modelling.IsBackboneScorerEnvSet(seed_rawmodel)):
-        modelling.SetupDefaultBackboneScorer(seed_rawmodel)
+    if not (IsBackboneScorerSet(seed_rawmodel) and\
+            IsBackboneScorerEnvSet(seed_rawmodel)):
+        SetupDefaultBackboneScorer(seed_rawmodel)
     
+    # assign default weights if not given
+    if scoring_weights == None:
+        scoring_weights = dict()
+        scoring_weights["hbond"] = 1
+        scoring_weights["clash"] = 0.08880086
+        scoring_weights["torsion"] = 0.69712072
+        scoring_weights["reduced"] = 2.13535565
+        scoring_weights["cbeta"] = 3.59177335
+        scoring_weights["cb_packing"] = 1.17280667
+
     # all chains get processed consecutively
     for ch_idx, ch_seqres in enumerate(seed_rawmodel.seqres):    
         
diff --git a/modelling/pymod/_pipeline.py b/modelling/pymod/_pipeline.py
index 9871d1a5..efc8c81c 100644
--- a/modelling/pymod/_pipeline.py
+++ b/modelling/pymod/_pipeline.py
@@ -8,6 +8,7 @@ from promod3 import loop, sidechain, core
 from _modelling import *
 from _closegaps import *
 from _ring_punches import *
+from _mtm import *
 # external
 import ost
 from ost import mol, conop
@@ -161,98 +162,6 @@ 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):
@@ -497,7 +406,8 @@ def CheckFinalModel(mhandle):
             ost.LogInfo("Stereo-chemical problem in sidechain " + \
                         "of residue " + str(res))
 
-def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
+def BuildFromRawModel(mhandle, alternative_mhandles = None, 
+                      use_amber_ff=False, extra_force_fields=list()):
     '''Build a model starting with a raw model (see :func:`BuildRawModel`).
 
     This function implements a recommended pipeline to generate complete models
@@ -559,10 +469,23 @@ def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     else:
         SetupDefaultBackboneScorer(mhandle)
 
+    # try to extend terminal coverage using alternative rawmodels
+    if alternative_mhandles != None:
+        # iterate extensions until nothing happens anymore
+        current_num_residues = mhandle.model.GetResidueCount()
+        new_num_residues = None
+        while current_num_residues != new_num_residues:
+          current_num_residues = mhandle.model.GetResidueCount()
+          ModelExtensions(mhandle, alternative_mhandles, 
+                          extension_strategy = "overlapping")
+          ModelExtensions(mhandle, alternative_mhandles,
+                          extension_strategy = "non-overlapping")
+          new_num_residues = mhandle.model.GetResidueCount()        
+
     # remove terminal gaps
     RemoveTerminalGaps(mhandle)
 
-    #close gaps
+    # close gaps
     CloseGaps(mhandle, merge_distance = merge_distance, 
               fragment_db = fragment_db, structure_db = structure_db, 
               torsion_sampler = torsion_sampler)
@@ -582,5 +505,5 @@ def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     return mhandle.model
 
 # these methods will be exported into module
-__all__ = ('BuildFromRawModel', 'CloseGaps','BuildSidechains', 
+__all__ = ('BuildFromRawModel', 'BuildSidechains', 
            'MinimizeModelEnergy', 'CheckFinalModel')
-- 
GitLab