implement subrotamer optimizer in C

sidechain/doc/CMakeLists.txt

@@ -9,6 +9,7 @@ rotamer_constructor.rst
 graph.rst
 disulfid.rst
 loading.rst
+subrotamer_optimizer.rst
 )
 
 add_doc_source(NAME sidechain RST ${SIDECHAIN_RST})

sidechain/doc/index.rst

@@ -44,6 +44,7 @@ Contents:
    graph
    disulfid
    loading
+   subrotamer_optimizer
 
 
 .. [krivov2009] Krivov GG, Shapovalov MV and Dunbrack RL Jr. (2009). Improved prediction of protein side-chain conformations with SCWRL4. Proteins.

sidechain/doc/subrotamer_optimizer.rst

+Subrotamer Optimization
+================================================================================
+
+.. currentmodule:: promod3.sidechain
+
+The idea of the flexible rotamer model is to have several subrotamers
+representing one single rotamer. 
+One subrotamer is the representative (the active subrotamer).
+Thats the one that gets inserted in the structure once the reconstruction has
+finished. This is not necessarily the optimal one. The SubrotamerOptimizer
+takes a list of flexible rotamers, converts every single flexible rotamer 
+in a rotamer group of rigid rotamers, solves the sidechain reconstruction
+problem and assigns the optimal subrotamers as active in the original 
+flexible rotamers. In terms of energies, the SubrotamerOptimizer expects
+the energy of the flexible rotamers towards the rigid frame already to be set.
+The internal energies can be controlled as parameters and a constant value is
+set for rotamers that are currently active or not. The reason for that is
+that you might want to prefer the currently active subrotamers if no
+pairwise energies to other rotamers are observed. 
+
+.. method:: SubrotamerOptimizer(rotamers, [active_internal_energy=-2.0, \
+                                inactive_internal_energy=0.0, \
+                                max_complexity=100000000, initial_epsilon=0.02)
+
+  Takes the **rotamers** of type :class:`FRMRotamer`, converts all their 
+  subrotamers into rigid rotamers, solves the sidechain reconstruction problem
+  and assigns the ideal subrotamers as active in the input **rotamers**.
+
+  :param rotamers:      The rotamers to be optimized               
+  :type rotamers:       :class:`list` of :class:`FRMRotamer`
+
+  :param active_internal_energy: Internal energy that gets assigned to all
+                                 currently active subrotamers
+  :type active_internal_energy:  :class:`float`
+
+  :param inactive_internal_energy: Internal energy that gets assigned to all
+                                   currently inactive subrotamers
+  :type inactive_internal_energy: :class:`float`
+
+  :param max_complexity: Max complexity of the internal :class:`RotamerGraph`
+  :type max_complexity: :class:`int`
+
+  :param initial_epsilon: Epsilon value controlling the pruning of 
+                          internal :class:`RotamerGraph`
+  :type initial_epsilon: :class:`float`

sidechain/pymod/CMakeLists.txt

@@ -12,6 +12,7 @@ export_rotamer_lib.cc
 export_sidechain_object_loader.cc
 export_sidechain_reconstructor.cc
 export_scwrl_rotamer_constructor.cc
+export_subrotamer_optimizer.cc
 wrap_sidechain.cc
 )
 

sidechain/pymod/_reconstruct_sidechains.py

@@ -350,35 +350,6 @@ def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_libra
     return disulfid_indices, disulfid_rotamers
 
 
-def RefineFRMRotamerGroups(solution, frm_rotamer_groups,
-                           graph_max_complexity=100000000 , 
-                           graph_initial_epsilon=0.02):
-    
-    rrm_rotamer_groups = list()
-    for i, rg in enumerate(frm_rotamer_groups):
-        frm_rotamer = rg[solution[i]]
-        rrm_rotamers = list()
-        original_active_subrotamer = frm_rotamer.GetActiveSubrotamer()
-        for j in range(frm_rotamer.GetNumSubrotamers()):
-            frm_rotamer.SetActiveSubrotamer(j)
-            new_rrm_rotamer = frm_rotamer.ToRRMRotamer()
-            if(j != original_active_subrotamer):
-                new_rrm_rotamer.SetInternalEnergy(0.0)
-            else:
-                new_rrm_rotamer.SetInternalEnergy(-0.5)
-            rrm_rotamers.append(new_rrm_rotamer)
-        new_rotamer_group = sidechain.RRMRotamerGroup(rrm_rotamers, 0)
-        rrm_rotamer_groups.append(new_rotamer_group)
-    graph = sidechain.RotamerGraph.CreateFromRRMList(rrm_rotamer_groups)
-
-    solution = graph.TreeSolve(max_complexity=graph_max_complexity,
-                              initial_epsilon=graph_initial_epsilon)[0]
-
-    return (solution, rrm_rotamer_groups)
-
-
-
-
 ###############################################################################
 
 def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
@@ -417,16 +388,10 @@ def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
     :param optimize_subrotamers: Only considered when **rotamer_model** 
                                  is "frm".
                                  If set to True, the FRM solution undergoes 
-                                 some postprocessing. The final FRM rotamers get 
-                                 turned into :class:`RRMRotamerGroup` objects 
-                                 and fed into a second run of graph solving to 
-                                 find the optimal subrotamers from the FRM 
-                                 model. This mainly improves the reconstruction 
-                                 performance of bulky sidechains such as 
-                                 PHE/TYR/TRP. Internal energies of the 
-                                 :class:`RRMRotamer` objects are set to 0.0,
-                                 -0.5 if they represent the active subrotamer 
-                                 in the :class:`FRMRotamer`.
+                                 some postprocessing by calling 
+                                 :func:`SubrotamerOptimizer` with default 
+                                 parametrization.
+    :type optimize_subrotamers: :class:`bool`
 
     :param graph_max_complexity: Max. complexity for
                                  :meth:`RotamerGraph.TreeSolve`.
@@ -520,10 +485,11 @@ def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
                                initial_epsilon=graph_initial_epsilon)[0]
 
     if use_frm and optimize_subrotamers:
-        solution, rotamer_groups = RefineFRMRotamerGroups(solution, 
-                                                          rotamer_groups,
-                                                          graph_max_complexity,
-                                                          graph_initial_epsilon)
+        rotamers_to_optimize = list()
+
+        for rot_group, sol in zip(rotamer_groups, solution): 
+            rotamers_to_optimize.append(rot_group[sol])
+        sidechain.SubrotamerOptimizer(rotamers_to_optimize)
 
     # update structure
     for i, rot_group, sol in zip(residues_with_rotamer_group, rotamer_groups,

sidechain/pymod/export_subrotamer_optimizer.cc

+#include <boost/python.hpp>
+#include <boost/python/iterator.hpp>
+#include <boost/python/register_ptr_to_python.hpp>
+
+#include <promod3/core/export_helper.hh>
+#include <promod3/sidechain/subrotamer_optimizer.hh>
+
+
+using namespace promod3;
+using namespace promod3::sidechain;
+using namespace boost::python;
+
+
+namespace{
+
+void WrapOptimizer(boost::python::list& rotamers,
+                   Real active_internal_energy,
+                   Real inactive_internal_energy,
+                   uint max_complexity,
+                   Real initial_epsilon) {
+
+  std::vector<FRMRotamerPtr> v_rotamers;
+  promod3::core::ConvertListToVector(rotamers, v_rotamers);
+
+  promod3::sidechain::SubrotamerOptimizer(v_rotamers, 
+                                          active_internal_energy,
+                                          inactive_internal_energy, 
+                                          max_complexity,
+                                          initial_epsilon);
+}
+
+}
+
+void export_SubrotamerOptimizer() {
+
+  def("SubrotamerOptimizer", &WrapOptimizer,(arg("rotamers"),
+                                             arg("active_internal_energy")=-2.0,
+                                             arg("inactive_internal_energy")=0.0,
+                                             arg("max_complexity")=100000000,
+                                             arg("initial_epsilon")=0.02));
+}
+
+
+
+
+
+

sidechain/pymod/wrap_sidechain.cc

@@ -13,6 +13,7 @@ void export_RotamerLib();
 void export_SidechainObjectLoader();
 void export_SidechainReconstructor();
 void export_SCWRLRotamerConstructor();
+void export_SubrotamerOptimizer();
 
 using namespace boost::python;
 
@@ -31,4 +32,5 @@ BOOST_PYTHON_MODULE(_sidechain)
   export_SidechainObjectLoader();
   export_SidechainReconstructor();
   export_SCWRLRotamerConstructor();
+  export_SubrotamerOptimizer();
 }

