Optimize subrotamers if Flexible rotamers participate in disulfid bonds

sidechain/doc/disulfid.rst

@@ -74,7 +74,7 @@ that has also been implemented here.
 
 
 .. method:: ResolveCysteins(rotamer_groups, ca_positions, cb_positions,\
-                            [score_threshold=45.0])
+                            [score_threshold=45.0, optimize_subrotamers=False])
 
   Tries to optimize disulfid bond network. In a first step, all disulfid bonds
   get detected using :func:`DisulfidScore`. If the value between two rotamers is
@@ -88,12 +88,20 @@ that has also been implemented here.
   :param cb_positions:   The CB positions of the according rotamers
   :param score_threshold: The score two rotamers must have to be considered
                           as a disulfid bond
+  :param optimize_subrotamers: If set to true and the input consists of flexible
+                               rotamer groups, the active subrotamers get 
+                               optimized. For every pair of rotamers
+                               participating in a disulfid bond, the subrotamers
+                               with best :func:`DisulfidScore` get activated in
+                               the input **rotamer_groups**. This has an effect
+                               when the rotamers get applied on residues.
 
   :type rotamer_groups:  :class:`list` of 
                          :class:`FRMRotamerGroup`/:class:`RRMRotamerGroup`
   :type ca_positions: :class:`list` of :class:`ost.geom.Vec3`
   :type cb_positions: :class:`list` of :class:`ost.geom.Vec3`
   :type score_threshold: :class:`float`
+  :type optimize_subrotamers: :class:`bool`
 
   :returns: A :class:`tuple` containing two :class:`list` objects with equal
             length. Both lists contain :class:`tuple` objects with two elements.

sidechain/pymod/_reconstruct_sidechains.py

@@ -207,23 +207,27 @@ def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
             incomplete_sidechains.append(idx)
 
 def _GetRotamerGroup(res_handle, rot_id, res_idx, rot_lib, rot_constructor,
-                     phi, psi, use_frm, bbdep):
+                     phi, psi, use_frm, bbdep, probability_cutoff = 0.98):
     if use_frm:
         if bbdep:
             return rot_constructor.ConstructFRMRotamerGroup(res_handle, rot_id,
                                                             res_idx, rot_lib,
-                                                            phi, psi)
+                                                            phi, psi,
+                                                            probability_cutoff)
         else:
             return rot_constructor.ConstructFRMRotamerGroup(res_handle, rot_id,
-                                                            res_idx, rot_lib)
+                                                            res_idx, rot_lib,
+                                                            probability_cutoff)
     else:
         if bbdep:
             return rot_constructor.ConstructRRMRotamerGroup(res_handle, rot_id,
                                                             res_idx, rot_lib,
-                                                            phi, psi)
+                                                            phi, psi,
+                                                            probability_cutoff)
         else:
             return rot_constructor.ConstructRRMRotamerGroup(res_handle, rot_id,
-                                                            res_idx, rot_lib)
+                                                            res_idx, rot_lib,
+                                                            probability_cutoff)
 
 def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_constructor,
                       phi_angles, psi_angles, use_frm, bbdep, frame_residues):
@@ -331,7 +335,8 @@ def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_libra
     bond_result, \
     rotamer_result = sidechain.ResolveCysteins(cystein_rot_groups, 
                                                cys_ca_positions,
-                                               cys_cb_positions, 45.0)
+                                               cys_cb_positions, 45.0,
+                                               True)
 
     # get CYS with disulfid bonds and the chosen rotamers
     disulfid_indices = list()

sidechain/pymod/export_disulfid.cc

