Skip to content
Snippets Groups Projects
Commit d60398e7 authored by Marco Biasini's avatar Marco Biasini
Browse files

fix BZDNG-261

Before falling back to alternative atom names as defined by the
components.cif dictionary, skim through all primary names. Some
residues in PDB files reuse alternative atom names as primary
names for other atoms. This resulted in wrong connectivity in,
for example, 1hiv.
parent c5e11cb1
Branches
Tags
No related merge requests found
...@@ -23,10 +23,17 @@ namespace ost { namespace conop { ...@@ -23,10 +23,17 @@ namespace ost { namespace conop {
int Compound::GetAtomSpecIndex(const String& name) const { int Compound::GetAtomSpecIndex(const String& name) const {
AtomSpecList::const_iterator i=atom_specs_.begin(); AtomSpecList::const_iterator i=atom_specs_.begin();
// BZDNG-261: first search all primary atom names before falling back to
// alternative names. There are some files where alternative atom names are
// used as primary names for other atoms
for (; i!=atom_specs_.end(); ++i) { for (; i!=atom_specs_.end(); ++i) {
if ((*i).name==name || (*i).alt_name==name) if ((*i).name==name)
return std::distance(atom_specs_.begin(), i);
}
for (; i!=atom_specs_.end(); ++i) {
if ((*i).alt_name==name)
return std::distance(atom_specs_.begin(), i); return std::distance(atom_specs_.begin(), i);
} }
return -1; return -1;
} }
}} }}
...@@ -238,18 +238,20 @@ void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh) ...@@ -238,18 +238,20 @@ void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh)
mol::AtomHandle a1=this->LocateAtom(atoms, bond.atom_one); mol::AtomHandle a1=this->LocateAtom(atoms, bond.atom_one);
mol::AtomHandle a2=this->LocateAtom(atoms, bond.atom_two); mol::AtomHandle a2=this->LocateAtom(atoms, bond.atom_two);
if (a1.IsValid() && a2.IsValid()) { if (a1.IsValid() && a2.IsValid()) {
if (this->GetBondFeasibilityCheck()==false) { if (this->GetBondFeasibilityCheck()==false) {
if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" || a2.GetElement()=="D")) { if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" ||
a2.GetElement()=="D")) {
continue; continue;
} }
e.Connect(a1, a2, bond.order); e.Connect(a1, a2, bond.order);
} else { } else {
if (IsBondFeasible(a1, a2)) { if (IsBondFeasible(a1, a2)) {
if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" || a2.GetElement()=="D")) { if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" ||
a2.GetElement()=="D")) {
continue; continue;
} }
e.Connect(a1, a2, bond.order); e.Connect(a1, a2, bond.order);
} }
} }
} }
} }
......
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
#include <ost/mol/mol.hh> #include <ost/mol/mol.hh>
#include <ost/platform.hh> #include <ost/platform.hh>
#include <ost/conop/rule_based_builder.hh> #include <ost/conop/rule_based_builder.hh>
#include <ost/conop/conop.hh>
#define BOOST_TEST_DYN_LINK #define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp> #include <boost/test/unit_test.hpp>
#include <boost/test/auto_unit_test.hpp> #include <boost/test/auto_unit_test.hpp>
...@@ -149,6 +149,139 @@ ResidueHandle make_defective_uracil2(ChainHandle chain) ...@@ -149,6 +149,139 @@ ResidueHandle make_defective_uracil2(ChainHandle chain)
return res; return res;
} }
ResidueHandle make_1zk(ChainHandle chain)
{
XCSEditor edi=chain.GetEntity().EditXCS();
ResidueHandle r=edi.AppendResidue(chain, "1ZK");
edi.InsertAtom(r, "C", geom::Vec3(3.946, 3.520, 10.861),"C");
edi.InsertAtom(r, "C1", geom::Vec3(6.308, 3.824, 8.036),"C");
edi.InsertAtom(r, "CA", geom::Vec3(4.856, 4.282, 9.931),"C");
edi.InsertAtom(r, "C2", geom::Vec3(6.770, 5.099, 8.019),"C");
edi.InsertAtom(r, "C3", geom::Vec3(7.805, 5.484, 7.142),"C");
edi.InsertAtom(r, "C4", geom::Vec3(8.363, 4.590, 6.284),"C");
edi.InsertAtom(r, "C5", geom::Vec3(8.469, 2.345, 5.318),"C");
edi.InsertAtom(r, "C6", geom::Vec3(7.986, 1.072, 5.322),"C");
edi.InsertAtom(r, "C4A", geom::Vec3(7.885, 3.277, 6.208),"C");
edi.InsertAtom(r, "C7", geom::Vec3(6.949, 0.682, 6.190),"C");
edi.InsertAtom(r, "C8", geom::Vec3(6.398, 1.548, 7.074),"C");
edi.InsertAtom(r, "C9", geom::Vec3(4.018, 0.269, 12.183),"C");
edi.InsertAtom(r, "C1A", geom::Vec3(6.824, 2.882, 7.117),"C");
edi.InsertAtom(r, "O", geom::Vec3(3.224, 4.213, 11.565),"O");
edi.InsertAtom(r, "O1", geom::Vec3(5.331, 3.476, 8.904),"O");
edi.InsertAtom(r, "O3", geom::Vec3(5.748, -4.044, 14.612),"O");
edi.InsertAtom(r, "N", geom::Vec3(3.855, 2.201, 10.790),"N");
edi.InsertAtom(r, "N1", geom::Vec3(4.814, 0.719, 13.173),"N");
edi.InsertAtom(r, "CA1", geom::Vec3(3.077, 1.346, 11.644),"C");
edi.InsertAtom(r, "O2", geom::Vec3(4.203, -0.810, 11.651),"O");
edi.InsertAtom(r, "O4", geom::Vec3(7.311, -5.667, 18.880),"O");
edi.InsertAtom(r, "CB", geom::Vec3(1.856, 0.712, 11.049),"C");
edi.InsertAtom(r, "CG", geom::Vec3(1.015, 1.845, 10.511),"C");
edi.InsertAtom(r, "ND1", geom::Vec3(1.467, 2.439, 9.321),"N");
edi.InsertAtom(r, "N2", geom::Vec3(6.478, -3.958, 16.751),"N");
edi.InsertAtom(r, "CD2", geom::Vec3(-0.105, 2.479, 10.887),"C");
edi.InsertAtom(r, "CE1", geom::Vec3(0.638, 3.428, 9.002),"C");
edi.InsertAtom(r, "NE2", geom::Vec3(-0.372, 3.461, 9.881),"N");
edi.InsertAtom(r, "N3", geom::Vec3(6.871, -7.136, 17.332),"N");
edi.InsertAtom(r, "CA2", geom::Vec3(5.881, -0.001, 13.808),"C");
edi.InsertAtom(r, "CB1", geom::Vec3(7.140, 0.860, 13.743),"C");
edi.InsertAtom(r, "CG1", geom::Vec3(7.503, 1.324, 12.299),"C");
edi.InsertAtom(r, "C21", geom::Vec3(5.126, -8.557, 18.179),"C");
edi.InsertAtom(r, "CD1", geom::Vec3(8.185, 0.142, 11.507),"C");
edi.InsertAtom(r, "CD21", geom::Vec3(8.420, 2.537, 12.325),"C");
edi.InsertAtom(r, "CE11", geom::Vec3(8.381, 0.689, 10.066),"C");
edi.InsertAtom(r, "CE2", geom::Vec3(8.907, 2.979, 10.922),"C");
edi.InsertAtom(r, "CZ", geom::Vec3(9.409, 1.807, 10.075),"C");
edi.InsertAtom(r, "CH", geom::Vec3(5.592, -0.511, 15.204),"C");
edi.InsertAtom(r, "OH", geom::Vec3(5.225, 0.377, 16.238),"O");
edi.InsertAtom(r, "CB11", geom::Vec3(4.426, -1.543, 15.170),"C");
edi.InsertAtom(r, "CA'", geom::Vec3(4.451, -2.730, 16.152),"C");
edi.InsertAtom(r, "CB'", geom::Vec3(3.124, -3.441, 16.281),"C");
edi.InsertAtom(r, "CG11", geom::Vec3(2.553, -3.986, 14.933),"C");
edi.InsertAtom(r, "C31", geom::Vec3(4.413, -7.811, 19.117),"C");
edi.InsertAtom(r, "CG2", geom::Vec3(3.204, -4.586, 17.345),"C");
edi.InsertAtom(r, "OB1", geom::Vec3(3.249, -0.875, 15.134),"O");
edi.InsertAtom(r, "CC", geom::Vec3(5.603, -3.655, 15.782),"C");
edi.InsertAtom(r, "CA3", geom::Vec3(7.592, -4.867, 16.603),"C");
edi.InsertAtom(r, "CD", geom::Vec3(7.274, -5.947, 17.691),"C");
edi.InsertAtom(r, "CB2", geom::Vec3(8.986, -4.351, 16.803),"C");
edi.InsertAtom(r, "CG12", geom::Vec3(9.488, -3.108, 16.016),"C");
edi.InsertAtom(r, "CG21", geom::Vec3(10.000, -5.461, 16.472),"C");
edi.InsertAtom(r, "CD11", geom::Vec3(9.099, -3.257, 14.571),"C");
edi.InsertAtom(r, "CM", geom::Vec3(6.587, -8.286, 18.106),"C");
edi.InsertAtom(r, "C41", geom::Vec3(3.045, -7.980, 19.287),"C");
edi.InsertAtom(r, "C51", geom::Vec3(2.423, -8.911, 18.456),"C");
edi.InsertAtom(r, "C61", geom::Vec3(3.164, -9.631, 17.518),"C");
edi.InsertAtom(r, "N11", geom::Vec3(4.497, -9.459, 17.386),"N");
return r;
}
void verify_1zk_connectivity(const ResidueHandle& r1)
{
BOOST_CHECK(BondExists(r1.FindAtom("C"), r1.FindAtom("CA")));
BOOST_CHECK(BondExists(r1.FindAtom("C"), r1.FindAtom("O")));
BOOST_CHECK(BondExists(r1.FindAtom("CA"), r1.FindAtom("O1")));
BOOST_CHECK(BondExists(r1.FindAtom("C1"), r1.FindAtom("C2")));
BOOST_CHECK(BondExists(r1.FindAtom("C1"), r1.FindAtom("C1A")));
BOOST_CHECK(BondExists(r1.FindAtom("C1"), r1.FindAtom("O1")));
BOOST_CHECK(BondExists(r1.FindAtom("C2"), r1.FindAtom("C3")));
BOOST_CHECK(BondExists(r1.FindAtom("C3"), r1.FindAtom("C4")));
BOOST_CHECK(BondExists(r1.FindAtom("C4"), r1.FindAtom("C4A")));
BOOST_CHECK(BondExists(r1.FindAtom("C4A"), r1.FindAtom("C5")));
BOOST_CHECK(BondExists(r1.FindAtom("C4A"), r1.FindAtom("C1A")));
BOOST_CHECK(BondExists(r1.FindAtom("C5"), r1.FindAtom("C6")));
BOOST_CHECK(BondExists(r1.FindAtom("C6"), r1.FindAtom("C7")));
BOOST_CHECK(BondExists(r1.FindAtom("C7"), r1.FindAtom("C8")));
BOOST_CHECK(BondExists(r1.FindAtom("C8"), r1.FindAtom("C1A")));
BOOST_CHECK(BondExists(r1.FindAtom("N"), r1.FindAtom("CA1")));
BOOST_CHECK(BondExists(r1.FindAtom("CA1"), r1.FindAtom("C9")));
BOOST_CHECK(BondExists(r1.FindAtom("CA1"), r1.FindAtom("CB")));
BOOST_CHECK(BondExists(r1.FindAtom("C9"), r1.FindAtom("O2")));
BOOST_CHECK(BondExists(r1.FindAtom("CB"), r1.FindAtom("CG")));
BOOST_CHECK(BondExists(r1.FindAtom("CG"), r1.FindAtom("ND1")));
BOOST_CHECK(BondExists(r1.FindAtom("CG"), r1.FindAtom("CD2")));
BOOST_CHECK(BondExists(r1.FindAtom("ND1"), r1.FindAtom("CE1")));
BOOST_CHECK(BondExists(r1.FindAtom("CD2"), r1.FindAtom("NE2")));
BOOST_CHECK(BondExists(r1.FindAtom("CE1"), r1.FindAtom("NE2")));
BOOST_CHECK(BondExists(r1.FindAtom("N1"), r1.FindAtom("CA2")));
BOOST_CHECK(BondExists(r1.FindAtom("CA2"), r1.FindAtom("CB1")));
BOOST_CHECK(BondExists(r1.FindAtom("CA2"), r1.FindAtom("CH")));
BOOST_CHECK(BondExists(r1.FindAtom("CB1"), r1.FindAtom("CG1")));
BOOST_CHECK(BondExists(r1.FindAtom("CG1"), r1.FindAtom("CD1")));
BOOST_CHECK(BondExists(r1.FindAtom("CG1"), r1.FindAtom("CD21")));
BOOST_CHECK(BondExists(r1.FindAtom("CD1"), r1.FindAtom("CE11")));
BOOST_CHECK(BondExists(r1.FindAtom("CD21"), r1.FindAtom("CE2")));
BOOST_CHECK(BondExists(r1.FindAtom("CE11"), r1.FindAtom("CZ")));
BOOST_CHECK(BondExists(r1.FindAtom("CE2"), r1.FindAtom("CZ")));
BOOST_CHECK(BondExists(r1.FindAtom("CH"), r1.FindAtom("OH")));
BOOST_CHECK(BondExists(r1.FindAtom("CH"), r1.FindAtom("CB11")));
BOOST_CHECK(BondExists(r1.FindAtom("CB11"), r1.FindAtom("CA'")));
BOOST_CHECK(BondExists(r1.FindAtom("CB11"), r1.FindAtom("OB1")));
BOOST_CHECK(BondExists(r1.FindAtom("CA'"), r1.FindAtom("CB'")));
BOOST_CHECK(BondExists(r1.FindAtom("CA'"), r1.FindAtom("CC")));
BOOST_CHECK(BondExists(r1.FindAtom("CB'"), r1.FindAtom("CG11")));
BOOST_CHECK(BondExists(r1.FindAtom("CB'"), r1.FindAtom("CG2")));
BOOST_CHECK(BondExists(r1.FindAtom("CC"), r1.FindAtom("O3")));
BOOST_CHECK(BondExists(r1.FindAtom("N2"), r1.FindAtom("CA3")));
BOOST_CHECK(BondExists(r1.FindAtom("CA3"), r1.FindAtom("CD")));
BOOST_CHECK(BondExists(r1.FindAtom("CA3"), r1.FindAtom("CB2")));
BOOST_CHECK(BondExists(r1.FindAtom("CD"), r1.FindAtom("O4")));
BOOST_CHECK(BondExists(r1.FindAtom("CB2"), r1.FindAtom("CG12")));
BOOST_CHECK(BondExists(r1.FindAtom("CB2"), r1.FindAtom("CG21")));
BOOST_CHECK(BondExists(r1.FindAtom("CG12"), r1.FindAtom("CD11")));
BOOST_CHECK(BondExists(r1.FindAtom("N3"), r1.FindAtom("CM")));
BOOST_CHECK(BondExists(r1.FindAtom("CM"), r1.FindAtom("C21")));
BOOST_CHECK(BondExists(r1.FindAtom("C21"), r1.FindAtom("C31")));
BOOST_CHECK(BondExists(r1.FindAtom("C21"), r1.FindAtom("N11")));
BOOST_CHECK(BondExists(r1.FindAtom("C31"), r1.FindAtom("C41")));
BOOST_CHECK(BondExists(r1.FindAtom("C41"), r1.FindAtom("C51")));
BOOST_CHECK(BondExists(r1.FindAtom("C51"), r1.FindAtom("C61")));
BOOST_CHECK(BondExists(r1.FindAtom("C61"), r1.FindAtom("N11")));
BOOST_CHECK(BondExists(r1.FindAtom("C"), r1.FindAtom("N")));
BOOST_CHECK(BondExists(r1.FindAtom("C9"), r1.FindAtom("N1")));
BOOST_CHECK(BondExists(r1.FindAtom("CC"), r1.FindAtom("N2")));
BOOST_CHECK(BondExists(r1.FindAtom("CD"), r1.FindAtom("N3")));
}
void verify_nucleotide_connectivity(const ResidueHandle& res) void verify_nucleotide_connectivity(const ResidueHandle& res)
{ {
BOOST_CHECK(BondExists(res.FindAtom("P"), BOOST_CHECK(BondExists(res.FindAtom("P"),
...@@ -297,4 +430,27 @@ BOOST_AUTO_TEST_CASE(nucleotide_based_connect) ...@@ -297,4 +430,27 @@ BOOST_AUTO_TEST_CASE(nucleotide_based_connect)
} }
BOOST_AUTO_TEST_CASE(rule_based_connect_1zk)
{
SetPrefixPath(getenv("OST_ROOT"));
String lib_path=GetSharedDataPath()+"/compounds.chemlib";
CompoundLibPtr compound_lib=CompoundLib::Load(lib_path);
if (!compound_lib) {
std::cout << "WARNING: skipping NUCLEOTIDE_BASED connect unit test. "
<< "Rule-based builder is required" << std::endl;
return;
}
boost::shared_ptr<RuleBasedBuilder> drb_builder(new RuleBasedBuilder(compound_lib));
Conopology::Instance().RegisterBuilder(drb_builder, "RBB");
Conopology::Instance().SetDefaultBuilder("RBB");
drb_builder->SetBondFeasibilityCheck(false);
EntityHandle e=CreateEntity();
ChainHandle c=e.EditXCS().InsertChain("A");
ResidueHandle r1=make_1zk(c);
Conopology::Instance().ConnectAll(drb_builder, e);
Conopology::Instance().SetDefaultBuilder("HEURISTIC");
verify_1zk_connectivity(r1);
}
BOOST_AUTO_TEST_SUITE_END( ); BOOST_AUTO_TEST_SUITE_END( );
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment