From 90c2fa894eb6ad0600acce449788bc585ccd0b03 Mon Sep 17 00:00:00 2001
From: Gerardo Tauriello <gerardo.tauriello@unibas.ch>
Date: Thu, 11 Feb 2016 17:58:11 +0100
Subject: [PATCH] Added check to disallow merging of insertions for
 frame-approach-scorer

---
 modelling/doc/index.rst           | 16 +++++++++++-
 modelling/pymod/_closegaps.py     | 19 +++++++++++---
 modelling/pymod/export_model.cc   | 16 +++++++++---
 modelling/src/model.cc            | 30 +++++++++++++++++++++
 modelling/src/model.hh            |  4 +++
 modelling/tests/test_modelling.py | 43 +++++++++++++++++++++++++++++++
 6 files changed, 120 insertions(+), 8 deletions(-)

diff --git a/modelling/doc/index.rst b/modelling/doc/index.rst
index a0e0a7b5..b7ae5aba 100644
--- a/modelling/doc/index.rst
+++ b/modelling/doc/index.rst
@@ -106,9 +106,23 @@ Closing Gaps
 
 .. autofunction:: CloseLargeDeletions
 
+.. function:: CountEnclosedGaps(mhandle, gap)
+              CountEnclosedInsertions(mhandle, gap)
+
+  Counts all gaps from `mhandle` which are fully enclosed by given `gap`.
+  This is either all gaps or only insertions.
+
+  :param mhandle: Modelling handle on which to apply change.
+  :type mhandle:  :class:`ModellingHandle`
+  :param gap:     Gap defining range in which gaps are to be removed.
+  :type gap:      :class:`StructuralGap`
+
+  :return: Number of gaps.
+  :rtype:  :class:`int`
+
 .. function:: ClearGaps(mhandle, gap)
 
-  Removes all gaps from mhandle which are fully enclosed by given gap.
+  Removes all gaps from `mhandle` which are fully enclosed by given `gap`.
 
   :param mhandle: Modelling handle on which to apply change.
   :type mhandle:  :class:`ModellingHandle`
diff --git a/modelling/pymod/_closegaps.py b/modelling/pymod/_closegaps.py
index 8cbf086b..55f3ba29 100644
--- a/modelling/pymod/_closegaps.py
+++ b/modelling/pymod/_closegaps.py
@@ -502,6 +502,8 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
         
         # keep copy of original gap
         gap_orig = mhandle.gaps[gap_idx].Copy()
+        if score_variant == 0:
+            gap_ins = CountEnclosedInsertions(mhandle, gap_orig)
 
         # terminal gaps are no loops, database needs 2 stems to work on
         if gap_orig.IsTerminal():
@@ -525,6 +527,9 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
                                           max_db_loop_len)
         # iteratively extend actual_gap until we have enough loops
         while actual_gap.length <= max_db_loop_len:
+            # check if we have enough
+            if found_loops >= max_loops_to_search:
+                break
             # extend the gap
             if not actual_extender.Extend():
                 break
@@ -532,10 +537,10 @@ def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db,
             if not actual_gap.before.IsValid() or \
                not actual_gap.after.IsValid():
                 break
