diff --git a/modules/io/pymod/wrap_io.cc b/modules/io/pymod/wrap_io.cc index 302720245b63d534d88aceb2bba9c4d64bba0464..8704a4657d5a7f152f96d339164976fe5bbeb0ce 100644 --- a/modules/io/pymod/wrap_io.cc +++ b/modules/io/pymod/wrap_io.cc @@ -33,6 +33,7 @@ using namespace boost::python; #include <ost/io/mol/entity_io_sdf_handler.hh> #include <ost/io/mol/pdb_reader.hh> #include <ost/io/mol/dcd_io.hh> +#include <ost/io/stereochemical_params_reader.hh> using namespace ost; using namespace ost::io; @@ -120,6 +121,15 @@ BOOST_PYTHON_MODULE(_ost_io) def("LoadMAE", &LoadMAE); def("LoadPQR", &LoadPQR); + class_<io::StereoChemicalParamsReader>("StereoChemicalParamsReader", + init<String>(arg("filename"))) + .def(init<>()) + .def("Read", &io::StereoChemicalParamsReader::Read, (arg("check")=false)) + .def_readwrite("filename", &io::StereoChemicalParamsReader::filename) + .def_readwrite("bond_table", &io::StereoChemicalParamsReader::bond_table) + .def_readwrite("angle_table", &io::StereoChemicalParamsReader::angle_table) + .def_readwrite("nonbonded_table", &io::StereoChemicalParamsReader::nonbonded_table); + export_pdb_io(); export_mmcif_io(); #if OST_IMG_ENABLED diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc index 9f6c7f850fcaa869be30b337b1c94c654e5d2927..fd2b9eb2d8fc97c886f0be1f8e133c3f03c2a2f3 100644 --- a/modules/mol/alg/pymod/wrap_mol_alg.cc +++ b/modules/mol/alg/pymod/wrap_mol_alg.cc @@ -50,6 +50,7 @@ namespace { std::pair<long int,long int> (*lddt_a)(const mol::EntityView&, const mol::alg::GlobalRDMap& , std::vector<Real>, int, const String&)=&mol::alg::LocalDistDiffTest; Real (*lddt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistDiffTest; +// Real (*lddt_d)(const mol::EntityView&, std::vector<mol::EntityView>&, const mol::alg::GlobalRDMap&, mol::alg::lDDTSettings&)=&mol::alg::LocalDistDiffTest; Real (*lddt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistDiffTest; std::pair<mol::EntityView,mol::alg::ClashingInfo> (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes; std::pair<mol::EntityView,mol::alg::ClashingInfo> (*fc_b)(const mol::EntityHandle&, const mol::alg::ClashingDistances&, bool)=&mol::alg::FilterClashes; @@ -58,6 +59,18 @@ std::pair<mol::EntityView,mol::alg::StereoChemistryInfo> (*csc_b)(const mol::Ent 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; + +Real local_dist_diff_test_wrapper(const mol::EntityView& v, list& ref_list, const mol::alg::GlobalRDMap& glob_dist_list, mol::alg::lDDTSettings& settings){ + 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 LocalDistDiffTest(v, ref_list_vector, glob_dist_list, settings); +} + 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__")()); @@ -100,6 +113,126 @@ ost::mol::alg::GlobalRDMap create_distance_list_from_multiple_references(const l return ost::mol::alg::CreateDistanceListFromMultipleReferences(ref_list_vector, cutoff_list_vector, sequence_separation, max_dist); } +object lDDTSettingsInitWrapper(tuple args, dict kwargs){ + + object self = args[0]; + args = tuple(args.slice(1,_)); + + Real bond_tolerance = 12; + if(kwargs.contains("bond_tolerance")){ + bond_tolerance = extract<bool>(kwargs["bond_tolerance"]); + kwargs["bond_tolerance"].del(); + } + + Real angle_tolerance = 12; + if(kwargs.contains("angle_tolerance")){ + angle_tolerance = extract<bool>(kwargs["angle_tolerance"]); + kwargs["angle_tolerance"].del(); + } + + Real radius = 15; + if(kwargs.contains("radius")){ + radius = extract<bool>(kwargs["radius"]); + kwargs["radius"].del(); + } + + int sequence_separation = 0; + if(kwargs.contains("sequence_separation")){ + sequence_separation = extract<bool>(kwargs["sequence_separation"]); + kwargs["sequence_separation"].del(); + } + + String sel = ""; + if(kwargs.contains("sel")){ + sel = extract<bool>(kwargs["sel"]); + kwargs["sel"].del(); + } + + String parameter_file_path = ""; + if(kwargs.contains("parameter_file_path")){ + parameter_file_path = extract<bool>(kwargs["parameter_file_path"]); + kwargs["parameter_file_path"].del(); + } + + bool structural_checks = true; + if(kwargs.contains("structural_checks")){ + structural_checks = extract<bool>(kwargs["structural_checks"]); + kwargs["structural_checks"].del(); + } + + bool consistency_checks = true; + if(kwargs.contains("consistency_checks")){ + consistency_checks = extract<bool>(kwargs["consistency_checks"]); + kwargs["consistency_checks"].del(); + } + + std::vector<Real> cutoffs; + if(kwargs.contains("cutoffs")){ + list cutoff_list = extract<list>(kwargs["cutoffs"]); + int cutoff_list_length = boost::python::extract<int>(cutoff_list.attr("__len__")()); + for (int i=0; i<cutoff_list_length; i++) { + cutoffs.push_back(boost::python::extract<Real>(cutoff_list[i])); + } + kwargs["cutoffs"].del(); + } else { + cutoffs.push_back(0.5); + cutoffs.push_back(1.0); + cutoffs.push_back(2.0); + cutoffs.push_back(4.0); + } + + String label = "localldt"; + if(kwargs.contains("label")){ + label = extract<bool>(kwargs["label"]); + kwargs["label"].del(); + } + + if(len(kwargs) > 0){ + std::stringstream ss; + ss << "Invalid keywords observed when setting up MolckSettings! "; + ss << "Or did you pass the same keyword twice? "; + ss << "Valid keywords are: bond_tolerance, angle_tolerance, radius, "; + ss << "sequence_separation, sel, parameter_file_path, structural_checks, "; + ss << "consistency_checks, cutoffs, label!"; + throw std::invalid_argument(ss.str()); + } + + return self.attr("__init__")(bond_tolerance, + angle_tolerance, + radius, + sequence_separation, + sel, + parameter_file_path, + structural_checks, + consistency_checks, + cutoffs, + label); +} + + +void clean_lddt_references_wrapper(const list& ref_list) +{ + 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::CleanlDDTReferences(ref_list_vector); +} + +ost::mol::alg::GlobalRDMap prepare_lddt_global_rdmap_wrapper(const list& ref_list, mol::alg::lDDTSettings settings) +{ + 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 mol::alg::PreparelDDTGlobalRDMap(ref_list_vector, settings); +} } @@ -120,6 +253,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg) def("LocalDistDiffTest", lddt_a, (arg("sequence_separation")=0,arg("local_lddt_property_string")="")); def("LocalDistDiffTest", lddt_c, (arg("local_lddt_property_string")="")); def("LocalDistDiffTest", lddt_b, (arg("ref_index")=0, arg("mdl_index")=1)); + def("LocalDistDiffTest", &local_dist_diff_test_wrapper, (arg("v"), arg("ref_list"), ("glob_dist_list"), arg("settings"))); 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)); @@ -161,13 +295,44 @@ BOOST_PYTHON_MODULE(_ost_mol_alg) .def("GetResidueName",&mol::alg::UniqueAtomIdentifier::GetResidueName) .def("GetAtomName",&mol::alg::UniqueAtomIdentifier::GetAtomName) .def("GetQualifiedAtomName",&mol::alg::UniqueAtomIdentifier::GetQualifiedAtomName) - ; + ; + + class_<mol::alg::lDDTSettings>("lDDTSettings", no_init) + .def("__init__", raw_function(lDDTSettingsInitWrapper)) + .def(init<Real, Real, Real, int, String, String, bool, bool, std::vector<Real>&, String>()) + .def("ToString", &mol::alg::lDDTSettings::ToString) + .def("SetStereoChemicalParamsPath", &mol::alg::lDDTSettings::SetStereoChemicalParamsPath) + .def("PrintParameters", &mol::alg::lDDTSettings::PrintParameters) + .def("__repr__", &mol::alg::lDDTSettings::ToString) + .def("__str__", &mol::alg::lDDTSettings::ToString) + .def_readwrite("bond_tolerance", &mol::alg::lDDTSettings::bond_tolerance) + .def_readwrite("angle_tolerance", &mol::alg::lDDTSettings::angle_tolerance) + .def_readwrite("radius", &mol::alg::lDDTSettings::radius) + .def_readwrite("sequence_separation", &mol::alg::lDDTSettings::sequence_separation) + .def_readwrite("sel", &mol::alg::lDDTSettings::sel) + .def_readwrite("parameter_file_path", &mol::alg::lDDTSettings::parameter_file_path) + .def_readwrite("structural_checks", &mol::alg::lDDTSettings::structural_checks) + .def_readwrite("consistency_checks", &mol::alg::lDDTSettings::consistency_checks) + .def_readwrite("cutoffs", &mol::alg::lDDTSettings::cutoffs) + .def_readwrite("label", &mol::alg::lDDTSettings::label); 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); + def("CleanlDDTReferences", &clean_lddt_references_wrapper); + def("PreparelDDTGlobalRDMap", + &prepare_lddt_global_rdmap_wrapper, + (arg("ref_list"), arg("settings"))); + def("CheckStructure", + &mol::alg::CheckStructure, + (arg("ent"), arg("bond_table"), arg("angle_table"), arg("nonbonded_table"), + arg("bond_tolerance"), arg("angle_tolerance"))); + def("PrintlDDTPerResidueStats", + &mol::alg::PrintlDDTPerResidueStats, + (arg("model"), arg("glob_dist_list"), arg("settings"))); + class_<mol::alg::PDBize>("PDBize", init<int>(arg("min_polymer_size")=10)) diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc index 2cff1c815d12dc2f591b4d6a7005f6afdd62cfef..357e766906e7ad97151830c9ec7eb926b50199af 100644 --- a/modules/mol/alg/src/lddt.cc +++ b/modules/mol/alg/src/lddt.cc @@ -157,6 +157,7 @@ int main (int argc, char **argv) // sets some default values for parameters String version = OST_VERSION_STRING; lDDTSettings settings; + settings.structural_checks=false; // creates the required loading profile IOProfile profile; // parses options @@ -278,7 +279,7 @@ int main (int argc, char **argv) ref_list.push_back(ref.CreateFullView()); } CleanlDDTReferences(ref_list); - glob_dist_list = PreparelDDTGlobalRDMap(ref_list,settings.cutoffs,settings.sequence_separation,settings.radius); + glob_dist_list = PreparelDDTGlobalRDMap(ref_list, settings); files.pop_back(); // prints out parameters used in the lddt calculation diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc index 79eb285a54de807ef8393302fac693c833d92efb..a019251b96492d54585a3c0fa81055d8fe15ac3b 100644 --- a/modules/mol/alg/src/local_dist_diff_test.cc +++ b/modules/mol/alg/src/local_dist_diff_test.cc @@ -5,12 +5,25 @@ #include <ost/platform.hh> #include "local_dist_diff_test.hh" #include <boost/concept_check.hpp> +#include <boost/filesystem/convenience.hpp> #include <ost/mol/alg/consistency_checks.hh> namespace ost { namespace mol { namespace alg { namespace { +// helper_function +String vector_to_string(std::vector<Real> vec) { + std::ostringstream str; + for (unsigned int i = 0; i < vec.size(); ++i) { + str << vec[i]; + if (i+1 != vec.size()) { + str << ", "; + } + } + return str.str(); +} + // helper function bool within_tolerance(Real mdl_dist, const std::pair<Real,Real>& values, Real tol) { @@ -382,22 +395,14 @@ lDDTSettings::lDDTSettings(): bond_tolerance(12.0), radius(15.0), sequence_separation(0), sel(""), - structural_checks(false), + structural_checks(true), consistency_checks(true), - print_stats(true), label("localldt") { cutoffs.push_back(0.5); cutoffs.push_back(1.0); cutoffs.push_back(2.0); cutoffs.push_back(4.0); - try { - parameter_file_path = ost::GetSharedDataPath() + "/stereo_chemical_props.txt"; - } catch (std::runtime_error& e) { - std::cerr << "WARNING: " << e.what(); - } - // if ((! boost::filesystem::exists(parameter_file_path)) || parameter_file_path == "") { - // std::cerr << "WARNING: Could not find stereo_chemical_props.txt at the provided location." << std::endl; - // } + SetStereoChemicalParamsPath(); } lDDTSettings::lDDTSettings(Real init_bond_tolerance, @@ -409,18 +414,42 @@ lDDTSettings::lDDTSettings(Real init_bond_tolerance, bool init_structural_checks, bool init_consistency_checks, std::vector<Real>& init_cutoffs, - bool init_print_stats, String init_label): bond_tolerance(init_bond_tolerance), angle_tolerance(init_angle_tolerance), radius(init_radius), sequence_separation(init_sequence_separation), sel(init_sel), - parameter_file_path(init_parameter_file_path), structural_checks(init_structural_checks), consistency_checks(init_consistency_checks), - print_stats(init_print_stats), - label(init_label) {} + cutoffs(init_cutoffs), + label(init_label) { + SetStereoChemicalParamsPath(init_parameter_file_path); +} + +void lDDTSettings::SetStereoChemicalParamsPath(const String& path) { + if (path == ""){ + try { + std::cerr << "Setting default stereochemical parameters file path." << std::endl; + parameter_file_path = ost::GetSharedDataPath() + "/stereo_chemical_props.txt"; + if (! structural_checks) { + std::cerr << "WARNING: Stereochemical parameters file is set but structural checks are disabled" << std::endl; + } + } catch (std::runtime_error& e) { + std::cerr << "WARNING: " << e.what() << std::endl; + std::cerr << "WARNING: Parameter file setting failed - structural check will be disabled" << std::endl; + structural_checks = false; + parameter_file_path = ""; + } + } else if (! boost::filesystem::exists(path)) { + std::cerr << "WARNING: Parameter file does not exist - structural check will be disabled" << std::endl; + structural_checks = false; + parameter_file_path = ""; + } else { + parameter_file_path = path; + structural_checks = true; + } +} std::string lDDTSettings::ToString() { std::ostringstream rep; @@ -429,9 +458,11 @@ std::string lDDTSettings::ToString() { } else { rep << "Stereo-chemical and steric clash checks: Off \n"; } - rep << "Inclusion Radius: " << radius << "\n"; + rep << "Inclusion Radius: " << radius << "\n"; rep << "Sequence separation: " << sequence_separation << "\n"; + rep << "Cutoffs: " << vector_to_string(cutoffs) << "\n"; + if (structural_checks) { rep << "Parameter filename: " + parameter_file_path + "\n"; rep << "Tolerance in stddevs for bonds: " << bond_tolerance << "\n"; @@ -586,7 +617,6 @@ Real LocalDistDiffTest(const EntityView& v, std::vector<EntityView>& ref_list, const GlobalRDMap& glob_dist_list, lDDTSettings& settings) { - for (std::vector<EntityView>::const_iterator ref_list_it = ref_list.begin(); ref_list_it != ref_list.end(); ++ref_list_it) { bool cons_check = ResidueNamesMatch(v,*ref_list_it,settings.consistency_checks); @@ -693,19 +723,17 @@ void CleanlDDTReferences(std::vector<EntityView>& ref_list){ } GlobalRDMap PreparelDDTGlobalRDMap(const std::vector<EntityView>& ref_list, - std::vector<Real>& cutoff_list, - int sequence_separation, - Real max_dist){ + lDDTSettings& settings){ GlobalRDMap glob_dist_list; if (ref_list.size()==1) { std::cout << "Multi-reference mode: Off" << std::endl; - glob_dist_list = CreateDistanceList(ref_list[0], max_dist); + glob_dist_list = CreateDistanceList(ref_list[0], settings.radius); } else { std::cout << "Multi-reference mode: On" << std::endl; glob_dist_list = CreateDistanceListFromMultipleReferences(ref_list, - cutoff_list, - sequence_separation, - max_dist); + settings.cutoffs, + settings.sequence_separation, + settings.radius); } return glob_dist_list; diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh index 7759447246922c98e777b72c4065049d2097edf5..f44d972a22fbc1398232926ed9ec6329b14c08b3 100644 --- a/modules/mol/alg/src/local_dist_diff_test.hh +++ b/modules/mol/alg/src/local_dist_diff_test.hh @@ -36,7 +36,6 @@ struct lDDTSettings { bool structural_checks; bool consistency_checks; std::vector<Real> cutoffs; - bool print_stats; String label; lDDTSettings(); @@ -49,8 +48,8 @@ struct lDDTSettings { bool init_structural_checks, bool init_consistency_checks, std::vector<Real>& init_cutoffs, - bool init_print_stats, String init_label); + void SetStereoChemicalParamsPath(const String& path=""); void PrintParameters(); std::string ToString(); }; @@ -180,9 +179,7 @@ void DLLEXPORT_OST_MOL_ALG CleanlDDTReferences(std::vector<EntityView>& ref_list // Prepare GlobalRDMap from reference list GlobalRDMap DLLEXPORT_OST_MOL_ALG PreparelDDTGlobalRDMap( const std::vector<EntityView>& ref_list, - std::vector<Real>& cutoff_list, - int sequence_separation, - Real max_dist); + lDDTSettings& settings); void DLLEXPORT_OST_MOL_ALG CheckStructure(EntityView& ent, StereoChemicalParams& bond_table,