diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 4f1f03d4ef8424beb299d419c577cce61eebb0a8..5f3e7f42ace1d12c81aabe0ab2bc53a25bb4b059 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -22,6 +22,7 @@
#include <ost/mol/alg/local_dist_test.hh>
#include <ost/mol/alg/superpose_frames.hh>
#include <ost/mol/alg/filter_clashes.hh>
+#include <ost/mol/alg/construct_cbeta.hh>
using namespace boost::python;
using namespace ost;
@@ -53,4 +54,6 @@ BOOST_PYTHON_MODULE(_mol_alg)
def("SuperposeFrames", &ost::mol::alg::SuperposeFrames,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
arg("end")=-1, arg("ref")=-1));
+// def("ConstructCBetas", one_arg, args("entity_handle"));
+ def("ConstructCBetas", &ost::mol::alg::ConstructCBetas, (arg("entity_handle"), arg("include_gly")=false));
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index d1c84aa186a4fb77a47bde198acf8adc1f472371..9ba04a5ee60442edc6ab9522b23222b46a8331f0 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -5,6 +5,7 @@ set(OST_MOL_ALG_HEADERS
local_dist_test.hh
superpose_frames.hh
filter_clashes.hh
+ construct_cbeta.hh
)
set(OST_MOL_ALG_SOURCES
@@ -13,6 +14,7 @@ set(OST_MOL_ALG_SOURCES
local_dist_test.cc
superpose_frames.cc
filter_clashes.cc
+ construct_cbeta.cc
)
set(MOL_ALG_DEPS mol seq)
diff --git a/modules/mol/alg/src/construct_cbeta.cc b/modules/mol/alg/src/construct_cbeta.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c7466ef1e46db1d8cff38e98f187021bf560234f
--- /dev/null
+++ b/modules/mol/alg/src/construct_cbeta.cc
@@ -0,0 +1,93 @@
+#include <cmath>
+
+#include <ost/mol/mol.hh>
+#include <ost/log.hh>
+#include "construct_cbeta.hh"
+
+namespace ost { namespace mol { namespace alg {
+
+namespace detail {
+
+
+
+class CBetaInserter:public EntityVisitor {
+public:
+ CBetaInserter(mol::EntityHandle ent);
+
+ virtual bool VisitResidue(const mol::ResidueHandle& res);
+private:
+ XCSEditor edi_;
+};
+CBetaInserter::CBetaInserter(mol::EntityHandle ent) {
+ edi_=ent.EditXCS(BUFFERED_EDIT);
+}
+
+bool CBetaInserter::VisitResidue(const mol::ResidueHandle& res)
+{
+ if (!res.IsPeptideLinking()) {
+ return false;
+ }
+ try {
+ if(res.GetOneLetterCode()=='G' and include_gly_==false) {
+ return false;
+ }
+ if (res.FindAtom("CB").IsValid()) {
+ return false;
+ }
+ geom::Vec3 cbeta_pos=CBetaPosition(res);
+ edi_.InsertAtom(res, "CB", cbeta_pos, "C");
+ }
+ catch(...) {
+ LOG_WARNING("residue " << res.GetQualifiedName()
+ << "doesn't have enough backbone atoms to reconstruct CBeta position");
+ return false;
+ }
+ return false;
+}
+
+} // namespace
+
+
+// Todo: Applying on an entity does not work since the handle is used and the constructed
+// C-betas are there not present in the resulting entity-reference -> use handle onyl or
+// ensure that C-betas are added to the view as well
+void ConstructCBetas(mol::EntityHandle& entity_handle, bool include_gly)
+{
+ include_gly_=include_gly;
+ detail::CBetaInserter cbi(entity_handle);
+ entity_handle.Apply(cbi);
+}
+
+
+geom::Vec3 CBetaPosition(const geom::Vec3& n_pos, const geom::Vec3& ca_pos,
+ const geom::Vec3& c_pos, Real bond_length) {
+
+ geom::Vec3 v1=geom::Normalize(ca_pos-n_pos);
+ geom::Vec3 v2=geom::Normalize(ca_pos-c_pos);
+ geom::Vec3 in_plane_v=geom::Normalize(v1+v2);
+ geom::Plane p(ca_pos ,n_pos, c_pos);
+ // rotate around vector perpendicular to p and in_plane_v
+ geom::Vec3 axis=geom::Normalize(geom::Cross(p.GetNormal(), in_plane_v));
+ geom::Mat3 rot_mat=geom::AxisRotation(axis, (-54/180.0)*M_PI);
+ return ca_pos+rot_mat*in_plane_v*bond_length;
+}
+
+
+geom::Vec3 CBetaPosition(const ResidueHandle& residue, Real bond_length)
+{
+ AtomHandle ca=residue.FindAtom("CA");
+ AtomHandle n=residue.FindAtom("N");
+ AtomHandle c=residue.FindAtom("C");
+ if (!(ca.IsValid() && c.IsValid() && n.IsValid())) {
+ std::stringstream ss;
+ ss << "residue " << residue.GetQualifiedName()
+ << "doesn't have enough backbone atoms to reconstruct CBeta position";
+ throw std::runtime_error(ss.str());
+ }
+ return CBetaPosition(n.GetPos(), ca.GetPos(), c.GetPos(), bond_length);
+}
+
+
+}}} // ns
+
+
diff --git a/modules/mol/alg/src/construct_cbeta.hh b/modules/mol/alg/src/construct_cbeta.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ad59bf8d5ffd50da772819fa5b6e7818120de8f5
--- /dev/null
+++ b/modules/mol/alg/src/construct_cbeta.hh
@@ -0,0 +1,42 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+#ifndef OST_MOL_ALG_CONSTRUCT_CBETA_HH
+#define OST_MOL_ALG_CONSTRUCT_CBETA_HH
+
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/alg/module_config.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+
+geom::Vec3 DLLEXPORT_OST_MOL_ALG CBetaPosition(const ResidueHandle& residue,
+ Real bond_length=1.5);
+
+geom::Vec3 DLLEXPORT_OST_MOL_ALG CBetaPosition(const geom::Vec3& n_pos,
+ const geom::Vec3& ca_pos,
+ const geom::Vec3& c_pos,
+ Real bond_length=1.5);
+
+void DLLEXPORT_OST_MOL_ALG ConstructCBetas(EntityHandle& entity_handle, bool include_gly=false);
+
+bool include_gly_;
+
+}}} // ns
+
+#endif // OST_MOL_ALG_CONSTRUCT_CBETA_HH
diff --git a/modules/qa/pymod/export_interaction.cc b/modules/qa/pymod/export_interaction.cc
index fe61394c0b819f8d32e96dbdc1c025cc9f174ac2..ad07cdfa3a23e0fcab697be91ac50b8bc6fc576e 100644
--- a/modules/qa/pymod/export_interaction.cc
+++ b/modules/qa/pymod/export_interaction.cc
@@ -34,9 +34,11 @@ F2 f2=&InteractionStatistics::GetCount;
-typedef float (AllAtomPotential::*oneArg)(mol::EntityView);
+typedef float (AllAtomPotential::*oneArg)(mol::EntityView,
+ std::string);
typedef float (AllAtomPotential::*twoArgs)(mol::EntityView,
- mol::EntityView);
+ mol::EntityView,
+ std::string);
typedef float (AllAtomPotential::*threeArgs)(atom::ChemType,
atom::ChemType, float);
@@ -251,6 +253,7 @@ void export_Interaction()
;
class_<AllAtomPotential>("AllAtomPotential", no_init)
.def("Create", &AllAtomPotential::Create).staticmethod("Create")
+ .def("Repair", &AllAtomPotential::Repair).staticmethod("Repair")
.def("LoadFromFile",
&AllAtomPotential::LoadFromFile).staticmethod("LoadFromFile")
.def("GetTotalEnergy", one_arg, args("view"))
diff --git a/modules/qa/src/all_atom_potential.cc b/modules/qa/src/all_atom_potential.cc
index 18a4ce5f01488a74a5d62d24553e5804f7938e10..d3af05a1c927523716b0aee63e09ee3162836157 100644
--- a/modules/qa/src/all_atom_potential.cc
+++ b/modules/qa/src/all_atom_potential.cc
@@ -37,18 +37,24 @@ class AllAtomPotentialCalculator : public mol::EntityVisitor {
public:
AllAtomPotentialCalculator(AllAtomPotential::AllAtomEnergies& energies,
AllAtomPotentialOpts& opts,
- mol::EntityView target_view):
+ mol::EntityView target_view,
+ int s):
energies_(energies),
options_(opts),
energy_(0.0),
target_view_(target_view),
- interaction_counts_(0)
+ interaction_counts_(0),
+ interaction_counts_vector_(std::vector<int>(s, 0)),
+ score_vector_(std::vector<float>(s, 0.0)),
+ pos_counter(0)
{
}
virtual bool VisitResidue(const mol::ResidueHandle& residue)
{
+ if(curr_aa_!=Xxx) pos_counter++;
curr_aa_=ResidueToAminoAcid(residue);
+// std::cout << residue_energy_ << " " << interaction_counts_ << std::endl;
return curr_aa_!=Xxx;
}
@@ -85,13 +91,16 @@ public:
continue;
}
energy_+=energies_.Get(type_a, type_b, d);
+ interaction_counts_vector_[pos_counter-1]++;
+ score_vector_[pos_counter-1]+=energies_.Get(type_a, type_b, d);
interaction_counts_++;
- }
+ }
}
}
}
return false;
}
+
float GetEnergy() const
{
return energy_;
@@ -102,13 +111,27 @@ public:
return interaction_counts_;
}
+ std::vector<float> GetResidueEnergyVector() const
+ {
+ return score_vector_;
+ }
+
+ std::vector<int> GetResidueCountsVector() const
+ {
+ return interaction_counts_vector_;
+ ;
+ }
+
private:
AllAtomPotential::AllAtomEnergies& energies_;
AllAtomPotentialOpts options_;
AminoAcid curr_aa_;
- float energy_;
+ float energy_;
mol::EntityView target_view_;
int interaction_counts_;
+ std::vector<int> interaction_counts_vector_;
+ std::vector<float> score_vector_;
+ int pos_counter;
};
}
@@ -154,22 +177,44 @@ AllAtomPotentialPtr AllAtomPotential::LoadFromFile(const String& filename)
}
float AllAtomPotential::GetTotalEnergy(mol::EntityView view,
- mol::EntityView target_view)
+ mol::EntityView target_view, std::string property_identifier="")
{
- AllAtomPotentialCalculator c(energies_, options_, target_view);
mol::EntityHandle e=view.GetHandle();
+ int res_count=view.GetResidueList().size();
+ AllAtomPotentialCalculator c(energies_, options_, target_view, res_count);
view.Apply(c);
interaction_counts_=c.GetEnergyCounts();
+
+ std::vector<float> rev=c.GetResidueEnergyVector();
+ std::vector<int> rcv=c.GetResidueCountsVector();
+
+ if(property_identifier!="") {
+ for(int i=0;i<res_count;++i) {
+ view.GetResidueList()[i].SetFloatProp(property_identifier, rev[i]);
+ view.GetResidueList()[i].SetIntProp(property_identifier+"_counts", rcv[i]);
+ }
+ }
return c.GetEnergy();
}
-float AllAtomPotential::GetTotalEnergy(mol::EntityView view)
+float AllAtomPotential::GetTotalEnergy(mol::EntityView view, std::string property_identifier="")
{
- AllAtomPotentialCalculator c(energies_, options_, view);
mol::EntityHandle e=view.GetHandle();
+ int res_count=view.GetResidueList().size();
+ AllAtomPotentialCalculator c(energies_, options_, view, res_count);
view.Apply(c);
interaction_counts_=c.GetEnergyCounts();
- return c.GetEnergy();
+ std::string a="";
+ std::vector<float> rev=c.GetResidueEnergyVector();
+ std::vector<int> rcv=c.GetResidueCountsVector();
+
+ if(property_identifier!="") {
+ for(int i=0;i<res_count;++i) {
+ view.GetResidueList()[i].SetFloatProp(property_identifier, rev[i]);
+ view.GetResidueList()[i].SetIntProp(property_identifier+"_counts", rcv[i]);
+ }
+ }
+ return c.GetEnergy();
}
void AllAtomPotential::SetSequenceSeparation(int seq_sep) {
@@ -208,7 +253,46 @@ void AllAtomPotentialOpts::Serialize(DS& ds)
ds & distance_bucket_size;
}
+void AllAtomPotential::Repair()
+{
+ // for C-apha (backbone only) models fill Calpha interactions with corresponding Cbeta values
+ typedef AllAtomPotential::AllAtomEnergies::IndexType Index;
+ int num=int((options_.upper_cutoff-
+ options_.lower_cutoff)/options_.distance_bucket_size);
+
+ for (int i=0; i<atom::UNKNOWN-1; ++i) {
+ for (int j=0; j<atom::UNKNOWN-1; ++j) {
+
+ for (int k=0; k<num; ++k) {
+ float e=energies_.Get(Index(i, j, k));
+ //check if Cbeta (has counts) and not Glycin-Calpha
+ if(i==0 or j==0){ //mysteriously needed for one case of j=0 (TODO: check later, perhaps in RepairCbetaStatistics())
+ continue;
+ }
+ if (e!= 0 and i==3 and energies_.Get(Index(i, j+1, k))==0) {
+ energies_.Set(Index(i, j-1, k), e);
+// std::cout << e << " ";
+ }
+ else if (e!= 0 and j==3 and energies_.Get(Index(i+1, j, k))==0) {
+ energies_.Set(Index(i-1, j, k), e);
+// std::cout << e << " ";
+ }
+ else {
+ if(e!= 0 and energies_.Get(Index(i, j+1, k))==0 and energies_.Get(Index(i+1, j, k))==0) {
+// std::cout << e << " ";
+
+ energies_.Set(Index(i-1, j, k), e);
+ energies_.Set(Index(i, j-1, k), e);
+ energies_.Set(Index(i-1, j-1, k), e);
+ }
+ }
+
+ }
+ }
+
+ }
+}
void AllAtomPotential::Fill(const InteractionStatisticsPtr& stats)
{
diff --git a/modules/qa/src/all_atom_potential.hh b/modules/qa/src/all_atom_potential.hh
index 0a1a99b7aeec1fa296645b6f7aeba81fd7a42ab5..9818ec7497732f2f7cccd153397b35b301290b5e 100644
--- a/modules/qa/src/all_atom_potential.hh
+++ b/modules/qa/src/all_atom_potential.hh
@@ -74,7 +74,7 @@ public:
void SaveToFile(const String& filename);
/// \brief calculate all-atom interaction score for whole entity
- float GetTotalEnergy(mol::EntityView view);
+ float GetTotalEnergy(mol::EntityView view, std::string property_identifier);
/// \brief extract energy of a specific interaction
/// (for plotting pseudo Lennard-Jones potential).
@@ -87,7 +87,7 @@ public:
/// Two entities need to be provided:
/// the atoms for which the energy should be derived and
/// the atoms with respect to which the energy should be calculted.
- float GetTotalEnergy(mol::EntityView view, mol::EntityView target_view);
+ float GetTotalEnergy(mol::EntityView view, mol::EntityView target_view, std::string property_identifier);
/// \brief retrieve total number of interactions (for normalisation)
int GetEnergyCounts() const { return interaction_counts_; }
@@ -99,6 +99,7 @@ public:
template <typename DS>
void Serialize(DS& ds);
public:
+ void Repair();
void Fill(const InteractionStatisticsPtr& stats);
/// parameters: atom type one, atom type two, distance
diff --git a/modules/qa/src/packing_potential.cc b/modules/qa/src/packing_potential.cc
index 10a7de4aeb20f1642f743915bc3cab7719d7062c..f46b0a6cf61c85772343903faf8d42d42bf28b4c 100644
--- a/modules/qa/src/packing_potential.cc
+++ b/modules/qa/src/packing_potential.cc
@@ -180,14 +180,18 @@ int PackingPotential::GetEnergyCounts()
bool PackingPotential::VisitAtom(const mol::AtomHandle& atom)
{
AminoAcid aa=ResidueToAminoAcid(atom.GetResidue());
- if (aa==Xxx)
+ if (aa==Xxx) {
+ atom.GetResidue().SetFloatProp("solvation_energy",0);
return false;
+ }
int count=0;
for (mol::EntityViewList::iterator i=views_.begin(),
e=views_.end(); i!=e; ++i) {
count+=i->FindWithin(atom.GetPos(), options_.cutoff).size();
}
- energy_+=this->GetPackingEnergy(aa, count);
+ float local_energy=this->GetPackingEnergy(aa, count);
+ atom.GetResidue().SetFloatProp("solvation_energy",local_energy);
+ energy_+=local_energy;
energy_counts_++;
return false;
}
diff --git a/modules/qa/src/torsion_potential.cc b/modules/qa/src/torsion_potential.cc
index 2ecf8c655f24de49220a00be1b09b24d1e14b856..1440f2d719221aa2d979de907cd4765cf48da199 100644
--- a/modules/qa/src/torsion_potential.cc
+++ b/modules/qa/src/torsion_potential.cc
@@ -41,13 +41,15 @@ namespace {
class TorsionEnergyCalc : public mol::EntityVisitor {
public:
TorsionEnergyCalc(const TorsionPotentialPtr& pot,
- TorsionPotentialOpts opts):
- pot_(pot), energy_(0.0), num_torsions_(0)
+ TorsionPotentialOpts opts, int length):
+ pot_(pot), energy_(0.0), num_torsions_(0), score_vector_(std::vector<float>(length, 0.0)),
+ pos_counter_(0)
{
}
virtual bool VisitResidue(const mol::ResidueHandle& rh)
{
+ pos_counter_++;
mol::ResidueHandle prev=rh.GetPrev();
mol::ResidueHandle next=rh.GetNext();
if (!(next && prev && next.IsPeptideLinking() &&
@@ -77,8 +79,11 @@ public:
// calculate position of the amino acid in the alphabet
- energy_+=pot_->GetTorsionEnergy(ca, prev_phi, prev_psi, central_phi,
+ float local_energy=pot_->GetTorsionEnergy(ca, prev_phi, prev_psi, central_phi,
central_psi, next_phi, next_psi);
+ energy_+=local_energy;
+
+ score_vector_[pos_counter_-1]=local_energy;
++num_torsions_;
return false;
}
@@ -91,13 +96,19 @@ public:
return num_torsions_;
}
+ std::vector<float> GetResidueEnergyVector() const {
+ return score_vector_;
+ }
+
private:
TorsionPotentialPtr pot_;
AminoAcid prev_;
AminoAcid center_;
mol::ResidueHandle cr_;
Real energy_;
- int num_torsions_;
+ int num_torsions_;
+ std::vector<float> score_vector_;
+ int pos_counter_;
};
}
@@ -211,17 +222,27 @@ void TorsionPotential::SaveToFile(const String& path)
Real TorsionPotential::GetTotalEnergy(mol::EntityHandle entity)
{
- TorsionEnergyCalc c(this->shared_from_this(), options_);
+ int res_count=entity.GetResidueList().size();
+ TorsionEnergyCalc c(this->shared_from_this(), options_, res_count);
entity.Apply(c);
num_torsions_ = c.GetEnergyCounts();
+ std::vector<float> rev=c.GetResidueEnergyVector();
+ for(int i=0;i<res_count;++i) {
+ entity.GetResidueList()[i].SetFloatProp("torsion_energy", rev[i]);
+ }
return c.GetEnergy();
}
Real TorsionPotential::GetTotalEnergy(mol::EntityView entity)
{
- TorsionEnergyCalc c(this->shared_from_this(), options_);
+ int res_count=entity.GetResidueList().size();
+ TorsionEnergyCalc c(this->shared_from_this(), options_, res_count);
entity.Apply(c);
num_torsions_ = c.GetEnergyCounts();
+ std::vector<float> rev=c.GetResidueEnergyVector();
+ for(int i=0;i<res_count;++i) {
+ entity.GetResidueList()[i].SetFloatProp("torsion_energy", rev[i]);
+ }
return c.GetEnergy();
}