@@ -31,7 +31,8 @@ Real WrapDisulfidTwo(FRMRotamerPtr rot_one, FRMRotamerPtr rot_two,
 boost::python::tuple WrapResolveCysteins(const boost::python::list& rot_groups,
                                          const boost::python::list& ca_pos,
                                          const boost::python::list& cb_pos,
-                                         Real score_threshold){
+                                         Real score_threshold, 
+                                         bool optimize_subrotamers){
 
   boost::python::list disulfid_return_list;
   boost::python::list rotamer_return_list;
@@ -53,13 +54,15 @@ boost::python::tuple WrapResolveCysteins(const boost::python::list& rot_groups,
   if(extractor.check()){
     std::vector<promod3::sidechain::FRMRotamerGroupPtr> v_rot_groups;
     promod3::core::ConvertListToVector(rot_groups, v_rot_groups);
-    ResolveCysteins(v_rot_groups, v_ca_pos, v_cb_pos, score_threshold, 
+    ResolveCysteins(v_rot_groups, v_ca_pos, v_cb_pos, score_threshold,
+                    optimize_subrotamers, 
                     disulfid_indices, rotamer_indices);
   }
   else{
     std::vector<promod3::sidechain::RRMRotamerGroupPtr> v_rot_groups;
     promod3::core::ConvertListToVector(rot_groups, v_rot_groups);
     ResolveCysteins(v_rot_groups, v_ca_pos, v_cb_pos, score_threshold,
+                    optimize_subrotamers,
                     disulfid_indices, rotamer_indices);
   }
 
@@ -86,15 +89,16 @@ void export_Disulfid(){
                                               arg("cb_pos_two"),
                                               arg("sg_pos_two")));
 
-  def("DisulfidScore",&WrapDisulfidOne,(arg("rotamer_one"),arg("rotamer_two"),
-                                        arg("ca_pos_one"),arg("cb_pos_one"),
-                                        arg("ca_pos_two"),arg("cb_pos_two")));
-  def("DisulfidScore",&WrapDisulfidTwo,(arg("rotamer_one"),arg("rotamer_two"),
-                                        arg("ca_pos_one"),arg("cb_pos_one"),
-                                        arg("ca_pos_two"),arg("cb_pos_two")));
+  def("DisulfidScore",&WrapDisulfidOne,(arg("rotamer_one"), arg("rotamer_two"),
+                                        arg("ca_pos_one"), arg("cb_pos_one"),
+                                        arg("ca_pos_two"), arg("cb_pos_two")));
+  def("DisulfidScore",&WrapDisulfidTwo,(arg("rotamer_one"), arg("rotamer_two"),
+                                        arg("ca_pos_one"), arg("cb_pos_one"),
+                                        arg("ca_pos_two"), arg("cb_pos_two")));
 
   def("ResolveCysteins", &WrapResolveCysteins, (arg("rotamer_groups"),
                                                 arg("ca_positions"),
                                                 arg("cb_positions"),
-                                                arg("score_threshold") = 45.0));
+                                                arg("score_threshold") = 45.0,
+                                                arg("optimize_subrotamers") = false));
 }
\ No newline at end of file

sidechain/src/disulfid.cc

@@ -190,7 +190,6 @@ void CysteineResolve(const std::vector<T>& rot_groups,
   }
 }
 
-
 }
 
 
@@ -225,6 +224,52 @@ Real DisulfidRawScore(const geom::Vec3& ca_pos_one,
 }
 
 