-
-            # check if we have enough
-            if found_loops >= max_loops_to_search:
-                break
+            # skip if we try to merge insertions with score_variant=0
+            if score_variant == 0 and \
+               CountEnclosedInsertions(mhandle, actual_gap) != gap_ins:
+                continue
 
             # get candidates for the current loop
             candidates = loop.LoopCandidates.FillFromDatabase(
@@ -643,6 +648,8 @@ def FillLoopsByMonteCarlo(mhandle, scorer, torsion_sampler, max_loops_to_search=
 
         # keep copy of original gap
         gap_orig = mhandle.gaps[gap_idx].Copy()
+        if score_variant == 0:
+            gap_ins = CountEnclosedInsertions(mhandle, gap_orig)
 
         # ignore terminal gaps
         if gap_orig.IsTerminal():
@@ -678,6 +685,10 @@ def FillLoopsByMonteCarlo(mhandle, scorer, torsion_sampler, max_loops_to_search=
             if not actual_gap.before.IsValid() or \
                not actual_gap.after.IsValid():
                 break
+            # skip if we try to merge insertions with score_variant=0
+            if score_variant == 0 and \
+               CountEnclosedInsertions(mhandle, actual_gap) != gap_ins:
+                continue
 
             # setup sampler, closer, scorer and cooler for MC
             try:
diff --git a/modelling/pymod/export_model.cc b/modelling/pymod/export_model.cc
index 28e3b743..5f457a27 100644
--- a/modelling/pymod/export_model.cc
+++ b/modelling/pymod/export_model.cc
@@ -12,6 +12,14 @@ ModellingHandle (*BuildRawModelHandleList)(const ost::seq::AlignmentList&,
                                            bool,
                                            const String&)=&BuildRawModel;
 
+namespace {
+int WrapCountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap) {
+  return CountEnclosedGaps(mhandle, gap, false);
+}
+int WrapCountEnclosedIns(ModellingHandle& mhandle, const StructuralGap& gap) {
+  return CountEnclosedGaps(mhandle, gap, true);
+}
+} // anon ns
 
 void export_model()
 {
@@ -20,9 +28,11 @@ void export_model()
     .def_readwrite("gaps", &ModellingHandle::gaps)
     .def_readwrite("seqres", &ModellingHandle::seqres)
   ;
-  def("ClearGaps",ClearGaps,(arg("gap")));
-  def("MergeGaps",MergeGaps,(arg("index")));
-  def("RemoveTerminalGaps",RemoveTerminalGaps);
+  def("ClearGaps", ClearGaps, (arg("mhandle"),arg("gap")));
+  def("CountEnclosedGaps", WrapCountEnclosedGaps, (arg("mhandle"),arg("gap")));
+  def("CountEnclosedInsertions", WrapCountEnclosedIns, (arg("mhandle"),arg("gap")));
+  def("MergeGaps", MergeGaps, (arg("mhandle"),arg("index")));
+  def("RemoveTerminalGaps", RemoveTerminalGaps, (arg("mhandle")));
   def("BuildRawModel", BuildRawModelHandle, 
       (arg("aln"),
        arg("include_ligands")=false,
diff --git a/modelling/src/model.cc b/modelling/src/model.cc
index c0965445..4c99d185 100644
--- a/modelling/src/model.cc
+++ b/modelling/src/model.cc
@@ -48,6 +48,36 @@ bool CheckBackboneAtoms(ResidueView res)
 
 }
 
+int CountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap,
+                      bool insertions_only)
+{
+  int num_gaps = 0;
+  
+  // extract info for this gap
+  String my_chain_name = gap.GetChainName();
+  ResNum my_b_res(0);
+  if (!gap.IsNTerminal()) my_b_res = gap.before.GetNumber();
+  ResNum my_a_res = my_b_res + gap.GetLength() + 1;
+
+  for (StructuralGapList::iterator i = mhandle.gaps.begin();
+       i != mhandle.gaps.end(); ++i) {
+    
+    // skip different chains
+    if (i->GetChainName() != my_chain_name) continue;
+
+    // get res_nums
+    ResNum b_res(0);
+    if (!i->IsNTerminal()) b_res = i->before.GetNumber();
+    ResNum a_res = b_res + i->GetLength() + 1;
+
+    // check full overlap
+    if (b_res >= my_b_res && a_res <= my_a_res) {
+      if (!insertions_only || i->GetLength() > 0) ++num_gaps;
+    }
+  }
+  return num_gaps;
+}
+
 int ClearGaps(ModellingHandle& mhandle, const StructuralGap& gap) {
   
   // return -1 if no more gaps
diff --git a/modelling/src/model.hh b/modelling/src/model.hh
index 8e7be7d2..6ee3516c 100644
--- a/modelling/src/model.hh
+++ b/modelling/src/model.hh
@@ -32,6 +32,10 @@ struct ModellingHandle {
   ost::seq::SequenceList seqres;
 };
 
+// see Python doc
+int CountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap,
+                      bool insertions_only=false);
+
 // see Python doc
 int ClearGaps(ModellingHandle& mhandle, const StructuralGap& gap);
 
diff --git a/modelling/tests/test_modelling.py b/modelling/tests/test_modelling.py
index 909fc2db..cfe41f7c 100644
--- a/modelling/tests/test_modelling.py
+++ b/modelling/tests/test_modelling.py
@@ -129,6 +129,49 @@ class ModellingTests(unittest.TestCase):
         with self.assertRaises(RuntimeError):
             modelling.MergeGaps(mhandle, 0)
 
+    def testCountEnclosedGaps(self):
+        tpl = io.LoadPDB('data/2dbs.pdb')
+        aln = seq.CreateAlignment(
+            seq.CreateSequence('trg', 'GATLNGFTVPAGNTLVLN---PD--KG---ATVTMAGA'),
+            seq.CreateSequence('tpl', '--NGG--TL--LI--PNGTYHFLGIQMKSNVHIRVE--'))
+        aln.AttachView(1, tpl.CreateFullView())
+        mhandle = modelling.BuildRawModel(aln)
+        # check
+        self.assertEqual(len(mhandle.gaps), 8)
+        for g in mhandle.gaps:
+          self.assertEqual(modelling.CountEnclosedGaps(mhandle, g), 1)
+          self.assertEqual(modelling.CountEnclosedInsertions(mhandle, g),
+                           1 if (g.length > 0) else 0)
+
+        mych = mhandle.model.chains[0]
+        # none
+        mygap = modelling.StructuralGap(mych.FindResidue(15),
+                                        mych.FindResidue(17), "L")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 0)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 0)
+        # extended singles
+        mygap = modelling.StructuralGap(mych.FindResidue(3),
+                                        mych.FindResidue(9), "LNGFT")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 1)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 1)
+        mygap = modelling.StructuralGap(mych.FindResidue(20),
+                                        mych.FindResidue(22), "K")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 1)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 0)
+        # doubles
+        mygap = modelling.StructuralGap(mych.FindResidue(3),
+                                        mych.FindResidue(12), "LNGFTVPA")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 2)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 2)
+        mygap = modelling.StructuralGap(mych.FindResidue(13),
+                                        mych.FindResidue(19), "TLVLN")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 2)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 1)
+        mygap = modelling.StructuralGap(mych.FindResidue(20),
+                                        mych.FindResidue(25), "KGAT")
+        self.assertEqual(modelling.CountEnclosedGaps(mhandle, mygap), 2)
+        self.assertEqual(modelling.CountEnclosedInsertions(mhandle, mygap), 0)
+
     def testClearGapsExceptions(self):
         tpl = io.LoadPDB('data/1mcg.pdb')
         aln = seq.CreateAlignment(seq.CreateSequence('trg', 'DDFAGTHN'),
-- 
GitLab