From 78606b846c4bd8e694ef97b72b8fb7f0e69e7e15 Mon Sep 17 00:00:00 2001
From: marco <marco@5a81b35b-ba03-0410-adc8-b2c5c5119f08>
Date: Tue, 21 Sep 2010 12:34:43 +0000
Subject: [PATCH] improved behaviour of rule-based builder
as soon as an atom is unknown, we switch back to a distance-based
connect mechanism for the full residue. This avoids problems with
residues coming from CHARMM/GROMACS and other force fields that do
not follow the naming convention of the PDB. To really become useful,
some additional changes to the compound library are required that
will be implemented in a next commit
git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2714 5a81b35b-ba03-0410-adc8-b2c5c5119f08
---
modules/conop/src/chemdict_parser.cc | 153 ++++++++++++++++++++++++
modules/conop/src/chemdict_parser.hh | 89 ++++++++++++++
modules/conop/src/conop.cc | 24 +++-
modules/conop/src/rule_based_builder.cc | 90 ++++++++++----
modules/conop/src/rule_based_builder.hh | 6 +-
5 files changed, 333 insertions(+), 29 deletions(-)
create mode 100644 modules/conop/src/chemdict_parser.cc
create mode 100644 modules/conop/src/chemdict_parser.hh
diff --git a/modules/conop/src/chemdict_parser.cc b/modules/conop/src/chemdict_parser.cc
new file mode 100644
index 000000000..2c7df33c0
--- /dev/null
+++ b/modules/conop/src/chemdict_parser.cc
@@ -0,0 +1,153 @@
+#include <ost/conop/chemdict_parser.hh>
+
+namespace ost { namespace conop {
+
+bool ChemdictParser::OnBeginData(const StringRef& data_name)
+{
+ compound_.reset(new Compound(data_name.str()));
+ if (last_!=data_name[0]) {
+ last_=data_name[0];
+ std::cout << last_ << std::flush;
+ }
+ atom_map_.clear();
+ insert_=true;
+ return true;
+}
+
+bool ChemdictParser::OnBeginLoop(const io::StarLoopDesc& header)
+{
+ if (header.GetCategory()=="chem_comp_atom") {
+ loop_type_=ATOM_SPEC;
+ indices_[ATOM_NAME]=header.GetIndex("atom_id");
+ indices_[ALT_ATOM_NAME]=header.GetIndex("alt_atom_id");
+ indices_[ELE]=header.GetIndex("type_symbol");
+ indices_[IS_LEAVING]=header.GetIndex("pdbx_leaving_atom_flag");
+ indices_[IS_AROMATIC]=header.GetIndex("pdbx_aromatic_flag");
+ indices_[ORDINAL]=header.GetIndex("pdbx_ordinal");
+ return true;
+ } else if (header.GetCategory()=="chem_comp_bond") {
+ loop_type_=BOND_SPEC;
+ indices_[ATOM_ID1]=header.GetIndex("atom_id_1");
+ indices_[ATOM_ID2]=header.GetIndex("atom_id_2");
+ indices_[BOND_ORDER]=header.GetIndex("value_order");
+ return true;
+ }
+ loop_type_=DONT_KNOW;
+ return false;
+}
+
+void ChemdictParser::OnDataRow(const io::StarLoopDesc& header,
+ const std::vector<StringRef>& columns)
+{
+ if (loop_type_==ATOM_SPEC) {
+ AtomSpec atom;
+ atom.is_leaving=columns[indices_[IS_LEAVING]][0]=='Y';
+ atom.name=columns[indices_[ATOM_NAME]].str();
+ atom.alt_name=columns[indices_[ALT_ATOM_NAME]].str();
+ atom.ordinal=columns[indices_[ORDINAL]].to_int().second-1;
+ atom.element=columns[indices_[ELE]].str();
+ atom.is_aromatic=columns[indices_[IS_AROMATIC]][0]=='Y';
+ compound_->AddAtom(atom);
+ atom_map_[atom.name]=atom.ordinal;
+ } else if (loop_type_==BOND_SPEC) {
+ BondSpec bond;
+ String atom_one=columns[indices_[ATOM_ID1]].str();
+ String atom_two=columns[indices_[ATOM_ID2]].str();
+ bond.atom_one=atom_map_[atom_one];
+ bond.atom_two=atom_map_[atom_two];
+ char b=columns[indices_[BOND_ORDER]][0];
+ bond.order=b=='D' ? 2 : (b=='S' ? 1 : 3);
+ if (bond.atom_one>bond.atom_two) {
+ std::swap(bond.atom_one, bond.atom_two);
+ }
+ compound_->AddBond(bond);
+ }
+}
+
+void ChemdictParser::OnDataItem(const io::StarDataItem& item)
+{
+ if (item.GetCategory()==StringRef("chem_comp", 9)) {
+ if (item.GetName()==StringRef("type", 4)) {
+ // convert type to uppercase
+ String type=item.GetValue().str();
+ for (String::iterator i=type.begin(), e=type.end(); i!=e; ++i) {
+ *i=toupper(*i);
+ }
+ // The type of water is set to "?". let's change it to water...
+ if (compound_->GetID()=="HOH") {
+ compound_->SetChemClass(mol::ChemClass(mol::ChemClass::Water));
+ } else {
+ std::map<String, mol::ChemClass>::iterator i=tm_.find(type);
+ if (i!=tm_.end()) {
+ compound_->SetChemClass(i->second);
+ } else {
+ std::cout << "unknown type '" << type << "' for compound "
+ << compound_->GetID() << std::endl;
+ }
+ }
+
+ } else if (item.GetName()==StringRef("formula", 7)) {
+ compound_->SetFormula(item.GetValue().str());
+ } else if (item.GetName()==StringRef("one_letter_code", 15)) {
+ if (item.GetValue().length()==1) {
+ compound_->SetOneLetterCode(item.GetValue()[0]);
+ }
+ } else if (item.GetName()==StringRef("pdbx_initial_date", 17)) {
+ compound_->SetCreationDate(Date::FromString(item.GetValue()));
+ } else if (item.GetName()==StringRef("pdbx_modified_date", 18)) {
+ compound_->SetModificationDate(Date::FromString(item.GetValue()));
+ }
+ } else if (item.GetName()==StringRef("atom_id", 7)) {
+ atom_.name=item.GetValue().str();
+ } else if (item.GetName()==StringRef("alt_atom_id", 11)) {
+ atom_.alt_name=item.GetValue().str();
+ } else if (item.GetName()==StringRef("type_symbol", 11)) {
+ atom_.element=item.GetValue().str();
+ } else if (item.GetName()==StringRef("pdbx_ordinal", 12)) {
+ atom_.ordinal=item.GetValue().to_int().second-1;
+ }
+}
+
+void ChemdictParser::OnEndData()
+{
+ if (insert_) {
+ if (compound_->GetAtomSpecs().empty()) {
+ compound_->AddAtom(atom_);
+ }
+ lib_->AddCompound(compound_);
+ }
+}
+
+std::map<String, mol::ChemClass> ChemdictParser::tm_=std::map<String, mol::ChemClass>();
+
+void ChemdictParser::InitTypeMap()
+{
+ if (!tm_.empty())
+ return;
+ tm_["L-PEPTIDE COOH CARBOXY TERMINUS"]=mol::ChemClass(mol::ChemClass::LPeptideLinking);
+ tm_["L-PEPTIDE NH3 AMINO TERMINUS"]=mol::ChemClass(mol::ChemClass::LPeptideLinking);
+ tm_["D-PEPTIDE NH3 AMINO TERMINUS"]=mol::ChemClass(mol::ChemClass::DPeptideLinking);
+ tm_["L-SACCHARIDE 1,4 AND 1,4 LINKING"]=mol::ChemClass(mol::ChemClass::LSaccharide);
+ tm_["D-SACCHARIDE 1,4 AND 1,4 LINKING"]=mol::ChemClass(mol::ChemClass::DSaccharide);
+ tm_["L-SACCHARIDE"]=mol::ChemClass(mol::ChemClass::LSaccharide);
+ tm_["D-SACCHARIDE"]=mol::ChemClass(mol::ChemClass::DSaccharide);
+ tm_["SACCHARIDE"]=mol::ChemClass(mol::ChemClass::Saccharide);
+ tm_["D-PEPTIDE LINKING"]=mol::ChemClass(mol::ChemClass::DPeptideLinking);
+ tm_["L-PEPTIDE LINKING"]=mol::ChemClass(mol::ChemClass::LPeptideLinking);
+ tm_["L-PEPTIDE-LINKING"]=mol::ChemClass(mol::ChemClass::LPeptideLinking);
+ tm_["DNA LINKING"]=mol::ChemClass(mol::ChemClass::DNALinking);
+ tm_["RNA LINKING"]=mol::ChemClass(mol::ChemClass::RNALinking);
+ tm_["L-DNA LINKING"]=mol::ChemClass(mol::ChemClass::DNALinking);
+ tm_["L-RNA LINKING"]=mol::ChemClass(mol::ChemClass::RNALinking);
+ tm_["R-DNA LINKING"]=mol::ChemClass(mol::ChemClass::DNALinking);
+ tm_["R-RNA LINKING"]=mol::ChemClass(mol::ChemClass::RNALinking);
+ tm_["DNA OH 3 PRIME TERMINUS"]=mol::ChemClass(mol::ChemClass::DNALinking);
+ tm_["PEPTIDE-LIKE"]=mol::ChemClass(mol::ChemClass::PeptideLinking);
+ tm_["PEPTIDE LINKING"]=mol::ChemClass(mol::ChemClass::PeptideLinking);
+ tm_["PEPTIDE-LINKING"]=mol::ChemClass(mol::ChemClass::PeptideLinking);
+ tm_["NON-POLYMER"]=mol::ChemClass(mol::ChemClass::NonPolymer);
+ tm_["RNA OH 3 PRIME TERMINUS"]=mol::ChemClass(mol::ChemClass::RNALinking);
+ tm_["?"]=mol::ChemClass(mol::ChemClass::Unknown);
+}
+
+}}
\ No newline at end of file
diff --git a/modules/conop/src/chemdict_parser.hh b/modules/conop/src/chemdict_parser.hh
new file mode 100644
index 000000000..4753a389a
--- /dev/null
+++ b/modules/conop/src/chemdict_parser.hh
@@ -0,0 +1,89 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+#ifndef OST_CONOP_CHEMDICT_PARSER_HH
+#define OST_CONOP_CHEMDICT_PARSER_HH
+
+/*
+ Author: Marco Biasini
+ */
+
+
+#include <ost/mol/mol.hh>
+#include <ost/io/mol/star_parser.hh>
+#include <ost/conop/compound_lib.hh>
+
+namespace ost { namespace conop {
+
+typedef enum {
+ ATOM_SPEC,
+ BOND_SPEC,
+ DONT_KNOW
+} LoopType;
+
+class DLLEXPORT_OST_CONOP ChemdictParser : public io::StarParser {
+public:
+ ChemdictParser(std::istream& stream):
+ io::StarParser(stream),
+ compound_(new Compound("UNK")), last_(0), loop_type_(DONT_KNOW)
+ {
+ this->InitTypeMap();
+ }
+
+ virtual bool OnBeginData(const StringRef& data_name);
+
+ virtual bool OnBeginLoop(const io::StarLoopDesc& header);
+
+ virtual void OnDataRow(const io::StarLoopDesc& header,
+ const std::vector<StringRef>& columns);
+
+ virtual void OnDataItem(const io::StarDataItem& item);
+
+ virtual void OnEndData();
+
+ void SetCompoundLib(const CompoundLibPtr& lib)
+ {
+ lib_=lib;
+ }
+private:
+ void InitTypeMap();
+ CompoundLibPtr lib_;
+ CompoundPtr compound_;
+ typedef enum {
+ ATOM_NAME=0,
+ ALT_ATOM_NAME=1,
+ IS_AROMATIC=2,
+ ORDINAL=3,
+ IS_LEAVING=4,
+ ELE=5,
+ STEREO_CONF=6,
+ ATOM_ID1=0,
+ ATOM_ID2=1,
+ BOND_ORDER=2
+ } PropIndex;
+ char last_;
+ int indices_[10];
+ bool insert_;
+ static std::map<String, mol::ChemClass> tm_;
+ std::map<String, int> atom_map_;
+ LoopType loop_type_;
+ AtomSpec atom_;
+};
+
+
+}}
diff --git a/modules/conop/src/conop.cc b/modules/conop/src/conop.cc
index 7d8b35867..9d34eba51 100644
--- a/modules/conop/src/conop.cc
+++ b/modules/conop/src/conop.cc
@@ -214,10 +214,13 @@ namespace {
class PropAssigner: public mol::EntityVisitor {
public:
PropAssigner(const BuilderP& builder): builder_(builder) {}
+
virtual bool VisitResidue(const mol::ResidueHandle& res)
{
String key=builder_->IdentifyResidue(res);
- builder_->CheckResidueCompleteness(res);
+ if (key=="UNK") {
+ unk_res_[res.GetKey()]+=1;
+ }
builder_->FillResidueProps(res);
return true;
}
@@ -227,8 +230,23 @@ public:
builder_->FillAtomProps(atom);
return false;
}
+
+ virtual void OnExit()
+ {
+ for (std::map<String, int>::iterator i=unk_res_.begin(),
+ e=unk_res_.end(); i!=e; ++i) {
+ if (i->second>1) {
+ LOG_WARNING("structure contains unknown residues with name " << i->first
+ << " ("<< i->second << "x)");
+ } else {
+ LOG_WARNING("structure contains unknown residue with name " << i->first);
+ }
+
+ }
+ }
private:
BuilderP builder_;
+ std::map<String, int> unk_res_;
};
@@ -247,8 +265,8 @@ public:
virtual bool VisitResidue(const mol::ResidueHandle& res) {
builder_->ConnectAtomsOfResidue(res);
- if(prev_) {
- builder_->ConnectResidueToPrev(res,prev_);
+ if (prev_) {
+ builder_->ConnectResidueToPrev(res,prev_);
}
prev_=res;
return false;
diff --git a/modules/conop/src/rule_based_builder.cc b/modules/conop/src/rule_based_builder.cc
index 3353dfb46..d79dbf75b 100644
--- a/modules/conop/src/rule_based_builder.cc
+++ b/modules/conop/src/rule_based_builder.cc
@@ -33,8 +33,9 @@ namespace ost { namespace conop {
void RuleBasedBuilder::CompleteAtoms(mol::ResidueHandle rh) {}
-void RuleBasedBuilder::CheckResidueCompleteness(const mol::ResidueHandle& rh) {
- LookupCompound(rh);
+void RuleBasedBuilder::CheckResidueCompleteness(const mol::ResidueHandle& rh)
+{
+ this->LookupCompound(rh);
if (!last_compound_) {
mol::AtomHandleList atoms=rh.GetAtomList();
for (mol::AtomHandleList::const_iterator i=atoms.begin(),
@@ -43,7 +44,7 @@ void RuleBasedBuilder::CheckResidueCompleteness(const mol::ResidueHandle& rh) {
}
return;
}
- ReorderAtoms(rh , last_compound_);
+ this->ReorderAtoms(rh , last_compound_);
AtomSpecList::const_iterator j=last_compound_->GetAtomSpecs().begin();
mol::AtomHandleList atoms=rh.GetAtomList();
mol::AtomHandleList::iterator i=atoms.begin();
@@ -53,12 +54,30 @@ void RuleBasedBuilder::CheckResidueCompleteness(const mol::ResidueHandle& rh) {
if ((*j).ordinal!=static_cast<int>((*i).Impl()->GetState())) {
this->OnMissingAtom(rh, (*j).name);
} else {
- this->FillAtomProps(*i, *j);
++i;
}
}
}
+bool RuleBasedBuilder::HasUnknownAtoms(mol::ResidueHandle res)
+{
+ this->LookupCompound(res);
+ if (!last_compound_) {
+ return true;
+ }
+ this->ReorderAtoms(res, last_compound_);
+ AtomSpecList::const_iterator j=last_compound_->GetAtomSpecs().begin();
+ mol::AtomHandleList atoms=res.GetAtomList();
+ mol::AtomHandleList::iterator i=atoms.begin();
+ for (mol::AtomHandleList::iterator
+ i=atoms.begin(), e=atoms.end(); i!=e; ++i) {
+ if ((*i).Impl()->GetState()==std::numeric_limits<unsigned int>::max()) {
+ return true;
+ }
+ }
+ return false;
+}
+
void RuleBasedBuilder::FillAtomProps(mol::AtomHandle atom, const AtomSpec& spec)
{
mol::AtomProp props=atom.GetAtomProps();
@@ -75,8 +94,9 @@ void RuleBasedBuilder::FillAtomProps(mol::AtomHandle atom, const AtomSpec& spec)
atom.SetAtomProps(props);
}
-void RuleBasedBuilder::FillResidueProps(mol::ResidueHandle residue) {
- LookupCompound(residue);
+void RuleBasedBuilder::FillResidueProps(mol::ResidueHandle residue)
+{
+ this->LookupCompound(residue);
if (!last_compound_)
return;
residue.SetChemClass(last_compound_->GetChemClass());
@@ -103,6 +123,9 @@ void RuleBasedBuilder::LookupCompound(const mol::ResidueHandle& rh)
void RuleBasedBuilder::ReorderAtoms(mol::ResidueHandle residue,
CompoundPtr compound)
{
+ if (last_residue_==residue) {
+ return;
+ }
mol::impl::ResidueImplPtr impl=residue.Impl();
mol::impl::AtomImplList::iterator i=impl->GetAtomList().begin();
for (; i!=impl->GetAtomList().end(); ++i) {
@@ -110,7 +133,7 @@ void RuleBasedBuilder::ReorderAtoms(mol::ResidueHandle residue,
int index=compound->GetAtomSpecIndex(atom->GetName());
if (index==-1) {
if (!this->OnUnknownAtom(mol::AtomHandle(atom))) {
- atom->SetState(std::numeric_limits<int>::max());
+ atom->SetState(std::numeric_limits<unsigned int>::max());
}
continue;
}
@@ -118,6 +141,14 @@ void RuleBasedBuilder::ReorderAtoms(mol::ResidueHandle residue,
}
std::sort(impl->GetAtomList().begin(), impl->GetAtomList().end(),
OrdinalComp());
+ last_residue_=residue;
+ unknown_atoms_=this->HasUnknownAtoms(residue);
+ if (unknown_atoms_) {
+ LOG_WARNING("residue " << residue << " doesn't look like a standard "
+ << residue.GetKey());
+ residue.SetChemClass(mol::ChemClass(mol::ChemClass::Unknown));
+ residue.SetOneLetterCode('?');
+ }
}
@@ -132,7 +163,9 @@ mol::ResidueKey RuleBasedBuilder::IdentifyResidue(const mol::ResidueHandle& rh)
}
}
-mol::AtomHandle RuleBasedBuilder::LocateAtom(const mol::AtomHandleList& ahl, int ordinal) {
+mol::AtomHandle RuleBasedBuilder::LocateAtom(const mol::AtomHandleList& ahl,
+ int ordinal)
+{
if (ahl.empty())
return mol::AtomHandle();
const mol::AtomHandle* r_it=&ahl.back();
@@ -148,14 +181,26 @@ mol::AtomHandle RuleBasedBuilder::LocateAtom(const mol::AtomHandleList& ahl, int
return not_found ? mol::AtomHandle() : *r_it;
}
-void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh) {
+
+inline void dist_connect(RuleBasedBuilder* b, mol::AtomHandleList atoms)
+{
+ for (mol::AtomHandleList::const_iterator i=atoms.begin(),
+ e=atoms.end(); i!=e; ++i) {
+ b->DistanceBasedConnect(*i);
+ }
+}
+void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh)
+{
+
LookupCompound(rh);
+
if (!last_compound_) {
- mol::AtomHandleList atoms=rh.GetAtomList();
- for (mol::AtomHandleList::const_iterator i=atoms.begin(),
- e=atoms.end(); i!=e; ++i) {
- this->DistanceBasedConnect(*i);
- }
+ dist_connect(this, rh.GetAtomList());
+ return;
+ }
+ this->ReorderAtoms(rh, last_compound_);
+ if (unknown_atoms_) {
+ dist_connect(this, rh.GetAtomList());
return;
}
mol::XCSEditor e=rh.GetEntity().RequestXCSEditor(mol::BUFFERED_EDIT);
@@ -166,19 +211,14 @@ void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh) {
mol::AtomHandle a1=this->LocateAtom(atoms, bond.atom_one);
mol::AtomHandle a2=this->LocateAtom(atoms, bond.atom_two);
if (a1.IsValid() && a2.IsValid() && this->IsBondFeasible(a1, a2)) {
- e.Connect(a1, a2, bond.order);
+ e.Connect(a1, a2, bond.order);
}
}
- for (mol::AtomHandleList::const_iterator i=atoms.begin(),
- e=atoms.end(); i!=e; ++i) {
- if (int((*i).Impl()->GetState())==std::numeric_limits<int>::max()) {
- this->DistanceBasedConnect(*i);
- }
- }
}
void RuleBasedBuilder::ConnectResidueToNext(mol::ResidueHandle rh,
- mol::ResidueHandle next) {
+ mol::ResidueHandle next)
+{
if (!next.IsValid()) {
return;
}
@@ -208,7 +248,8 @@ void RuleBasedBuilder::ConnectResidueToNext(mol::ResidueHandle rh,
}
-void RuleBasedBuilder::AssignTorsions(mol::ChainHandle chain) {
+void RuleBasedBuilder::AssignTorsions(mol::ChainHandle chain)
+{
if (chain.GetResidueCount()==0)
return;
std::vector<mol::ResidueHandle> rlist = chain.GetResidueList();
@@ -219,7 +260,8 @@ void RuleBasedBuilder::AssignTorsions(mol::ChainHandle chain) {
}
}
-void RuleBasedBuilder::AssignTorsionsToResidue(mol::ResidueHandle residue) {
+void RuleBasedBuilder::AssignTorsionsToResidue(mol::ResidueHandle residue)
+{
/// The only components having named torsions are the standard set of amino
/// acids, plus some of compounds derived from them such as selenium
/// methionine. Things are simplified a lot by only storing the torsions
diff --git a/modules/conop/src/rule_based_builder.hh b/modules/conop/src/rule_based_builder.hh
index 0ab392b24..d04f42272 100644
--- a/modules/conop/src/rule_based_builder.hh
+++ b/modules/conop/src/rule_based_builder.hh
@@ -108,14 +108,16 @@ public:
/// \brief Set residue properties such as chemical class
virtual void FillResidueProps(mol::ResidueHandle residue);
-
-
+ /// \brief whether the residue has unknown atoms
+ bool HasUnknownAtoms(mol::ResidueHandle res);
/// \brief Check whether the residue has all required atoms. This does not
/// include hydrogens and leaving atoms such as the terminal OXT.
virtual bool IsResidueComplete(const mol::ResidueHandle& residue);
private:
CompoundLibPtr compound_lib_;
CompoundPtr last_compound_;
+ mol::ResidueHandle last_residue_;
+ bool unknown_atoms_;
void LookupCompound(const mol::ResidueHandle& rh);
/// Change internal order of atoms in residue to the order given by compound
void ReorderAtoms(mol::ResidueHandle residue, CompoundPtr compound);
--
GitLab