added LocalDistTest

modules/index.rst

@@ -11,6 +11,7 @@ OpenStructure documentation
   base/base
   geom/geom
   mol/base/mol
+  mol/alg/molalg
   conop/conop
   img/base/img
   img/alg/alg

modules/mol/alg/molalg.rst

+:mod:`mol.alg <ost.mol.alg>` -- Algorithms for Structures
+================================================================================
+
+.. module:: ost.mol.alg
+   :synopsis: Algorithms operating on molecular structures
+
+.. function:: LocalDistTest(model, reference, tolerance, radius)
+  
+  This function calculates the agreement of local contacts between the model and 
+  the reference structure. The overlap is a number between zero and one, where 
+  one indicates complete agreement, zero indicates no agreement at all. This 
+  score is similar to the GDT, but does not require any superposition of the 
+  model and the reference.
+  
+  The distance of atom pairs in the reference structure that are closer than a 
+  certain distance (radius) to each other is compared to the distance found in 
+  the model. If the difference between these two distances is smaller than a 
+  threshold value (tolerance), the model and the reference agree. Missing atoms 
+  in the model are treated disagreement and thus lower the score.
+  
+  For residue with symmetric sidechains (GLU, ASP, ARG, VAL, PHE, TYR), the 
+  naming of the atoms is ambigous. For these residues, the overlap of both 
+  possible solutions to the fixed atoms, that is, everything that is not 
+  ambigous is calculated. The solution that gives higher overlap is then used to 
+  calculate the actual overlap score.
\ No newline at end of file

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -22,6 +22,8 @@
  */
 #include <boost/python.hpp>
 #include <ost/config.hh>
+#include <ost/mol/alg/local_dist_test.hh>
+
 using namespace boost::python;
 
 void export_svdSuperPose();
@@ -36,4 +38,6 @@ BOOST_PYTHON_MODULE(_mol_alg)
   #if OST_IMG_ENABLED
   export_entity_to_density();
   #endif
+  
+  def("LocalDistTest", &ost::mol::alg::LocalDistTest);
 }

modules/mol/alg/src/CMakeLists.txt

@@ -2,11 +2,13 @@ set(OST_MOL_ALG_HEADERS
   svd_superpose.hh
   module_config.hh
   sec_structure_segments.hh
+  local_dist_test.hh
 )
 
 set(OST_MOL_ALG_SOURCES
   svd_superpose.cc
   sec_structure_segments.cc
+  local_dist_test.cc
 )
 
 set(MOL_ALG_DEPS mol)

modules/mol/alg/src/local_dist_test.cc

+#include <ost/mol/mol.hh>
+#include "local_dist_test.hh"
+
+namespace ost { namespace mol { namespace alg {
+
+namespace {
+
+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;
+}
+
+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=="ARG") {
+    return (aname=="NH1" || aname=="NH2");
+  }  
+  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> overlap(0.0, 0.0);
+  AtomViewList ref_atoms=ref_res.GetAtomList();
+
+  ResidueView mdl_res=mdl_chain.FindResidue(ref_res.GetNumber());
+  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();
+    AtomViewList 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;
+        }
+        AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(), 
+                                             aj->GetName());
+        overlap.second+=1.0;
+        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;
+        }
+        continue;
+      }
+      if (aj->GetResidue().GetNumber()>ref_res.GetNumber()) {
+        if (aj->GetResidue().GetNumber()<=ref_res.GetNumber()) {
+          continue;
+        }
+        AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(), 
+                                             aj->GetName());
+        overlap.second+=1.0;
+        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;
+        }
+      }
+    }      
+  }
+  return overlap;
+}
+
+}
+
+Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
+                   Real cutoff, Real max_dist)
+{
+  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 
+  // fixed atoms and take the solution that gives bigger scores.
+  XCSEditor edi=ref.GetHandle().RequestXCSEditor(BUFFERED_EDIT);
+  for (ResidueViewList::iterator
+       i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
+    const String rname=i->GetName();
+    ResidueView mdl_res=mdl_chain.FindResidue(i->GetNumber());
+    if (!mdl_res) {
+     continue;
+    }     
+    if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || 
+         rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+      continue;
+    }
+
+    std::pair<Real, Real> ov1=calc_overlap(*i, ref, mdl_chain, 
+                                          cutoff, max_dist, true, 
+                                          false);
+    std::pair<Real, Real> ov2=calc_overlap(*i, ref, mdl_chain, 
+                                          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 (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, 
+                                            max_dist, false, false);
+     total_ov.first+=ov1.first;
+     total_ov.second+=ov1.second;
+  }
+  return total_ov.first/total_ov.second;
+}
+
+}}}

modules/mol/alg/src/local_dist_test.hh

+#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>
+
+namespace ost { namespace mol { namespace alg {
+  
+Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
+                                         const EntityView& ref,
+                                         Real cutoff, Real max_dist);
+}}}
+
+#endif
+