diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 2989eb5da746c61f06f2c072f60b84c13ba98b12..81bdfe50c701abd142799a069331602b7e6cb44a 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -18,7 +18,7 @@
//------------------------------------------------------------------------------
#include <boost/python.hpp>
-#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
+#include <boost/python/suite/indexing/map_indexing_suite.hpp>
#include <ost/config.hh>
#include <ost/mol/alg/local_dist_test.hh>
#include <ost/mol/alg/superpose_frames.hh>
@@ -36,7 +36,7 @@ void export_entity_to_density();
namespace {
-Real (*ldt_a)(const mol::EntityView&, const mol::alg::GlobalDistanceList& , Real, const String&)=&mol::alg::LocalDistTest;
+Real (*ldt_a)(const mol::EntityView&, const mol::alg::GlobalRDMap& , std::vector<Real>, const String&)=&mol::alg::LocalDistTest;
Real (*ldt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistTest;
Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes;
@@ -70,8 +70,19 @@ ost::mol::alg::ClashingDistances fill_clashing_distances_wrapper (const list& st
return ost::mol::alg::FillClashingDistances(stereo_chemical_props_file_vector,min_default_distance,min_distance_tolerance);
}
+ost::mol::alg::GlobalRDMap create_distance_list_from_multiple_references(const list& ref_list, Real cutoff, Real max_dist)
+{
+ int ref_list_length = boost::python::extract<int>(ref_list.attr("__len__")());
+ std::vector<ost::mol::EntityView> ref_list_vector(ref_list_length);
+
+ for (int i=0; i<ref_list_length; i++) {
+ ref_list_vector[i] = boost::python::extract<ost::mol::EntityView>(ref_list[i]);
+ }
+
+ return ost::mol::alg::CreateDistanceListFromMultipleReferences(ref_list_vector, cutoff, max_dist);
}
+}
BOOST_PYTHON_MODULE(_ost_mol_alg)
@@ -92,6 +103,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("CheckStereoChemistry", csc_b, (arg("ent"), arg("bonds"), arg("angles"), arg("bond_tolerance"), arg("angle_tolerance"), arg("always_remove_bb")=false));
def("LDTHA",&mol::alg::LDTHA);
def("CreateDistanceList",&mol::alg::CreateDistanceList);
+ def("CreateDistanceListFromMultipleReferences",&create_distance_list_from_multiple_references);
def("SuperposeFrames", superpose_frames1,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
@@ -127,23 +139,18 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
;
- class_<mol::alg::ReferenceDistance> ("ReferenceDistance", init <const mol::alg::UniqueAtomIdentifier&,const mol::alg::UniqueAtomIdentifier&, Real, Real>())
- .def("GetFirstAtom",&mol::alg::ReferenceDistance::GetFirstAtom)
- .def("GetSecondAtom",&mol::alg::ReferenceDistance::GetSecondAtom)
- .def("GetMinDistance",&mol::alg::ReferenceDistance::GetMinDistance)
- .def("GetMaxDistance",&mol::alg::ReferenceDistance::GetMaxDistance)
- ;
-
- class_<std::vector<mol::alg::ReferenceDistance> >("ResidueDistanceList")
- .def(vector_indexing_suite<std::vector<mol::alg::ReferenceDistance > >())
+ class_<mol::alg::ResidueRDMap>("ResidueRDMap")
+ .def(map_indexing_suite<mol::alg::ResidueRDMap>())
;
- class_<std::vector<mol::alg::ResidueDistanceList> >("GlobalDistanceList")
- .def(vector_indexing_suite<std::vector<mol::alg::ResidueDistanceList > >())
+ class_<mol::alg::GlobalRDMap>("GlobalRDMap")
+ .def(map_indexing_suite<mol::alg::GlobalRDMap>())
;
def("FillClashingDistances",&fill_clashing_distances_wrapper);
def("FillStereoChemicalParams",&fill_stereochemical_params_wrapper);
def("IsStandardResidue",&mol::alg::IsStandardResidue);
+ def("PrintGlobalRDMap",&mol::alg::PrintGlobalRDMap);
+ def("PrintResidueRDMap",&mol::alg::PrintResidueRDMap);
}
diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc
index 709adeb9909ecf5dd0b7187155ed7e74272b0113..23f11fd5fe714b66768687039a46a0ec2c1dc583 100644
--- a/modules/mol/alg/src/ldt.cc
+++ b/modules/mol/alg/src/ldt.cc
@@ -77,7 +77,7 @@ void usage()
std::cerr << " -e print version" << std::endl;
}
-std::pair<int,int> compute_coverage (const EntityView& v,const GlobalDistanceList& glob_dist_list)
+std::pair<int,int> compute_coverage (const EntityView& v,const GlobalRDMap& glob_dist_list)
{
int second=0;
int first=0;
@@ -85,15 +85,12 @@ std::pair<int,int> compute_coverage (const EntityView& v,const GlobalDistanceLis
return std::make_pair<int,int>(0,1);
}
ChainView vchain=v.GetChainList()[0];
- for (std::vector<ResidueDistanceList>::const_iterator i=glob_dist_list.begin();i!=glob_dist_list.end();++i)
+ for (GlobalRDMap::const_iterator i=glob_dist_list.begin();i!=glob_dist_list.end();++i)
{
- ResNum rnum = (*i)[0].GetFirstAtom().GetResNum();
- String rname = (*i)[0].GetFirstAtom().GetResidueName();
- if (IsStandardResidue(rname)) {
- second++;
- if (vchain.FindResidue(rnum)) {
- first++;
- }
+ ResNum rnum = (*i).first;
+ second++;
+ if (vchain.FindResidue(rnum)) {
+ first++;
}
}
return std::make_pair<int,int>(first,second);
@@ -200,7 +197,7 @@ int main (int argc, char **argv)
}
files.pop_back();
EntityView ref_view=ref.Select(sel);
- GlobalDistanceList glob_dist_list = CreateDistanceList(ref_view,radius);
+ GlobalRDMap glob_dist_list = CreateDistanceList(ref_view,radius);
std::cout << "Verbosity level: " << verbosity_level << std::endl;
if (structural_checks) {
std::cout << "Stereo-chemical and steric clash checks: On " << std::endl;
@@ -296,22 +293,18 @@ int main (int argc, char **argv)
std::cout << "Global LDT score: 0.0" << std::endl;
return 0;
}
- Real ldt=LDTHA(v, glob_dist_list);
+ Real ldt=LDTHA(v, glob_dist_list);
std::cout << "Global LDT score: " << ldt << std::endl;
std::cout << "Local LDT Score:" << std::endl;
std::cout << "Chain\tResName\tResNum\tScore" << std::endl;
- String labels[]={"localldt0.5","localldt1","localldt2","ldtlocal4"};
for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
ResidueView ritv = *rit;
- Real ldt_local_sum = 0;
- if (ritv.HasProp("localldt0.5")) {
- for (int n=0; n<4; ++n) {
- ldt_local_sum+=ritv.GetFloatProp(labels[n]);
- }
- ldt_local_sum/=4.0;
- std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << ldt_local_sum << std::endl;
+ Real ldt_local = 0;
+ if (ritv.HasProp("localldt")) {
+ ldt_local=ritv.GetFloatProp("localldt");
}
+ std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << ldt_local << std::endl;
}
std::cout << std::endl;
}
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
index 89da8c9ae2255c3caa2793e0534761a8d6592288..9f508280d6cc5c26e615dee5dab8d77780fe1873 100644
--- a/modules/mol/alg/src/local_dist_test.cc
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -51,10 +51,10 @@ bool swappable(const String& rname, const String& aname)
return false;
}
-std::pair<bool,Real> within_tolerance(Real mdl_dist, const ReferenceDistance& ref_dist, Real tol)
+std::pair<bool,Real> within_tolerance(Real mdl_dist, std::pair<Real,Real>& values, Real tol)
{
- Real min = ref_dist.GetMaxDistance();
- Real max = ref_dist.GetMaxDistance();
+ Real min = values.first;
+ Real max = values.second;
bool within_tol = false;
Real difference = 0;
if (mdl_dist>=min && mdl_dist <=max) {
@@ -75,18 +75,18 @@ std::pair<bool,Real> within_tolerance(Real mdl_dist, const ReferenceDistance& re
}
-std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list,
+std::pair<Real, Real> calc_overlap1(const ResidueRDMap& res_distance_list, const ResNum& rnum,
ChainView mdl_chain,
- Real tol, bool only_fixed,
+ std::vector<Real>& tol_list, bool only_fixed,
bool swap,std::vector<std::pair<Real, Real> >& overlap_list, bool log )
{
std::pair<Real, Real> overlap(0.0, 0.0);
-
- ResidueView mdl_res=mdl_chain.FindResidue(res_distance_list[0].GetFirstAtom().GetResNum());
- for (ResidueDistanceList::const_iterator ai=res_distance_list.begin(),
- ae=res_distance_list.end(); ai!=ae; ++ai) {
- UniqueAtomIdentifier first_atom=ai->GetFirstAtom();
- UniqueAtomIdentifier second_atom=ai->GetSecondAtom();
+ ResidueView mdl_res=mdl_chain.FindResidue(rnum);
+ for (ResidueRDMap::const_iterator ai=res_distance_list.begin(); ai!=res_distance_list.end(); ++ai) {
+ UAtomIdentifiers uais = (*ai).first;
+ std::pair <Real,Real> values = (*ai).second;
+ UniqueAtomIdentifier first_atom=uais.first;
+ UniqueAtomIdentifier second_atom=uais.second;
String name=swap ? swapped_name(first_atom.GetAtomName()) : first_atom.GetAtomName();
AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
@@ -94,58 +94,42 @@ std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list
if (swappable(second_atom.GetResidueName(), second_atom.GetAtomName())) {
continue;
}
- AtomView av2=mdl_chain.FindAtom(second_atom.GetResNum(),second_atom.GetAtomName());
- overlap.second+=1.0;
- if (!(av1 && av2)) {
- continue;
- }
- Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
- ReferenceDistance ref_dist=*ai;
- if (within_tolerance(mdl_dist,ref_dist,tol).first) {
+ }
+ AtomView av2=mdl_chain.FindAtom(second_atom.GetResNum(),second_atom.GetAtomName());
+ overlap.second+=tol_list.size();
+ if (av1) {
+ overlap_list[av1.GetResidue().GetIndex()].second+=1.0;
+ }
+ if (av2) {
+ overlap_list[av2.GetResidue().GetIndex()].second+=1.0;
+ }
+ if (!(av1 && av2)) {
+ continue;
+ }
+ Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+ for (std::vector<Real>::const_reverse_iterator tol_list_it=tol_list.rbegin();tol_list_it!=tol_list.rend();++tol_list_it) {
+ Real tol = * tol_list_it;
+ if (within_tolerance(mdl_dist,values,tol).first) {
if (log) {
LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
<< " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << ref_dist.GetMinDistance() << " " << ref_dist.GetMaxDistance() << " " << tol << " " << "PASS")
+ << mdl_dist << " " << values.first << " " << values.second << " " << tol << " " << "PASS")
}
overlap.first+=1;
+ overlap_list[av1.GetResidue().GetIndex()].first+=1;
+ overlap_list[av2.GetResidue().GetIndex()].first+=1;
} else {
if (log) {
LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
<< " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "FAIL")
+ << mdl_dist << " " << values.first << " " << values.second << " " << tol << " " << "FAIL")
}
+ break;
}
- continue;
- }
- AtomView av2=mdl_chain.FindAtom(second_atom.GetResNum(), second_atom.GetAtomName());
- overlap.second+=1.0;
- if (av1) {
- overlap_list[av1.GetResidue().GetIndex()].second+=1.0;
- }
- if (av2) {
- overlap_list[av2.GetResidue().GetIndex()].second+=1.0;
- }
- if (!(av1 && av2)) {
- continue;
- }
-
- Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
- ReferenceDistance ref_dist=*ai;
- if (within_tolerance(mdl_dist,ref_dist,tol).first) {
- LOG_INFO("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
- << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "PASS")
- overlap.first+=1;
- overlap_list[av1.GetResidue().GetIndex()].first+=1.0;
- overlap_list[av2.GetResidue().GetIndex()].first+=1.0;
- } else {
- LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
- << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "FAIL")
- }
- }
+ }
+ }
return overlap;
-}
+}
std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
@@ -232,8 +216,145 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
return overlap;
}
+void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, std::vector<Real> cutoff_list, std::vector<std::pair<Real, Real> > overlap_list)
+{
+ ChainView mdl_chain=mdl.GetChainList()[0];
+ XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
+ for (GlobalRDMap::const_iterator i=glob_dist_list.begin(); i!=glob_dist_list.end(); ++i) {
+ ResidueRDMap rdl = (*i).second;
+ ResNum rnum = (*i).first;
+ if (rdl.size()==0) {
+ continue;
+ }
+ ResidueView mdl_res=mdl_chain.FindResidue(rnum);
+ if (!mdl_res) {
+ continue;
+ }
+ String rname = mdl_res.GetName();
+ if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+ continue;
+ }
+ std::pair<Real, Real> ov1=calc_overlap1(rdl, rnum,mdl_chain,
+ cutoff_list, true,
+ false, overlap_list,false);
+
+ std::pair<Real, Real> ov2=calc_overlap1(rdl, rnum, mdl_chain,
+ cutoff_list, true,
+ true, overlap_list,false);
+
+ if (ov1.first/ov1.second<ov2.first/ov2.second) {
+ AtomViewList atoms=mdl_res.GetAtomList();
+ for (AtomViewList::iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
+ if (swappable(rname, j->GetName())) {
+ edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
+ }
+ }
+ }
+ }
+}
+
+void update_existence_map (ExistenceMap& ex_map, const EntityView& ev, int ref_counter)
+{
+ AtomViewList ev_atom=ev.GetAtomList();
+ for (AtomViewList::iterator ev_atom_it=ev_atom.begin(); ev_atom_it!=ev_atom.end();++ev_atom_it) {
+ AtomView ev_atom = (*ev_atom_it);
+ UniqueAtomIdentifier uai (ev_atom.GetResidue().GetChain().GetName(),ev_atom.GetResidue().GetNumber(),ev_atom.GetResidue().GetName(),ev_atom.GetName());
+ ExistenceMap::const_iterator find_uai_ci = ex_map.find(uai);
+ int uai_value = 0;
+ if (find_uai_ci!=ex_map.end()) {
+ uai_value=(*find_uai_ci).second;
+ }
+ uai_value |= 1 << (ref_counter);
+ ex_map[uai]=uai_value;
+ }
+}
+
+int in_existence_map(const ExistenceMap& ex_map, const UniqueAtomIdentifier& uai)
+{
+ ExistenceMap::const_iterator find_uai_ci = ex_map.find(uai);
+ int return_value = 0;
+ if (find_uai_ci!=ex_map.end()) {
+ return_value=(*find_uai_ci).second;
+ }
+ return return_value;
+}
+
+void merge_distance_lists(GlobalRDMap& ref_dist_map, const GlobalRDMap& new_dist_map, ExistenceMap& ex_map, const EntityView& ref,int ref_counter)
+{
+ // iterate over the residues in the ref_dist_map
+ for (GlobalRDMap::iterator ref_dist_map_it=ref_dist_map.begin();ref_dist_map_it!=ref_dist_map.end();++ref_dist_map_it) {
+ ResidueRDMap ref_res_map = (*ref_dist_map_it).second;
+ ResNum ref_resnum = (*ref_dist_map_it).first;
+ GlobalRDMap::const_iterator find_new_res_ci = new_dist_map.find(ref_resnum);
+ //if the residue is found in new_dist_map,
+ if (find_new_res_ci != new_dist_map.end()) {
+ ResidueRDMap found_res_new = (*find_new_res_ci).second;
+ //iterate over the the reference distances in the ResidueDistanceMap
+ for (ResidueRDMap::iterator ref_res_map_it = ref_res_map.begin(); ref_res_map_it!=ref_res_map.end();++ref_res_map_it) {
+ UAtomIdentifiers ref_rd = (*ref_res_map_it).first;
+ std::pair<Real,Real> ref_minmax = (*ref_res_map_it).second;
+ ResidueRDMap::const_iterator find_new_rd_ci = found_res_new.find(ref_rd);
+ // if you find the distance in the residue new, udate min and max
+ if (find_new_rd_ci != found_res_new.end()) {
+ std::pair <Real,Real> new_minmax = (*find_new_rd_ci).second;
+ Real min = ref_minmax.first;
+ Real max = ref_minmax.second;
+ if (new_minmax.first < min) min = new_minmax.first;
+ if (new_minmax.second > max) max = new_minmax.second;
+ ref_res_map[ref_rd] = std::make_pair<Real,Real>(min,max);
+ } else {
+ // if you don't find it in the residue new, check that it is not missing because it is too long
+ UniqueAtomIdentifier first_atom_to_find = ref_rd.first;
+ UniqueAtomIdentifier second_atom_to_find = ref_rd.second;
+ // if both atoms are there, remove the distance from the ref_dist_map,
+ if ((ref.FindAtom(first_atom_to_find.GetChainName(),first_atom_to_find.GetResNum(),first_atom_to_find.GetAtomName()).IsValid() &&
+ ref.FindAtom(second_atom_to_find.GetChainName(),second_atom_to_find.GetResNum(),second_atom_to_find.GetAtomName()).IsValid()) ) {
+ ref_res_map.erase(ref_res_map_it);
+ }
+ }
+ }
+ // now iterate over the the new reference distances in residue new
+ for (ResidueRDMap::const_iterator new_res_map_it = found_res_new.begin(); new_res_map_it!=found_res_new.end();++new_res_map_it) {
+ UAtomIdentifiers new_rd = (*new_res_map_it).first;
+ std::pair<Real,Real> new_minmax = (*new_res_map_it).second;
+ ResidueRDMap::const_iterator find_ref_rd_ci = ref_res_map.find(new_rd);
+ // if the distance is found in the residue ref,
+ // it has been taken care of before. If not
+ if (find_ref_rd_ci==ref_res_map.end()) {
+ UniqueAtomIdentifier first_atom_to_find = new_rd.first;
+ UniqueAtomIdentifier second_atom_to_find = new_rd.second;
+ // check that there isn't a structure already processed where both atoms are in
+ // if there is none, add the distance to the residue ref map
+ if (!(in_existence_map(ex_map,first_atom_to_find) & in_existence_map(ex_map,second_atom_to_find))) {
+ ref_res_map[new_rd]= new_minmax;
+ }
+ }
+ }
+ }
+ // if the residue was not found in the new list, it means that it is
+ // absent in the new structure, no new information
+ }
+
+ // now iterate over the residues in the new_list
+ for (GlobalRDMap::const_iterator new_dist_map_it=new_dist_map.begin();new_dist_map_it!=new_dist_map.end();++new_dist_map_it) {
+ ResidueRDMap new_res_map = (*new_dist_map_it).second;
+ ResNum new_resnum = (*new_dist_map_it).first;
+ GlobalRDMap::const_iterator find_ref_res_ci = ref_dist_map.find(new_resnum);
+ // if the residue is found in new_dist_map, it has been taken care before,
+ // if not, add it to the res_dist_map:
+ if (find_ref_res_ci == ref_dist_map.end()) {
+ ResidueRDMap found_res_ref = (*find_ref_res_ci).second;
+ ref_dist_map[new_resnum] = found_res_ref;
+ }
+ }
+
+ // finally, update the existence map
+ update_existence_map (ex_map,ref,ref_counter);
+}
+
}
+
bool IsStandardResidue(String rn)
{
String upper_rn=rn;
@@ -263,6 +384,7 @@ bool IsStandardResidue(String rn)
return false;
}
+
bool UniqueAtomIdentifier::operator==(const UniqueAtomIdentifier& rhs) const
{
if (chain_ == rhs.GetChainName() &&
@@ -274,51 +396,47 @@ bool UniqueAtomIdentifier::operator==(const UniqueAtomIdentifier& rhs) const
return false;
}
-bool ReferenceDistance::IsValid() const
-{
- if (mind_ == -1.0 && maxd_ == -1.0) {
- return false;
- }
- return true;
-}
-
-void ReferenceDistance::Print() const
-{
- if (this->IsValid() == true) {
- std::cout << first_atom_.GetChainName() << " " << first_atom_.GetResNum() << " " << first_atom_.GetResidueName() << " " << first_atom_.GetAtomName() << " " <<
- second_atom_.GetChainName() << " " << second_atom_.GetResNum() << " " << second_atom_.GetResidueName() << " " << second_atom_.GetAtomName() << " " <<
- mind_ << " " << maxd_ << std::endl;
- } else {
- std::cout << first_atom_.GetChainName() << " " << first_atom_.GetResNum() << " " << first_atom_.GetResidueName() << " Placeholder" << std::endl;
- }
-}
-
-bool ReferenceDistance::operator==(const ReferenceDistance& rhs) const
+bool UniqueAtomIdentifier::operator<(const UniqueAtomIdentifier& rhs) const
{
- if (first_atom_ == rhs.GetFirstAtom() &&
- second_atom_ == rhs.GetSecondAtom() &&
- mind_ == rhs.GetMinDistance() &&
- maxd_ == rhs.GetMaxDistance() ) {
+ if (chain_ < rhs.GetChainName()) {
return true;
+ } else if (chain_ > rhs.GetChainName()) {
+ return false;
+ } else {
+ if (residue_ < rhs.GetResNum()) {
+ return true;
+ } else if (residue_ > rhs.GetResNum()) {
+ return false;
+ } else {
+ if (residue_name_ < rhs.GetResidueName()) {
+ return true;
+ } else if (residue_name_ > rhs.GetResidueName()) {
+ return false;
+ } else {
+ if (atom_ < rhs.GetAtomName()) {
+ return true;
+ } else if (atom_ > rhs.GetAtomName()) {
+ return false;
+ } else {
+ return false;
+ }
+ }
+ }
}
- return false;
-}
-
-
+}
-
-GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist)
+GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist)
{
- GlobalDistanceList dist_list;
- ResidueViewList ref_residues=ref.GetChainList()[0].GetResidueList();
+ GlobalRDMap dist_list;
+ ResidueViewList ref_residues=ref.GetChainList()[0].GetResidueList();
for (ResidueViewList::iterator i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
ResidueView rview = (*i);
if (IsStandardResidue(rview.GetName())) {
- ResidueDistanceList res_dist_list;
+ ResidueRDMap res_dist_list;
+ ResNum rnum = rview.GetNumber();
AtomViewList ref_atoms=(*i).GetAtomList();
AtomViewList within;
if (max_dist<0){
- dist_list.push_back(res_dist_list);
within=ref.GetAtomList();
}
for (AtomViewList::iterator ai=ref_atoms.begin(), ae=ref_atoms.end(); ai!=ae; ++ai) {
@@ -335,25 +453,60 @@ GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist)
}
if (aj->GetResidue().GetNumber()>i->GetNumber()) {
Real dist=geom::Length(ai->GetPos()-aj->GetPos());
- ReferenceDistance ref_dist(first_atom,second_atom,dist,dist);
- res_dist_list.push_back(ref_dist);
+ UAtomIdentifiers atoms = std::make_pair<UniqueAtomIdentifier,UniqueAtomIdentifier>(first_atom,second_atom);
+ std::pair<Real,Real> values = std::make_pair<Real,Real>(dist,dist);
+ res_dist_list[atoms]=values;
}
}
- }
- if (res_dist_list.size()==0) {
- UniqueAtomIdentifier current_residue_fake_atom(rview.GetChain().GetName(),rview.GetNumber(),rview.GetName(),"CA");
- ReferenceDistance fake_ref_distance(current_residue_fake_atom,current_residue_fake_atom,-1.0,-1.0);
- res_dist_list.push_back(fake_ref_distance);
- }
- dist_list.push_back(res_dist_list);
+ }
+ dist_list[rnum]=res_dist_list;
}
- }
+ }
return dist_list;
}
+GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list, std::vector<Real> cutoff_list, Real max_dist)
+{
+ int ref_counter=0;
+ ExistenceMap ex_map;
+ GlobalRDMap glob_dist_list = CreateDistanceList(ref_list[0],max_dist);
+ update_existence_map (ex_map,ref_list[0],ref_counter);
+ ref_counter++;
+ for (std::vector<EntityView>::const_iterator ref_list_it=ref_list.begin()+1;ref_list_it!=ref_list.end();++ref_list_it) {
+ EntityView ref = *ref_list_it;
+ std::vector<std::pair<Real, Real> > overlap_list(ref.GetResidueCount(), std::pair<Real, Real>(0.0, 0.0));
+ check_and_swap(glob_dist_list,ref,cutoff_list,overlap_list);
+ GlobalRDMap new_dist_list=CreateDistanceList(ref,max_dist);
+ merge_distance_lists(glob_dist_list,new_dist_list,ex_map,ref,ref_counter);
+ ref_counter++;
+ }
+ return glob_dist_list;
+}
-Real LocalDistTest(const EntityView& mdl, const GlobalDistanceList& glob_dist_list,
- Real cutoff, const String& local_ldt_property_string)
+void PrintResidueRDMap(const ResidueRDMap& res_dist_list)
+{
+ for (ResidueRDMap::const_iterator res_dist_list_it = res_dist_list.begin();res_dist_list_it!=res_dist_list.end();++res_dist_list_it) {
+ UAtomIdentifiers uais = (*res_dist_list_it).first;
+ std::pair<Real,Real> minmax = (*res_dist_list_it).second;
+ std::cout << uais.first.GetChainName() << " " << uais.first.GetResNum() << " " << uais.first.GetResidueName() << " " << uais.first.GetAtomName() << " " <<
+ uais.second.GetChainName() << " " << uais.second.GetResNum() << " " << uais.second.GetResidueName() << " " << uais.second.GetAtomName() << " " <<
+ minmax.first << " " << minmax.second << std::endl;
+ }
+}
+
+
+void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list){
+ for (GlobalRDMap::const_iterator glob_dist_list_it = glob_dist_list.begin();glob_dist_list_it!=glob_dist_list.end();++glob_dist_list_it) {
+ if ((*glob_dist_list_it).second.size()!=0) {
+ PrintResidueRDMap((*glob_dist_list_it).second);
+ }
+ }
+}
+
+
+
+Real LocalDistTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list,
+ std::vector<Real> cutoff_list, const String& local_ldt_property_string)
{
if (!mdl.GetResidueCount()) {
LOG_WARNING("model structures doesn't contain any residues");
@@ -364,60 +517,25 @@ Real LocalDistTest(const EntityView& mdl, const GlobalDistanceList& glob_dist_li
return 0.0;
}
std::vector<std::pair<Real, Real> > overlap_list(mdl.GetResidueCount(), std::pair<Real, Real>(0.0, 0.0));
+ check_and_swap(glob_dist_list,mdl,cutoff_list,overlap_list);
ChainView mdl_chain=mdl.GetChainList()[0];
- // Residues with symmetric side-chains require special treatment as there are
- // two possible ways to name the atoms. We calculate the overlap with the
- // fixed atoms and take the solution that gives bigger scores.
- XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
- for (GlobalDistanceList::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
- ResidueDistanceList rdl = *i;
- if (!rdl[0].IsValid()) {
- continue;
- }
-
- ResNum rnum = rdl[0].GetFirstAtom().GetResNum() ;
- String rname=rdl[0].GetFirstAtom().GetResidueName();
- ResidueView mdl_res=mdl_chain.FindResidue(rnum);
- if (!mdl_res) {
- continue;
- }
- if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || rname=="PHE" || rname=="LYS" || rname=="ARG")) {
- continue;
- }
- std::pair<Real, Real> ov1=calc_overlap1(*i, mdl_chain,
- cutoff, true,
- false, overlap_list,false);
-
- std::pair<Real, Real> ov2=calc_overlap1(*i, mdl_chain,
- cutoff, true,
- true, overlap_list,false);
- if (ov1.first/ov1.second<ov2.first/ov2.second) {
- AtomViewList atoms=mdl_res.GetAtomList();
- for (AtomViewList::iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
- if (swappable(rname, j->GetName())) {
- edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
- }
- }
- }
- }
overlap_list.clear();
std::pair<Real, Real> total_ov(0.0, 0.0);
- for (GlobalDistanceList::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
- ResidueDistanceList rdl = *i;
- if (rdl[0].IsValid()) {
- std::pair<Real, Real> ov1=calc_overlap1(*i, mdl_chain, cutoff,
+ for (GlobalRDMap::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
+ ResidueRDMap rdl = (*i).second;
+ ResNum rn = (*i).first;
+ if (rdl.size()!=0) {
+ std::pair<Real, Real> ov1=calc_overlap1(rdl, rn, mdl_chain, cutoff_list,
false, false, overlap_list,true);
total_ov.first+=ov1.first;
total_ov.second+=ov1.second;
}
if(local_ldt_property_string!="") {
- ResNum rn = rdl[0].GetFirstAtom().GetResNum();
ResidueView mdlr=mdl_chain.FindResidue(rn);
if (mdlr.IsValid()) {
int mdl_res_index =mdlr.GetIndex();
Real local_ldt=overlap_list[mdl_res_index].first/(overlap_list[mdl_res_index].second ? overlap_list[mdl_res_index].second : 1);
- ResidueView res_to_wr = mdl_chain.FindResidue((*i)[0].GetFirstAtom().GetResNum());
- res_to_wr.SetFloatProp(local_ldt_property_string, local_ldt);
+ mdlr.SetFloatProp(local_ldt_property_string, local_ldt);
}
}
}
@@ -427,8 +545,10 @@ Real LocalDistTest(const EntityView& mdl, const GlobalDistanceList& glob_dist_li
Real LocalDistTest(const EntityView& mdl, const EntityView& target, Real cutoff, Real max_dist, const String& local_ldt_property_string)
{
- GlobalDistanceList glob_dist_list = CreateDistanceList(target,max_dist);
- return LocalDistTest(mdl, glob_dist_list, cutoff, local_ldt_property_string);
+ std::vector<Real> cutoffs;
+ cutoffs.push_back(cutoff);
+ GlobalRDMap glob_dist_list = CreateDistanceList(target,max_dist);
+ return LocalDistTest(mdl, glob_dist_list, cutoffs, local_ldt_property_string);
}
@@ -484,16 +604,41 @@ Real LocalDistTest(const ost::seq::AlignmentHandle& aln,
return 0.0;
}
-Real LDTHA(EntityView& v, const GlobalDistanceList& global_dist_list)
+Real LDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
{
-
- Real cutoffs[]={0.5,1,2,4};
- String labels[]={"localldt0.5","localldt1","localldt2","ldtlocal4"};
- Real ldt=0.0;
- for (int n=0; n<4; ++n) {
- ldt+=alg::LocalDistTest(v, global_dist_list, cutoffs[n], labels[n]);
- }
- ldt/=4.0;
+ std::vector<Real> cutoffs;
+ cutoffs.push_back(0.5);
+ cutoffs.push_back(1.0);
+ cutoffs.push_back(2.0);
+ cutoffs.push_back(4.0);
+ String label="localldt";
+ Real ldt=alg::LocalDistTest(v, global_dist_list, cutoffs, label);
+ return ldt;
+}
+
+Real OldStyleLDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
+{
+ Real ldt =0;
+ std::vector<Real> cutoffs05;
+ cutoffs05.push_back(0.5);
+ Real ldt05=alg::LocalDistTest(v, global_dist_list, cutoffs05, "localldt0.5");
+ std::vector<Real> cutoffs1;
+ cutoffs1.push_back(1.0);
+ Real ldt1=alg::LocalDistTest(v, global_dist_list, cutoffs1, "localldt1");
+ std::vector<Real> cutoffs2;
+ cutoffs2.push_back(2.0);
+ Real ldt2=alg::LocalDistTest(v, global_dist_list, cutoffs2, "localldt2");
+ std::vector<Real> cutoffs4;
+ cutoffs4.push_back(4.0);
+ Real ldt4=alg::LocalDistTest(v, global_dist_list, cutoffs4, "localldt4");
+ ldt = (ldt05+ldt1+ldt2+ldt4)/4.0;
+ for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
+ ResidueView ritv = *rit;
+ if (ritv.HasProp("localldt0.5")) {
+ Real ldt_local = (ritv.GetFloatProp("localldt0.5")+ritv.GetFloatProp("localldt1")+ritv.GetFloatProp("localldt2")+ritv.GetFloatProp("localldt4"))/4.0;
+ ritv.SetFloatProp("localldt",ldt_local);
+ }
+ }
return ldt;
}
diff --git a/modules/mol/alg/src/local_dist_test.hh b/modules/mol/alg/src/local_dist_test.hh
index 9c174f9698bb37ccd7e23f16febf77e1928624c4..8a26e1c7dbf8ef67f50a8fcfce326cd45abc16ce 100644
--- a/modules/mol/alg/src/local_dist_test.hh
+++ b/modules/mol/alg/src/local_dist_test.hh
@@ -31,12 +31,16 @@ class UniqueAtomIdentifier
public:
UniqueAtomIdentifier(const String& chain,const ResNum& residue,const String& residue_name, const String& atom): chain_(chain),residue_(residue),residue_name_(residue_name),atom_(atom) {}
-
+
+ // to make the compiler happy (boost python map suite)
+ UniqueAtomIdentifier(): chain_(""),residue_(ResNum(1)),residue_name_(""),atom_("") {}
+
String GetChainName() const { return chain_; }
ResNum GetResNum() const { return residue_; }
String GetResidueName() const { return residue_name_; }
String GetAtomName() const { return atom_; }
bool operator==(const UniqueAtomIdentifier& rhs) const;
+ bool operator<(const UniqueAtomIdentifier& rhs) const;
private:
@@ -46,40 +50,19 @@ private:
String atom_;
};
-class ReferenceDistance
-{
-
-public:
-
- ReferenceDistance(const UniqueAtomIdentifier& first_atom, const UniqueAtomIdentifier& second_atom, Real mind, Real maxd):
- first_atom_(first_atom),second_atom_(second_atom),mind_(mind),maxd_(maxd) {}
- UniqueAtomIdentifier GetFirstAtom() const {return first_atom_;}
- UniqueAtomIdentifier GetSecondAtom() const {return second_atom_;}
- Real GetMinDistance() const {return mind_ ;}
- Real GetMaxDistance() const {return maxd_ ;}
- bool IsValid() const;
- void Print() const;
- bool operator==(const ReferenceDistance& rhs) const;
-
-private:
-
- UniqueAtomIdentifier first_atom_;
- UniqueAtomIdentifier second_atom_;
- Real mind_;
- Real maxd_;
-};
+typedef std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier> UAtomIdentifiers;
+typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<float,float> > ResidueRDMap;
+typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
+typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
-typedef std::vector<ReferenceDistance> ResidueDistanceList;
-typedef std::vector<ResidueDistanceList> GlobalDistanceList;
-
Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
- const GlobalDistanceList& dist_list,
- Real cutoff,
+ const GlobalRDMap& dist_list,
+ std::vector<Real> cutoff_list,
const String& local_ldt_property_string="");
Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
const EntityView& target,
- Real cutoff,
+ Real cutoff_list,
Real max_dist,
const String& local_ldt_property_string="");
@@ -88,12 +71,15 @@ Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
int ref_index=0, int mdl_index=1);
-Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView& v, const GlobalDistanceList& global_dist_list);
-
-GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist);
-bool IsStandardResidue(String rn);
+Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView& v, const GlobalRDMap& global_dist_list);
+Real DLLEXPORT_OST_MOL_ALG OldStyleLDTHA(EntityView& v, const GlobalRDMap& global_dist_list);
+GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist);
+GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list,Real cutoff, Real max_dist);
+void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list);
+void PrintResidueRDMap(const ResidueRDMap& res_dist_list);
+bool IsStandardResidue(String rn);
}}}
diff --git a/modules/mol/base/src/residue_prop.hh b/modules/mol/base/src/residue_prop.hh
index e86eef5f11b16ba4eed55153e0cb38901faedd1e..4a411cca6e27250fdd23d4b538ffea6064b8d106 100644
--- a/modules/mol/base/src/residue_prop.hh
+++ b/modules/mol/base/src/residue_prop.hh
@@ -34,6 +34,12 @@ class DLLEXPORT ResNum: private
boost::unit_steppable<ResNum> > > > >
{
public:
+
+ // needed to wrap certain map classes
+ ResNum():
+ num_(1),alt_('\0')
+ {}
+
ResNum(int n):
num_(n), alt_('\0')
{ }