From ae3c5c0d480641fad81c3e4cb2ba58b5a130db37 Mon Sep 17 00:00:00 2001
From: Marco Biasini <marco.biasini@unibas.ch>
Date: Thu, 26 Jul 2012 12:17:51 +0200
Subject: [PATCH] added RuleBasedConop
RuleBased conop is an experimental new interface for processing
structures in the conop module. The class exposes a single method
Process with takes care of the complete conop procedure. There are
still a few methods missing, they will be added shortly.
---
CMakeLists.txt | 15 ++
modules/conop/pymod/CMakeLists.txt | 1 +
modules/conop/pymod/export_rule_based.cc | 12 ++
modules/conop/pymod/wrap_conop.cc | 2 +
modules/conop/src/CMakeLists.txt | 2 +
modules/conop/src/rule_based.cc | 210 +++++++++++++++++++++++
modules/conop/src/rule_based.hh | 32 ++++
modules/conop/src/rule_based_builder.cc | 1 +
8 files changed, 275 insertions(+)
create mode 100644 modules/conop/pymod/export_rule_based.cc
create mode 100644 modules/conop/src/rule_based.cc
create mode 100644 modules/conop/src/rule_based.hh
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 91c10973b..fba8fd543 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -189,6 +189,21 @@ if (ENABLE_INFO)
find_package(Qt4 4.5.0 REQUIRED)
endif()
+if (OPTIMIZE)
+ if (CXX_FLAGS)
+ set (CMAKE_CXX_FLAGS_RELEASE ${CXX_FLAGS})
+ endif()
+ if (C_FLAGS)
+ set (CMAKE_C_FLAGS_RELEASE ${C_FLAGS})
+ endif()
+else()
+ if (CXX_FLAGS)
+ set(CMAKE_CXX_FLAGS_DEBUG ${CXX_FLAGS})
+ endif()
+ if (C_FLAGS)
+ set(CMAKE_C_FLAGS_DEBUG ${C_FLAGS})
+ endif()
+endif()
if (ENABLE_GFX)
if(USE_MESA)
find_package(Mesa REQUIRED)
diff --git a/modules/conop/pymod/CMakeLists.txt b/modules/conop/pymod/CMakeLists.txt
index c1bf89d0e..2ceff8be3 100644
--- a/modules/conop/pymod/CMakeLists.txt
+++ b/modules/conop/pymod/CMakeLists.txt
@@ -4,6 +4,7 @@ set(OST_CONOP_PYMOD_SOURCES
export_compound.cc
export_amino_acids.cc
export_conop.cc
+ export_rule_based.cc
export_non_standard.cc
export_ring_finder.cc
)
diff --git a/modules/conop/pymod/export_rule_based.cc b/modules/conop/pymod/export_rule_based.cc
new file mode 100644
index 000000000..74112149b
--- /dev/null
+++ b/modules/conop/pymod/export_rule_based.cc
@@ -0,0 +1,12 @@
+#include <boost/python.hpp>
+#include <ost/conop/rule_based.hh>
+
+using namespace boost::python;
+using namespace ost::conop;
+
+void export_rule_based()
+{
+ class_<RuleBasedConop>("RuleBasedConop", init<CompoundLibPtr>())
+ .def("Process", &RuleBasedConop::Process)
+ ;
+}
diff --git a/modules/conop/pymod/wrap_conop.cc b/modules/conop/pymod/wrap_conop.cc
index 6573bf41a..99ceec316 100644
--- a/modules/conop/pymod/wrap_conop.cc
+++ b/modules/conop/pymod/wrap_conop.cc
@@ -27,10 +27,12 @@ void export_RingFinder();
void export_AminoAcids();
void export_NonStandard();
+void export_rule_based();
BOOST_PYTHON_MODULE(_ost_conop)
{
export_Builder();
export_Conop();
+ export_rule_based();
export_Compound();
export_RingFinder();
export_AminoAcids();
diff --git a/modules/conop/src/CMakeLists.txt b/modules/conop/src/CMakeLists.txt
index 8f493e83e..55eb837f9 100644
--- a/modules/conop/src/CMakeLists.txt
+++ b/modules/conop/src/CMakeLists.txt
@@ -9,6 +9,7 @@ model_check.hh
compound.hh
compound_lib.hh
module_config.hh
+rule_based.hh
nonstandard.hh
rule_based_builder.hh
ring_finder.hh
@@ -19,6 +20,7 @@ builder.cc
amino_acids.cc
conop.cc
diag.cc
+rule_based.cc
model_check.cc
heuristic_builder.cc
compound.cc
diff --git a/modules/conop/src/rule_based.cc b/modules/conop/src/rule_based.cc
new file mode 100644
index 000000000..941ecade4
--- /dev/null
+++ b/modules/conop/src/rule_based.cc
@@ -0,0 +1,210 @@
+#include <ost/log.hh>
+#include <ost/mol/xcs_editor.hh>
+#include <ost/mol/bond_handle.hh>
+#include <ost/mol/impl/residue_impl.hh>
+#include <ost/mol/impl/atom_impl.hh>
+#include <ost/mol/residue_handle.hh>
+#include "rule_based.hh"
+namespace ost { namespace conop {
+
+struct OrdinalAtomComp {
+ bool operator()(const mol::impl::AtomImplPtr& a,
+ const mol::impl::AtomImplPtr& b) const {
+ return a->GetState()<b->GetState();
+ }
+};
+
+void RuleBasedConop::Process(mol::EntityHandle ent)
+{
+ mol::ChainHandleList chains=ent.GetChainList();
+ for (mol::ChainHandleList::iterator
+ i = chains.begin(), e = chains.end(); i!= e; ++i) {
+ mol::ChainHandle& chain = *i;
+ mol::ResidueHandleList residues = chain.GetResidueList();
+ mol::ResidueHandle prev;
+ for (mol::ResidueHandleList::iterator
+ j = residues.begin(), e2 = residues.end(); j != e2; ++j) {
+ mol::ResidueHandle residue = *j;
+ CompoundPtr compound = lib_->FindCompound(residue.GetName(), Compound::PDB);
+ if (!compound) {
+ // process unknown residue...
+ continue;
+ }
+ this->ReorderAtoms(residue, compound);
+ bool unks = this->HasUnknownAtoms(residue, compound);
+ if (unks) {
+ LOG_WARNING("residue " << residue << " doesn't look like a standard "
+ << residue.GetKey() << " (" << compound->GetFormula() << ")");
+ residue.SetChemClass(mol::ChemClass(mol::ChemClass::UNKNOWN));
+ residue.SetChemType(mol::ChemType(mol::ChemType::UNKNOWN));
+ residue.SetOneLetterCode('?');
+ continue;
+ }
+ this->FillResidueProps(residue, compound);
+ this->ConnectAtomsOfResidue(residue, compound);
+ this->ConnectResidues(prev, residue);
+ prev = residue;
+ }
+ }
+}
+
+void RuleBasedConop::ConnectResidues(mol::ResidueHandle rh,
+ mol::ResidueHandle next)
+{
+ if (!next.IsValid() || !rh.IsValid()) {
+ return;
+ }
+
+ mol::XCSEditor e=rh.GetEntity().EditXCS(mol::BUFFERED_EDIT);
+
+ // check if both of the residues are able to form a peptide bond.
+ if (rh.IsPeptideLinking() && next.IsPeptideLinking()) {
+ // If we have an OXT then there is no peptide bond connecting the two
+ // residues.
+ if (rh.FindAtom("OXT"))
+ return;
+ mol::AtomHandle c=rh.FindAtom("C");
+ mol::AtomHandle n=next.FindAtom("N");
+ // Give subclasses a chance to give us their opinions on the feasibility of
+ // the peptide bond.
+ if (c.IsValid() && n.IsValid() && this->IsBondFeasible(c, n)) {
+ e.Connect(c, n, 1);
+ rh.SetIsProtein(true);
+ next.SetIsProtein(true);
+ }
+ } else if (rh.IsNucleotideLinking() && next.IsNucleotideLinking()) {
+ mol::AtomHandle c=rh.FindAtom("O3'");
+ mol::AtomHandle n=next.FindAtom("P");
+ if (c.IsValid() && n.IsValid() && this->IsBondFeasible(c, n)) {
+ e.Connect(c, n, 1);
+ }
+ }
+}
+
+void RuleBasedConop::FillResidueProps(mol::ResidueHandle residue,
+ CompoundPtr compound)
+{
+ residue.SetChemClass(compound->GetChemClass());
+ residue.SetChemType(compound->GetChemType());
+ residue.SetOneLetterCode(compound->GetOneLetterCode());
+}
+
+mol::AtomHandle RuleBasedConop::LocateAtom(const mol::AtomHandleList& ahl,
+ int ordinal)
+{
+ if (ahl.empty())
+ return mol::AtomHandle();
+ const mol::AtomHandle* r_it=&ahl.back();
+ if (static_cast<int>(ahl.size())>ordinal) {
+ r_it=&ahl.front()+ordinal;
+ }
+ while ((r_it>=&ahl.front()) &&
+ (static_cast<int>(r_it->Impl()->GetState())>ordinal)) {
+ --r_it;
+ }
+ bool not_found=(r_it<&ahl.front() ||
+ static_cast<int>(r_it->Impl()->GetState())!=ordinal);
+ return not_found ? mol::AtomHandle() : *r_it;
+}
+
+
+void RuleBasedConop::ConnectAtomsOfResidue(mol::ResidueHandle rh, CompoundPtr compound)
+{
+
+ //if (!compound) {
+ // dist_connect(this, rh.GetAtomList());
+ // return;
+ //}
+ //if (unknown_atoms_) {
+ // dist_connect(this, rh.GetAtomList());
+ // return;
+ //}
+ mol::XCSEditor e=rh.GetEntity().EditXCS(mol::BUFFERED_EDIT);
+ BondSpecList::const_iterator j=compound->GetBondSpecs().begin();
+ mol::AtomHandleList atoms=rh.GetAtomList();
+ for(; j!=compound->GetBondSpecs().end(); ++j) {
+ const BondSpec& bond=*j;
+ mol::AtomHandle a1=this->LocateAtom(atoms, bond.atom_one);
+ mol::AtomHandle a2=this->LocateAtom(atoms, bond.atom_two);
+ if (a1.IsValid() && a2.IsValid()) {
+ if (!check_bond_feasibility_) {
+ if (strict_hydrogens_ && (a1.GetElement()=="H" ||
+ a2.GetElement()=="D")) {
+ continue;
+ }
+ e.Connect(a1, a2, bond.order);
+ } else {
+ if (this->IsBondFeasible(a1, a2)) {
+ if (strict_hydrogens_ && (a1.GetElement()=="H" ||
+ a2.GetElement()=="D")) {
+ continue;
+ }
+ e.Connect(a1, a2, bond.order);
+ }
+ }
+ }
+ }
+ //for (mol::AtomHandleList::iterator i=atoms.begin(), e=atoms.end(); i!=e; ++i) {
+ // if (((*i).GetElement()=="H" || (*i).GetElement()=="D") &&
+ // (*i).GetBondCount()==0) {
+ // this->DistanceBasedConnect(*i);
+ // }
+ // }
+}
+
+bool RuleBasedConop::IsBondFeasible(const mol::AtomHandle& atom_a,
+ const mol::AtomHandle& atom_b)
+{
+ Real radii=0.0;
+ if (atom_a.GetRadius()>0.0) {
+ radii=atom_a.GetRadius();
+ } else {
+ return false;
+ }
+ if (atom_b.GetRadius()>0.0) {
+ radii+=atom_b.GetRadius();
+ } else {
+ return false;
+ }
+ Real len=geom::Length2(atom_a.GetPos()-atom_b.GetPos());
+ Real lower_bound=radii*radii*0.0625;
+ Real upper_bound=lower_bound*6.0;
+ return (len<=upper_bound && len>=lower_bound);
+}
+void RuleBasedConop::ReorderAtoms(mol::ResidueHandle residue,
+ CompoundPtr compound)
+{
+ mol::impl::ResidueImplPtr impl=residue.Impl();
+ mol::impl::AtomImplList::iterator i=impl->GetAtomList().begin();
+ for (; i!=impl->GetAtomList().end(); ++i) {
+ mol::impl::AtomImplPtr atom=*i;
+ atom->SetState(std::numeric_limits<unsigned int>::max());
+ int index=compound->GetAtomSpecIndex(atom->GetName());
+ if (index==-1) {
+ atom->SetState(std::numeric_limits<unsigned int>::max());
+ continue;
+ }
+ atom->SetState((compound->GetAtomSpecs())[index].ordinal);
+ }
+ std::sort(impl->GetAtomList().begin(), impl->GetAtomList().end(),
+ OrdinalAtomComp());
+}
+
+bool RuleBasedConop::HasUnknownAtoms(mol::ResidueHandle res, CompoundPtr compound)
+{
+ AtomSpecList::const_iterator j=compound->GetAtomSpecs().begin();
+ mol::AtomHandleList atoms=res.GetAtomList();
+ mol::AtomHandleList::iterator i=atoms.begin();
+ for (mol::AtomHandleList::iterator
+ i=atoms.begin(), e=atoms.end(); i!=e; ++i) {
+ if ((*i).Impl()->GetState()==std::numeric_limits<unsigned int>::max()) {
+ if (((*i).GetElement()=="H" || (*i).GetElement()=="D") &&
+ !strict_hydrogens_) {
+ continue;
+ }
+ return true;
+ }
+ }
+ return false;
+}
+}}
diff --git a/modules/conop/src/rule_based.hh b/modules/conop/src/rule_based.hh
new file mode 100644
index 000000000..2a8c95b8d
--- /dev/null
+++ b/modules/conop/src/rule_based.hh
@@ -0,0 +1,32 @@
+#ifndef OST_CONOP_RULE_BASED_HH
+#define OST_CONOP_RULE_BASED_HH
+
+#include <ost/mol/entity_handle.hh>
+#include "compound_lib.hh"
+#include "diag.hh"
+namespace ost { namespace conop {
+
+
+class DLLEXPORT_OST_CONOP RuleBasedConop {
+public:
+ RuleBasedConop(CompoundLibPtr compound_lib): lib_(compound_lib),
+ strict_hydrogens_(false), check_bond_feasibility_(false) {}
+
+ void Process(mol::EntityHandle ent);
+
+private:
+ bool HasUnknownAtoms(mol::ResidueHandle residue, CompoundPtr compound);
+ void ReorderAtoms(mol::ResidueHandle residue, CompoundPtr compound);
+ void FillResidueProps(mol::ResidueHandle residue, CompoundPtr compound);
+ void ConnectAtomsOfResidue(mol::ResidueHandle residue, CompoundPtr compound);
+ void ConnectResidues(mol::ResidueHandle residue, mol::ResidueHandle next);
+ bool IsBondFeasible(const mol::AtomHandle&, const mol::AtomHandle&);
+ mol::AtomHandle LocateAtom(const mol::AtomHandleList&, int ordinal);
+ CompoundLibPtr lib_;
+ bool strict_hydrogens_;
+ bool check_bond_feasibility_;
+};
+
+}}
+#endif
+
diff --git a/modules/conop/src/rule_based_builder.cc b/modules/conop/src/rule_based_builder.cc
index 529fbcbc6..53410e52d 100644
--- a/modules/conop/src/rule_based_builder.cc
+++ b/modules/conop/src/rule_based_builder.cc
@@ -216,6 +216,7 @@ inline void dist_connect(RuleBasedBuilder* b, mol::AtomHandleList atoms)
b->DistanceBasedConnect(*i);
}
}
+
void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh)
{
--
GitLab