From 9f40b1e587140ccca0e683ed7a103759ccb7e9b6 Mon Sep 17 00:00:00 2001
From: Marco Biasini <marco.biasini@unibas.ch>
Date: Sat, 30 Oct 2010 14:46:36 +0200
Subject: [PATCH] added Builder::GuessChemClass

So far, the method only guesses whether a residue is a peptide based
on the availability of the backbone atoms N, C, CA, O with correct
connectivity.
---
 modules/conop/src/builder.cc        | 19 +++++++++++++++++++
 modules/conop/src/builder.hh        |  4 +++-
 modules/conop/src/conop.cc          | 16 +++++++++++++++-
 modules/conop/tests/test_builder.cc | 19 +++++++++++++++++++
 4 files changed, 56 insertions(+), 2 deletions(-)

diff --git a/modules/conop/src/builder.cc b/modules/conop/src/builder.cc
index 60bf273a5..b50048c5e 100644
--- a/modules/conop/src/builder.cc
+++ b/modules/conop/src/builder.cc
@@ -163,6 +163,25 @@ bool Builder::AreResiduesConsecutive(const mol::ResidueHandle& r1,
          r2.GetNumber().GetInsCode()==r1.GetNumber().NextInsertionCode().GetInsCode();
 }
 
+void Builder::GuessChemClass(mol::ResidueHandle res)
+{
+  // try peptide
+  res.SetChemClass(mol::ChemClass());
+  mol::AtomHandle ca=res.FindAtom("CA");
+  if (!ca.IsValid() || ca.GetElement()!="C") return;
+  mol::AtomHandle n=res.FindAtom("N");
+  if (!n.IsValid() || n.GetElement()!="N") return;
+  mol::AtomHandle c=res.FindAtom("C");
+  if (!c.IsValid() || c.GetElement()!="C") return;
+  mol::AtomHandle o=res.FindAtom("O");
+  if (!o.IsValid() || o.GetElement()!="O") return;
+  if (this->IsBondFeasible(n, ca) && this->IsBondFeasible(ca, c) &&
+      this->IsBondFeasible(c, o)) {
+    res.SetChemClass(mol::ChemClass(mol::ChemClass::PeptideLinking));
+  }
+}
+
+
 void Builder::AssignBackBoneTorsionsToResidue(mol::ResidueHandle res)
 {
 
diff --git a/modules/conop/src/builder.hh b/modules/conop/src/builder.hh
index 5312d80e9..3946780d1 100644
--- a/modules/conop/src/builder.hh
+++ b/modules/conop/src/builder.hh
@@ -137,7 +137,9 @@ public:
   /// \brief whether the r1 and r2 have consecutive residue numbers
   static bool AreResiduesConsecutive(const mol::ResidueHandle& r1, 
                                      const mol::ResidueHandle& r2);
-
+  /// \brief guess and assign chemical class of residue based on atoms and
+  ///     connectivity
+  void GuessChemClass(mol::ResidueHandle res);
   /// |brief Connect \p atom with all atoms for whith IsBondFeasible() and 
   ///    AreResiduesConsecutive() returns true
   void DistanceBasedConnect(mol::AtomHandle atom);
diff --git a/modules/conop/src/conop.cc b/modules/conop/src/conop.cc
index f546f5204..76b45d7bf 100644
--- a/modules/conop/src/conop.cc
+++ b/modules/conop/src/conop.cc
@@ -213,7 +213,19 @@ private:
   std::map<String, int>  unk_res_;  
 };
 
-
+class ChemClassAssigner : public mol::EntityVisitor {
+public:
+  ChemClassAssigner(BuilderP b): builder(b) { }
+  virtual bool VisitResidue(const mol::ResidueHandle& res)
+  {
+    if (res.GetChemClass()==mol::ChemClass()) {
+      builder->GuessChemClass(res);
+    }
+    return false;
+  }
+private:
+  BuilderP builder;
+};
 class Connector: public mol::EntityVisitor
 {
 public:
@@ -265,6 +277,8 @@ void Conopology::ConnectAll(const BuilderP& b, mol::EntityHandle eh, int flag)
   mol::XCSEditor xcs_e=eh.EditXCS(mol::BUFFERED_EDIT);
   PropAssigner a(b);
   eh.Apply(a);
+  ChemClassAssigner cca(b);
+  eh.Apply(cca);
   LOG_DEBUG("Conopology: ConnectAll: connecting all bonds");
   Connector connector(b);
   eh.Apply(connector);
diff --git a/modules/conop/tests/test_builder.cc b/modules/conop/tests/test_builder.cc
index 4ba860bb1..9a4ff6aa6 100644
--- a/modules/conop/tests/test_builder.cc
+++ b/modules/conop/tests/test_builder.cc
@@ -18,6 +18,7 @@
 //------------------------------------------------------------------------------
 #include <ost/mol/mol.hh>
 #include <ost/conop/builder.hh>
+#include <ost/conop/heuristic_builder.hh>
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 
@@ -29,6 +30,24 @@ using namespace ost::mol;
 
 BOOST_AUTO_TEST_SUITE( conop )
 
+BOOST_AUTO_TEST_CASE(test_guess_chem_class)
+{
+  HeuristicBuilder h;
+  EntityHandle ent=mol::CreateEntity();
+  XCSEditor edi=ent.EditXCS();
+  ChainHandle chain=edi.InsertChain("A");
+  ResidueHandle res=edi.AppendResidue(chain, "DUMMY");
+  AtomHandle n=edi.InsertAtom(res, "N", geom::Vec3(1, 0, 0), "N");
+  AtomHandle ca=edi.InsertAtom(res, "CA", geom::Vec3(2, 0, 0), "C");
+  AtomHandle c=edi.InsertAtom(res, "C", geom::Vec3(3, 0, 0), "C");
+  AtomHandle o=edi.InsertAtom(res, "O", geom::Vec3(4, 0, 0), "O");
+  h.GuessChemClass(res);
+  BOOST_CHECK(res.IsPeptideLinking());
+  res.SetChemClass(ChemClass());
+  edi.SetAtomPos(n, geom::Vec3(-1,0,0));
+  h.GuessChemClass(res);
+  BOOST_CHECK(!res.IsPeptideLinking());
+}
 
 BOOST_AUTO_TEST_CASE( test_builder )
 {
-- 
GitLab