From f18f0cb0f1e6c8d32d1e2d782c5192b9a6e10b60 Mon Sep 17 00:00:00 2001
From: Marco Biasini <mvbiasini@gmail.com>
Date: Sat, 13 Jul 2013 18:27:46 +0200
Subject: [PATCH] fix quack-typing for residues with element-less atoms
---
modules/conop/src/heuristic.cc | 9 +++--
modules/conop/src/processor.cc | 6 ++--
modules/conop/src/processor.hh | 2 +-
modules/conop/src/rule_based.cc | 7 ++--
modules/conop/src/rule_based.hh | 1 +
modules/conop/tests/test_heuristic_conop.cc | 40 +++++++++++++++++++++
6 files changed, 57 insertions(+), 8 deletions(-)
diff --git a/modules/conop/src/heuristic.cc b/modules/conop/src/heuristic.cc
index cecf6dcb2..8c5277206 100644
--- a/modules/conop/src/heuristic.cc
+++ b/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);
diff --git a/modules/conop/src/processor.cc b/modules/conop/src/processor.cc
index d8f5c2698..4ce93eba6 100644
--- a/modules/conop/src/processor.cc
+++ b/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";
}
diff --git a/modules/conop/src/processor.hh b/modules/conop/src/processor.hh
index 6001b5cbf..d66f6ccbc 100644
--- a/modules/conop/src/processor.hh
+++ b/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);
diff --git a/modules/conop/src/rule_based.cc b/modules/conop/src/rule_based.cc
index 2a18ddb20..c4c870005 100644
--- a/modules/conop/src/rule_based.cc
+++ b/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()));
}
}
diff --git a/modules/conop/src/rule_based.hh b/modules/conop/src/rule_based.hh
index 4275d7e40..3b91771a2 100644
--- a/modules/conop/src/rule_based.hh
+++ b/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,
diff --git a/modules/conop/tests/test_heuristic_conop.cc b/modules/conop/tests/test_heuristic_conop.cc
index 494339316..61155c644 100644
--- a/modules/conop/tests/test_heuristic_conop.cc
+++ b/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();
--
GitLab