diff --git a/sidechain/doc/disulfid.rst b/sidechain/doc/disulfid.rst
index 6bf6078b9dc311a000be25d2fbf9c408d20aa781..5b330287975708bdf84fa008b20d2e63315f4895 100644
--- a/sidechain/doc/disulfid.rst
+++ b/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.
diff --git a/sidechain/pymod/_reconstruct_sidechains.py b/sidechain/pymod/_reconstruct_sidechains.py
index b9e13dbd689f1200263fc2a923ca38f04f798ccf..ac41d155fe72ce701872f10305f3b00641b2a9fa 100644
--- a/sidechain/pymod/_reconstruct_sidechains.py
+++ b/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()
diff --git a/sidechain/pymod/export_disulfid.cc b/sidechain/pymod/export_disulfid.cc
index b9000e9a21aac42c1a1810b0f3961f066ea6111b..72ea6c3c91febe9133ecdfde62ecb92951c63875 100644
--- a/sidechain/pymod/export_disulfid.cc
+++ b/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
diff --git a/sidechain/src/disulfid.cc b/sidechain/src/disulfid.cc
index f7837bafefa8782953b8a8bacb8260a49832d9fd..a0e1d3c07b4cc0b0c94994825d489f614be5f1fd 100644
--- a/sidechain/src/disulfid.cc
+++ b/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]);
+ }
+ }
}
diff --git a/sidechain/src/disulfid.hh b/sidechain/src/disulfid.hh
index 40613ca192d58f655d319de66b68803a2df30890..84f98aaed96b86342031e10b9f30ad1cafc4016d 100644
--- a/sidechain/src/disulfid.hh
+++ b/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);
diff --git a/sidechain/src/sidechain_reconstructor.cc b/sidechain/src/sidechain_reconstructor.cc
index 71cb3ccf7f3cd89d81bb23801ca0679afa579429..92ff3157cf4f71301ee5d6dac8af8f9d7f612c7a 100644
--- a/sidechain/src/sidechain_reconstructor.cc
+++ b/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);