From 0d7ea45a848484cbf7fb95effb6b83df8d5b9271 Mon Sep 17 00:00:00 2001
From: Valerio Mariani <valerio.mariani@unibas.ch>
Date: Fri, 24 Feb 2012 16:24:17 +0100
Subject: [PATCH] Refactored and wrapped ldt-related fuctions
---
modules/mol/alg/pymod/wrap_mol_alg.cc | 61 ++++++-
modules/mol/alg/src/filter_clashes.cc | 174 ++++++++++++++++++-
modules/mol/alg/src/filter_clashes.hh | 38 +++--
modules/mol/alg/src/ldt.cc | 228 +++++--------------------
modules/mol/alg/src/local_dist_test.cc | 14 +-
modules/mol/alg/src/local_dist_test.hh | 4 +
6 files changed, 303 insertions(+), 216 deletions(-)
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 55de61740..52b7a703f 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -38,8 +38,35 @@ Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real, co
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;
mol::EntityView (*fc_b)(const mol::EntityHandle&, const mol::alg::ClashingDistances&, bool)=&mol::alg::FilterClashes;
+mol::EntityView (*csc_a)(const mol::EntityView&, const mol::alg::StereoChemicalParams&, const mol::alg::StereoChemicalParams&, Real, Real, bool)=&mol::alg::CheckStereoChemistry;
+mol::EntityView (*csc_b)(const mol::EntityHandle&, const mol::alg::StereoChemicalParams&, const mol::alg::StereoChemicalParams&, Real, Real, bool)=&mol::alg::CheckStereoChemistry;
mol::CoordGroupHandle (*superpose_frames1)(mol::CoordGroupHandle&, mol::EntityView&, int, int, int)=&mol::alg::SuperposeFrames;
mol::CoordGroupHandle (*superpose_frames2)(mol::CoordGroupHandle&, mol::EntityView&, mol::EntityView&, int, int)=&mol::alg::SuperposeFrames;
+
+ost::mol::alg::StereoChemicalParams fill_stereochemical_params_wrapper (const String& header, const list& stereo_chemical_props_file)
+{
+ int stereo_chemical_props_file_length = boost::python::extract<int>(stereo_chemical_props_file.attr("__len__")());
+ std::vector<String> stereo_chemical_props_file_vector;
+
+ for (int i=0; i<stereo_chemical_props_file_length; i++) {
+ String extracted_string = boost::python::extract<String>(stereo_chemical_props_file[i]);
+ stereo_chemical_props_file_vector.push_back(extracted_string);
+ }
+ return ost::mol::alg::FillStereoChemicalParams(header,stereo_chemical_props_file_vector);
+}
+
+ost::mol::alg::ClashingDistances fill_clashing_distances_wrapper (const list& stereo_chemical_props_file, Real min_default_distance, Real min_distance_tolerance)
+{
+ int stereo_chemical_props_file_length = boost::python::extract<int>(stereo_chemical_props_file.attr("__len__")());
+ std::vector<String> stereo_chemical_props_file_vector(stereo_chemical_props_file_length);
+
+ for (int i=0; i<stereo_chemical_props_file_length; i++) {
+ stereo_chemical_props_file_vector[i] = boost::python::extract<char const*>(stereo_chemical_props_file[i]);
+ }
+
+ return ost::mol::alg::FillClashingDistances(stereo_chemical_props_file_vector,min_default_distance,min_distance_tolerance);
+}
+
}
@@ -54,12 +81,40 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("LocalDistTest", ldt_a, (arg("local_ldt_property_string")=""));
def("LocalDistTest", ldt_b, (arg("ref_index")=0, arg("mdl_index")=1));
- def("FilterClashes", fc_a, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
- def("FilterClashes", fc_b, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
+ def("FilterClashes", fc_a, (arg("ent"), arg("clashing_distances"), arg("always_remove_bb")=false));
+ def("FilterClashes", fc_b, (arg("ent"), arg("clashing_distances"), arg("always_remove_bb")=false));
+ def("CheckStereoChemistry", csc_a, (arg("ent"), arg("bonds"), arg("angles"), arg("bond_tolerance"), arg("angle_tolerance"), arg("always_remove_bb")=false));
+ 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("SuperposeFrames", superpose_frames1,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0,
arg("end")=-1, arg("ref")=-1));
def("SuperposeFrames", superpose_frames2,
(arg("source"), arg("sel")=ost::mol::EntityView(), arg("ref_view")=ost::mol::EntityView(),arg("begin")=0,
arg("end")=-1));
-}
+
+
+ class_<mol::alg::ClashingDistances> ("ClashingDistances" ,init<Real,Real>())
+ .def("SetClashingDistance",&mol::alg::ClashingDistances::SetClashingDistance)
+ .def("GetClashingDistance",&mol::alg::ClashingDistances::GetClashingDistance)
+ .def("GetMaxAdjustedDistance",&mol::alg::ClashingDistances::GetMaxAdjustedDistance)
+ .def("IsEmpty",&mol::alg::ClashingDistances::IsEmpty)
+
+ .def("PrintAllDistances",&mol::alg::ClashingDistances::PrintAllDistances)
+ ;
+
+ class_<mol::alg::StereoChemicalParams> ("StereoChemicalParams" ,init<>())
+ .def("SetParam",&mol::alg::StereoChemicalParams::SetParam)
+ .def("GetParam",&mol::alg::StereoChemicalParams::GetParam)
+ .def("ContainsParam",&mol::alg::StereoChemicalParams::ContainsParam)
+ .def("IsEmpty",&mol::alg::StereoChemicalParams::IsEmpty)
+
+ .def("PrintAllParameters",&mol::alg::StereoChemicalParams::PrintAllParameters)
+
+ ;
+
+ def("FillClashingDistances",&fill_clashing_distances_wrapper);
+ def("FillStereoChemicalParams",&fill_stereochemical_params_wrapper);
+
+}
\ No newline at end of file
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index 8bb40bf66..e84943581 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -109,6 +109,14 @@ Real ClashingDistances::GetMaxAdjustedDistance() const
return max_adjusted_distance;
}
+bool ClashingDistances::IsEmpty() const
+{
+ if (min_distance_.size()==0) {
+ return true;
+ }
+ return false;
+}
+
void StereoChemicalParams::SetParam(const String& param, const String& residue, Real value, Real st_dev)
{
std::pair<String,String> key = std::make_pair<String,String>(param,residue);
@@ -144,6 +152,158 @@ void StereoChemicalParams::PrintAllParameters() const
}
};
+bool StereoChemicalParams::IsEmpty() const
+{
+ if (params_.size()==0) {
+ return true;
+ }
+ return false;
+}
+
+StereoChemicalParams FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file)
+{
+ StereoChemicalParams table;
+ bool found=false;
+ std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
+ while (line_iter!=stereo_chemical_props_file.end()) {
+ if ((*line_iter).length()!=0 && (*line_iter).length()!=1) {
+ StringRef line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> line_str_vec = line_string_ref.split();
+ if (line_str_vec[0].str()==header) {
+ found=true;
+ line_iter++;
+ while ((*line_iter)[0]!='-') {
+ if ((*line_iter)[0]!='#') {
+ StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
+ if (second_line_str_vec.size()!=4) {
+ std::cout << "The number of elements in one of the lines is wrong" << std::endl;
+ return StereoChemicalParams();
+ }
+ StringRef item = second_line_str_vec[0];
+ String res = second_line_str_vec[1].str();
+ std::pair<bool,float> parse_value = second_line_str_vec[2].to_float();
+ std::pair<bool,float> parse_stddev = second_line_str_vec[3].to_float();
+ Real value,stddev;
+ if (parse_value.first==true) {
+ value=static_cast<Real>(parse_value.second);
+ } else {
+ std::cout << "One of the values in the third column is not a number" << std::endl;
+ return StereoChemicalParams();
+ };
+ if (parse_stddev.first==true) {
+ stddev=static_cast<Real>(parse_stddev.second);
+ } else {
+ std::cout << "One of the values in the fourth column is not a number" << std::endl;
+ return StereoChemicalParams();
+ };
+ std::vector<StringRef> split_item = item.split('-');
+ String rearranged_item;
+ if (split_item.size() == 2) {
+ String atom1 = split_item[0].str();
+ String atom2 = split_item[1].str();
+ if (atom2 < atom1) {
+ std::stringstream srearr;
+ srearr << atom2 << "-" << atom1;
+ rearranged_item=srearr.str();
+ } else {
+ rearranged_item = item.str();
+ }
+ } else if (split_item.size() == 3) {
+ String atom1 = split_item[0].str();
+ String atom = split_item[1].str();
+ String atom2 = split_item[2].str();
+ if (atom2 < atom1) {
+ std::stringstream srearr;
+ srearr << atom2 << "-" << atom << "-" << atom1;
+ rearranged_item=srearr.str();
+ } else {
+ rearranged_item = item.str();
+ }
+ } else {
+ std::cout << "One of the strings describing the parameter has the wrong format" << std::endl;
+ return StereoChemicalParams();
+ }
+ table.SetParam(rearranged_item,res,value,stddev);
+ line_iter++;
+ }
+ }
+ }
+ }
+ line_iter++;
+ }
+ if (found==false) {
+ std::cout << "Could not find the relevant section in the stereo-chemical parameter file" << std::endl;
+ return StereoChemicalParams();
+ };
+ return table;
+};
+
+ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_props_file, Real min_default_distance, Real min_distance_tolerance)
+{
+ ClashingDistances table(min_default_distance,min_distance_tolerance);
+ bool found=false;
+ std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
+ while (line_iter!=stereo_chemical_props_file.end()) {
+ if ((*line_iter).length()!=0 && (*line_iter).length()!=1) {
+ StringRef line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> line_str_vec = line_string_ref.split();
+ if (line_str_vec[0].str()=="Non-bonded") {
+ found=true;
+ line_iter++;
+ while ((*line_iter)[0]!='-') {
+ if ((*line_iter)[0]!='#') {
+ StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
+ if (second_line_str_vec.size()!=3) {
+ std::cout << "The number of elements in one of the lines is wrong" << std::endl;
+ return ClashingDistances(min_default_distance,min_distance_tolerance);
+ }
+ String item = second_line_str_vec[0].str();
+
+ std::pair<bool,float> parse_value = second_line_str_vec[1].to_float();
+ std::pair<bool,float> parse_stddev = second_line_str_vec[2].to_float();
+ Real value,stddev;
+ if (parse_value.first==true) {
+ value=static_cast<Real>(parse_value.second);
+ } else {
+ std::cout << "One of the distance values is not a number" << std::endl;
+ return ClashingDistances(min_default_distance,min_distance_tolerance);
+ };
+ if (parse_stddev.first==true) {
+ stddev=static_cast<Real>(parse_stddev.second);
+ } else {
+ std::cout << "One of the tolerance values is not a number" << std::endl;
+ return ClashingDistances(min_default_distance,min_distance_tolerance);
+ }
+ StringRef itemsr(item.data(),item.length());
+ std::vector<StringRef> eles = itemsr.split('-');
+ if (itemsr.size() != 3) {
+ std::cout << "One of the strings describing the interacting atoms has the wrong format" << std::endl;
+ return ClashingDistances(min_default_distance,min_distance_tolerance);
+ }
+ String ele1=eles[0].str();
+ String ele2=eles[1].str();
+ if (ele2 < ele1) {
+ table.SetClashingDistance(ele2,ele1,value,stddev);
+ } else {
+ table.SetClashingDistance(ele1,ele2,value,stddev);
+ }
+ line_iter++;
+ }
+ }
+ }
+ }
+ line_iter++;
+ }
+ if (found==false) {
+ std::cout << "Could not find the relevant section in the stereo-chemical parameter file" << std::endl;
+ return ClashingDistances(min_default_distance,min_distance_tolerance);
+ }
+ return table;
+}
+
+
EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb)
{
LOG_INFO("Checking stereo-chemistry")
@@ -187,7 +347,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
String bond_str = bond_string(atom,other_atom);
std::pair<Real,Real> length_stddev = bond_table.GetParam(bond_str,residue_str);
Real ref_length = length_stddev.first;
- Real ref_stddev = length_stddev.second;
+ Real ref_stddev = length_stddev.second;
Real min_length = ref_length - bond_tolerance*ref_stddev;
Real max_length = ref_length + bond_tolerance*ref_stddev;
Real zscore = (blength - ref_length)/ref_stddev;
@@ -219,7 +379,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
continue;
}
if (atom1.GetResidue()!=res.GetHandle()) {
- continue;
+ continue;
}
for (BondHandeList::const_iterator bond_iter2=bonds.begin(); bond_iter2!=bonds.end(); ++bond_iter2) {
BondHandle bond2=*bond_iter2;
@@ -232,7 +392,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
continue;
}
if (atom2.GetResidue()!=res.GetHandle()) {
- continue;
+ continue;
}
if (atom1.GetHashCode() > atom2.GetHashCode()) {
Real awidth;
@@ -245,7 +405,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
String angle_str = angle_string(atom1,atom,atom2);
std::pair<Real,Real> width_stddev = angle_table.GetParam(angle_str,residue_str);
Real ref_width = width_stddev.first;
- Real ref_stddev = width_stddev.second;
+ Real ref_stddev = width_stddev.second;
Real min_width = ref_width - angle_tolerance*ref_stddev;
Real max_width = ref_width + angle_tolerance*ref_stddev;
Real zscore = (awidth - ref_width)/ref_stddev;
@@ -344,9 +504,9 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
Real d=geom::Length2(atom.GetPos()-atom2.GetPos());
std::pair <Real,Real> distance_tolerance=min_distances.GetClashingDistance(ele1, ele2);
- Real distance=distance_tolerance.first;
- Real tolerance=distance_tolerance.second;
- Real threshold=distance-tolerance;
+ Real distance=distance_tolerance.first;
+ Real tolerance=distance_tolerance.second;
+ Real threshold=distance-tolerance;
if (d<threshold*threshold) {
LOG_INFO(atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetName() << " " << atom2.GetResidue().GetName() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "FAIL")
diff --git a/modules/mol/alg/src/filter_clashes.hh b/modules/mol/alg/src/filter_clashes.hh
index 0e41cbe17..5d64fc230 100644
--- a/modules/mol/alg/src/filter_clashes.hh
+++ b/modules/mol/alg/src/filter_clashes.hh
@@ -28,10 +28,11 @@ class ClashingDistances
{
public:
- ClashingDistances(Real default_dist, Real tolerance): default_min_distance_(default_dist), default_min_distance_tolerance_(tolerance) {}
+ ClashingDistances(Real default_dist, Real tolerance): default_min_distance_(default_dist), default_min_distance_tolerance_(tolerance), valid_flag_(true) {}
void SetClashingDistance(const String& ele1,const String& ele2, Real min_distance, Real tolerance);
std::pair<Real,Real> GetClashingDistance(const String& ele1,const String& ele2) const;
- Real GetMaxAdjustedDistance() const;
+ Real GetMaxAdjustedDistance() const;
+ bool IsEmpty() const;
//DEBUG
void PrintAllDistances() const;
@@ -41,6 +42,7 @@ private:
std::map <String,std::pair<float,float> > min_distance_;
Real default_min_distance_;
Real default_min_distance_tolerance_;
+ bool valid_flag_;
};
@@ -51,6 +53,7 @@ public:
void SetParam(const String& param, const String& residue, Real value, Real st_dev);
std::pair<Real,Real> GetParam(const String& element,const String& residue) const;
bool ContainsParam(const String& param,const String& residue) const;
+ bool IsEmpty() const;
//DEBUG
void PrintAllParameters() const;
@@ -58,32 +61,31 @@ public:
private:
std::map<std::pair<String,String>,std::pair<float,float> > params_;
-
-};
-
-
-
-
+
+};
+
+ClashingDistances DLLEXPORT_OST_MOL_ALG FillClashingDistances(std::vector<String>& stereo_chemical_props_file, Real min_default_distance, Real min_distance_tolerance);
+StereoChemicalParams DLLEXPORT_OST_MOL_ALG FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file);
+
EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityView& ent,
const ClashingDistances& min_distances, bool always_remove_bb=false);
EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityHandle& ent,
const ClashingDistances& min_distances, bool always_remove_bb=false);
-
EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityView& ent,
const StereoChemicalParams& bond_table,
- const StereoChemicalParams& angle_table,
- Real bond_tolerance,
- Real angle_tolerance,
- bool always_remove_bb=false);
+ const StereoChemicalParams& angle_table,
+ Real bond_tolerance,
+ Real angle_tolerance,
+ bool always_remove_bb=false);
EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityHandle& ent,
- const StereoChemicalParams& bond_table,
- const StereoChemicalParams& angle_table,
- Real bond_tolerance,
- Real angle_tolerance,
- bool always_remove_bb=false);
+ const StereoChemicalParams& bond_table,
+ const StereoChemicalParams& angle_table,
+ Real bond_tolerance,
+ Real angle_tolerance,
+ bool always_remove_bb=false);
}}}
diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc
index ec2f0fd07..a89e294b5 100644
--- a/modules/mol/alg/src/ldt.cc
+++ b/modules/mol/alg/src/ldt.cc
@@ -40,137 +40,6 @@ using namespace ost::mol;
using namespace ost::mol::alg;
namespace po=boost::program_options;
-
-
-void FillStereoChemicalParams(const String& header, StereoChemicalParams& table, std::vector<String>& stereo_chemical_props_file)
-{
- bool found=false;
- std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
- while (line_iter!=stereo_chemical_props_file.end()) {
- if ((*line_iter).length()!=0) {
- StringRef line_string_ref(line_iter->data(),(*line_iter).length());
- std::vector<StringRef> line_str_vec = line_string_ref.split();
- if (line_str_vec[0].str()==header) {
- found=true;
- line_iter++;
- while ((*line_iter)[0]!='-') {
- if ((*line_iter)[0]!='#') {
- StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
- std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
- if (second_line_str_vec.size()!=4) {
- throw Error("The number of elements in one of the lines is wrong");
- }
- StringRef item = second_line_str_vec[0];
- String res = second_line_str_vec[1].str();
- std::pair<bool,float> parse_value = second_line_str_vec[2].to_float();
- std::pair<bool,float> parse_stddev = second_line_str_vec[3].to_float();
- Real value,stddev;
- if (parse_value.first==true) {
- value=static_cast<Real>(parse_value.second);
- } else {
- throw Error("One of the values in the third column is not a number");
- };
- if (parse_stddev.first==true) {
- stddev=static_cast<Real>(parse_stddev.second);
- } else {
- throw Error("One of the values in the fourth column is not a number");
- };
- std::vector<StringRef> split_item = item.split('-');
- String rearranged_item;
- if (split_item.size() == 2) {
- String atom1 = split_item[0].str();
- String atom2 = split_item[1].str();
- if (atom2 < atom1) {
- std::stringstream srearr;
- srearr << atom2 << "-" << atom1;
- rearranged_item=srearr.str();
- } else {
- rearranged_item = item.str();
- }
- } else if (split_item.size() == 3) {
- String atom1 = split_item[0].str();
- String atom = split_item[1].str();
- String atom2 = split_item[2].str();
- if (atom2 < atom1) {
- std::stringstream srearr;
- srearr << atom2 << "-" << atom << "-" << atom1;
- rearranged_item=srearr.str();
- } else {
- rearranged_item = item.str();
- }
- } else {
- throw Error("One of the strings describing the parameter has the wrong format");
- }
- table.SetParam(rearranged_item,res,value,stddev);
- line_iter++;
- }
- }
- }
- }
- line_iter++;
- }
- if (found==false) {
- throw Error("Could not find the relevant section in the stereo-chemical parameter file");
- };
-};
-
-void FillClashingDistances(ClashingDistances& table, std::vector<String>& stereo_chemical_props_file)
-{
- bool found=false;
- std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
- while (line_iter!=stereo_chemical_props_file.end()) {
- if ((*line_iter).length()!=0) {
- StringRef line_string_ref(line_iter->data(),(*line_iter).length());
- std::vector<StringRef> line_str_vec = line_string_ref.split();
- if (line_str_vec[0].str()=="Non-bonded") {
- found=true;
- line_iter++;
- while ((*line_iter)[0]!='-') {
- if ((*line_iter)[0]!='#') {
- StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
- std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
- if (second_line_str_vec.size()!=3) {
- throw Error("The number of elements in one of the lines is wrong");
- }
- String item = second_line_str_vec[0].str();
-
- std::pair<bool,float> parse_value = second_line_str_vec[1].to_float();
- std::pair<bool,float> parse_stddev = second_line_str_vec[2].to_float();
- Real value,stddev;
- if (parse_value.first==true) {
- value=static_cast<Real>(parse_value.second);
- } else {
- throw Error("One of the distance values is not a number");
- };
- if (parse_stddev.first==true) {
- stddev=static_cast<Real>(parse_stddev.second);
- } else {
- throw Error("One of the tolerance values is not a number");
- }
- StringRef itemsr(item.data(),item.length());
- std::vector<StringRef> eles = itemsr.split('-');
- if (itemsr.size() != 3) {
- throw Error("One of the strings describing the interacting atoms has the wrong format");
- }
- String ele1=eles[0].str();
- String ele2=eles[1].str();
- if (ele2 < ele1) {
- table.SetClashingDistance(ele2,ele1,value,stddev);
- } else {
- table.SetClashingDistance(ele1,ele2,value,stddev);
- }
- line_iter++;
- }
- }
- }
- }
- line_iter++;
- }
- if (found==false) {
- throw Error("Could not find the relevant section in the stereo-chemical parameter file");
- }
-}
-
EntityHandle load(const String& file, const IOProfile& profile)
{
try {
@@ -213,6 +82,7 @@ int main (int argc, char **argv)
Real min_distance_tolerance = 0.0;
Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0;
+ Real radius=8.0;
IOProfile profile;
profile.bond_feasibility_check=false;
@@ -222,16 +92,16 @@ int main (int argc, char **argv)
po::options_description desc("Options");
desc.add_options()
("calpha,c", "consider only calpha atoms")
- ("sel,s", po::value<String>(&sel)->default_value(""), "selection for reference")
+ ("sel,s", po::value<String>(&sel)->default_value(""), "selection performed on reference structure")
("tolerant,t", "fault tolerant mode")
- ("structural-checks,f", "structural checks")
+ ("structural-checks,f", "perform stereo-chemical and clash checks")
("version,e", "version")
- ("parameter-file,p", po::value<String>(), "stereo chemical parameter file")
+ ("parameter-file,p", po::value<String>(), "stereo-chemical parameter file")
("verbosity,v", po::value<int>(), "verbosity level")
- ("tolerance_bonds,b", po::value<Real>(), "tolerance in stddev for bonds")
- ("tolerance_angles,a", po::value<Real>(), "tolerance in stddev for angles")
+ ("bond_tolerance,b", po::value<Real>(), "tolerance in stddev for bonds")
+ ("angle_tolerance,a", po::value<Real>(), "tolerance in stddev for angles")
("default_clash,m", po::value<Real>(), "clashing distance for unknown atom types")
- ("files", po::value< std::vector<String> >(), "input file")
+ ("files", po::value< std::vector<String> >(), "input file(s)")
;
po::positional_options_description p;
p.add("files", -1);
@@ -282,12 +152,11 @@ int main (int argc, char **argv)
}
}
-
- if (vm.count("tolerance_bonds")) {
- bond_tolerance=vm["tolerance_bonds"].as<Real>();
+ if (vm.count("bond_tolerance")) {
+ bond_tolerance=vm["bond_tolerance"].as<Real>();
}
- if (vm.count("tolerance_angles")) {
- angle_tolerance=vm["tolerance_angles"].as<Real>();
+ if (vm.count("angle_tolerance")) {
+ angle_tolerance=vm["angle_tolerance"].as<Real>();
}
if (vm.count("default_clash")) {
min_default_distance=vm["default_clash"].as<Real>();
@@ -300,40 +169,37 @@ int main (int argc, char **argv)
}
files.pop_back();
EntityView ref_view=ref.Select(sel);
- std::cout << "#Parameter filename: " << parameter_filename << std::endl;
- std::cout << "#Verbosity level: " << verbosity_level << std::endl;
+ std::cout << "Verbosity level: " << verbosity_level << std::endl;
if (structural_checks) {
- std::cout << "#Stereo-chemical and steric clash checks: On " << std::endl;
+ std::cout << "Stereo-chemical and steric clash checks: On " << std::endl;
} else {
- std::cout << "#Stereo-chemical and steric clash checks: Off " << std::endl;
+ std::cout << "Stereo-chemical and steric clash checks: Off " << std::endl;
}
if (structural_checks) {
- std::cout << "#Tolerance in stddevs for bonds: " << bond_tolerance << std::endl;
- std::cout << "#Tolerance in stddevs for angles: " << angle_tolerance << std::endl;
- std::cout << "#Clashing distance for unknown atom types: " << min_default_distance << std::endl;
- LOG_INFO("#Log entries format:");
- LOG_INFO("#BOND INFO FORMAT:" << " " << "Chain" << " " << "Residue" << " " << "ResNum" << " " << "Bond" << " " << "Min" << " " << "Max" << " " << "Observed" << " " << "Z-score" << " " << "Status");
- LOG_INFO("#ANGLE INFO FORMAT:" << " " << "Chain" << " " << "Residue" << " " << "ResNum" << " " << "Angle" << " " << "Min" << " " << "Max" << " " << "Observed" << " " << "Z-score" << " " << "Status");
- LOG_INFO("#CLASH INFO FORMAT:" << " " << "Chain1" << " " << "Residue1" << " " << "ResNum1" << " " << "Atom1" << " " << "Chain2" << " " << "Residue2" << " " << "ResNum2" << " " << "Atom2" << " " << "Min" << " " << "Observed" << " " << "Difference" << " " << "Status");
- LOG_INFO("#LDT INFO FORMAT:" << " " << "Chain1" << " " << "Residue1" << " " << "ResNum1" << " " << "Atom1" << " " << "Chain2" << " " << "Residue2" << " " << "ResNum2" << " " << "Atom2" << " " << "Min" << " " << "ModelDist" << " " << "TargetDist" << " " << "Difference" << " " << "Tolerance"
- << " " << "Status");
+ std::cout << "Parameter filename: " << parameter_filename << std::endl;
+ std::cout << "Tolerance in stddevs for bonds: " << bond_tolerance << std::endl;
+ std::cout << "Tolerance in stddevs for angles: " << angle_tolerance << std::endl;
+ std::cout << "Clashing distance for unknown atom types: " << min_default_distance << std::endl;
+ LOG_INFO("Log entries format:");
+ LOG_INFO("BOND INFO FORMAT: Chain Residue ResNum Bond Min Max Observed Z-score Status");
+ LOG_INFO("ANGLE INFO FORMAT: Chain Residue ResNum Angle Min Max Observed Z-score Status");
+ LOG_INFO("CLASH INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 Min Observed Difference Status");
+ LOG_INFO("LDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 Min ModelDist TargetDist Difference Tolerance Status");
}
for (size_t i=0; i<files.size(); ++i) {
EntityHandle model=load(files[i], profile);
if (!model) {
if (!profile.fault_tolerant) {
- exit(-1);
+ exit(-1);
}
continue;
}
- StereoChemicalParams bond_table, angle_table;
- ClashingDistances nonbonded_table(min_default_distance,min_distance_tolerance);
EntityView v=model.CreateFullView();
if (structural_checks) {
boost::filesystem::path loc(parameter_filename);
boost::filesystem::ifstream infile(loc);
if (!infile) {
- std::cout << "Couldn't find " << parameter_filename << std::endl;
+ std::cout << "Could not find " << parameter_filename << std::endl;
exit(-1);
}
std::vector<String> stereo_chemical_props;
@@ -343,47 +209,35 @@ int main (int argc, char **argv)
std::stringstream line_stream(line);
stereo_chemical_props.push_back(line);
}
- try {
- FillStereoChemicalParams("Bond",bond_table,stereo_chemical_props);
- }
- catch (Error& e) {
- std::cout << "Error reading 'Bond Lengths' section of the stereo-chemical parameter file." << std::endl;
- std::cout << e.what() << std::endl;
+ StereoChemicalParams bond_table = FillStereoChemicalParams("Bond",stereo_chemical_props);
+ if (bond_table.IsEmpty()) {
+ std::cout << "Error reading the Bond section of the stereo-chemical parameter file." << std::endl;
exit(-1);
}
- try {
- FillStereoChemicalParams("Angle",angle_table,stereo_chemical_props);
- }
- catch (Error& e) {
- std::cout << "Error reading 'Angles' section of the stereo-chemical parameter file." << std::endl;
- std::cout << e.what() << std::endl;
+ StereoChemicalParams angle_table = FillStereoChemicalParams("Angle",stereo_chemical_props);
+ if (angle_table.IsEmpty()) {
+ std::cout << "Error reading the Angles section of the stereo-chemical parameter file." << std::endl;
exit(-1);
}
- try {
- FillClashingDistances(nonbonded_table,stereo_chemical_props);
- }
- catch (Error& e) {
- std::cout << "Error reading 'Non-bonded Distances' section of the stereo-chemical parameter file." << std::endl;
- std::cout << e.what() << std::endl;
+ ClashingDistances nonbonded_table = FillClashingDistances(stereo_chemical_props,min_default_distance,min_distance_tolerance);
+
+ if (nonbonded_table.IsEmpty()) {
+ std::cout << "Error reading the Clashing section of the stereo-chemical parameter file." << std::endl;
exit(-1);
}
v=alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance);
v=alg::FilterClashes(v,nonbonded_table);
}
- float cutoffs[]={0.5,1,2,4};
- String labels[]={"localldt0.5","localldt1","localldt2","ldtlocal4"};
- float ldt=0.0;
- for (int n=0; n<4; ++n) {
- ldt+=alg::LocalDistTest(v, ref_view, cutoffs[n], 8.0,labels[n]);
- }
- ldt/=4.0;
- std::cout << "#File: " << files[i] << std::endl;
+
+ Real ldt=LDTHA(v, ref_view, radius);
+
+ std::cout << "File: " << files[i] << std::endl;
std::cout << "Global LDT score: " << ldt << std::endl;
std::cout << "Local LDT Score:" << std::endl;
- std::cout << "#Chain\tResName\tResNum\tScore" << 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;
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
index b77fb9048..f9a467103 100644
--- a/modules/mol/alg/src/local_dist_test.cc
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -334,4 +334,16 @@ Real LocalDistTest(const ost::seq::AlignmentHandle& aln,
return 0.0;
}
-}}}
+Real LDTHA(EntityView&v, const EntityView& ref_view, Real radius)
+{
+ 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, ref_view, cutoffs[n], radius,labels[n]);
+ }
+ ldt/=4.0;
+ return ldt;
+}
+
+}}}
\ No newline at end of file
diff --git a/modules/mol/alg/src/local_dist_test.hh b/modules/mol/alg/src/local_dist_test.hh
index 76ef23803..c6dd4d7e0 100644
--- a/modules/mol/alg/src/local_dist_test.hh
+++ b/modules/mol/alg/src/local_dist_test.hh
@@ -34,6 +34,10 @@ Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
Real cutoff, Real max_dist,
int ref_index=0, int mdl_index=1);
+
+Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView&v, const EntityView& ref_view, Real radius);
+
+
}}}
#endif
--
GitLab