fix quack-typing for residues with element-less atoms

modules/conop/src/heuristic.cc

@@ -32,14 +32,16 @@ namespace ost { namespace conop {
 
 void HeuristicProcessor::ProcessUnkResidue(DiagnosticsPtr diags, 
                                            mol::ResidueHandle res) const {
-  res.SetChemClass(GuessChemClass(res));
   mol::AtomHandleList atoms = res.GetAtomList(); 
   for (mol::AtomHandleList::iterator 
        i = atoms.begin(), e = atoms.end(); i != e; ++i) {
     mol::AtomHandle a = *i;
     if (!Conopology::Instance().IsValidElement(a.GetElement()))
-      a.SetElement(GuessAtomElement(a.GetName(), a.IsHetAtom()));
+      a.SetElement(GuessAtomElement(a.GetName(), atoms.size()));
   }
+
+  res.SetChemClass(GuessChemClass(res));
+
   if (!this->GetConnect()) { return; }
   for (mol::AtomHandleList::iterator 
        i = atoms.begin(), e = atoms.end(); i != e; ++i) {
@@ -64,6 +66,9 @@ void HeuristicProcessor::DoProcess(DiagnosticsPtr diags,
       if (!compound) {
         // process unknown residue...
         this->ProcessUnkResidue(diags, residue);
+        if (this->GetConnectAminoAcids()) 
+          this->ConnectResidues(prev, residue);
+        prev = residue;
         continue;
       }
       this->ReorderAtoms(residue, compound, true);  

modules/conop/src/processor.cc

@@ -115,6 +115,7 @@ void AssignBackboneTorsions(mol::ResidueHandle prev,
     }
   }
 } 
+
 mol::ChemClass GuessChemClass(mol::ResidueHandle res)
 {
   mol::AtomHandle ca=res.FindAtom("CA");
@@ -131,7 +132,8 @@ mol::ChemClass GuessChemClass(mol::ResidueHandle res)
   }
   return mol::ChemClass();
 }
-String GuessAtomElement(const String& aname, bool hetatm)
+
+String GuessAtomElement(const String& aname, int atom_count)
 {
   static String l1[] = {
     "H","C","N","O","P","S","K"
@@ -156,7 +158,7 @@ String GuessAtomElement(const String& aname, bool hetatm)
   if(ele[0]=='H') {
     return "H";
   }
-  if (hetatm==false) {
+  if (atom_count > 1) {
     if (ele=="CA" || ele=="CB") {
       return "C";
     }

modules/conop/src/processor.hh

@@ -126,7 +126,7 @@ String DLLEXPORT_OST_CONOP StringFromConopAction(ConopAction action);
 
 
 /// \brief guess element of atom based on name and hetatm flag
-String DLLEXPORT_OST_CONOP GuessAtomElement(const String& atom_name, bool hetatm);
+String DLLEXPORT_OST_CONOP GuessAtomElement(const String& atom_name, int atom_count);
 
 /// \brief guess chemclass based on atoms of residue
 mol::ChemClass DLLEXPORT_OST_CONOP GuessChemClass(mol::ResidueHandle res);

modules/conop/src/rule_based.cc

@@ -63,7 +63,7 @@ void RuleBasedProcessor::DoProcess(DiagnosticsPtr diags,
         mol::AtomHandleList unk_atoms;
         unk_atoms = GetUnknownAtomsOfResidue(residue, compound, 
                                              this->GetStrictHydrogens());
-        this->ProcessUnkAtoms(diags, unk_atoms, atoms_to_connect);
+        this->ProcessUnkAtoms(diags, residue, unk_atoms, atoms_to_connect);
         residue.SetChemClass(mol::ChemClass(mol::ChemClass::UNKNOWN));
         residue.SetChemType(mol::ChemType(mol::ChemType::UNKNOWN));
         residue.SetOneLetterCode('?');
@@ -102,7 +102,7 @@ void RuleBasedProcessor::ProcessUnkResidue(DiagnosticsPtr diags,
          e = atoms.end(); i !=e; ++i) {
       remaining_atoms.push_back(*i);
       if (!Conopology::Instance().IsValidElement(i->GetElement()))
-        i->SetElement(GuessAtomElement(i->GetName(), i->IsHetAtom()));
+        i->SetElement(GuessAtomElement(i->GetName(), atoms.size()));
     }
   }
   // Don't do anything if treatment is set to SILENT
@@ -129,6 +129,7 @@ void RuleBasedProcessor::ProcessUnkResidue(DiagnosticsPtr diags,
 }
 
 void RuleBasedProcessor::ProcessUnkAtoms(DiagnosticsPtr diags,
+                                         mol::ResidueHandle res,
                                          mol::AtomHandleList unks,
                                          mol::AtomHandleList& remaining_atoms) const
 {
@@ -138,7 +139,7 @@ void RuleBasedProcessor::ProcessUnkAtoms(DiagnosticsPtr diags,
          e = unks.end(); i !=e; ++i) {
       remaining_atoms.push_back(*i);
       if (!Conopology::Instance().IsValidElement(i->GetElement()))
-        i->SetElement(GuessAtomElement(i->GetName(), i->IsHetAtom()));
+        i->SetElement(GuessAtomElement(i->GetName(), res.GetAtomCount()));
     }
   }
 

modules/conop/src/rule_based.hh

@@ -85,6 +85,7 @@ protected:
                          mol::ResidueHandle res, 
                          mol::AtomHandleList& remaining) const;
   void ProcessUnkAtoms(DiagnosticsPtr diags,
+                       mol::ResidueHandle res,
                        mol::AtomHandleList unks,
                        mol::AtomHandleList& remaining) const;
   virtual void DoProcess(DiagnosticsPtr diags, 

modules/conop/tests/test_heuristic_conop.cc

@@ -19,6 +19,7 @@
 
 #include <ost/log.hh>
 #include <ost/mol/mol.hh>
+#include <ost/mol/builder.hh>
 #include <ost/conop/heuristic.hh>
 
 #define BOOST_TEST_DYN_LINK
@@ -140,5 +141,44 @@ BOOST_AUTO_TEST_CASE(does_assign_torsions) {
   BOOST_CHECK(a2.GetPsiTorsion().IsValid());
 }
 
+
+BOOST_AUTO_TEST_CASE(quack_types_unknown_residues) {
+  
+  mol::EntityHandle e = Builder()
+    .Chain("A")
+      .Residue("GLY")
+        .Atom("N", geom::Vec3(-8.22, 35.20, 22.39))
+        .Atom("CA", geom::Vec3(-8.28, 36.36, 21.49))
+        .Atom("C", geom::Vec3(-8.59, 35.93, 20.06))
+        .Atom("O", geom::Vec3(-7.88, 36.30, 19.12))
+        .Atom("CB", geom::Vec3(-6.96, 37.11, 21.53))
+      .Residue("SEP")
+        .Atom("N", geom::Vec3(-9.66, 35.16, 19.91))
+        .Atom("CA", geom::Vec3(-10.08, 34.66, 18.61))
+        .Atom("CB", geom::Vec3(-10.90, 33.38, 18.79))
+        .Atom("OG", geom::Vec3(-12.22, 33.66, 19.22))
+        .Atom("C", geom::Vec3(-10.84, 35.69, 17.77))
+        .Atom("O", geom::Vec3(-10.99, 35.54, 16.56))
+        .Atom("P", geom::Vec3(-12.43, 33.20, 20.70))
+        .Atom("O1P", geom::Vec3(-11.55, 34.04, 21.65))
+        .Atom("O2P", geom::Vec3(-12.05, 31.71, 20.86))
+        .Atom("O3P", geom::Vec3(-13.92, 33.41, 21.07))
+      .Residue("GLY")
+        .Atom("N", geom::Vec3(-11.31, 36.75, 18.43))
+        .Atom("CA", geom::Vec3(-12.04, 37.81, 17.75))
+        .Atom("C", geom::Vec3(-12.34, 38.94, 18.73))
+        .Atom("O", geom::Vec3(-13.45, 39.03, 19.26))
+  ;
+  HeuristicProcessor proc;
+  proc.Process(e);
+  BOOST_CHECK(e.FindResidue("A", 1).IsPeptideLinking());
+  BOOST_CHECK(e.FindResidue("A", 2).IsPeptideLinking());
+  BOOST_CHECK(e.FindResidue("A", 3).IsPeptideLinking());
+
+  BOOST_CHECK(mol::BondExists(e.FindAtom("A", 1, "C"), e.FindAtom("A", 2, "N")));
+  BOOST_CHECK(mol::BondExists(e.FindAtom("A", 2, "C"), e.FindAtom("A", 3, "N")));
+
+}
+
 BOOST_AUTO_TEST_SUITE_END();