diff --git a/modules/conop/pymod/export_builder.cc b/modules/conop/pymod/export_builder.cc
index 4bd468f36b23180c726f661f34a89227ee7168c0..28608ab6ee36e49f6ad539d56b210be6c5cd9cb2 100644
--- a/modules/conop/pymod/export_builder.cc
+++ b/modules/conop/pymod/export_builder.cc
@@ -36,6 +36,8 @@ void export_Builder() {
.add_property("dialect", &Builder::GetDialect, &Builder::SetDialect)
.add_property("strict_hydrogens", &Builder::GetStrictHydrogenMode,
&Builder::SetStrictHydrogenMode)
+ .add_property("feasibility_check", &Builder::GetBondFeasibilityCheck,
+ &Builder::SetBondFeasibilityCheck)
.def("GetDialect", &Builder::GetDialect)
.def("SetDialect", &Builder::SetDialect)
.def("CompleteAtoms", &Builder::CompleteAtoms)
@@ -48,6 +50,8 @@ void export_Builder() {
.def("AssignTorsionsToResidue", &Builder::AssignTorsionsToResidue)
.def("FillAtomProps", &Builder::FillAtomProps)
.def("IsResidueComplete", &Builder::IsResidueComplete)
+ .def("SetBondFeasibilityFlag", &Builder::SetBondFeasibilityCheck)
+ .def("GetBondFeasibilityFlag", &Builder::GetBondFeasibilityCheck)
;
class_<HeuristicBuilder, bases<Builder> >("HeuristicBuilder", init<>())
@@ -56,5 +60,6 @@ void export_Builder() {
init<const CompoundLibPtr&>())
.add_property("compound_lib", &RuleBasedBuilder::GetCompoundLib)
.def("GetUnknownAtoms", &RuleBasedBuilder::GetUnknownAtoms)
+
;
}
diff --git a/modules/conop/src/builder.hh b/modules/conop/src/builder.hh
index f912aba702cd075768da77c1c837f558877f3fda..4b915681a51d2d47c6d36fe449424f5b189de772 100644
--- a/modules/conop/src/builder.hh
+++ b/modules/conop/src/builder.hh
@@ -54,7 +54,7 @@ typedef enum {
class DLLEXPORT_OST_CONOP Builder {
public:
- Builder(): dialect_(PDB_DIALECT), strict_(false) { }
+ Builder(): dialect_(PDB_DIALECT), strict_(false), bond_feasibility_check_(true) { }
virtual ~Builder();
/// \brief add any missing atoms to the residue based on its key,
@@ -143,9 +143,17 @@ public:
/// |brief Connect \p atom with all atoms for whith IsBondFeasible() and
/// AreResiduesConsecutive() returns true
void DistanceBasedConnect(mol::AtomHandle atom);
+
+ /// \brief Set bond feasibility check flag
+ void SetBondFeasibilityCheck(bool b_feas_flag) { bond_feasibility_check_ = b_feas_flag; }
+
+ /// \brief Get bond feasibility check flag
+ bool GetBondFeasibilityCheck() const { return bond_feasibility_check_; }
+
private:
Dialect dialect_;
bool strict_;
+ bool bond_feasibility_check_;
};
diff --git a/modules/conop/src/heuristic_builder.cc b/modules/conop/src/heuristic_builder.cc
index 89b8ad8e1720754b0941e8b1a7d8977e252c14fc..d0a82e6601fd207de856e87b482b039e506657d0 100644
--- a/modules/conop/src/heuristic_builder.cc
+++ b/modules/conop/src/heuristic_builder.cc
@@ -209,21 +209,37 @@ void HeuristicBuilder::ConnectivityFromAtomNames(const mol::ResidueHandle& res,
<< it2->GetName() << ") in connectivity table of "
<< res.GetKey() << "... ");
int conn=centry.Check(it1->GetName(),it2->GetName());
- if (conn==1 && this->IsBondFeasible(*it1, *it2)) {
- LOG_TRACE( "found");
- editor.Connect(*it1,*it2);
- } else if(conn==2 && this->IsBondFeasible(*it2, *it1)) {
- LOG_TRACE( "found (reversed)");
- editor.Connect(*it2,*it1);
+ if (conn==1) {
+ if (this->GetBondFeasibilityCheck()==false) {
+ LOG_TRACE( "found");
+ editor.Connect(*it1,*it2);
+ } else {
+ if (this->IsBondFeasible(*it1, *it2)) {
+ LOG_TRACE( "found");
+ editor.Connect(*it1,*it2);
+ } else {
+ LOG_TRACE( "not found");
+ }
+ }
+ } else if (conn==2) {
+ if (this->GetBondFeasibilityCheck()==false) {
+ LOG_TRACE( "found (reversed)");
+ editor.Connect(*it2,*it1);
+ } else {
+ if(this->IsBondFeasible(*it2, *it1)) {
+ LOG_TRACE( "found (reversed)");
+ editor.Connect(*it2,*it1);
+ } else {
+ LOG_TRACE( "not found");
+ }
+ }
} else {
- LOG_TRACE( "not found");
- }
+ unknown_atoms.push_back(*it1);
+ }
}
- } else {
- unknown_atoms.push_back(*it1);
}
}
-}
+}
void HeuristicBuilder::ConnectAtomsOfResidue(mol::ResidueHandle res)
{
diff --git a/modules/conop/src/rule_based_builder.cc b/modules/conop/src/rule_based_builder.cc
index 35709c287d61181fa0beb5c5920a48ce7e8326b3..33d33c667faebb47593eebe53f28285c448b1956 100644
--- a/modules/conop/src/rule_based_builder.cc
+++ b/modules/conop/src/rule_based_builder.cc
@@ -237,12 +237,20 @@ void RuleBasedBuilder::ConnectAtomsOfResidue(mol::ResidueHandle rh)
const BondSpec& bond=*j;
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)) {
- if (this->GetStrictHydrogenMode() &&
- (a1.GetElement()=="H" || a2.GetElement()=="D")) {
- continue;
+ if (a1.IsValid() && a2.IsValid()) {
+ if (this->GetBondFeasibilityCheck()==false) {
+ if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" || a2.GetElement()=="D")) {
+ continue;
+ }
+ e.Connect(a1, a2, bond.order);
+ } else {
+ if (IsBondFeasible(a1, a2)) {
+ if (this->GetStrictHydrogenMode() && (a1.GetElement()=="H" || a2.GetElement()=="D")) {
+ continue;
+ }
+ e.Connect(a1, a2, bond.order);
+ }
}
- e.Connect(a1, a2, bond.order);
}
}
for (mol::AtomHandleList::iterator i=atoms.begin(), e=atoms.end(); i!=e; ++i) {
diff --git a/modules/conop/src/rule_based_builder.hh b/modules/conop/src/rule_based_builder.hh
index 9ecc45c101cbbb2c6777ff7fb985bb2e9e921db1..71a6fd03453ba9c418a6e39d8c0e011161f13be6 100644
--- a/modules/conop/src/rule_based_builder.hh
+++ b/modules/conop/src/rule_based_builder.hh
@@ -109,7 +109,7 @@ public:
/// by ignoring it or by inserting a dummy atom.
virtual void OnMissingAtom(const mol::ResidueHandle& residue,
const String& atom_name) { }
-
+
/// \brief Fill in missing information based on atom name.
virtual void FillAtomProps(mol::AtomHandle atom, const AtomSpec& spec);
@@ -123,6 +123,7 @@ public:
virtual bool IsResidueComplete(const mol::ResidueHandle& residue);
CompoundLibPtr GetCompoundLib() const { return compound_lib_; }
+
private:
CompoundLibPtr compound_lib_;
CompoundPtr last_compound_;
@@ -139,7 +140,6 @@ private:
void AssignBackBoneTorsionsToResidue(mol::ResidueHandle residue);
-
};
typedef boost::shared_ptr<RuleBasedBuilder> RuleBasedBuilderPtr;
diff --git a/modules/conop/tests/test_heuristic_builder.cc b/modules/conop/tests/test_heuristic_builder.cc
index 4220f3110ca1af0ecec5acc183dba6d1ad5b5f11..0e3a811ae922d0d8375fe98ab96163765f5766f4 100644
--- a/modules/conop/tests/test_heuristic_builder.cc
+++ b/modules/conop/tests/test_heuristic_builder.cc
@@ -65,6 +65,21 @@ ResidueHandle make_leu(ChainHandle chain)
return res;
}
+ResidueHandle make_defective_leu(ChainHandle chain)
+{
+ XCSEditor e=chain.GetEntity().EditXCS();
+ ResidueHandle res=e.AppendResidue(chain, "LEU");
+
+ e.InsertAtom(res, "N", geom::Vec3(19.003,32.473,60.366));
+ e.InsertAtom(res, "CA", geom::Vec3(18.330,32.402,61.664));
+ e.InsertAtom(res, "C", geom::Vec3(17.884,33.787,62.117));
+ e.InsertAtom(res, "O", geom::Vec3(17.853,34.091,63.308));
+ e.InsertAtom(res, "CB", geom::Vec3(19.269,31.793,102.710));
+ e.InsertAtom(res, "CG", geom::Vec3(19.695,30.340,62.501));
+ e.InsertAtom(res, "CD1", geom::Vec3(20.585,29.897,63.648));
+ e.InsertAtom(res, "CD2", geom::Vec3(18.461,29.459,62.420));
+ return res;
+}
void verify_connectivity_x(const ResidueHandle& res)
{
@@ -154,12 +169,25 @@ BOOST_AUTO_TEST_CASE(name_based_connect)
heuristic_builder.FillAtomProps(*i);
}
+ EntityHandle de=CreateEntity();
+ ChainHandle dc=de.EditXCS().InsertChain("A");
+ ResidueHandle dile=make_defective_leu(dc);
+ HeuristicBuilder dheuristic_builder;
+ dheuristic_builder.SetBondFeasibilityCheck(false);
+ for (AtomHandleIter i=de.AtomsBegin(),x=de.AtomsEnd(); i!=x; ++i) {
+ dheuristic_builder.FillAtomProps(*i);
+ }
+
BOOST_MESSAGE("running distance based checks on arginine");
heuristic_builder.ConnectAtomsOfResidue(arg);
verify_connectivity(arg);
BOOST_MESSAGE("running distance based checks on leu");
heuristic_builder.ConnectAtomsOfResidue(ile);
verify_connectivity(ile);
+
+ BOOST_MESSAGE("running distance based checks on defective leu");
+ dheuristic_builder.ConnectAtomsOfResidue(dile);
+ verify_connectivity(dile);
}
BOOST_AUTO_TEST_CASE(test_assign_torsions){
diff --git a/modules/conop/tests/test_rule_based_builder.cc b/modules/conop/tests/test_rule_based_builder.cc
index 0bc08d4047027a9d7895a6ae6e824b4358155ff0..50e34860a606912aa978fc3429dc1b1f87c2656a 100644
--- a/modules/conop/tests/test_rule_based_builder.cc
+++ b/modules/conop/tests/test_rule_based_builder.cc
@@ -119,6 +119,35 @@ ResidueHandle make_uracil2(ChainHandle chain)
return res;
}
+ResidueHandle make_defective_uracil2(ChainHandle chain)
+{
+ XCSEditor e=chain.GetEntity().EditXCS();
+ ResidueHandle res = e.AppendResidue(chain, "U");
+
+ e.InsertAtom(res, "P", geom::Vec3(16.249, 28.254, 44.759));
+ e.InsertAtom(res, "OP1", geom::Vec3(16.654, 28.055, 46.181));
+ e.InsertAtom(res, "OP2", geom::Vec3(17.115, 27.731, 43.656));
+ e.InsertAtom(res, "O5'", geom::Vec3(14.826, 27.572, 44.574));
+ e.InsertAtom(res, "C5'", geom::Vec3(14.711, 26.206, 44.216));
+ e.InsertAtom(res, "C4'", geom::Vec3(13.279, 25.889, 43.887));
+ e.InsertAtom(res, "O4'", geom::Vec3(12.832, 26.793, 42.838));
+ e.InsertAtom(res, "C3'", geom::Vec3(13.155, 24.434, 43.329));
+ e.InsertAtom(res, "O3'", geom::Vec3(12.269, 23.625, 44.098));
+ e.InsertAtom(res, "C2'", geom::Vec3(12.871, 24.595, 41.875));
+ e.InsertAtom(res, "O2'", geom::Vec3(11.811, 23.752, 41.462));
+ e.InsertAtom(res, "C1'", geom::Vec3(12.424, 26.056, 41.694));
+ e.InsertAtom(res, "N1", geom::Vec3(13.030, 26.692, 40.497));
+ e.InsertAtom(res, "C2", geom::Vec3(12.517, 26.365, 39.228));
+ e.InsertAtom(res, "O2", geom::Vec3(11.579, 25.594, 39.068));
+ e.InsertAtom(res, "N3", geom::Vec3(13.141, 26.987, 38.161));
+ e.InsertAtom(res, "C4", geom::Vec3(14.197, 27.888, 38.210));
+ e.InsertAtom(res, "O4", geom::Vec3(14.627, 28.368, 37.156));
+ e.InsertAtom(res, "C5", geom::Vec3(14.671, 28.189, 39.542));
+ e.InsertAtom(res, "C6", geom::Vec3(14.087, 27.597, 80.612));
+
+ return res;
+}
+
void verify_nucleotide_connectivity(const ResidueHandle& res)
{
BOOST_CHECK(BondExists(res.FindAtom("P"),
@@ -217,16 +246,28 @@ BOOST_AUTO_TEST_CASE(nucleotide_based_connect)
}
RuleBasedBuilder rb_builder = RuleBasedBuilder(compound_lib);
+ RuleBasedBuilder drb_builder = RuleBasedBuilder(compound_lib);
+ drb_builder.SetBondFeasibilityCheck(false);
+
EntityHandle e=CreateEntity();
ChainHandle c=e.EditXCS().InsertChain("A");
ResidueHandle c0=make_cytosine(c);
ResidueHandle u1=make_uracil1(c);
ResidueHandle u2=make_uracil2(c);
+
+ EntityHandle de=CreateEntity();
+ ChainHandle dc=de.EditXCS().InsertChain("A");
+ ResidueHandle du2=make_defective_uracil2(dc);
+
for (AtomHandleIter i=e.AtomsBegin(),x=e.AtomsEnd(); i!=x; ++i) {
rb_builder.FillAtomProps(*i);
}
+ for (AtomHandleIter i=de.AtomsBegin(),x=de.AtomsEnd(); i!=x; ++i) {
+ drb_builder.FillAtomProps(*i);
+ }
+
// running positive test
BOOST_MESSAGE("running distance based checks on cytosine");
rb_builder.ConnectAtomsOfResidue(c0);
@@ -247,6 +288,12 @@ BOOST_AUTO_TEST_CASE(nucleotide_based_connect)
BOOST_MESSAGE("connecting cytosine to second uracil");
rb_builder.ConnectResidueToNext(c0, u2);
verify_nucleotide_nolink(c0, u2);
+
+ // running positive test
+ BOOST_MESSAGE("running distance based checks on defective uracil");
+ drb_builder.ConnectAtomsOfResidue(du2);
+ verify_nucleotide_connectivity(du2);
+
}
BOOST_AUTO_TEST_SUITE_END( )
diff --git a/modules/io/pymod/__init__.py b/modules/io/pymod/__init__.py
index 143fc2fe22545ec56667ba17602310b8626bce64..ab19f54eedafd79ff57b017e4bc94b79011dea86 100644
--- a/modules/io/pymod/__init__.py
+++ b/modules/io/pymod/__init__.py
@@ -86,7 +86,7 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
fault_tolerant=None, load_multi=False, quack_mode=None,
join_spread_atom_records=None, calpha_only=None,
profile='DEFAULT', remote=False, dialect=None,
- strict_hydrogens=None, seqres=False):
+ strict_hydrogens=None, seqres=False, bond_feasibility=None):
"""
Load PDB file from disk and return one or more entities. Several options
allow to customize the exact behaviour of the PDB import. For more information
@@ -147,6 +147,7 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
prof.quack_mode=_override(prof.quack_mode, quack_mode)
prof.strict_hydrogens=_override(prof.strict_hydrogens, strict_hydrogens)
prof.fault_tolerant=_override(prof.fault_tolerant, fault_tolerant)
+ prof.bond_feasibility=_override(prof.bond_feasibility, bond_feasibility)
prof.join_spread_atom_records=_override(prof.join_spread_atom_records,
join_spread_atom_records)
@@ -164,6 +165,7 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
elif prof.dialect=='CHARMM':
builder.dialect=conop.CHARMM_DIALECT
builder.strict_hydrogens=prof.strict_hydrogens
+ builder.bond_feasibility=prof.bond_feasibility
reader=PDBReader(filename, prof)
reader.read_seqres=seqres
try:
diff --git a/modules/io/pymod/export_pdb_io.cc b/modules/io/pymod/export_pdb_io.cc
index e431239e2e7d591d508f05d01464be380f468b49..6f7566906d43c81d8b1921168fadbf551c9f7a11 100644
--- a/modules/io/pymod/export_pdb_io.cc
+++ b/modules/io/pymod/export_pdb_io.cc
@@ -36,13 +36,14 @@ void (PDBWriter::*write_b)(const mol::EntityView&)=&PDBWriter::Write;
void export_pdb_io()
{
class_<IOProfile>("IOProfile",
- init<String,bool,bool,bool,bool,bool,bool>((arg("dialect")="PDB",
- arg("strict_hydrogens")=false,
- arg("quack_mode")=false,
- arg("fault_tolerant")=false,
- arg("join_spread_atom_records")=false,
- arg("no_hetatms")=false,
- arg("calpha_only")=false)))
+ init<String,bool,bool,bool,bool,bool,bool,bool>((arg("dialect")="PDB",
+ arg("strict_hydrogens")=false,
+ arg("quack_mode")=false,
+ arg("fault_tolerant")=false,
+ arg("join_spread_atom_records")=false,
+ arg("no_hetatms")=false,
+ arg("calpha_only")=false,
+ arg("bond_feasibility")=true)))
.def_readwrite("dialect", &IOProfile::dialect)
.def_readwrite("fault_tolerant", &IOProfile::fault_tolerant)
.def_readwrite("quack_mode", &IOProfile::quack_mode)
@@ -50,6 +51,7 @@ void export_pdb_io()
.def_readwrite("no_hetatms", &IOProfile::no_hetatms)
.def_readwrite("calpha_only", &IOProfile::calpha_only)
.def_readwrite("join_spread_atom_records", &IOProfile::join_spread_atom_records)
+ .def_readwrite("bond_feasibility", &IOProfile::bond_feasibilty)
.def("Copy", &IOProfile::Copy)
.def(self_ns::str(self))
;
diff --git a/modules/io/src/mol/io_profile.hh b/modules/io/src/mol/io_profile.hh
index 91fe27b7793350e10d29d922cf0ca501a1697da5..a51572636e43cc1da175929b00d94e5c290b855e 100644
--- a/modules/io/src/mol/io_profile.hh
+++ b/modules/io/src/mol/io_profile.hh
@@ -27,13 +27,13 @@ namespace ost { namespace io {
struct DLLEXPORT IOProfile {
public:
- IOProfile(String d, bool sh, bool qm, bool ft, bool js, bool nh, bool co):
+ IOProfile(String d, bool sh, bool qm, bool ft, bool js, bool nh, bool co, bool bf):
dialect(d), strict_hydrogens(sh), quack_mode(qm), fault_tolerant(ft),
- join_spread_atom_records(js), no_hetatms(nh), calpha_only(co)
+ join_spread_atom_records(js), no_hetatms(nh), calpha_only(co), bond_feasibilty(bf)
{ }
IOProfile(): dialect("PDB"), strict_hydrogens(true), quack_mode(false),
fault_tolerant(false), join_spread_atom_records(false), no_hetatms(false),
- calpha_only(false)
+ calpha_only(false), bond_feasibilty(true)
{ }
virtual ~IOProfile() { }
@@ -44,11 +44,12 @@ public:
bool join_spread_atom_records;
bool no_hetatms;
bool calpha_only;
+ bool bond_feasibilty;
IOProfile Copy()
{
return IOProfile(dialect, strict_hydrogens, quack_mode, fault_tolerant,
- join_spread_atom_records, no_hetatms, calpha_only);
+ join_spread_atom_records, no_hetatms, calpha_only, bond_feasibilty);
}
virtual void PostImport(mol::EntityHandle ent) { }
};
@@ -61,7 +62,8 @@ inline std::ostream& operator<<(std::ostream& stream, const IOProfile& p)
<< (p.join_spread_atom_records ? "True" : "False") << ", no_hetatms="
<< (p.no_hetatms ? "True" : "False") << ", calpha_only="
<< (p.calpha_only ? "True" : "False") << ", fault_tolerant="
- << (p.fault_tolerant ? "True" : "False") << ")";
+ << (p.fault_tolerant ? "True" : "False") << ", bond_feasibilty="
+ << (p.bond_feasibilty ? "True" : "False") << ")";
return stream;
}
diff --git a/modules/io/tests/test_io_pdb.cc b/modules/io/tests/test_io_pdb.cc
index 5bbfbf699a9684f0eee761fbe04f42698f77daa2..48ba5667caf3a5dfe279f5d2852c52693bf3e02e 100644
--- a/modules/io/tests/test_io_pdb.cc
+++ b/modules/io/tests/test_io_pdb.cc
@@ -836,7 +836,7 @@ BOOST_AUTO_TEST_CASE(charmm_rname)
{
PDBWriter writer(String("testfiles/pdb/charmm_rname-out.pdb"),
IOProfile("CHARMM", true, false, false,
- false, false, false));
+ false, false, false, true));
mol::EntityHandle ent=mol::CreateEntity();
mol::XCSEditor edi=ent.EditXCS();
@@ -856,7 +856,7 @@ BOOST_AUTO_TEST_CASE(charmm_longcname)
{
PDBWriter writer(String("testfiles/pdb/charmm_longcname-out.pdb"),
IOProfile("CHARMM", true, false, false,
- false, false, false));
+ false, false, false, true));
mol::EntityHandle ent=mol::CreateEntity();
mol::XCSEditor edi=ent.EditXCS();
@@ -876,7 +876,7 @@ BOOST_AUTO_TEST_CASE(write_charmm_ter)
{
PDBWriter writer(String("testfiles/pdb/charmm_ter-out.pdb"),
IOProfile("CHARMM", true, false, false,
- false, false, false));
+ false, false, false, true));
mol::EntityHandle ent=mol::CreateEntity();
mol::XCSEditor edi=ent.EditXCS();