diff --git a/sidechain/doc/CMakeLists.txt b/sidechain/doc/CMakeLists.txt
index b38d2f0a696ad677f23cc4361c77cf83ce9cfe90..09dbdca6be55cf02fa7774e38edd120219130931 100644
--- a/sidechain/doc/CMakeLists.txt
+++ b/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})
diff --git a/sidechain/doc/index.rst b/sidechain/doc/index.rst
index a1ba3806d3fff16687944b8a4067be5401680f96..69da3d9180dc58002e0e340b08ad2d7d482d9b89 100644
--- a/sidechain/doc/index.rst
+++ b/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.
diff --git a/sidechain/doc/subrotamer_optimizer.rst b/sidechain/doc/subrotamer_optimizer.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e4124466394f2bdc0248aa40056a590f5b556dcb
--- /dev/null
+++ b/sidechain/doc/subrotamer_optimizer.rst
@@ -0,0 +1,45 @@
+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`
diff --git a/sidechain/pymod/CMakeLists.txt b/sidechain/pymod/CMakeLists.txt
index afc814bac04e1dbe04614a30e51cf5be5ef65cbf..3a745384883693f4c46e412d53182e2549b84501 100644
--- a/sidechain/pymod/CMakeLists.txt
+++ b/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
)
diff --git a/sidechain/pymod/_reconstruct_sidechains.py b/sidechain/pymod/_reconstruct_sidechains.py
index eacf17c86804661c502f8340f3cba8df63b9f129..2b1601a98a380f3f293fc2c48f7e3dada105662e 100644
--- a/sidechain/pymod/_reconstruct_sidechains.py
+++ b/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,
diff --git a/sidechain/pymod/export_subrotamer_optimizer.cc b/sidechain/pymod/export_subrotamer_optimizer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1b168786700daca626ac00b535921e04c6808663
--- /dev/null
+++ b/sidechain/pymod/export_subrotamer_optimizer.cc
@@ -0,0 +1,47 @@
+#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));
+}
+
+
+
+
+
+
diff --git a/sidechain/pymod/wrap_sidechain.cc b/sidechain/pymod/wrap_sidechain.cc
index 0b5547b1c38bbd0f0fb4c18a053455a01bb8a580..5e3f80d7a2b0d6834f4c742a24392b3a15359183 100644
--- a/sidechain/pymod/wrap_sidechain.cc
+++ b/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();
}
diff --git a/sidechain/src/CMakeLists.txt b/sidechain/src/CMakeLists.txt
index 8175397ca8ba5a012fd47888ec739fd865987f6b..34a9a875640fc05ecd20e24c7bb0067f078de118 100644
--- a/sidechain/src/CMakeLists.txt
+++ b/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}
diff --git a/sidechain/src/subrotamer_optimizer.cc b/sidechain/src/subrotamer_optimizer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..64067573a88c06f741216c3faf6ef8f7254bd788
--- /dev/null
+++ b/sidechain/src/subrotamer_optimizer.cc
@@ -0,0 +1,54 @@
+#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
diff --git a/sidechain/src/subrotamer_optimizer.hh b/sidechain/src/subrotamer_optimizer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..594d3ffd7aa87b07d95569aa2729b6b27022c28b
--- /dev/null
+++ b/sidechain/src/subrotamer_optimizer.hh
@@ -0,0 +1,18 @@
+#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