diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 67f4fbf9e13454d105b0f69844297b06a1508309..91a9a0904041a90bf502a8c4ecb00e7cbe99c7c0 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -17,14 +17,12 @@
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
-/*
- * Author Juergen Haas
- */
#include <boost/python.hpp>
#include <ost/config.hh>
#include <ost/mol/alg/local_dist_test.hh>
#include <ost/mol/alg/superpose_frames.hh>
using namespace boost::python;
+using namespace ost;
void export_svdSuperPose();
@@ -32,6 +30,12 @@ void export_svdSuperPose();
void export_entity_to_density();
#endif
+namespace {
+
+Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real)=&mol::alg::LocalDistTest;
+Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
+}
+
BOOST_PYTHON_MODULE(_mol_alg)
{
export_svdSuperPose();
@@ -39,8 +43,9 @@ BOOST_PYTHON_MODULE(_mol_alg)
export_entity_to_density();
#endif
- def("LocalDistTest", &ost::mol::alg::LocalDistTest);
+ def("LocalDistTest", ldt_a);
+ def("LocalDistTest", ldt_b, (arg("ref_index")=0, arg("mdl_index")=1));
def("SuperposeFrames", &ost::mol::alg::SuperposeFrames,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
- arg("end")=-1, arg("ref")=-1));
+ arg("end")=-1, arg("ref")=-1));
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index fef02349e2e5031cad63583693692230bd181664..15e93c0e72ade2450464aaa55551828073dc1ce3 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -13,7 +13,7 @@ set(OST_MOL_ALG_SOURCES
superpose_frames.cc
)
-set(MOL_ALG_DEPS mol)
+set(MOL_ALG_DEPS mol seq)
if (ENABLE_IMG)
set(OST_MOL_ALG_HEADERS
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
index 2d97aa915587805c6729327d22320e7db25c7dc3..de78df711c03dcfaf619de243c193b0bd70576f8 100644
--- a/modules/mol/alg/src/local_dist_test.cc
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -49,11 +49,11 @@ bool swappable(const String& rname, const String& aname)
return false;
}
-std::pair<Real, Real> calc_overlap(ResidueView ref_res,
- EntityView ref,
- ChainView mdl_chain,
- Real tol, Real max_dist, bool only_fixed,
- bool swap)
+std::pair<Real, Real> calc_overlap1(ResidueView ref_res,
+ EntityView ref,
+ ChainView mdl_chain,
+ Real tol, Real max_dist, bool only_fixed,
+ bool swap)
{
std::pair<Real, Real> overlap(0.0, 0.0);
AtomViewList ref_atoms=ref_res.GetAtomList();
@@ -82,7 +82,7 @@ std::pair<Real, Real> calc_overlap(ResidueView ref_res,
continue;
}
AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(),
- aj->GetName());
+ aj->GetName());
overlap.second+=1.0;
if (!(av1 && av2)) {
continue;
@@ -112,20 +112,101 @@ std::pair<Real, Real> calc_overlap(ResidueView ref_res,
return overlap;
}
+std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
+ const seq::ConstSequenceHandle& mdl_seq,
+ int pos, Real tol, Real max_dist,
+ bool only_fixed, bool swap)
+{
+ std::pair<Real, Real> overlap(0.0, 0.0);
+ EntityView ref=ref_seq.GetAttachedView();
+ ResidueView ref_res=ref_seq.GetResidue(pos);
+ if (!ref_res.IsValid()) {
+ return std::pair<Real,Real>(0.0, 0.0);
+ }
+ AtomViewList ref_atoms=ref_res.GetAtomList();
+ ResidueView mdl_res=mdl_seq.GetResidue(pos);
+ AtomViewList within;
+ if (max_dist<0) {
+ within=ref.GetAtomList();
+ }
+ 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();
+ AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
+ if (max_dist>=0){
+ within=ref.FindWithin(ai->GetPos(), max_dist);
+ }
+ for (AtomViewList::iterator aj=within.begin(),
+ ae2=within.end(); aj!=ae2; ++aj) {
+ if (aj->GetElement()=="H") { continue; }
+ if (only_fixed) {
+ if (aj->GetResidue().GetNumber()==ref_res.GetNumber()) {
+ continue;
+ }
+ if (swappable(aj->GetResidue().GetName(), aj->GetName())) {
+ continue;
+ }
+ overlap.second+=1.0;
+ // map from residue index to position in alignment
+ try {
+ int aln_pos=ref_seq.GetPos(aj->GetResidue().GetIndex());
+ ResidueView r2=mdl_seq.GetResidue(aln_pos);
+ if (!r2.IsValid()) {
+ continue;
+ }
+ AtomView av2=r2.FindAtom(aj->GetName());
+ if (!(av1 && av2)) {
+ continue;
+ }
+ Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+ Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
+ if (std::abs(mdl_dist-ref_dist)<tol) {
+ overlap.first+=1;
+ }
+ } catch(...) { }
+ continue;
+ } else {
+ if (aj->GetResidue().GetNumber()>ref_res.GetNumber()) {
+ overlap.second+=1.0;
+ try {
+ int aln_pos=ref_seq.GetPos(aj->GetResidue().GetIndex());
+ ResidueView r2=mdl_seq.GetResidue(aln_pos);
+ if (!r2.IsValid()) {
+ continue;
+ }
+ AtomView av2=r2.FindAtom(aj->GetName());
+ if (!(av1 && av2)) {
+ continue;
+ }
+ Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+ Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
+ if (std::abs(mdl_dist-ref_dist)<tol) {
+ overlap.first+=1;
+ }
+ } catch (...) { }
+ }
+ }
+ }
+ }
+ return overlap;
+}
+
}
Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
Real cutoff, Real max_dist)
{
- ResidueViewList ref_residues=ref.GetResidueList();
if (!mdl.GetResidueCount()) {
- LOG_WARNING("model structures doesn't contain any residues");
- return 0.0;
+ LOG_WARNING("model structures doesn't contain any residues");
+ return 0.0;
}
- if (ref_residues.empty()) {
- LOG_WARNING("reference structures doesn't contain any residues");
- return 0.0;
+
+ if (!ref.GetResidueCount()) {
+ LOG_WARNING("reference structures doesn't contain any residues");
+ return 0.0;
}
+ ResidueViewList ref_residues=ref.GetResidueList();
ChainView mdl_chain=mdl.GetChainList()[0];
// Residues with symmetric side-chains require special treatment as there are
// two possible ways to name the atoms. We calculate the overlap with the
@@ -143,10 +224,10 @@ Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
continue;
}
- std::pair<Real, Real> ov1=calc_overlap(*i, ref, mdl_chain,
+ std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain,
cutoff, max_dist, true,
false);
- std::pair<Real, Real> ov2=calc_overlap(*i, ref, mdl_chain,
+ std::pair<Real, Real> ov2=calc_overlap1(*i, ref, mdl_chain,
cutoff, max_dist, true,
true);
if (ov1.first/ov1.second<ov2.first/ov2.second) {
@@ -157,18 +238,68 @@ Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
}
}
-
}
}
std::pair<Real, Real> total_ov(0.0, 0.0);
for (ResidueViewList::iterator
i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
- std::pair<Real, Real> ov1=calc_overlap(*i, ref, mdl_chain, cutoff,
+ std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain, cutoff,
max_dist, false, false);
total_ov.first+=ov1.first;
total_ov.second+=ov1.second;
}
- return total_ov.first/total_ov.second;
+ return total_ov.first/(total_ov.second ? total_ov.second : 1);
+}
+
+Real LocalDistTest(const ost::seq::AlignmentHandle& aln,
+ Real cutoff, Real max_dist, int ref_index, int mdl_index)
+{
+ seq::ConstSequenceHandle ref_seq=aln.GetSequence(ref_index);
+ seq::ConstSequenceHandle mdl_seq=aln.GetSequence(mdl_index);
+ if (!ref_seq.HasAttachedView()) {
+ LOG_ERROR("reference sequence doesn't have a view attached.");
+ return 0.0;
+ }
+ if (!mdl_seq.HasAttachedView()) {
+ LOG_ERROR("model sequence doesn't have a view attached");
+ return 0.0;
+ }
+ XCSEditor edi=ref_seq.GetAttachedView().GetHandle().EditXCS(BUFFERED_EDIT);
+ for (int i=0; i<aln.GetLength(); ++i) {
+ ResidueView mdl_res=mdl_seq.GetResidue(i);
+ if (!mdl_res) {
+ continue;
+ }
+ String rname=mdl_res.GetName();
+ if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" ||
+ rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+ continue;
+ }
+ std::pair<Real, Real> ov1=calc_overlap2(ref_seq, mdl_seq, i,
+ cutoff, max_dist, true,
+ false);
+ std::pair<Real, Real> ov2=calc_overlap2(ref_seq, mdl_seq, i,
+ cutoff, max_dist, true,
+ true);
+ if (ov1.first/ov1.second<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()));
+ }
+ }
+ }
+ }
+ std::pair<Real, Real> total_ov(0.0, 0.0);
+ for (int i=0; i<aln.GetLength(); ++i) {
+ std::pair<Real, Real> ov1=calc_overlap2(ref_seq, mdl_seq, i, cutoff,
+ max_dist, false, false);
+ total_ov.first+=ov1.first;
+ total_ov.second+=ov1.second;
+ }
+ return total_ov.first/(total_ov.second ? total_ov.second : 1);
+ return 0.0;
}
}}}
diff --git a/modules/mol/alg/src/local_dist_test.hh b/modules/mol/alg/src/local_dist_test.hh
index f0ba3ebea20a6138a0cf14f2a34d24839abe3ef4..ab7588cdf84d401aaee052bcfe2b390da8174b15 100644
--- a/modules/mol/alg/src/local_dist_test.hh
+++ b/modules/mol/alg/src/local_dist_test.hh
@@ -3,12 +3,17 @@
#include <ost/mol/entity_view.hh>
#include <ost/mol/alg/module_config.hh>
+#include <ost/seq/alignment_handle.hh>
namespace ost { namespace mol { namespace alg {
Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
const EntityView& ref,
Real cutoff, Real max_dist);
+
+Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
+ Real cutoff, Real max_dist,
+ int ref_index=0, int mdl_index=1);
}}}
#endif