+namespace{
+
+// this is only in here because it requires the DisulfidRawScore function..
+void SetOptimalSubrotamer(promod3::sidechain::FRMRotamerPtr rot_one, 
+                          promod3::sidechain::FRMRotamerPtr rot_two,
+                          const geom::Vec3& ca_pos_one, 
+                          const geom::Vec3& cb_pos_one,
+                          const geom::Vec3& ca_pos_two, 
+                          const geom::Vec3& cb_pos_two) {
+
+  std::vector<geom::Vec3> sg_pos_one, sg_pos_two;
+  ExtractCYSSG(*rot_one, *rot_two, sg_pos_one, sg_pos_two);
+
+  if(sg_pos_one.size() != rot_one->subrotamer_size() || 
+     sg_pos_two.size() != rot_two->subrotamer_size() ){
+    throw promod3::Error("Expect cysteins to have exactly one gamma sulfur "
+                         "per subrotamer!");
+  }
+
+  Real min_score = std::numeric_limits<Real>::max();
+  uint min_i = 0;
+  uint min_j = 0;
+
+  for(uint i = 0; i < sg_pos_one.size(); ++i){
+    for(uint j = 0; j < sg_pos_two.size(); ++j){
+
+      Real score = DisulfidRawScore(ca_pos_one,cb_pos_one,
+                                    sg_pos_one[i], ca_pos_two,
+                                    cb_pos_two,sg_pos_two[j]);
+      if(score < min_score){
+        min_score = score;
+        min_i = i;
+        min_j = j;
+      }
+    }
+  }
+
+  if(min_score < std::numeric_limits<Real>::max()){
+    rot_one->SetActiveSubrotamer(min_i);
+    rot_two->SetActiveSubrotamer(min_j);
+  }
+}
+
+}
+
+
 Real DisulfidScore(RRMRotamerPtr rot_one, RRMRotamerPtr rot_two, 
                    const geom::Vec3& ca_pos_one, const geom::Vec3& cb_pos_one, 
                    const geom::Vec3& ca_pos_two, const geom::Vec3& cb_pos_two){
@@ -273,13 +318,13 @@ Real DisulfidScore(FRMRotamerPtr rot_one, FRMRotamerPtr rot_two,
   return min_raw_score + self_energy;
 }
 
-// template functionility of following functions has been hidden in
+// template functionality of following functions has been hidden in
 // the source file to avoid bloating the header...
 void ResolveCysteins(
       const std::vector<RRMRotamerGroupPtr>& rot_groups,
       const std::vector<geom::Vec3>& ca_pos,
       const std::vector<geom::Vec3>& cb_pos,
-      Real disulfid_score_thresh,
+      Real disulfid_score_thresh, bool optimize_subrotamers,
       std::vector<std::pair<uint, uint> >& disulfid_indices,
       std::vector<std::pair<uint, uint> >& rotamer_indices){
 
@@ -292,13 +337,29 @@ void ResolveCysteins(
 void ResolveCysteins(const std::vector<FRMRotamerGroupPtr>& rot_groups,
                      const std::vector<geom::Vec3>& ca_pos,
                      const std::vector<geom::Vec3>& cb_pos,
-                     Real disulfid_score_thresh,
+                     Real disulfid_score_thresh, bool optimize_rotamers,
                      std::vector<std::pair<uint, uint> >& disulfid_indices,
                      std::vector<std::pair<uint, uint> >& rotamer_indices){
 
   CysteineResolve<FRMRotamerGroupPtr>(rot_groups, ca_pos, cb_pos,
                                       disulfid_score_thresh, disulfid_indices,
                                       rotamer_indices);
+
+  if(optimize_rotamers){
+
+    for(uint i = 0; i < disulfid_indices.size(); ++i){
+
+      const std::pair<uint, uint>& group_pair = disulfid_indices[i];
+      const std::pair<uint, uint>& rotamer_pair = rotamer_indices[i];
+
+      SetOptimalSubrotamer((*rot_groups[group_pair.first])[rotamer_pair.first],
+                           (*rot_groups[group_pair.second])[rotamer_pair.second],
+                           ca_pos[group_pair.first], 
+                           cb_pos[group_pair.first],
+                           ca_pos[group_pair.second], 
+                           cb_pos[group_pair.second]);
+    }
+  }
 }
 
 

sidechain/src/disulfid.hh

@@ -28,14 +28,14 @@ Real DisulfidScore(FRMRotamerPtr rot_one, FRMRotamerPtr rot_two,
  void ResolveCysteins(const std::vector<FRMRotamerGroupPtr>& rot_groups,
                       const std::vector<geom::Vec3>& ca_pos,
                       const std::vector<geom::Vec3>& cb_pos,
-                      Real disulfid_score_thresh,
+                      Real disulfid_score_thresh, bool optimize_subrotamers,
                       std::vector<std::pair<uint, uint> >& disulfid_indices,
                       std::vector<std::pair<uint, uint> >& rotamer_indices);
 
  void ResolveCysteins(const std::vector<RRMRotamerGroupPtr>& rot_groups,
                       const std::vector<geom::Vec3>& ca_pos,
                       const std::vector<geom::Vec3>& cb_pos,
-                      Real disulfid_score_thresh,
+                      Real disulfid_score_thresh, bool optimize_subrotamers,
                       std::vector<std::pair<uint, uint> >& disulfid_indices,
                       std::vector<std::pair<uint, uint> >& rotamer_indices);
 

sidechain/src/sidechain_reconstructor.cc

@@ -241,7 +241,8 @@ template<typename RotamerGroup> void SidechainReconstructor::BuildDisulfids_(
   std::vector<std::pair<uint, uint> > rot_solution;
 
   ResolveCysteins(cys_rotamer_groups, ca_pos, cb_pos,
-                  disulfid_score_thresh_, bond_solution, rot_solution);
+                  disulfid_score_thresh_, true,
+                  bond_solution, rot_solution);
 
   std::vector<int> is_bonded(cys_rotamer_groups.size(), 0);