From 81485e55096770b8bb16c88d93f5aa669a86f12a Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Sat, 19 Aug 2017 17:25:05 +0200
Subject: [PATCH] Properly handle keep_sidechains flag in case of disulfid
bonds
The ResolveCysteins function is used to optimize the disulfid bond
network. The function uses RotamerGroupPtr as input.
To allow fixed cysteins also being part of the optimization,
the rotamer groups get newly constructed containing only one
rotamer with the same SG position as in the input structure.
---
modelling/pymod/_reconstruct_sidechains.py | 40 ++++++++++++++++++----
modelling/src/sidechain_reconstructor.cc | 40 ++++++++++++++++++----
modelling/src/sidechain_reconstructor.hh | 1 -
3 files changed, 66 insertions(+), 15 deletions(-)
diff --git a/modelling/pymod/_reconstruct_sidechains.py b/modelling/pymod/_reconstruct_sidechains.py
index aaf1a0eb..26d31567 100644
--- a/modelling/pymod/_reconstruct_sidechains.py
+++ b/modelling/pymod/_reconstruct_sidechains.py
@@ -297,8 +297,9 @@ def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_constructor,
return rotamer_groups, residues_with_rotamer_group
-def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_library,
- use_frm, bbdep, rotamer_constructor, phi_angles, psi_angles):
+def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
+ res_list, rotamer_library, use_frm, bbdep,
+ rotamer_constructor, phi_angles, psi_angles):
'''Get disulfid bridges for CYS and according rotamers.
CYS are identified by by items in cystein_indices (into res_list).
Returns: disulfid_indices: list of res. index in bridge,
@@ -322,10 +323,34 @@ def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_libra
continue
cys_ca_positions.append(ca.GetPos())
cys_cb_positions.append(cb.GetPos())
- # get RotamerGroup in case of disulfids, we do FRM in any case
- rot_group = _GetRotamerGroup(r.handle, sidechain.CYD, i, rotamer_library,
- rotamer_constructor, phi_angles[i],
- psi_angles[i], True, bbdep)
+
+ # get RotamerGroup in case of cysteins, we do FRM in any case...
+ # If we're suposed to keep the existing sidechains and the processed
+ # cystein contains all required atoms, we manually construct an FRM
+ # rotamer to still allow the rotamer to enter the disulfid bridge
+ # resolving algorithm.
+ rot_group = None
+ if keep_sidechains:
+ sg = r.FindAtom("SG")
+ if sg.IsValid():
+ particle = sidechain.Particle(sidechain.SidechainParticle.SParticle,
+ sg.GetPos(),
+ -0.19, "SG")
+ particle_list = [particle]
+ # The temperature and self internal_e_prefactor parameter have
+ # been copied from the SCWRLRotamerConstructor. Hacky...
+ frm_rotamer = sidechain.FRMRotamer(particle_list,
+ 1.69, 1.0, 4.07)
+ frm_rotamer.AddSubrotamerDefinition([0])
+ frm_rotamer.SetInternalEnergy(0.0)
+ rot_group = sidechain.FRMRotamerGroup([frm_rotamer], i)
+
+ if rot_group == None:
+ rot_group = _GetRotamerGroup(r.handle, sidechain.CYD, i,
+ rotamer_library,
+ rotamer_constructor,
+ phi_angles[i],
+ psi_angles[i], True, bbdep)
rot_group.AddFrameEnergy(frame)
cystein_rot_groups.append(rot_group)
indices_with_rot_groups.append(i)
@@ -455,7 +480,8 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
if len(cystein_indices) > 0:
# get disulfid bridges and according rotamers
disulfid_indices, disulfid_rotamers = \
- _GetDisulfidBridges(frame_residues, cystein_indices, prot.residues,
+ _GetDisulfidBridges(frame_residues, keep_sidechains,
+ cystein_indices, prot.residues,
rotamer_library, use_frm, bbdep,
rotamer_constructor, phi_angles, psi_angles)
# update frame_residues, incomplete_sidechains
diff --git a/modelling/src/sidechain_reconstructor.cc b/modelling/src/sidechain_reconstructor.cc
index 8ab263f3..5db32178 100644
--- a/modelling/src/sidechain_reconstructor.cc
+++ b/modelling/src/sidechain_reconstructor.cc
@@ -201,17 +201,19 @@ void SidechainReconstructor::SetEnvData_(loop::AllAtomEnv& env) {
attached_environment_ = true;
}
-// NOTE: it's ok to have private template functions...
-template<typename RotamerGroup> void SidechainReconstructor::BuildDisulfids_(
+void SidechainReconstructor::BuildDisulfids_(
SidechainReconstructionDataPtr res,
std::vector<uint>& cys_indices,
std::vector<sidechain::FrameResiduePtr>& frame_residues,
std::vector<bool>& has_sidechain) const {
// setup
std::vector<uint>& res_indices = res->env_pos->res_indices;
- typedef boost::shared_ptr<RotamerGroup> RotamerGroupPtr;
- typedef typename RotamerGroup::value_type RotamerPtr;
+ typedef sidechain::FRMRotamer Rotamer;
+ typedef sidechain::FRMRotamerPtr RotamerPtr;
+ typedef sidechain::FRMRotamerGroup RotamerGroup;
+ typedef sidechain::FRMRotamerGroupPtr RotamerGroupPtr;
+
loop::AllAtomPositions& out_pos = *(res->env_pos->all_pos);
// collect rotamers and data required to resolve the disulfid bridges
@@ -232,6 +234,32 @@ template<typename RotamerGroup> void SidechainReconstructor::BuildDisulfids_(
RotamerGroupPtr p;
env_->GetCydRotamerGroup(p, global_idx);
if (p) {
+ if(keep_sidechains_) {
+ sidechain::FrameResiduePtr f_res = env_->GetScFrameResidue(global_idx);
+ if(f_res) {
+ // We want to keep existing sidechains (keep_sidechains_ == true)
+ // and everything seems to be present since the sidechain frame
+ // residue exists. We therefore replace the extracted FRMRotamerGroup
+ // (p) with position data from the frame.
+ // The reason for still creating a FRMRotamerGroup is that it can
+ // enter the ResolveCystein evaluation.
+ std::vector<sidechain::Particle> particles;
+ particles.push_back(sidechain::Particle((*f_res)[0].GetParticleType(),
+ (*f_res)[0].GetPos(),
+ (*f_res)[0].GetCharge(),
+ (*f_res)[0].GetName()));
+ RotamerPtr r(new Rotamer(particles, (*p)[0]->GetTemperature(), 1.0,
+ (*p)[0]->GetInternalEnergyPrefactor()));
+ std::vector<int> subrotamer_definition(1,0);
+ r->AddSubrotamerDefinition(subrotamer_definition);
+ r->SetInternalEnergy(0.0);
+ std::vector<RotamerPtr> rotamers;
+ rotamers.push_back(r);
+ RotamerGroupPtr new_rot_group(new RotamerGroup(rotamers, global_idx));
+ p = new_rot_group;
+ }
+ }
+
p->SetFrameEnergy(frame);
cys_rotamer_groups.push_back(p);
ca_pos.push_back(all_pos_->GetPos(global_idx, loop::BB_CA_INDEX));
@@ -380,9 +408,7 @@ void SidechainReconstructor::SolveSystem_(
// handle cysteins
res->disulfid_bridges.clear();
if (build_disulfids_) {
- // we use the FRM rotamer model in any case for reconstructing disulfids
- BuildDisulfids_<sidechain::FRMRotamerGroup>(res, cys_indices,
- frame_residues, has_sidechain);
+ BuildDisulfids_(res, cys_indices, frame_residues, has_sidechain);
}
// collect rotamers
diff --git a/modelling/src/sidechain_reconstructor.hh b/modelling/src/sidechain_reconstructor.hh
index 565b830a..94e05496 100644
--- a/modelling/src/sidechain_reconstructor.hh
+++ b/modelling/src/sidechain_reconstructor.hh
@@ -84,7 +84,6 @@ private:
// -> pairs of cys_indices[i] for bridges added to res->disulfid_bridges
// -> unhandled CYS with known SC are added to frame_residues
// -> has_sidechain is updated as needed (i.e. when SC has been dealt with)
- template<typename RotamerGroup>
void BuildDisulfids_(SidechainReconstructionDataPtr res,
std::vector<uint>& cys_indices,
std::vector<sidechain::FrameResiduePtr>& frame_residues,
--
GitLab