diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 7b403f02d410525986db487beec6397c678c4fc3..1a02874c599e05c66af011a60f737207b26844c0 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -21,10 +21,13 @@
#include <boost/python/suite/indexing/map_indexing_suite.hpp>
#include <ost/config.hh>
#include <ost/mol/alg/local_dist_diff_test.hh>
+#include <ost/mol/alg/distance_test_common.hh>
#include <ost/mol/alg/superpose_frames.hh>
#include <ost/mol/alg/filter_clashes.hh>
#include <ost/mol/alg/consistency_checks.hh>
#include <ost/export_helper/pair_to_tuple_conv.hh>
+#include <ost/export_helper/vec_to_list_conv.hh>
+
using namespace boost::python;
using namespace ost;
@@ -91,6 +94,7 @@ ost::mol::alg::GlobalRDMap create_distance_list_from_multiple_references(const l
return ost::mol::alg::CreateDistanceListFromMultipleReferences(ref_list_vector, cutoff_list_vector, sequence_separation, max_dist);
}
+
}
@@ -117,7 +121,9 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("LDDTHA",&mol::alg::LDDTHA, (arg("sequence_separation")=0));
def("CreateDistanceList",&mol::alg::CreateDistanceList);
def("CreateDistanceListFromMultipleReferences",&create_distance_list_from_multiple_references);
-
+
+
+
def("SuperposeFrames", superpose_frames1,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
arg("end")=-1, arg("ref")=-1));
@@ -150,8 +156,8 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
.def("GetResNum",&mol::alg::UniqueAtomIdentifier::GetResNum)
.def("GetResidueName",&mol::alg::UniqueAtomIdentifier::GetResidueName)
.def("GetAtomName",&mol::alg::UniqueAtomIdentifier::GetAtomName)
+ .def("GetQualifiedAtomName",&mol::alg::UniqueAtomIdentifier::GetQualifiedAtomName)
;
-
class_<mol::alg::ResidueRDMap>("ResidueRDMap")
.def(map_indexing_suite<mol::alg::ResidueRDMap>())
@@ -177,14 +183,47 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
.def("GetCount",&mol::alg::BondLengthInfo::GetCount)
;
+ class_<mol::alg::ClashEvent> ("ClashEvent" ,init <>())
+ .def(init<const mol::alg::UniqueAtomIdentifier&,const mol::alg::UniqueAtomIdentifier&, Real, Real>())
+ .def("GetFirstAtom",&mol::alg::ClashEvent::GetFirstAtom)
+ .def("GetSecondAtom",&mol::alg::ClashEvent::GetSecondAtom)
+ .def("GetModelDistance",&mol::alg::ClashEvent::GetModelDistance)
+ .def("GetAdjustedReferenceDistance",&mol::alg::ClashEvent::GetAdjustedReferenceDistance)
+ ;
+
+ class_<mol::alg::StereoChemicalBondViolation> ("StereoChemicalBondViolation" ,init <>())
+ .def(init<const mol::alg::UniqueAtomIdentifier&,const mol::alg::UniqueAtomIdentifier&,
+ Real, const std::pair<Real, Real>& >())
+ .def("GetFirstAtom",&mol::alg::StereoChemicalBondViolation::GetFirstAtom)
+ .def("GetSecondAtom",&mol::alg::StereoChemicalBondViolation::GetSecondAtom)
+ .def("GetModelValue",&mol::alg::StereoChemicalBondViolation::GetModelValue)
+ .def("GetAllowedRange",&mol::alg::StereoChemicalBondViolation::GetAllowedRange)
+ ;
+
+ class_<mol::alg::StereoChemicalAngleViolation> ("StereoChemicalAngleViolation" ,init <>())
+ .def(init<const mol::alg::UniqueAtomIdentifier&,const mol::alg::UniqueAtomIdentifier&,
+ const mol::alg::UniqueAtomIdentifier&, Real, const std::pair<Real, Real>& >())
+ .def("GetFirstAtom",&mol::alg::StereoChemicalAngleViolation::GetFirstAtom)
+ .def("GetSecondAtom",&mol::alg::StereoChemicalAngleViolation::GetSecondAtom)
+ .def("GetThirdAtom",&mol::alg::StereoChemicalAngleViolation::GetThirdAtom)
+ .def("GetModelValue",&mol::alg::StereoChemicalAngleViolation::GetModelValue)
+ .def("GetAllowedRange",&mol::alg::StereoChemicalAngleViolation::GetAllowedRange)
+ ;
+
class_<mol::alg::ClashingInfo> ("ClashingInfo" ,init <>())
- .def(init<int,Real>())
+ .def(init<int,Real,const std::vector<mol::alg::ClashEvent> >())
.def("GetClashCount",&mol::alg::ClashingInfo::GetClashCount)
.def("GetAverageOffset",&mol::alg::ClashingInfo::GetAverageOffset)
+ .def("GetClashList",&mol::alg::ClashingInfo::GetClashList)
;
+ to_python_converter<std::pair<mol::EntityView,mol::alg::ClashingInfo>,
+ PairToTupleConverter<mol::EntityView, mol::alg::ClashingInfo> >();
+
class_<mol::alg::StereoChemistryInfo> ("StereoChemistryInfo" ,init <>())
- .def(init<Real,int,int,Real,int,int, const std::map<String,mol::alg::BondLengthInfo>&>())
+ .def(init<Real,int,int,Real,int,int, const std::map<String,mol::alg::BondLengthInfo>&,
+ const std::vector<mol::alg::StereoChemicalBondViolation>&,
+ const std::vector<mol::alg::StereoChemicalAngleViolation>& >())
.def("GetAvgZscoreBonds",&mol::alg::StereoChemistryInfo::GetAvgZscoreBonds)
.def("GetBadBondCount",&mol::alg::StereoChemistryInfo::GetBadBondCount)
.def("GetBondCount",&mol::alg::StereoChemistryInfo::GetBondCount)
@@ -192,8 +231,21 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
.def("GetBadAngleCount",&mol::alg::StereoChemistryInfo::GetBadAngleCount)
.def("GetAngleCount",&mol::alg::StereoChemistryInfo::GetAngleCount)
.def("GetAvgBondLengthInfo",&mol::alg::StereoChemistryInfo::GetAvgBondLengthInfo)
+ .def("GetBondViolationList",&mol::alg::StereoChemistryInfo::GetBondViolationList)
+ .def("GetAngleViolationList",&mol::alg::StereoChemistryInfo::GetAngleViolationList)
;
+ to_python_converter<std::pair<mol::EntityView,mol::alg::StereoChemistryInfo>,
+ PairToTupleConverter<mol::EntityView, mol::alg::StereoChemistryInfo> >();
+
+
+ to_python_converter<std::vector<mol::alg::ClashEvent>,
+ VectorToListConverter<mol::alg::ClashEvent> >();
+
+ to_python_converter<std::vector<mol::alg::StereoChemicalBondViolation>,
+ VectorToListConverter<mol::alg::StereoChemicalBondViolation> >();
+ to_python_converter<std::vector<mol::alg::StereoChemicalAngleViolation>,
+ VectorToListConverter<mol::alg::StereoChemicalAngleViolation> >();
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index eca9aa6e892d9a8bd463be00b1e2f01fc3e59125..68a7a9104b20e60818479a816f47ea95a6094d7e 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -3,6 +3,7 @@ set(OST_MOL_ALG_HEADERS
module_config.hh
sec_structure_segments.hh
local_dist_diff_test.hh
+ distance_test_common.hh
superpose_frames.hh
filter_clashes.hh
construct_cbeta.hh
@@ -17,6 +18,7 @@ set(OST_MOL_ALG_SOURCES
clash_score.cc
sec_structure_segments.cc
local_dist_diff_test.cc
+ distance_test_common.cc
superpose_frames.cc
filter_clashes.cc
construct_cbeta.cc
diff --git a/modules/mol/alg/src/distance_test_common.cc b/modules/mol/alg/src/distance_test_common.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2835205fd7ed8418f6727f84e1adcf521c461dac
--- /dev/null
+++ b/modules/mol/alg/src/distance_test_common.cc
@@ -0,0 +1,81 @@
+#include <ost/mol/alg/distance_test_common.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+
+String UniqueAtomIdentifier::GetQualifiedAtomName() const
+{
+ std::stringstream output_string_s;
+ output_string_s << chain_ << "." << residue_name_ << residue_ << "." << atom_;
+ return (output_string_s.str());
+}
+
+bool UniqueAtomIdentifier::operator==(const UniqueAtomIdentifier& rhs) const
+{
+ return chain_==rhs.chain_ && residue_==rhs.residue_ && atom_==rhs.atom_;
+}
+
+bool UniqueAtomIdentifier::operator<(const UniqueAtomIdentifier& rhs) const
+{
+ int cc=chain_.compare(rhs.chain_);
+ if (cc) {
+ return cc<0;
+ }
+ if (residue_<rhs.residue_) {
+ return true;
+ } else if (residue_>rhs.residue_) {
+ return false;
+ }
+ return atom_.compare(rhs.atom_)<0;
+}
+
+// helper function
+String SwappedName(const String& name)
+{
+ if (name=="OE1") return "OE2";
+ if (name=="OE2") return "OE1";
+
+ if (name=="OD1") return "OD2";
+ if (name=="OD2") return "OD1";
+
+ if (name=="CG1") return "CG2";
+ if (name=="CG2") return "CG1";
+
+ if (name=="CE1") return "CE2";
+ if (name=="CE2") return "CE1";
+
+ if (name=="CD1") return "CD2";
+ if (name=="CD2") return "CD1";
+
+ if (name=="NH1") return "NH2";
+ if (name=="NH2") return "NH1";
+
+ return name;
+}
+
+// helper function
+bool Swappable(const String& rname, const String& aname)
+{
+ if (rname=="GLU") {
+ return (aname=="OE1" || aname=="OE2");
+ }
+ if (rname=="ASP") {
+ return (aname=="OD1" || aname=="OD2");
+ }
+ if (rname=="VAL") {
+
+ return (aname=="CG1" || aname=="CG2");
+ }
+ if (rname=="TYR" || rname=="PHE") {
+ return (aname=="CD1" || aname=="CD2" || aname=="CE1" || aname=="CE2");
+ }
+ if (rname=="LEU") {
+ return (aname=="CD1" || aname=="CD2");
+ }
+ if (rname=="ARG") {
+ return (aname=="NH1" || aname=="NH2");
+ }
+ return false;
+}
+
+}}}
diff --git a/modules/mol/alg/src/distance_test_common.hh b/modules/mol/alg/src/distance_test_common.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ca76a33b5e4c2d06c10b11e15301acdc855146a1
--- /dev/null
+++ b/modules/mol/alg/src/distance_test_common.hh
@@ -0,0 +1,92 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 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_DISTANCE_TEST_COMMON_HH
+#define OST_MOL_ALG_DISTANCE_TEST_COMMON_HH
+
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/alg/module_config.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+String DLLEXPORT_OST_MOL_ALG SwappedName(const String& name);
+bool DLLEXPORT_OST_MOL_ALG Swappable(const String& rname, const String& aname);
+
+/// \brief Contains the infomation needed to uniquely identify an atom in a structure
+///
+/// Used by the the Local Distance Difference Test classes and functions
+class DLLEXPORT_OST_MOL_ALG UniqueAtomIdentifier
+{
+
+public:
+ /// \brief Constructor with all the relevant information
+ UniqueAtomIdentifier(const String& chain,const ResNum& residue,const String& residue_name, const String& atom): chain_(chain),residue_(residue),residue_name_(residue_name),atom_(atom) {}
+
+ // to make the compiler happy (boost python map suite)
+ UniqueAtomIdentifier(): chain_(""),residue_(ResNum(1)),residue_name_(""),atom_("") {}
+
+ /// \brief Returns the name of the chain to which the atom belongs, as a String
+ String GetChainName() const { return chain_; }
+
+ /// \brief Returns the ResNum of the residue to which the atom belongs
+ ResNum GetResNum() const { return residue_; }
+
+ /// \brief Returns the name of the residue to which the atom belongs, as a String
+ String GetResidueName() const { return residue_name_; }
+
+ /// \brief Returns the name of the atom, as a String
+ String GetAtomName() const { return atom_; }
+
+ /// \brief Returns the qualified name of the atom, as a String
+ String GetQualifiedAtomName() const ;
+
+ // required because UniqueAtomIdentifier is used as a key for a std::map
+ bool operator==(const UniqueAtomIdentifier& rhs) const ;
+
+ // required because UniqueAtomIdentifier is used as a key for a std::map
+ bool operator<(const UniqueAtomIdentifier& rhs) const;
+
+private:
+
+ String chain_;
+ ResNum residue_;
+ String residue_name_;
+ String atom_;
+};
+
+// typedef used to make the code cleaner
+typedef std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier> UAtomIdentifiers;
+
+/// \brief Residue distance list.
+///
+/// Container for all the interatomic distances that are checked in a Local Distance Difference Test
+/// and are originating from a single specific residue
+typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<Real,Real> > ResidueRDMap;
+
+/// \brief Global distance list.
+///
+/// Container for all the residue-based interatomic distance lists that are checked in a Local Distance Difference Test
+/// and belong to the same structure
+typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
+
+// used by the multi-reference distance-list generator function
+typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
+
+}}}
+
+#endif
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index e544ce95b724bd0d01372bd18696b7d1ee1d0e3f..875db0aa178c32269868638d82e3b2c273380922 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -82,6 +82,22 @@ String angle_string(const ost::mol::AtomHandle& atom1, const ost::mol::AtomView&
namespace ost { namespace mol { namespace alg {
+std::vector<ClashEvent> ClashingInfo::GetClashList() const
+{
+ std::vector<ClashEvent> return_list;
+ std::vector<UAtomIdentifiers> seen_list;
+ for (std::vector<ClashEvent>::const_iterator cl_it=clash_list_.begin();cl_it!=clash_list_.end();++cl_it) {
+ UniqueAtomIdentifier atom1 = cl_it->GetFirstAtom();
+ UniqueAtomIdentifier atom2 = cl_it->GetSecondAtom();
+ UAtomIdentifiers check = std::make_pair<UniqueAtomIdentifier,UniqueAtomIdentifier>(atom2,atom1);
+ if (std::find(seen_list.begin(),seen_list.end(),check)==seen_list.end()) {
+ return_list.push_back(*cl_it);
+ seen_list.push_back(std::make_pair<UniqueAtomIdentifier,UniqueAtomIdentifier>(atom1,atom2));
+ }
+ }
+ return return_list;
+}
+
void ClashingDistances::SetClashingDistance(const String& ele1,const String& ele2, Real min_distance, Real tolerance)
{
std::stringstream stkey;
@@ -331,6 +347,8 @@ ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_pro
std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb)
{
+ std::vector<StereoChemicalBondViolation> bond_violation_list;
+ std::vector<StereoChemicalAngleViolation> angle_violation_list;
Real running_sum_zscore_bonds=0.0;
Real running_sum_zscore_angles=0.0;
int bond_count = 0;
@@ -382,6 +400,10 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
if (blength < min_length || blength > max_length) {
LOG_INFO("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "FAIL")
bad_bond_count++;
+ UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
+ UniqueAtomIdentifier other_atom_ui(other_atom.GetResidue().GetChain().GetName(),other_atom.GetResidue().GetNumber(),other_atom.GetResidue().GetName(),other_atom.GetName());
+ StereoChemicalBondViolation bond_v(atom_ui,other_atom_ui,blength,std::make_pair<Real,Real>(min_length,max_length));
+ bond_violation_list.push_back(bond_v);
remove_sc=true;
if (always_remove_bb==true) {
remove_bb=true;
@@ -446,7 +468,12 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
Real zscore = (awidth - ref_width)/ref_stddev;
if (awidth < min_width || awidth > max_width) {
LOG_INFO("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "FAIL")
- bad_angle_count++;
+ bad_angle_count++;
+ UniqueAtomIdentifier atom1_ui(atom1.GetResidue().GetChain().GetName(),atom1.GetResidue().GetNumber(),atom1.GetResidue().GetName(),atom1.GetName());
+ UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
+ UniqueAtomIdentifier atom2_ui(atom2.GetResidue().GetChain().GetName(),atom2.GetResidue().GetNumber(),atom2.GetResidue().GetName(),atom2.GetName());
+ StereoChemicalAngleViolation angle_v(atom1_ui,atom_ui,atom2_ui,awidth,std::make_pair<Real,Real>(min_width,max_width));
+ angle_violation_list.push_back(angle_v);
remove_sc=true;
if (always_remove_bb==true) {
remove_bb=true;
@@ -499,7 +526,8 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
avg_bond_length_info[key]=bond_info;
}
StereoChemistryInfo info(avg_zscore_bonds, bad_bond_count, bond_count,avg_zscore_angles,
- bad_angle_count, angle_count,avg_bond_length_info);
+ bad_angle_count, angle_count,avg_bond_length_info,
+ bond_violation_list,angle_violation_list);
filtered.AddAllInclusiveBonds();
return std::make_pair<EntityView,StereoChemistryInfo>(filtered,info);
}
@@ -513,6 +541,7 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityHandl
std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb)
{
+ std::vector<ClashEvent> clash_list;
int distance_count = 0;
int bad_distance_count = 0;
Real average_offset_sum = 0.0;
@@ -564,6 +593,10 @@ std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const Cl
if (d<threshold*threshold) {
LOG_INFO(atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetName() << " " << atom2.GetResidue().GetNumber() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "FAIL")
bad_distance_count++;
+ UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
+ UniqueAtomIdentifier atom2_ui(atom2.GetResidue().GetChain().GetName(),atom2.GetResidue().GetNumber(),atom2.GetResidue().GetName(),atom2.GetName());
+ ClashEvent clash_event(atom_ui,atom2_ui,sqrt(d),threshold);
+ clash_list.push_back(clash_event);
average_offset_sum+=sqrt(d)-threshold;
remove_sc=true;
if (always_remove_bb==true) {
@@ -604,7 +637,7 @@ std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const Cl
if (bad_distance_count!=0) {
average_offset = average_offset_sum / static_cast<Real>(bad_distance_count);
}
- ClashingInfo info(bad_distance_count,average_offset);
+ ClashingInfo info(bad_distance_count,average_offset,clash_list);
filtered.AddAllInclusiveBonds();
return std::make_pair<EntityView,ClashingInfo>(filtered,info);
}
diff --git a/modules/mol/alg/src/filter_clashes.hh b/modules/mol/alg/src/filter_clashes.hh
index d8e1c4c2ea7978536e7e4ad8073bf0e08e1b87a3..a11bd5d8a842581e4f11a049fff4a395d1bad538 100644
--- a/modules/mol/alg/src/filter_clashes.hh
+++ b/modules/mol/alg/src/filter_clashes.hh
@@ -21,7 +21,7 @@
#include <ost/mol/entity_view.hh>
#include <ost/mol/alg/module_config.hh>
-
+#include <ost/mol/alg/distance_test_common.hh>
namespace ost { namespace mol { namespace alg {
class BondLengthInfo
@@ -42,20 +42,81 @@ private:
int count_;
};
+class ClashEvent
+{
+public:
+ ClashEvent(): atom1_(UniqueAtomIdentifier()),atom2_(UniqueAtomIdentifier()),mdl_dist_(0.0),adjusted_ref_dist_(0.0) {}
+ ClashEvent(const UniqueAtomIdentifier& atom1,const UniqueAtomIdentifier& atom2, Real mdl_dist, Real adjusted_ref_dist ):
+ atom1_(atom1),atom2_(atom2),mdl_dist_(mdl_dist),adjusted_ref_dist_(adjusted_ref_dist) {}
+ UniqueAtomIdentifier GetFirstAtom() const { return atom1_; }
+ UniqueAtomIdentifier GetSecondAtom() const { return atom2_; }
+ Real GetModelDistance() const { return mdl_dist_; }
+ Real GetAdjustedReferenceDistance() const { return adjusted_ref_dist_; }
+private:
+ UniqueAtomIdentifier atom1_;
+ UniqueAtomIdentifier atom2_;
+ Real mdl_dist_;
+ Real adjusted_ref_dist_;
+};
class ClashingInfo
{
public:
- ClashingInfo(): clash_count_(0), average_offset_ (0) {}
- ClashingInfo (int clash_count, Real average_offset):
- clash_count_(clash_count), average_offset_ (average_offset) {}
- int GetClashCount() {return clash_count_;}
- Real GetAverageOffset() {return average_offset_;}
+ ClashingInfo(): clash_count_(0), average_offset_ (0), clash_list_(std::vector<ClashEvent>()) {}
+ ClashingInfo (int clash_count, Real average_offset, const std::vector<ClashEvent>& clash_list):
+ clash_count_(clash_count), average_offset_ (average_offset),clash_list_(clash_list) {}
+ int GetClashCount() const {return clash_count_/2.0;}
+ Real GetAverageOffset() const {return average_offset_;}
+ std::vector<ClashEvent> GetClashList() const;
private:
int clash_count_;
Real average_offset_;
+ std::vector<ClashEvent> clash_list_;
+};
+
+class StereoChemicalBondViolation
+{
+public:
+ StereoChemicalBondViolation():
+ atom1_(UniqueAtomIdentifier()),atom2_(UniqueAtomIdentifier()),mdl_value_(0.0),allowed_range_(std::pair<Real,Real>(0.0,0.0)) {}
+ StereoChemicalBondViolation(const UniqueAtomIdentifier& atom1,
+ const UniqueAtomIdentifier& atom2,
+ Real mdl_value, std::pair<Real,Real> allowed_range ):
+ atom1_(atom1),atom2_(atom2),mdl_value_(mdl_value),allowed_range_(allowed_range) {}
+ UniqueAtomIdentifier GetFirstAtom() const { return atom1_; }
+ UniqueAtomIdentifier GetSecondAtom() const { return atom2_; }
+ Real GetModelValue() const { return mdl_value_; }
+ std::pair<Real,Real> GetAllowedRange() const { return allowed_range_; }
+private:
+ UniqueAtomIdentifier atom1_;
+ UniqueAtomIdentifier atom2_;
+ Real mdl_value_;
+ std::pair<Real,Real> allowed_range_;
+};
+
+class StereoChemicalAngleViolation
+{
+public:
+ StereoChemicalAngleViolation():
+ atom1_(UniqueAtomIdentifier()),atom2_(UniqueAtomIdentifier()),atom3_(UniqueAtomIdentifier()),mdl_value_(0.0),allowed_range_(std::pair<Real,Real>(0.0,0.0)) {}
+ StereoChemicalAngleViolation(const UniqueAtomIdentifier& atom1,
+ const UniqueAtomIdentifier& atom2,
+ const UniqueAtomIdentifier& atom3,
+ Real mdl_value, std::pair<Real,Real> allowed_range ):
+ atom1_(atom1),atom2_(atom2),atom3_(atom3),mdl_value_(mdl_value),allowed_range_(allowed_range) {}
+ UniqueAtomIdentifier GetFirstAtom() const { return atom1_; }
+ UniqueAtomIdentifier GetSecondAtom() const { return atom2_; }
+ UniqueAtomIdentifier GetThirdAtom() const { return atom3_; }
+ Real GetModelValue() const { return mdl_value_; }
+ std::pair<Real,Real> GetAllowedRange() const { return allowed_range_; }
+private:
+ UniqueAtomIdentifier atom1_;
+ UniqueAtomIdentifier atom2_;
+ UniqueAtomIdentifier atom3_;
+ Real mdl_value_;
+ std::pair<Real,Real> allowed_range_;
};
class StereoChemistryInfo
@@ -68,24 +129,33 @@ public:
avg_zscore_angles_(0),
bad_angle_count_(0),
angle_count_(0),
- avg_bond_length_info_(std::map<String,BondLengthInfo>()) {}
+ avg_bond_length_info_(std::map<String,BondLengthInfo>()),
+ bond_violation_list_(std::vector<StereoChemicalBondViolation>()),
+ angle_violation_list_(std::vector<StereoChemicalAngleViolation>()) {}
StereoChemistryInfo(Real avg_zscore_bonds, int bad_bond_count, int bond_count,
Real avg_zscore_angles, int bad_angle_count, int angle_count,
- const std::map<String,BondLengthInfo>& avg_bond_length_info):
+ const std::map<String,BondLengthInfo>& avg_bond_length_info,
+ const std::vector<StereoChemicalBondViolation>& bond_violation_list,
+ const std::vector<StereoChemicalAngleViolation>& angle_violation_list):
avg_zscore_bonds_(avg_zscore_bonds),
bad_bond_count_(bad_bond_count),
bond_count_(bond_count),
avg_zscore_angles_(avg_zscore_angles),
bad_angle_count_(bad_angle_count),
angle_count_(angle_count),
- avg_bond_length_info_(avg_bond_length_info) {}
- Real GetAvgZscoreBonds() {return avg_zscore_bonds_;}
- int GetBadBondCount() {return bad_bond_count_;}
- int GetBondCount() {return bond_count_;}
- Real GetAvgZscoreAngles() {return avg_zscore_angles_;}
- int GetBadAngleCount() {return bad_angle_count_;}
- int GetAngleCount() {return angle_count_;}
+ avg_bond_length_info_(avg_bond_length_info),
+ bond_violation_list_(bond_violation_list),
+ angle_violation_list_(angle_violation_list) {}
+ Real GetAvgZscoreBonds() const {return avg_zscore_bonds_;}
+ int GetBadBondCount() const {return bad_bond_count_;}
+ int GetBondCount() const {return bond_count_;}
+ Real GetAvgZscoreAngles() const {return avg_zscore_angles_;}
+ int GetBadAngleCount() const {return bad_angle_count_;}
+ int GetAngleCount() const {return angle_count_;}
std::map<String,BondLengthInfo> GetAvgBondLengthInfo() {return avg_bond_length_info_;}
+ std::vector<StereoChemicalBondViolation> GetBondViolationList() { return bond_violation_list_; }
+ std::vector<StereoChemicalAngleViolation> GetAngleViolationList() { return angle_violation_list_; }
+
private:
Real avg_zscore_bonds_;
@@ -95,6 +165,8 @@ private:
int bad_angle_count_;
int angle_count_;
std::map<String,BondLengthInfo> avg_bond_length_info_;
+ std::vector<StereoChemicalBondViolation> bond_violation_list_;
+ std::vector<StereoChemicalAngleViolation> angle_violation_list_;
};
diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc
index 8f936411652c58d407fe70706f17c09e6d087a05..2bf01896c77f29d48306946ad850027d02083a89 100644
--- a/modules/mol/alg/src/local_dist_diff_test.cc
+++ b/modules/mol/alg/src/local_dist_diff_test.cc
@@ -7,55 +7,6 @@ namespace ost { namespace mol { namespace alg {
namespace {
-// helper function
-String swapped_name(const String& name)
-{
- if (name=="OE1") return "OE2";
- if (name=="OE2") return "OE1";
-
- if (name=="OD1") return "OD2";
- if (name=="OD2") return "OD1";
-
- if (name=="CG1") return "CG2";
- if (name=="CG2") return "CG1";
-
- if (name=="CE1") return "CE2";
- if (name=="CE2") return "CE1";
-
- if (name=="CD1") return "CD2";
- if (name=="CD2") return "CD1";
-
- if (name=="NH1") return "NH2";
- if (name=="NH2") return "NH1";
-
- return name;
-}
-
-// helper function
-bool swappable(const String& rname, const String& aname)
-{
- if (rname=="GLU") {
- return (aname=="OE1" || aname=="OE2");
- }
- if (rname=="ASP") {
- return (aname=="OD1" || aname=="OD2");
- }
- if (rname=="VAL") {
-
- return (aname=="CG1" || aname=="CG2");
- }
- if (rname=="TYR" || rname=="PHE") {
- return (aname=="CD1" || aname=="CD2" || aname=="CE1" || aname=="CE2");
- }
- if (rname=="LEU") {
- return (aname=="CD1" || aname=="CD2");
- }
- if (rname=="ARG") {
- return (aname=="NH1" || aname=="NH2");
- }
- return false;
-}
-
// helper function
bool within_tolerance(Real mdl_dist, const std::pair<Real,Real>& values, Real tol)
{
@@ -76,14 +27,14 @@ std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_lis
const std::pair <Real,Real>& values = ai->second;
const UniqueAtomIdentifier& first_atom=uais.first;
const UniqueAtomIdentifier& second_atom=uais.second;
- String name=swap ? swapped_name(first_atom.GetAtomName()) : first_atom.GetAtomName();
+ String name=swap ? SwappedName(first_atom.GetAtomName()) : first_atom.GetAtomName();
AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
if (only_fixed) {
if (std::abs(first_atom.GetResNum().GetNum()-second_atom.GetResNum().GetNum())<=sequence_separation) {
continue;
}
- if (swappable(second_atom.GetResidueName(), second_atom.GetAtomName())) {
+ if (Swappable(second_atom.GetResidueName(), second_atom.GetAtomName())) {
continue;
}
}
@@ -157,7 +108,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
for (AtomViewList::iterator ai=ref_atoms.begin(),
ae=ref_atoms.end(); ai!=ae; ++ai) {
if (ai->GetElement()=="H") { continue; }
- String name=swap ? swapped_name(ai->GetName()) : ai->GetName();
+ String name=swap ? SwappedName(ai->GetName()) : ai->GetName();
AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
if (max_dist>=0){
within=ref.FindWithin(ai->GetPos(), max_dist);
@@ -172,7 +123,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
if (aj->GetResidue().GetNumber()==ref_res.GetNumber()) {
continue;
}
- if (swappable(aj->GetResidue().GetName(), aj->GetName())) {
+ if (Swappable(aj->GetResidue().GetName(), aj->GetName())) {
continue;
}
overlap.second+=1.0;
@@ -254,8 +205,8 @@ void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, st
static_cast<Real>(ov2.first)/ov2.second) {
AtomViewList atoms=mdl_res.GetAtomList();
for (AtomViewList::iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
- if (swappable(rname, j->GetName())) {
- edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
+ if (Swappable(rname, j->GetName())) {
+ edi.RenameAtom(j->GetHandle(), SwappedName(j->GetName()));
}
}
}
@@ -557,8 +508,8 @@ Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln,
AtomViewList atoms=mdl_res.GetAtomList();
for (AtomViewList::iterator j=atoms.begin(),
e2=atoms.end(); j!=e2; ++j) {
- if (swappable(rname, j->GetName())) {
- edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
+ if (Swappable(rname, j->GetName())) {
+ edi.RenameAtom(j->GetHandle(), SwappedName(j->GetName()));
}
}
}
diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh
index 9803f0a723d3c77ae8b1332ef5ab09e444ae5ff1..6f0d2c20f7530e8ba98967612c7d7906cfc547b7 100644
--- a/modules/mol/alg/src/local_dist_diff_test.hh
+++ b/modules/mol/alg/src/local_dist_diff_test.hh
@@ -19,81 +19,12 @@
#ifndef OST_MOL_ALG_LOCAL_DIST_TEST_HH
#define OST_MOL_ALG_LOCAL_DIST_TEST_HH
-#include <ost/mol/entity_view.hh>
#include <ost/mol/alg/module_config.hh>
#include <ost/seq/alignment_handle.hh>
+#include <ost/mol/alg/distance_test_common.hh>
namespace ost { namespace mol { namespace alg {
-/// \brief Contains the infomation needed to uniquely identify an atom in a structure
-///
-/// Used by the the Local Distance Difference Test classes and functions
-class DLLEXPORT_OST_MOL_ALG UniqueAtomIdentifier
-{
-
-public:
- /// \brief Constructor with all the relevant information
- UniqueAtomIdentifier(const String& chain,const ResNum& residue,const String& residue_name, const String& atom): chain_(chain),residue_(residue),residue_name_(residue_name),atom_(atom) {}
-
- // to make the compiler happy (boost python map suite)
- UniqueAtomIdentifier(): chain_(""),residue_(ResNum(1)),residue_name_(""),atom_("") {}
-
- /// \brief Returns the name of the chain to which the atom belongs, as a String
- String GetChainName() const { return chain_; }
-
- /// \brief Returns the ResNum of the residue to which the atom belongs
- ResNum GetResNum() const { return residue_; }
-
- /// \brief Returns the name of the residue to which the atom belongs, as a String
- String GetResidueName() const { return residue_name_; }
-
- /// \brief Returns the name of the atom, as a String
- String GetAtomName() const { return atom_; }
-
- // required because UniqueAtomIdentifier is used as a key for a std::map
- bool operator==(const UniqueAtomIdentifier& rhs) const {
- return chain_==rhs.chain_ && residue_==rhs.residue_ && atom_==rhs.atom_;
- }
-
- // required because UniqueAtomIdentifier is used as a key for a std::map
- bool operator<(const UniqueAtomIdentifier& rhs) const {
- int cc=chain_.compare(rhs.chain_);
- if (cc) {
- return cc<0;
- }
- if (residue_<rhs.residue_) {
- return true;
- } else if (residue_>rhs.residue_) {
- return false;
- }
- return atom_.compare(rhs.atom_)<0;
- }
-private:
-
- String chain_;
- ResNum residue_;
- String residue_name_;
- String atom_;
-};
-
-// typedef used to make the code cleaner
-typedef std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier> UAtomIdentifiers;
-
-/// \brief Residue distance list.
-///
-/// Container for all the interatomic distances that are checked in a Local Distance Difference Test
-/// and are originating from a single specific residue
-typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<Real,Real> > ResidueRDMap;
-
-/// \brief Global distance list.
-///
-/// Container for all the residue-based interatomic distance lists that are checked in a Local Distance Difference Test
-/// and belong to the same structure
-typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
-
-// used by the multi-reference distance-list generator function
-typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
-
/// \brief Calculates number of distances conserved in a model, given a list of distances to check and a model
///
/// Calculates the two values needed to determine the Local Distance Difference Test for a given model, i.e.