Add SCWRL3 scoring function for sidechain modelling

This should be an example of how to add a scoring function

actions/doc/index.rst

@@ -190,3 +190,8 @@ Several flags control the modelling behaviour:
 
   Dont do subrotamer optimization if flexible rotamer model is used
 
+.. option:: -f, --energy_function
+
+  The energy function to be used. Default is SCWRL4, can be any function
+  supported by :meth:`promod3.modelling.ReconstructSidechains`.
+

actions/pm-build-sidechains

@@ -54,6 +54,8 @@ parser.add_argument('-i', '--backbone-independent', action='store_true',
 parser.add_argument('-s', '--no-subrotamer-optimization', action='store_true',
                     help='Dont do subrotamer optimization if flexible ' +
                          'rotamer model is used.')
+parser.add_argument('-f', '--energy_function', action = 'store', 
+                    dest = "energy_function", default = "SCWRL4")
 
 opts = parser.Parse()
 
@@ -80,7 +82,8 @@ modelling.ReconstructSidechains(prot, keep_sidechains=opts.keep_sidechains,
                                 build_disulfids=opts.no_disulfids==False, 
                                 rotamer_model=rotamer_model, 
                                 rotamer_library=lib, 
-                                optimize_subrotamers=opt_subrot)
+                                optimize_subrotamers=opt_subrot,
+                                energy_function = opts.energy_function)
 
 # output
 ost.PopVerbosityLevel()

modelling/pymod/_reconstruct_sidechains.py

@@ -402,7 +402,8 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
                           rotamer_model="frm", consider_ligands=True, 
                           rotamer_library=None, optimize_subrotamers=True,
                           graph_max_complexity=100000000, 
-                          graph_initial_epsilon=0.02):
+                          graph_initial_epsilon=0.02, 
+                          energy_function = "SCWRL4"):
     '''Reconstruct sidechains for the given structure.
 
     :param ent: Structure for sidechain reconstruction. Note, that the sidechain
@@ -445,6 +446,9 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
     :param graph_intial_epsilon: Initial epsilon for
                                  :meth:`RotamerGraph.TreeSolve`.
     :type graph_intial_epsilon:  :class:`float`
+    :param energy_function:      What energy function to use can be any in 
+                                 ["SCWRL4", "SCWRL3"]
+    :type energy_function:       :class:`str`
     '''
     prof_name = 'modelling::ReconstructSidechains'
     prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
@@ -463,8 +467,16 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
     if type(rotamer_library) is sidechain.BBDepRotamerLib:
         bbdep = True
 
-    rotamer_constructor = sidechain.SCWRL4RotamerConstructor(False)
-    
+    rotamer_constructor = None
+
+    if energy_function == "SCWRL4":
+        rotamer_constructor = sidechain.SCWRL4RotamerConstructor(False)
+    if energy_function == "SCWRL3":
+        rotamer_constructor = sidechain.SCWRL3RotamerConstructor(False)
+
+    if rotamer_constructor == None:
+        raise RuntimeError("Invalid Energy function to reconstruct sidechains!")
+
     # take out ligand chain and any non-peptides
     prot = ent.Select("peptide=true and cname!='_'")
     

sidechain/pymod/export_rotamer_constructor.cc

@@ -21,6 +21,7 @@
 #include <promod3/sidechain/rotamer.hh>
 #include <promod3/sidechain/rotamer_group.hh>
 #include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl3_rotamer_constructor.hh>
 
 using namespace boost::python;
 using namespace promod3::sidechain;
@@ -34,6 +35,11 @@ SCWRL4RotamerConstructorPtr WrapSCWRL4RotamerConstructorInit(bool cb_in_sidechai
   return p;
 }
 
+SCWRL3RotamerConstructorPtr WrapSCWRL3RotamerConstructorInit(bool cb_in_sidechain) {
+  SCWRL3RotamerConstructorPtr p(new SCWRL3RotamerConstructor(cb_in_sidechain));
+  return p;
+}
+
 RRMRotamerGroupPtr WrapRRMGroup_res(RotamerConstructorPtr constructor,
                                     const ost::mol::ResidueHandle& res,
                                     RotamerID id,
@@ -444,4 +450,10 @@ void export_RotamerConstructor(){
                                                                     arg("comp_lib")))
   ;
 
+
+  class_<SCWRL3RotamerConstructor, SCWRL3RotamerConstructorPtr, 
+         bases<RotamerConstructor> >("SCWRL3RotamerConstructor", no_init)
+    .def("__init__", boost::python::make_constructor(&WrapSCWRL3RotamerConstructorInit))
+  ;
+
 }

sidechain/src/CMakeLists.txt

@@ -7,6 +7,7 @@ set(SIDECHAIN_HEADERS
   particle_scoring_base.hh
   particle_scoring.hh
   scwrl4_particle_scoring.hh
+  scwrl3_particle_scoring.hh
   rotamer.hh
   rotamer_cruncher.hh
   rotamer_density.hh
@@ -20,6 +21,7 @@ set(SIDECHAIN_HEADERS
   rotamer_constructor.hh
   rotamer_lookup.hh
   scwrl4_rotamer_constructor.hh
+  scwrl3_rotamer_constructor.hh
   subrotamer_optimizer.hh
 )
 
@@ -29,6 +31,7 @@ set(SIDECHAIN_SOURCES
   frame.cc
   particle_scoring.cc
   scwrl4_particle_scoring.cc
+  scwrl3_particle_scoring.cc
   rotamer.cc
   rotamer_cruncher.cc
   rotamer_density.cc
@@ -42,6 +45,7 @@ set(SIDECHAIN_SOURCES
   rotamer_constructor.cc
   rotamer_lookup.cc
   scwrl4_rotamer_constructor.cc
+  scwrl3_rotamer_constructor.cc
   subrotamer_optimizer.cc
 )
 

sidechain/src/particle_scoring.cc

@@ -16,6 +16,7 @@
 #include <promod3/core/message.hh>
 #include <promod3/sidechain/particle_scoring.hh>
 #include <promod3/sidechain/scwrl4_particle_scoring.hh>
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
 
 namespace promod3{ namespace sidechain{
 
@@ -31,6 +32,10 @@ Real PairwiseParticleScore(PScoringParam* p_one, PScoringParam* p_two) {
       return SCWRL4PairwiseScore(reinterpret_cast<SCWRL4Param*>(p_one),
                                  reinterpret_cast<SCWRL4Param*>(p_two));
     }
+    case SCWRL3: {
+      return SCWRL3PairwiseScore(reinterpret_cast<SCWRL3Param*>(p_one),
+                                 reinterpret_cast<SCWRL3Param*>(p_two));
+    }
     default: {
       throw promod3::Error("Don't know what scoring function to call between "
                            "particles... maybe a lazy developer forgot to "
@@ -85,6 +90,16 @@ void PScoringParamContainer::PairwiseScores(const PScoringParamContainer& other,
       }
       break;
     }
+    case SCWRL3: {
+      scores.resize(idx.size());
+      for(uint i = 0; i < idx.size(); ++i) {
+        scores[i] = 
+        SCWRL3PairwiseScore(reinterpret_cast<SCWRL3Param*>(params_[idx[i].first]),
+                            reinterpret_cast<SCWRL3Param*>(other.params_[idx[i].second]));
+      }
+      break;
+
+    }
     default: {
       throw promod3::Error("Don't know what scoring function to call between "
                            "particles... maybe a lazy developer forgot to "

sidechain/src/particle_scoring_base.hh

@@ -25,6 +25,7 @@ namespace promod3{ namespace sidechain{
 // the available particle scoring functions
 enum PScoringFunction{
 SCWRL4,
+SCWRL3,
 INVALID_PSCORING_FUNCTION
 };
 

sidechain/src/scwrl3_particle_scoring.cc

+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
+
+
+namespace promod3 { namespace sidechain {
+
+void SCWRL3Param::ApplyTransform(const geom::Mat4& t) {
+
+  // get transformed pos
+  Real a = t(0,0)*pos_[0] + t(0,1)*pos_[1] + t(0,2)*pos_[2] + t(0,3);
+  Real b = t(1,0)*pos_[0] + t(1,1)*pos_[1] + t(1,2)*pos_[2] + t(1,3);
+  Real c = t(2,0)*pos_[0] + t(2,1)*pos_[1] + t(2,2)*pos_[2] + t(2,3);
+  pos_ = geom::Vec3(a, b, c);
+}
+
+
+SCWRL3Param* SCWRL3Param::Copy() const {
+  SCWRL3Param* return_ptr = new SCWRL3Param(pos_, radius_);
+  return return_ptr;
+}
+
+
+bool SCWRL3Param::EqualTo(PScoringParam* other) const {
+
+  if(other == NULL) return false;
+  if(other->GetScoringFunction() != SCWRL3) return false;
+
+  // as the other scoring function is also SCWRL3, we can assume that
+  // the following dynamic cast is fine...
+  SCWRL3Param* p = dynamic_cast<SCWRL3Param*>(other);
+
+  return pos_ == p->pos_ &&
+         radius_ == p->radius_;
+}
+
+
+Real SCWRL3PairwiseScore(SCWRL3Param* p_one, SCWRL3Param* p_two) {
+
+  Real d = geom::Distance(p_one->pos_, p_two->pos_);
+  Real R = p_one->radius_ + p_two->radius_;
+  Real e = 0.0;
+  if(d < Real(0.8254) * R) {
+    e = 10.0;
+  } else if(d <= R) {
+    e = 57.273 * (1 - d/R);
+  }
+  return e;
+}
+
+}}//ns

sidechain/src/scwrl3_particle_scoring.hh

+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef PROMOD3_SCWRL3_PARTICLE_SCORING_HH
+#define PROMOD3_SCWRL3_PARTICLE_SCORING_HH
+
+#include <ost/geom/mat4.hh>
+#include <ost/geom/vecmat3_op.hh>
+#include <promod3/sidechain/particle_scoring_base.hh>
+
+namespace promod3{ namespace sidechain{
+
+struct SCWRL3Param : public PScoringParam {
+
+  SCWRL3Param(const geom::Vec3& pos, Real r): pos_(pos), radius_(r) { }
+
+  virtual ~SCWRL3Param() { }
+
+  virtual const geom::Vec3& GetPos() const { return pos_; }
+
+  virtual Real GetCollisionDistance() const { return radius_; }
+
+  virtual void ApplyTransform(const geom::Mat4& t);
+
+  virtual SCWRL3Param* Copy() const;
+
+  virtual bool EqualTo(PScoringParam* other) const;
+
+  virtual PScoringFunction GetScoringFunction() const { return SCWRL3; }
+
+  geom::Vec3 pos_;
+  Real radius_;
+};
+
+
+Real SCWRL3PairwiseScore(SCWRL3Param* p_one, SCWRL3Param* p_two);
+
+}} //ns
+
+#endif

sidechain/src/scwrl3_rotamer_constructor.cc

+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#include <promod3/sidechain/scwrl3_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
+#include <promod3/core/runtime_profiling.hh>
+#include <promod3/core/message.hh>
+#include <promod3/core/runtime_profiling.hh>
+#include <ost/conop/conop.hh>
+
+namespace promod3 { namespace sidechain {
+
+struct SCWRL3RotamerParam : public RotamerLookupParam {
+
+  SCWRL3RotamerParam() {
+
+    // overwrite default internal_e_prefactor from base class
+    internal_e_prefactor[ost::conop::ARG] = 3.0;
+    internal_e_prefactor[ost::conop::ASN] = 3.0;
+    internal_e_prefactor[ost::conop::ASP] = 3.0;
+    internal_e_prefactor[ost::conop::GLN] = 3.0;
+    internal_e_prefactor[ost::conop::GLU] = 3.0;
+    internal_e_prefactor[ost::conop::LYS] = 3.0;
+    internal_e_prefactor[ost::conop::SER] = 3.0;
+    internal_e_prefactor[ost::conop::CYS] = 3.0;
+    internal_e_prefactor[ost::conop::MET] = 3.0;
+    internal_e_prefactor[ost::conop::TRP] = 3.0;
+    internal_e_prefactor[ost::conop::TYR] = 3.0;
+    internal_e_prefactor[ost::conop::THR] = 3.0;
+    internal_e_prefactor[ost::conop::VAL] = 3.0;
+    internal_e_prefactor[ost::conop::ILE] = 3.0;
+    internal_e_prefactor[ost::conop::LEU] = 3.0;
+    internal_e_prefactor[ost::conop::PRO] = 3.0;
+    internal_e_prefactor[ost::conop::HIS] = 3.0;
+    internal_e_prefactor[ost::conop::PHE] = 3.0;
+  }  
+
+  static const SCWRL3RotamerParam& Instance() {
+    static SCWRL3RotamerParam scwrl3_param;
+    return scwrl3_param;
+  }
+};
+
+
+SCWRL3RotamerConstructor::SCWRL3RotamerConstructor(bool cb_in_sidechain): 
+                          RotamerConstructor(cb_in_sidechain, HEAVY_ATOM_MODE,
+                                             SCWRL3RotamerParam::Instance()) { }
+
+void SCWRL3RotamerConstructor::AssignInternalEnergies(RRMRotamerGroupPtr group,
+                                                      RotamerID id, 
+                                                      uint residue_index,
+                                                      Real phi, Real psi,
+                                                      bool n_ter, bool c_ter) {
+
+  core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
+          "SCWRL3RotamerConstructor::AssignInternalEnergies", 2);
+
+  Real max_p = group->GetMaxP();
+  for(uint i = 0; i < group->size(); ++i){
+    RRMRotamerPtr r = (*group)[i];
+    Real internal_e_prefactor = r->GetInternalEnergyPrefactor();
+    Real probability = r->GetProbability();
+    Real e = -internal_e_prefactor * std::log(probability/max_p);
+    r->SetInternalEnergy(e);
+  }
+}
+
+
+void SCWRL3RotamerConstructor::AssignInternalEnergies(FRMRotamerGroupPtr group,
+                                                      RotamerID id, 
+                                                      uint residue_index,
+                                                      Real phi, Real psi,
+                                                      bool n_ter, bool c_ter) {
+
+  core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
+          "SCWRL3RotamerConstructor::AssignInternalEnergies", 2);
+
+  Real max_p = group->GetMaxP();
+  for(uint i = 0; i < group->size(); ++i){
+    FRMRotamerPtr r = (*group)[i];
+    Real internal_e_prefactor = r->GetInternalEnergyPrefactor();
+    Real probability = r->GetProbability();
+    Real e = -internal_e_prefactor * std::log(probability/max_p);
+    r->SetInternalEnergy(e);
+  }
+}
+
+
+void SCWRL3RotamerConstructor::ParametrizeParticle(int at_idx, 
+                                                   bool is_hydrogen, 
+                                                   Particle& particle) {
+
+  if(is_hydrogen) {
+    throw promod3::Error("Expect no hydrogen in SCWRL3RotamerConstructor!");
+  }
+
+  geom::Vec3 pos;
+  Real radius = 0.0;
+
+  if(at_idx == -2) {
+    pos = terminal_o_pos_;
+    radius = 1.3;
+  } else if (at_idx == -1){
+    pos = terminal_oxt_pos_;
+    radius = 1.3;
+  } else {
+    pos = pos_buffer_->GetPos(id_, at_idx);
+    ost::conop::AminoAcid aa = RotIDToAA(id_);
+    promod3::loop::AminoAcidAtom aaa =
+    promod3::loop::AminoAcidLookup::GetInstance().GetAAA(aa, at_idx);
+    String ele = promod3::loop::AminoAcidLookup::GetInstance().GetElement(aaa);
+    if(ele == "C") {
+      radius = 1.6;
+    } else if(ele == "N") {
+      radius = 1.3;
+    } else if(ele == "O") {
+      radius = 1.3;
+    } else if (ele == "S") {
+      radius = 1.7;
+    }
+  }
+
+  SCWRL3Param* p = new SCWRL3Param(pos, radius);
+  particle.SetSParam(p);
+}
+
+}} //ns

sidechain/src/scwrl3_rotamer_constructor.hh

+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#ifndef PROMOD3_SCWRL3_ROTAMER_CONSTRUCTOR_HH
+#define PROMOD3_SCWRL3_ROTAMER_CONSTRUCTOR_HH
+
+#include <promod3/sidechain/rotamer_constructor.hh>
+
+namespace promod3 { namespace sidechain {
+
+class SCWRL3RotamerConstructor;
+typedef boost::shared_ptr<SCWRL3RotamerConstructor> SCWRL3RotamerConstructorPtr;
+
+class SCWRL3RotamerConstructor : public RotamerConstructor{
+
+public:
+
+  SCWRL3RotamerConstructor(bool cb_in_sidechain = false);
+
+  virtual ~SCWRL3RotamerConstructor() { }
+
+  // Assign internal energies to rotamer groups
+  virtual void AssignInternalEnergies(RRMRotamerGroupPtr group,
+                                      RotamerID id,
+                                      uint residue_index,
+                                      Real phi = -1.0472, 
+                                      Real psi =  -0.7854,
+                                      bool n_ter = false,
+                                      bool c_ter = false);
+
+  virtual void AssignInternalEnergies(FRMRotamerGroupPtr group,
+                                      RotamerID id,
+                                      uint residue_index,
+                                      Real phi = -1.0472, 
+                                      Real psi =  -0.7854,
+                                      bool n_ter = false,
+                                      bool c_ter = false);
+
+private:
+
+  void ParametrizeParticle(int atom_idx, bool is_hydrogen, Particle& p);
+};
+
+}} // ns
+
+#endif