sidechain/src/CMakeLists.txt

@@ -17,6 +17,7 @@ set(SIDECHAIN_HEADERS
   sidechain_object_loader.hh
   sidechain_reconstructor.hh
   scwrl_rotamer_constructor.hh
+  subrotamer_optimizer.hh
 )
 
 set(SIDECHAIN_SOURCES
@@ -38,6 +39,7 @@ set(SIDECHAIN_SOURCES
   sidechain_object_loader.cc
   sidechain_reconstructor.cc
   scwrl_rotamer_constructor.cc
+  subrotamer_optimizer.cc
 )
 
 module(NAME sidechain HEADERS ${SIDECHAIN_HEADERS} SOURCES ${SIDECHAIN_SOURCES}

sidechain/src/subrotamer_optimizer.cc

+#include <promod3/sidechain/subrotamer_optimizer.hh>
+
+namespace promod3{ namespace sidechain{
+
+void SubrotamerOptimizer(std::vector<FRMRotamerPtr>& rotamers,
+                         Real active_internal_energy,
+                         Real inactive_internal_energy,
+                         uint max_complexity,
+                         Real initial_epsilon) {
+
+  // buildup a list of RRMRotamer groups
+  std::vector<RRMRotamerGroupPtr> rrm_rotamer_groups;
+
+  for(uint rot_idx = 0; rot_idx < rotamers.size(); ++rot_idx) {
+    
+    // extract all subrotamers of current FRMRotamer as RRMRotamers
+    FRMRotamerPtr frm_rotamer = rotamers[rot_idx];
+    std::vector<RRMRotamerPtr> rrm_rotamers;
+    uint active_idx = frm_rotamer->GetActiveSubrotamer();
+
+    for(uint subrot_idx = 0; subrot_idx < frm_rotamer->subrotamer_size(); 
+        ++subrot_idx) {
+      
+      frm_rotamer->SetActiveSubrotamer(subrot_idx);
+      RRMRotamerPtr new_rrm_rotamer = frm_rotamer->ToRRMRotamer();
+
+      if(subrot_idx == active_idx) {
+        new_rrm_rotamer->SetInternalEnergy(active_internal_energy);
+      }   
+      else {
+        new_rrm_rotamer->SetInternalEnergy(inactive_internal_energy);
+      }
+      rrm_rotamers.push_back(new_rrm_rotamer);
+    }
+    
+    RRMRotamerGroupPtr rrm_rotamer_group(new RRMRotamerGroup(rrm_rotamers, 0));
+    rrm_rotamer_groups.push_back(rrm_rotamer_group);
+  }
+
+  // generate a graph from the RRMRotamerGroups to find the optimal subrotamers
+  RotamerGraphPtr graph = RotamerGraph::CreateFromRRMList(rrm_rotamer_groups);
+
+  std::pair<std::vector<int>, Real> full_solution = 
+  graph->TreeSolve(max_complexity, initial_epsilon);
+
+  const std::vector<int>& solution = full_solution.first;
+
+  // activate the appropriate subrotamers
+  for(uint rot_idx = 0; rot_idx < solution.size(); ++rot_idx) {
+     rotamers[rot_idx]->SetActiveSubrotamer(solution[rot_idx]); 
+  }
+}
+
+}} // ns

sidechain/src/subrotamer_optimizer.hh

+#ifndef PROMOD3_SUBROTAMER_OPTIMIZER_HH
+#define PROMOD3_SUBROTAMER_OPTIMIZER_HH
+
+#include <vector>
+
+#include <promod3/sidechain/rotamer_graph.hh>
+
+
+namespace promod3{ namespace sidechain{
+
+void SubrotamerOptimizer(std::vector<FRMRotamerPtr>& rotamers,
+                         Real active_internal_energy = -0.5,
+                         Real inactive_internal_energy = 0.0,
+                         uint max_complexity = 100000000,
+                         Real initial_epsilon = 0.02);
+}} // ns
+
+#endif