Skip to content
Snippets Groups Projects
Commit 7f49c161 authored by Marco Biasini's avatar Marco Biasini
Browse files

added residue level (reduced) statistical potential

parent cdf8cccb
No related branches found
No related tags found
No related merge requests found
......@@ -2,6 +2,7 @@ set(OST_QA_PYMOD_SOURCES
export_interaction.cc
export_torsion.cc
export_packing.cc
export_reduced.cc
wrap_qa.cc
export_clash.cc
)
......
//------------------------------------------------------------------------------
// 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
//------------------------------------------------------------------------------
#include <boost/python.hpp>
#include <ost/mol/mol.hh>
#include <ost/qa/reduced_statistics.hh>
#include <ost/qa/reduced_potential.hh>
using namespace ost::qa;
using namespace boost::python;
uint64_t (ReducedStatistics::*get_count)(AminoAcid,AminoAcid,int,int)=&ReducedStatistics::GetCount;
void export_Reduced()
{
class_<ReducedStatOptions>("ReducedStatOptions", init<>())
.def(init<Real,Real,uint,uint,uint>((arg("l_cutoff"), arg("d_cutoff"),
arg("num_ang_bins"), arg("num_dist_bins"), arg("sequence_sep"))))
.def_readwrite("lower_cutoff", &ReducedStatOptions::lower_cutoff)
.def_readwrite("upper_cutoff", &ReducedStatOptions::upper_cutoff)
.def_readwrite("num_angular_bins", &ReducedStatOptions::num_angular_bins)
.def_readwrite("num_dist_bins", &ReducedStatOptions::num_dist_bins)
.def_readwrite("sequence_sep", &ReducedStatOptions::sequence_sep)
;
class_<ReducedStatistics, ReducedStatisticsPtr>("ReducedStatistics", no_init)
.def(init<Real,Real,uint,uint,uint>((arg("l_cutoff"), arg("u_cutoff"),
arg("num_ang_bins"), arg("num_dist_bins"), arg("sequence_sep"))))
.def("Extract", &ReducedStatistics::Extract)
.add_property("options", make_function(&ReducedStatistics::GetOptions,
return_value_policy<copy_const_reference>()))
.def("Save", &ReducedStatistics::Save)
.def("Load", &ReducedStatistics::Load).staticmethod("Load")
.def("GetCount", get_count)
;
class_<ReducedPotential, ReducedPotentialPtr>("ReducedPotential", no_init)
.def("Load", &ReducedPotential::Load).staticmethod("Load")
.def("Create", &ReducedPotential::Create).staticmethod("Create")
.def("Save", &ReducedPotential::Save)
.def("GetEnergy", &ReducedPotential::GetEnergy)
.def("GetTotalEnergy", &ReducedPotential::GetTotalEnergy,
(arg("ent"), arg("normalize")=true))
.add_property("options", make_function(&ReducedPotential::GetOptions,
return_value_policy<copy_const_reference>()))
;
}
\ No newline at end of file
......@@ -84,6 +84,7 @@ void export_Torsion()
.value("His", His)
.value("Phe", Phe)
.value("Trp", Trp)
.export_values()
;
class_<AminoAcidSet>("AminoAcidSet", init<>())
......
......@@ -22,10 +22,12 @@ void export_Torsion();
void export_Interaction();
void export_Packing();
void export_Clash();
void export_Reduced();
BOOST_PYTHON_MODULE(_qa)
{
export_Torsion();
export_Interaction();
export_Packing();
export_Clash();
export_Reduced();
}
......@@ -13,10 +13,14 @@ packing_statistics.hh
atom_types.hh
packing_potential.hh
clash_score.hh
reduced_statistics.hh
reduced_potential.hh
module_config.hh
)
set(OST_QA_SOURCES
reduced_statistics.cc
reduced_potential.cc
torsion_statistics.cc
torsion_potential.cc
solis_torsion_statistics.cc
......@@ -28,8 +32,10 @@ atom_types.cc
all_atom_potential.cc
packing_potential.cc
clash_score.cc
impl/reduced_impl.cc
)
module(NAME qa SOURCES ${OST_QA_SOURCES}
HEADERS ${OST_QA_HEADERS} DEPENDS_ON io)
HEADERS reduced_impl.hh IN_DIR impl ${OST_QA_HEADERS} DEPENDS_ON io)
#include <ost/log.hh>
#include "reduced_impl.hh"
namespace ost { namespace qa { namespace impl {
bool ReducedPotentialImpl::VisitResidue(const mol::ResidueHandle& res)
{
if (!res.IsPeptideLinking()) {
return false;
}
AminoAcid aa_one=OneLetterCodeToAminoAcid(res.GetOneLetterCode());
if (aa_one==Xxx) {
return false;
}
geom::Vec3 ca_pos_one;
geom::Vec3 cb_pos_one;
uint index=res.GetIndex();
if (!this->GetCAlphaCBetaPos(res, ca_pos_one, cb_pos_one)) {
return false;
}
mol::AtomHandleList within=ent_.FindWithin(ca_pos_one,
opts_.upper_cutoff-0.00001);
for (mol::AtomHandleList::const_iterator
i=within.begin(), e=within.end(); i!=e; ++i) {
mol::ResidueHandle res_two=i->GetResidue();
if (res_two.GetIndex()-index<opts_.sequence_sep) {
continue;
}
if (i->GetName()!="CA" || !res_two.IsPeptideLinking()) {
continue;
}
AminoAcid aa_two=OneLetterCodeToAminoAcid(res_two.GetOneLetterCode());
if (aa_two==Xxx) {
continue;
}
geom::Vec3 ca_pos_two;
geom::Vec3 cb_pos_two;
if (!this->GetCAlphaCBetaPos(res_two, ca_pos_two, cb_pos_two)) {
continue;
}
Real dist=geom::Length(ca_pos_one-ca_pos_two);
if (dist<opts_.lower_cutoff) {
continue;
}
Real angle=geom::Angle(cb_pos_one-ca_pos_one, cb_pos_two-ca_pos_two);
this->OnInteraction(aa_one, aa_two, dist, angle);
}
return false;
}
bool ReducedPotentialImpl::GetCAlphaCBetaPos(const mol::ResidueHandle& res,
geom::Vec3& ca_pos,
geom::Vec3& cb_pos)
{
const static Real bond_length=1.5;
mol::AtomHandle ca=res.FindAtom("CA");
if (!ca.IsValid()) {
return false;
}
ca_pos=ca.GetPos();
mol::AtomHandle cb=res.FindAtom("CB");
if (cb.IsValid()) {
cb_pos=cb.GetPos();
return true;
}
mol::AtomHandle n=res.FindAtom("N");
mol::AtomHandle c=res.FindAtom("C");
if (!(ca.IsValid() && c.IsValid() && n.IsValid())) {
LOG_WARNING("residue " << res.GetQualifiedName()
<< " doesn't have enough atoms to reconstruct Cbeta position");
return false;
}
geom::Vec3 v1=geom::Normalize(ca.GetPos()-n.GetPos());
geom::Vec3 v2=geom::Normalize(ca.GetPos()-c.GetPos());
geom::Vec3 in_plane_v=geom::Normalize(v1+v2);
geom::Plane p(ca.GetPos() ,n.GetPos(), c.GetPos());
// 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);
cb_pos=ca.GetPos()+rot_mat*in_plane_v*bond_length;
return true;
}
}}}
//------------------------------------------------------------------------------
// 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_QA_IMPL_REDUCED_HH
#define OST_QA_IMPL_REDUCED_HH
#include <ost/mol/mol.hh>
#include <ost/qa/module_config.hh>
#include <ost/qa/amino_acids.hh>
#include <ost/qa/reduced_statistics.hh>
namespace ost { namespace qa { namespace impl {
class DLLEXPORT_OST_QA ReducedPotentialImpl : public mol::EntityVisitor {
public:
ReducedPotentialImpl(const ReducedStatOptions& opts, mol::EntityHandle ent):
opts_(opts), ent_(ent)
{ }
virtual bool VisitResidue(const mol::ResidueHandle& res);
bool GetCAlphaCBetaPos(const mol::ResidueHandle& res, geom::Vec3& ca_pos,
geom::Vec3& cb_pos);
virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two,
Real dist, Real angle)=0;
private:
ReducedStatOptions opts_;
mol::EntityHandle ent_;
};
}}}
#endif
//------------------------------------------------------------------------------
// 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
//------------------------------------------------------------------------------
#include <ost/io/binary_data_source.hh>
#include <ost/io/binary_data_sink.hh>
#include <ost/io/container_serialization.hh>
#include "reduced_potential.hh"
#include "impl/reduced_impl.hh"
namespace ost { namespace qa {
namespace {
class ReducedEnergiesCalc : public impl::ReducedPotentialImpl {
public:
ReducedEnergiesCalc(const ReducedStatOptions& opts, ReducedEnergies& energies,
mol::EntityHandle ent, bool norm):
impl::ReducedPotentialImpl(opts, ent),
energies_(energies), energy_(0.0), norm_(norm), count_(0)
{ }
virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two,
Real dist, Real angle)
{
//std::cout << aa_one << " " << aa_two << " " << dist << " " << angle << std::endl;
energy_+=energies_.Get(aa_one, aa_two, dist, angle);
count_+=1;
}
Real GetEnergy() const
{
if (norm_) {
return count_ > 0 ? energy_ /count_ : 0.0;
}
return energy_;
}
private:
ReducedEnergies& energies_;
Real energy_;
bool norm_;
uint64_t count_;
};
}
ReducedPotentialPtr ReducedPotential::Load(const String& filename)
{
std::ifstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSource ds(stream);
ReducedPotentialPtr p(new ReducedPotential);
ds >> *p.get();
return p;
}
ReducedPotentialPtr ReducedPotential::Create(const ReducedStatisticsPtr& stat)
{
ReducedPotentialPtr pot(new ReducedPotential);
pot->opts_=stat->GetOptions();
pot->energies_=ReducedEnergies(0.0, IntegralClassifier(20, 0),
IntegralClassifier(20, 0),
ContinuousClassifier(pot->opts_.num_dist_bins,
pot->opts_.lower_cutoff,
pot->opts_.upper_cutoff),
ContinuousClassifier(pot->opts_.num_angular_bins,
0.0, M_PI));
pot->InitFromStats(stat);
return pot;
}
void ReducedPotential::InitFromStats(const ReducedStatisticsPtr& stats)
{
typedef ReducedEnergies::IndexType Index;
uint64_t total_count=stats->GetTotalCount();
for (size_t i=0; i<Xxx; ++i) {
for (size_t j=0; j<Xxx; ++j) {
uint64_t t2=stats->GetCount(AminoAcid(i), AminoAcid(j));
for (size_t k=0; k<opts_.num_dist_bins; ++k) {
for (size_t l=0; l<opts_.num_angular_bins; ++l) {
uint64_t t3=stats->GetCount(AminoAcid(i), AminoAcid(j), k, l);
// could be moved outside of the loop and cached...
uint64_t t4=stats->GetCount(k, l);
float f1=t2 > 0 ?(float(t3)/float(t2)) : 0.0;
float f2=total_count > 0 ? (float(t4)/float(total_count)) : 0.0;
float d = f2>1e-10 ? f1/f2: 0.0;
float e=0.582*(log(1+0.02*t2)-log(1+0.02*t2*d));
energies_.Set(Index(i, j, k, l), e);
}
}
}
}
}
void ReducedPotential::Save(const String& filename)
{
std::ofstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSink ds(stream);
ds << *this;
}
template <typename DS>
void ReducedPotential::Serialize(DS& ds)
{
ds & opts_;
ds & energies_;
}
Real ReducedPotential::GetTotalEnergy(ost::mol::EntityHandle ent, bool norm)
{
ReducedEnergiesCalc calc(opts_, energies_, ent, norm);
ent.Apply(calc);
return calc.GetEnergy();
}
}}
//------------------------------------------------------------------------------
// 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_QA_REDUCED_POTENTIAL_HH
#define OST_QA_REDUCED_POTENTIAL_HH
#include <ost/qa/reduced_statistics.hh>
namespace ost { namespace qa {
class ReducedPotential;
typedef boost::shared_ptr<ReducedPotential> ReducedPotentialPtr;
typedef MultiClassifier<float, int, int, float, float> ReducedEnergies;
/// \brief reduced model statistical potential operating on the residue level
class DLLEXPORT_OST_QA ReducedPotential {
public:
static ReducedPotentialPtr Create(const ReducedStatisticsPtr& stats);
static ReducedPotentialPtr Load(const String& filename);
void Save(const String& filename);
template <typename DS>
void Serialize(DS& ds);
const ReducedStatOptions& GetOptions() const { return opts_; }
Real GetTotalEnergy(ost::mol::EntityHandle ent, bool normalize=true);
Real GetEnergy(AminoAcid aa_one, AminoAcid aa_two,
Real ca_dist, Real angle) const
{
if (aa_one==Xxx || aa_two==Xxx) { return 0.0; }
if (ca_dist<opts_.lower_cutoff || ca_dist>=opts_.upper_cutoff) {
return 0.0;
}
if (angle<0.0 || angle>M_PI) { return 0.0; }
return energies_.Get(aa_one, aa_two, ca_dist, angle);
}
private:
void InitFromStats(const ReducedStatisticsPtr& stats);
ReducedPotential() { }
ReducedStatOptions opts_;
ReducedEnergies energies_;
};
}}
#endif
//------------------------------------------------------------------------------
// 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
//------------------------------------------------------------------------------
#include <ost/log.hh>
#include <ost/mol/mol.hh>
#include "amino_acids.hh"
#include "reduced_statistics.hh"
#include "impl/reduced_impl.hh"
#include <ost/io/binary_data_source.hh>
#include <ost/io/binary_data_sink.hh>
#include <ost/io/container_serialization.hh>
namespace ost { namespace qa {
namespace {
class ReducedStatExtractor : public impl::ReducedPotentialImpl {
public:
ReducedStatExtractor(const ReducedStatOptions& opts, ReducedHistogram& histo,
mol::EntityHandle ent):
impl::ReducedPotentialImpl(opts, ent),
histo_(histo)
{ }
virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two,
Real dist, Real angle)
{
histo_.Add(1, aa_one, aa_two, dist, angle);
}
private:
ReducedHistogram& histo_;
};
}
uint64_t ReducedStatistics::GetTotalCount() const
{
typedef ReducedHistogram::IndexType Index;
uint64_t count=0;
for (size_t i=0; i<Xxx; ++i) {
for (size_t j=0; j<Xxx; ++j) {
for (size_t k=0; k<opts_.num_dist_bins; ++k) {
for (size_t l=0; l<opts_.num_angular_bins; ++l) {
count+=histo_.Get(Index(i, j, k, l));
}
}
}
}
return count;
}
uint64_t ReducedStatistics::GetCount(AminoAcid aa_one, AminoAcid aa_two) const
{
typedef ReducedHistogram::IndexType Index;
uint64_t count=0;
for (size_t k=0; k<opts_.num_dist_bins; ++k) {
for (size_t l=0; l<opts_.num_angular_bins; ++l) {
count+=histo_.Get(Index(aa_one, aa_two, k, l));
}
}
return count;
}
uint64_t ReducedStatistics::GetCount(int dist_bin, int ang_bin) const
{
typedef ReducedHistogram::IndexType Index;
uint64_t count=0;
for (size_t i=0; i<Xxx; ++i) {
for (size_t j=0; j<Xxx; ++j) {
count+=histo_.Get(Index(i, j, dist_bin, ang_bin));
}
}
return count;
}
void ReducedStatistics::Save(const String& filename)
{
std::ofstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSink ds(stream);
ds << *this;
}
ReducedStatisticsPtr ReducedStatistics::Load(const String& filename)
{
std::ifstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSource ds(stream);
ReducedStatisticsPtr p(new ReducedStatistics);
ds >> *p.get();
return p;
}
template <typename DS>
void ReducedStatistics::Serialize(DS& ds)
{
ds & opts_;
ds & histo_;
}
void ReducedStatistics::Extract(mol::EntityHandle ent)
{
if (ent.GetChainCount()!=1) {
std::stringstream ss;
ss << "Expected exactly one chain, but entity has "
<< ent.GetChainCount() << "chains";
throw std::runtime_error(ss.str());
}
ReducedStatExtractor extractor(opts_, histo_, ent);
ent.Apply(extractor);
}
}}
\ No newline at end of file
//------------------------------------------------------------------------------
// 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_QA_REDUCED_STATISTICS_HH
#define OST_QA_REDUCED_STATISTICS_HH
#include <ost/qa/histogram.hh>
#include <ost/qa/amino_acids.hh>
namespace ost {
namespace mol {
class EntityHandle;
}
namespace qa {
/// \brief aggregates the option for the reduced model statistics
struct DLLEXPORT_OST_QA ReducedStatOptions {
ReducedStatOptions(): lower_cutoff(0), upper_cutoff(0),
num_angular_bins(0), num_dist_bins(0),
sequence_sep(0)
{ }
ReducedStatOptions(Real l_cutoff, Real u_cutoff, uint n_ang_bins,
uint n_dist_bins, uint ssep):
lower_cutoff(l_cutoff), upper_cutoff(u_cutoff),
num_angular_bins(n_ang_bins), num_dist_bins(n_dist_bins),
sequence_sep(ssep)
{ }
Real lower_cutoff;
Real upper_cutoff;
uint num_angular_bins;
uint num_dist_bins;
uint sequence_sep;
template <typename DS>
void Serialize(DS& ds)
{
ds & lower_cutoff;
ds & upper_cutoff;
ds & num_angular_bins;
ds & num_dist_bins;
ds & sequence_sep;
}
};
class ReducedStatistics;
typedef boost::shared_ptr<ReducedStatistics> ReducedStatisticsPtr;
// parametrized as first amino acid, second amino acid,
// calpha-calpha distance, angle
typedef Histogram<int, int, float, float> ReducedHistogram;
class DLLEXPORT_OST_QA ReducedStatistics {
public:
ReducedStatistics(Real l_cutoff, Real u_cutoff, uint num_ang_bins,
uint num_dist_bins, uint ssep):
opts_(l_cutoff, u_cutoff, num_ang_bins, num_dist_bins, ssep),
histo_(IntegralClassifier(20, 0), IntegralClassifier(20, 0),
ContinuousClassifier(num_dist_bins, l_cutoff, u_cutoff),
ContinuousClassifier(num_ang_bins, 0.0, M_PI))
{ }
const ReducedStatOptions& GetOptions() const { return opts_; }
/// \brief extract the statistics from the given entity
void Extract(mol::EntityHandle ent);
void Save(const String& filename);
static ReducedStatisticsPtr Load(const String& filename);
/// \internal
template <typename DS>
void Serialize(DS& ds);
uint64_t GetTotalCount() const;
uint64_t GetCount(AminoAcid aa_one, AminoAcid aa_two) const;
uint64_t GetCount(AminoAcid aa_one, AminoAcid aa_two, int dist_bin,
int ang_bin)
{
return histo_.Get(ReducedHistogram::IndexType(aa_one, aa_two, dist_bin,
ang_bin));
}
uint64_t GetCount(int dist_bin, int ang_bin) const;
private:
ReducedStatistics(): opts_(), histo_() {}
ReducedStatOptions opts_;
ReducedHistogram histo_;
};
}}
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment