diff --git a/modules/io/pymod/wrap_io.cc b/modules/io/pymod/wrap_io.cc
index 8704a4657d5a7f152f96d339164976fe5bbeb0ce..4b64e61891f710fc4f4aecf4fa9d4eb44add1942 100644
--- a/modules/io/pymod/wrap_io.cc
+++ b/modules/io/pymod/wrap_io.cc
@@ -68,6 +68,11 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_handle_ov,
BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_view_ov,
save_ent_view, 2, 3)
+ost::mol::alg::StereoChemicalProps (*read_props_a)(String filename, bool check) = &ReadStereoChemicalPropsFile;
+ost::mol::alg::StereoChemicalProps (*read_props_b)(bool check) = &ReadStereoChemicalPropsFile;
+
+
+Real (*lddt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistDiffTest;
}
void export_pdb_io();
@@ -121,14 +126,11 @@ 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);
+ def("ReadStereoChemicalPropsFile", read_props_a,
+ (arg("filename"), arg("check")=true));
+ def("ReadStereoChemicalPropsFile", read_props_b,
+ (arg("check")=true));
+ def("GetStereoChemicalPropsFile", &GetStereoChemicalPropsFile);
export_pdb_io();
export_mmcif_io();
diff --git a/modules/io/src/mol/stereochemical_params_reader.cc b/modules/io/src/mol/stereochemical_params_reader.cc
index f0317a3ea19ae8dd78b675a4c08576555e5e659a..0201caf573168a860ae10190cef51ec526238ed2 100644
--- a/modules/io/src/mol/stereochemical_params_reader.cc
+++ b/modules/io/src/mol/stereochemical_params_reader.cc
@@ -18,47 +18,50 @@
//------------------------------------------------------------------------------
#include <boost/filesystem/fstream.hpp>
#include <ost/platform.hh>
-#include <ost/mol/alg/filter_clashes.hh>
#include <ost/io/stereochemical_params_reader.hh>
-namespace {
- std::vector<String> ReadStereoChemicalFile(const String& filename){
- boost::filesystem::path loc(filename);
- boost::filesystem::ifstream infile(loc);
- if (!infile) {
- std::stringstream serr;
- serr << "Could not find " << filename;
- throw ost::Error(serr.str());
- }
- std::vector<String> stereo_chemical_props;
- String line;
- while (std::getline(infile, line))
- {
- std::stringstream line_stream(line);
- stereo_chemical_props.push_back(line);
- }
-
- return stereo_chemical_props;
- }
-}
-
namespace ost { namespace io {
-StereoChemicalParamsReader::StereoChemicalParamsReader() {
- filename = ost::GetSharedDataPath() + "/stereo_chemical_props.txt";
-}
+ost::mol::alg::StereoChemicalProps ReadStereoChemicalPropsFile(String filename, bool check){
+ boost::filesystem::path loc(filename);
+ boost::filesystem::ifstream infile(loc);
+ if (!infile) {
+ std::stringstream serr;
+ serr << "Could not find parameter file '" << filename << "'";
+ throw ost::Error(serr.str());
+ }
+ std::vector<String> stereo_chemical_props;
+ String line;
+ while (std::getline(infile, line))
+ {
+ std::stringstream line_stream(line);
+ stereo_chemical_props.push_back(line);
+ }
-StereoChemicalParamsReader::StereoChemicalParamsReader(const String& init_filename): filename(init_filename) {}
+ ost::mol::alg::StereoChemicalParams bond_table;
+ ost::mol::alg::StereoChemicalParams angle_table;
+ ost::mol::alg::ClashingDistances nonbonded_table;
-void StereoChemicalParamsReader::Read(bool check) {
- std::vector<String> stereo_chemical_props = ReadStereoChemicalFile(filename);
- // Bonds
bond_table = ost::mol::alg::FillStereoChemicalParams("Bond", stereo_chemical_props, check);
// Angles
angle_table = ost::mol::alg::FillStereoChemicalParams("Angle", stereo_chemical_props, check);
// Not bonded
nonbonded_table = ost::mol::alg::FillClashingDistances(stereo_chemical_props, check);
+
+ return ost::mol::alg::StereoChemicalProps(bond_table, angle_table, nonbonded_table);
+}
+
+ost::mol::alg::StereoChemicalProps ReadStereoChemicalPropsFile(bool check) {
+ String filename = GetStereoChemicalPropsFile();
+ return ReadStereoChemicalPropsFile(filename, check);
}
+
+String GetStereoChemicalPropsFile() {
+ String filename;
+ filename = ost::GetSharedDataPath() + "/stereo_chemical_props.txt";
+ return filename;
+}
+
}} // ns
diff --git a/modules/io/src/mol/stereochemical_params_reader.hh b/modules/io/src/mol/stereochemical_params_reader.hh
index f8990f55268de150d6f17a28160f53e695781bfa..a5d33acce7888bdcd88b74ec0aecd38a120eb9d6 100644
--- a/modules/io/src/mol/stereochemical_params_reader.hh
+++ b/modules/io/src/mol/stereochemical_params_reader.hh
@@ -20,20 +20,13 @@
#define OST_IO_STEREOCHEMICAL_PARAMS_READER_H
#include <ost/io/module_config.hh>
-#include <ost/mol/alg/filter_clashes.hh>
+#include <ost/mol/alg/local_dist_diff_test.hh>
namespace ost { namespace io {
-struct StereoChemicalParamsReader {
- String filename;
- ost::mol::alg::StereoChemicalParams bond_table;
- ost::mol::alg::StereoChemicalParams angle_table;
- ost::mol::alg::ClashingDistances nonbonded_table;
-
- StereoChemicalParamsReader();
- StereoChemicalParamsReader(const String& filename);
- void Read(bool check=false);
-};
+ost::mol::alg::StereoChemicalProps ReadStereoChemicalPropsFile(String filename, bool check=false);
+ost::mol::alg::StereoChemicalProps ReadStereoChemicalPropsFile(bool check=false);
+String GetStereoChemicalPropsFile();
}} // ns
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index ac7471758bd2ab4a88230d9d7a2aa7a20324563f..d4a7c208c0d518bb08141a4b2f5cd34c788c802b 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -20,6 +20,7 @@
#include <boost/python.hpp>
#include <boost/python/raw_function.hpp>
#include <boost/python/suite/indexing/map_indexing_suite.hpp>
+#include <ost/log.hh>
#include <ost/config.hh>
#include <ost/mol/alg/local_dist_diff_test.hh>
#include <ost/mol/alg/distance_test_common.hh>
@@ -151,12 +152,6 @@ object lDDTSettingsInitWrapper(tuple args, dict kwargs){
kwargs["sel"].del();
}
- String parameter_file_path = "";
- if(kwargs.contains("parameter_file_path")){
- parameter_file_path = extract<String>(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"]);
@@ -192,7 +187,7 @@ object lDDTSettingsInitWrapper(tuple args, dict kwargs){
if(len(kwargs) > 0){
std::stringstream ss;
- ss << "Invalid keywords observed when setting up MolckSettings! ";
+ ss << "Invalid keywords observed when setting up lDDTSettings! ";
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, ";
@@ -205,13 +200,65 @@ object lDDTSettingsInitWrapper(tuple args, dict kwargs){
radius,
sequence_separation,
sel,
- parameter_file_path,
structural_checks,
consistency_checks,
cutoffs,
label);
}
+object lDDTScorerInitWrapper(tuple args, dict kwargs){
+
+ object self = args[0];
+ args = tuple(args.slice(1, len(args)));
+
+ std::vector<ost::mol::EntityHandle> reference_list_vector;
+ if(kwargs.contains("references")){
+ list reference_list = boost::python::extract<list>(kwargs["references"]);
+ int reference_list_length = boost::python::extract<int>(reference_list.attr("__len__")());
+ for (int i=0; i<reference_list_length; i++) {
+ reference_list_vector.push_back(boost::python::extract<ost::mol::EntityHandle>(reference_list[i]));
+ }
+ kwargs["references"].del();
+ } else {
+ throw std::invalid_argument("'references' argument is required");
+ }
+
+ ost::mol::EntityHandle model;
+ if(kwargs.contains("model")){
+ model = boost::python::extract<ost::mol::EntityHandle>(kwargs["model"]);
+ kwargs["model"].del();
+ } else {
+ throw std::invalid_argument("'model' argument is required");
+ }
+
+ ost::mol::alg::lDDTSettings settings;
+ if(kwargs.contains("settings")){
+ settings = extract<ost::mol::alg::lDDTSettings>(kwargs["settings"]);
+ kwargs["settings"].del();
+ } else {
+ throw std::invalid_argument("'settings' argument is required");
+ }
+
+ ost::mol::alg::StereoChemicalProps stereochemical_params;
+ if(kwargs.contains("stereochemical_params")){
+ stereochemical_params = extract<ost::mol::alg::StereoChemicalProps>(kwargs["stereochemical_params"]);
+ kwargs["stereochemical_params"].del();
+ }
+
+ if(len(kwargs) > 0){
+ std::stringstream ss;
+ ss << "Invalid keywords observed when setting up lDDTScorer! ";
+ ss << "Or did you pass the same keyword twice? ";
+ ss << "Valid keywords are: settings and stereochemical_params!";
+ throw std::invalid_argument(ss.str());
+ }
+
+ return self.attr("__init__")(reference_list_vector,
+ model,
+ settings,
+ stereochemical_params);
+}
+
void clean_lddt_references_wrapper(const list& reference_list)
{
@@ -225,7 +272,10 @@ void clean_lddt_references_wrapper(const list& reference_list)
return ost::mol::alg::CleanlDDTReferences(reference_list_vector);
}
-ost::mol::alg::GlobalRDMap prepare_lddt_global_rdmap_wrapper(const list& reference_list, mol::alg::lDDTSettings settings)
+ost::mol::alg::GlobalRDMap prepare_lddt_global_rdmap_wrapper(const list& reference_list,
+ list& cutoff_list,
+ int sequence_separation,
+ Real max_dist)
{
int reference_list_length = boost::python::extract<int>(reference_list.attr("__len__")());
std::vector<ost::mol::EntityView> reference_list_vector(reference_list_length);
@@ -233,14 +283,22 @@ ost::mol::alg::GlobalRDMap prepare_lddt_global_rdmap_wrapper(const list& referen
for (int i=0; i<reference_list_length; i++) {
reference_list_vector[i] = boost::python::extract<ost::mol::EntityView>(reference_list[i]);
}
+
+ int cutoff_list_length = boost::python::extract<int>(cutoff_list.attr("__len__")());
+ std::vector<Real> cutoff_list_vector(cutoff_list_length);
+
+ for (int i=0; i<cutoff_list_length; i++) {
+ cutoff_list_vector[i] = boost::python::extract<Real>(cutoff_list[i]);
+ }
- return mol::alg::PreparelDDTGlobalRDMap(reference_list_vector, settings);
+ return mol::alg::PreparelDDTGlobalRDMap(reference_list_vector, cutoff_list_vector, sequence_separation, max_dist);
}
list get_lddt_per_residue_stats_wrapper(mol::EntityHandle& model,
ost::mol::alg::GlobalRDMap& distance_list,
- ost::mol::alg::lDDTSettings& settings) {
- std::vector<mol::alg::lDDTLocalScore> scores = GetlDDTPerResidueStats(model, distance_list, settings);
+ bool structural_checks,
+ String label) {
+ std::vector<mol::alg::lDDTLocalScore> scores = GetlDDTPerResidueStats(model, distance_list, structural_checks, label);
list local_scores_list;
for (std::vector<mol::alg::lDDTLocalScore>::const_iterator sit = scores.begin(); sit != scores.end(); ++sit) {
local_scores_list.append(*sit);
@@ -248,7 +306,16 @@ list get_lddt_per_residue_stats_wrapper(mol::EntityHandle& model,
return local_scores_list;
}
-void print_lddt_per_residue_stats_wrapper(list& scores, mol::alg::lDDTSettings settings){
+list get_local_scores_wrapper(mol::alg::lDDTScorer& scorer) {
+ std::vector<mol::alg::lDDTLocalScore> scores = scorer.GetLocalScores();
+ list local_scores_list;
+ for (std::vector<mol::alg::lDDTLocalScore>::const_iterator sit = scores.begin(); sit != scores.end(); ++sit) {
+ local_scores_list.append(*sit);
+ }
+ return local_scores_list;
+}
+
+void print_lddt_per_residue_stats_wrapper(list& scores, bool structural_checks, int cutoffs_size){
int scores_length = boost::python::extract<int>(scores.attr("__len__")());
std::vector<mol::alg::lDDTLocalScore> scores_vector(scores_length);
@@ -256,7 +323,7 @@ void print_lddt_per_residue_stats_wrapper(list& scores, mol::alg::lDDTSettings s
scores_vector[i] = boost::python::extract<mol::alg::lDDTLocalScore>(scores[i]);
}
- return mol::alg::PrintlDDTPerResidueStats(scores_vector, settings);
+ return mol::alg::PrintlDDTPerResidueStats(scores_vector, structural_checks, cutoffs_size);
}
}
@@ -326,9 +393,8 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
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(init<Real, Real, Real, int, 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)
@@ -337,7 +403,6 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
.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)
@@ -346,6 +411,8 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
class_<mol::alg::lDDTLocalScore>("lDDTLocalScore", init<String, String, int, String, String, Real, int, int>())
.def("ToString", &mol::alg::lDDTLocalScore::ToString)
.def("GetHeader", &mol::alg::lDDTLocalScore::GetHeader)
+ .def("__str__", &mol::alg::lDDTLocalScore::Repr)
+ .def("__repr__", &mol::alg::lDDTLocalScore::Repr)
.def_readwrite("cname", &mol::alg::lDDTLocalScore::cname)
.def_readwrite("rname", &mol::alg::lDDTLocalScore::rname)
.def_readwrite("rnum", &mol::alg::lDDTLocalScore::rnum)
@@ -354,6 +421,24 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
.def_readwrite("local_lddt", &mol::alg::lDDTLocalScore::local_lddt)
.def_readwrite("conserved_dist", &mol::alg::lDDTLocalScore::conserved_dist)
.def_readwrite("total_dist", &mol::alg::lDDTLocalScore::total_dist);
+
+ class_<mol::alg::lDDTScorer>("lDDTScorer", no_init)
+ .def("__init__", raw_function(lDDTScorerInitWrapper))
+ .def(init<std::vector<mol::EntityHandle>&, mol::EntityHandle&, mol::alg::lDDTSettings&, mol::alg::StereoChemicalProps&>())
+ .add_property("global_score", &mol::alg::lDDTScorer::GetGlobalScore)
+ .add_property("conserved_residues", &mol::alg::lDDTScorer::GetNumConservedResidues)
+ .add_property("total_residues", &mol::alg::lDDTScorer::GetNumTotalResidues)
+ .def("PrintPerResidueStats", &mol::alg::lDDTScorer::PrintPerResidueStats)
+ .add_property("local_scores", &get_local_scores_wrapper)
+ .add_property("is_valid", &mol::alg::lDDTScorer::IsValid);
+
+ class_<mol::alg::StereoChemicalProps>("StereoChemicalProps",
+ init<mol::alg::StereoChemicalParams&,
+ mol::alg::StereoChemicalParams&,
+ mol::alg::ClashingDistances&>())
+ .def_readwrite("bond_table", &mol::alg::StereoChemicalProps::bond_table)
+ .def_readwrite("angle_table", &mol::alg::StereoChemicalProps::angle_table)
+ .def_readwrite("nonbonded_table", &mol::alg::StereoChemicalProps::nonbonded_table);
def("FillClashingDistances",&fill_clashing_distances_wrapper);
def("FillStereoChemicalParams",&fill_stereochemical_params_wrapper);
@@ -363,17 +448,17 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
def("CleanlDDTReferences", &clean_lddt_references_wrapper);
def("PreparelDDTGlobalRDMap",
&prepare_lddt_global_rdmap_wrapper,
- (arg("reference_list"), arg("settings")));
+ (arg("reference_list"), arg("cutoffs"), arg("sequence_separation"), arg("radius")));
def("CheckStructure",
&mol::alg::CheckStructure,
(arg("ent"), arg("bond_table"), arg("angle_table"), arg("nonbonded_table"),
arg("bond_tolerance"), arg("angle_tolerance")));
def("GetlDDTPerResidueStats",
&get_lddt_per_residue_stats_wrapper,
- (arg("model"), arg("distance_list"), arg("settings")));
+ (arg("model"), arg("distance_list"), arg("structural_checks"), arg("label")));
def("PrintlDDTPerResidueStats",
&print_lddt_per_residue_stats_wrapper,
- (arg("scores"), arg("settings")));
+ (arg("scores"), arg("structural_checks"), arg("cutoff_list_length")));
class_<mol::alg::PDBize>("PDBize",
diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index 60156b86b4515c773994cd3031252915525c07fe..e9117757b0696ca982561ca8f82dd5aad328f491 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;
+ String parameter_file_path;
settings.structural_checks=false;
// creates the required loading profile
IOProfile profile;
@@ -224,7 +225,7 @@ int main (int argc, char **argv)
}
if (vm.count("parameter-file")) {
- settings.parameter_file_path=vm["parameter-file"].as<String>();
+ parameter_file_path=vm["parameter-file"].as<String>();
} else if (settings.structural_checks==true) {
std::cout << "Please specify a stereo-chemical parameter file" << std::endl;
exit(-1);
@@ -290,7 +291,10 @@ int main (int argc, char **argv)
}
}
CleanlDDTReferences(ref_list);
- glob_dist_list = PreparelDDTGlobalRDMap(ref_list, settings);
+ glob_dist_list = PreparelDDTGlobalRDMap(ref_list,
+ settings.cutoffs,
+ settings.sequence_separation,
+ settings.radius);
files.pop_back();
// prints out parameters used in the lddt calculation
@@ -326,9 +330,9 @@ int main (int argc, char **argv)
std::cout << "File: " << files[i] << std::endl;
if (settings.structural_checks) {
- StereoChemicalParamsReader stereochemical_params(settings.parameter_file_path);
+ StereoChemicalProps stereochemical_params;
try {
- stereochemical_params.Read(true);
+ stereochemical_params = ost::io::ReadStereoChemicalPropsFile(parameter_file_path, true);
} catch (std::exception& e) {
std::cout << e.what() << std::endl;
exit(-1);
@@ -352,8 +356,8 @@ int main (int argc, char **argv)
// prints the residue-by-residue statistics
std::vector<lDDTLocalScore> local_scores;
- local_scores = GetlDDTPerResidueStats(model, glob_dist_list, settings);
- PrintlDDTPerResidueStats(local_scores, settings);
+ local_scores = GetlDDTPerResidueStats(model, glob_dist_list, settings.structural_checks, settings.label);
+ PrintlDDTPerResidueStats(local_scores, settings.structural_checks, settings.cutoffs.size());
}
return 0;
}
diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc
index 9c988d1005c1ca2614a48fbf74e7258a07caf135..d02f41a32e2d89a7ba1c2255b6e4d4c0033557c2 100644
--- a/modules/mol/alg/src/local_dist_diff_test.cc
+++ b/modules/mol/alg/src/local_dist_diff_test.cc
@@ -7,6 +7,7 @@
#include <boost/concept_check.hpp>
#include <boost/filesystem/convenience.hpp>
#include <ost/mol/alg/consistency_checks.hh>
+#include <ost/io/stereochemical_params_reader.hh>
namespace ost { namespace mol { namespace alg {
@@ -388,7 +389,19 @@ bool IsStandardResidue(String rn)
return true;
}
return false;
-}
+}
+
+StereoChemicalProps::StereoChemicalProps():
+ is_valid(false) {}
+
+StereoChemicalProps::StereoChemicalProps(
+ ost::mol::alg::StereoChemicalParams& init_bond_table,
+ ost::mol::alg::StereoChemicalParams& init_angle_table,
+ ost::mol::alg::ClashingDistances& init_nonbonded_table):
+ bond_table(init_bond_table),
+ angle_table(init_angle_table),
+ nonbonded_table(init_nonbonded_table),
+ is_valid(true) {}
lDDTSettings::lDDTSettings(): bond_tolerance(12.0),
angle_tolerance(12.0),
@@ -402,7 +415,6 @@ lDDTSettings::lDDTSettings(): bond_tolerance(12.0),
cutoffs.push_back(1.0);
cutoffs.push_back(2.0);
cutoffs.push_back(4.0);
- SetStereoChemicalParamsPath();
}
lDDTSettings::lDDTSettings(Real init_bond_tolerance,
@@ -410,7 +422,6 @@ lDDTSettings::lDDTSettings(Real init_bond_tolerance,
Real init_radius,
int init_sequence_separation,
String init_sel,
- String init_parameter_file_path,
bool init_structural_checks,
bool init_consistency_checks,
std::vector<Real>& init_cutoffs,
@@ -423,33 +434,7 @@ lDDTSettings::lDDTSettings(Real init_bond_tolerance,
structural_checks(init_structural_checks),
consistency_checks(init_consistency_checks),
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;
- }
-}
+ label(init_label) {}
std::string lDDTSettings::ToString() {
std::ostringstream rep;
@@ -464,7 +449,6 @@ std::string lDDTSettings::ToString() {
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";
rep << "Tolerance in stddevs for angles: " << angle_tolerance << "\n";
}
@@ -517,6 +501,12 @@ String lDDTLocalScore::ToString(bool structural_checks) const {
return outstr.str();
}
+String lDDTLocalScore::Repr() const {
+ std::stringstream outstr;
+ outstr << "<lDDTLocalScore " << cname << "." << rname << "." << rnum << ">";
+ return outstr.str();
+}
+
String lDDTLocalScore::GetHeader(bool structural_checks, int cutoffs_length) {
std::stringstream outstr;
if (structural_checks) {
@@ -527,6 +517,205 @@ String lDDTLocalScore::GetHeader(bool structural_checks, int cutoffs_length) {
return outstr.str();
}
+
+lDDTScorer::lDDTScorer(std::vector<EntityHandle>& init_references,
+ ost::mol::EntityHandle& init_model,
+ lDDTSettings& init_settings):
+ references(init_references),
+ model(init_model),
+ settings(init_settings),
+ _score_calculated(false),
+ _score_valid(false),
+ _has_local_scores(false),
+ _num_con_res(-1),
+ _num_tot_res(-1),
+ _global_score(-1.0)
+ {
+ model_view = model.GetChainList()[0].Select("peptide=true");
+ _PrepareReferences();
+ _PrepareGlobalRDMap();
+ if (settings.structural_checks) {
+ if (!stereochemical_params.is_valid) {
+ throw std::runtime_error(
+ "Structural checks are ON. Please provide stereochemical_params to enable structural checks.");
+ }
+ try {
+ CheckStructure(model_view,
+ stereochemical_params.bond_table,
+ stereochemical_params.angle_table,
+ stereochemical_params.nonbonded_table,
+ settings.bond_tolerance,
+ settings.angle_tolerance);
+ } catch (std::exception& e) {
+ std::cout << e.what() << std::endl;
+ std::stringstream errstr;
+ errstr << "Structure check failed: " << e.what();
+ throw std::runtime_error(errstr.str());
+ }
+ }
+ }
+
+lDDTScorer::lDDTScorer(std::vector<EntityHandle>& init_references,
+ ost::mol::EntityHandle& init_model,
+ lDDTSettings& init_settings,
+ StereoChemicalProps& init_stereochemical_params):
+ references(init_references),
+ model(init_model),
+ settings(init_settings),
+ stereochemical_params(init_stereochemical_params),
+ _score_calculated(false),
+ _score_valid(false),
+ _has_local_scores(false),
+ _num_con_res(-1),
+ _num_tot_res(-1),
+ _global_score(-1.0) {
+ model_view = model.GetChainList()[0].Select("peptide=true");
+ _PrepareReferences();
+ _PrepareGlobalRDMap();
+ if (settings.structural_checks) {
+ if (!stereochemical_params.is_valid) {
+ throw std::runtime_error(
+ "Please provide stereochemical_params to enable structural_checks.");
+ }
+ try {
+ CheckStructure(model_view,
+ stereochemical_params.bond_table,
+ stereochemical_params.angle_table,
+ stereochemical_params.nonbonded_table,
+ settings.bond_tolerance,
+ settings.angle_tolerance);
+ } catch (std::exception& e) {
+ std::cout << e.what() << std::endl;
+ std::stringstream errstr;
+ errstr << "Structure check failed: " << e.what();
+ throw std::runtime_error(errstr.str());
+ }
+ }
+ }
+
+Real lDDTScorer::GetGlobalScore(){
+ if (!_score_calculated) {
+ _ComputelDDT();
+ }
+ return _global_score;
+}
+
+int lDDTScorer::GetNumConservedResidues(){
+ if (!_score_calculated) {
+ _ComputelDDT();
+ }
+ return _num_con_res;
+}
+
+int lDDTScorer::GetNumTotalResidues(){
+ if (!_score_calculated) {
+ _ComputelDDT();
+ }
+ return _num_tot_res;
+}
+
+std::vector<lDDTLocalScore> lDDTScorer::GetLocalScores(){
+ if (!_has_local_scores) {
+ _GetLocallDDT();
+ }
+ return _local_scores;
+}
+
+void lDDTScorer::PrintPerResidueStats(){
+ if (!_has_local_scores) {
+ _GetLocallDDT();
+ }
+ PrintlDDTPerResidueStats(_local_scores,
+ settings.structural_checks,
+ settings.cutoffs.size());
+}
+
+void lDDTScorer::_PrepareReferences(){
+ for (unsigned int i = 0; i < references.size(); i++) {
+ if (settings.sel != ""){
+ std::cout << "Performing \"" << settings.sel << "\" selection on reference " << references[i].GetName() << std::endl;
+ try {
+ references_view.push_back(references[i].Select(settings.sel));
+ } catch (const ost::mol::QueryError& e) {
+ std::stringstream errstr;
+ errstr << "Provided selection argument failed." << std::endl << e.GetFormattedMessage();
+ throw std::runtime_error(errstr.str());
+ }
+ }
+ else {
+ references_view.push_back(references[i].CreateFullView());
+ }
+ }
+
+ CleanlDDTReferences(references_view);
+}
+
+void lDDTScorer::_PrepareGlobalRDMap(){
+ glob_dist_list = PreparelDDTGlobalRDMap(references_view,
+ settings.cutoffs,
+ settings.sequence_separation,
+ settings.radius);
+}
+
+bool lDDTScorer::IsValid(){
+ return _score_valid;
+}
+
+void lDDTScorer::_CheckConsistency(){
+ for (std::vector<EntityView>::const_iterator ref_list_it = references_view.begin();
+ ref_list_it != references_view.end(); ++ref_list_it) {
+ bool cons_check = ResidueNamesMatch(model_view, *ref_list_it, settings.consistency_checks);
+ if (cons_check == false) {
+ if (settings.consistency_checks == true) {
+ throw std::runtime_error("Residue names in model and in reference structure(s) are inconsistent.");
+ } else {
+ LOG_WARNING("Residue names in model and in reference structure(s) are inconsistent.");
+ }
+ }
+ }
+}
+
+void lDDTScorer::_ComputelDDT(){
+ std::pair<int,int> cov = ComputeCoverage(model_view,glob_dist_list);
+ if (cov.second == -1) {
+ std::cout << "Coverage: 0 (0 out of 0 residues)" << std::endl;
+ } else {
+ std::cout << "Coverage: " << (float(cov.first)/float(cov.second)) << " (" << cov.first << " out of " << cov.second << " residues)" << std::endl;
+ }
+
+ if (cov.first == 0) {
+ std::cout << "Global LDDT score: 0.0" << std::endl;
+ _num_tot_res = 0;
+ _num_con_res = 0;
+ _global_score = 0.0;
+ _score_calculated = true;
+ _score_valid = false;
+ }
+
+ std::pair<int,int> total_ov=alg::LocalDistDiffTest(model_view, glob_dist_list, settings.cutoffs, settings.sequence_separation, settings.label);
+ Real lddt = static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
+ std::cout << "Global LDDT score: " << std::setprecision(4) << lddt << std::endl;
+ std::cout << "(" << std::fixed << total_ov.first << " conserved distances out of " << total_ov.second
+ << " checked, over " << settings.cutoffs.size() << " thresholds)" << std::endl;
+ _num_tot_res = total_ov.second ? total_ov.second : 1;
+ _num_con_res = total_ov.first;
+ _global_score = lddt;
+ _score_calculated = true;
+ _score_valid = true;
+}
+
+void lDDTScorer::_GetLocallDDT(){
+ if (!_score_calculated){
+ _ComputelDDT();
+ }
+ _local_scores = GetlDDTPerResidueStats(model,
+ glob_dist_list,
+ settings.structural_checks,
+ settings.label);
+ _has_local_scores = true;
+}
+
+
GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist)
{
GlobalRDMap dist_list;
@@ -772,17 +961,19 @@ void CleanlDDTReferences(std::vector<EntityView>& ref_list){
}
GlobalRDMap PreparelDDTGlobalRDMap(const std::vector<EntityView>& ref_list,
- lDDTSettings& settings){
+ std::vector<Real>& cutoff_list,
+ int sequence_separation,
+ Real max_dist){
GlobalRDMap glob_dist_list;
if (ref_list.size()==1) {
std::cout << "Multi-reference mode: Off" << std::endl;
- glob_dist_list = CreateDistanceList(ref_list[0], settings.radius);
+ glob_dist_list = CreateDistanceList(ref_list[0], max_dist);
} else {
std::cout << "Multi-reference mode: On" << std::endl;
glob_dist_list = CreateDistanceListFromMultipleReferences(ref_list,
- settings.cutoffs,
- settings.sequence_separation,
- settings.radius);
+ cutoff_list,
+ sequence_separation,
+ max_dist);
}
return glob_dist_list;
@@ -831,7 +1022,10 @@ void CheckStructure(EntityView& ent,
std::cout << "Distances shorter than tolerance are on average shorter by: " << std::fixed << std::setprecision(5) << clash_info.GetAverageOffset() << std::endl;
}
-std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model, GlobalRDMap& glob_dist_list, lDDTSettings& settings){
+std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model,
+ GlobalRDMap& glob_dist_list,
+ bool structural_checks,
+ String label){
std::vector<lDDTLocalScore> scores;
EntityView outv = model.GetChainList()[0].Select("peptide=true");
for (ChainViewList::const_iterator ci = outv.GetChainList().begin(),
@@ -844,7 +1038,7 @@ std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model, GlobalRD
bool assessed = false;
String assessed_string="No";
String quality_problems_string;
- if (settings.structural_checks) {
+ if (structural_checks) {
quality_problems_string="No";
} else {
quality_problems_string="NA";
@@ -867,15 +1061,15 @@ std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model, GlobalRD
}
if (assessed==true) {
- if (ritv.HasProp(settings.label)) {
- lddt_local=ritv.GetFloatProp(settings.label);
+ if (ritv.HasProp(label)) {
+ lddt_local=ritv.GetFloatProp(label);
std::stringstream stkeylddt;
stkeylddt << std::fixed << std::setprecision(4) << lddt_local;
lddt_local_string=stkeylddt.str();
- conserved_dist=ritv.GetIntProp(settings.label+"_conserved");
- total_dist=ritv.GetIntProp(settings.label+"_total");
+ conserved_dist=ritv.GetIntProp(label+"_conserved");
+ total_dist=ritv.GetIntProp(label+"_total");
} else {
- std::cout << settings.label << std::endl;
+ std::cout << label << std::endl;
lddt_local = 0;
lddt_local_string="0.0000";
conserved_dist = 0;
@@ -898,8 +1092,8 @@ std::vector<lDDTLocalScore> GetlDDTPerResidueStats(EntityHandle& model, GlobalRD
return scores;
}
-void PrintlDDTPerResidueStats(std::vector<lDDTLocalScore>& scores, lDDTSettings& settings){
- if (settings.structural_checks) {
+void PrintlDDTPerResidueStats(std::vector<lDDTLocalScore>& scores, bool structural_checks, int cutoffs_length){
+ if (structural_checks) {
std::cout << "Local LDDT Scores:" << std::endl;
std::cout << "(A 'Yes' in the 'Quality Problems' column stands for problems" << std::endl;
std::cout << "in the side-chain of a residue, while a 'Yes+' for problems" << std::endl;
@@ -907,9 +1101,9 @@ void PrintlDDTPerResidueStats(std::vector<lDDTLocalScore>& scores, lDDTSettings&
} else {
std::cout << "Local LDDT Scores:" << std::endl;
}
- std::cout << lDDTLocalScore::GetHeader(settings.structural_checks, settings.cutoffs.size()) << std::endl;
+ std::cout << lDDTLocalScore::GetHeader(structural_checks, cutoffs_length) << std::endl;
for (std::vector<lDDTLocalScore>::const_iterator sit = scores.begin(); sit != scores.end(); ++sit) {
- std::cout << sit->ToString(settings.structural_checks) << std::endl;
+ std::cout << sit->ToString(structural_checks) << std::endl;
}
std::cout << std::endl;
}
diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh
index 21aee5a1d7a035ae21450b05c84212bc103a1b9d..ebcf71579a838422b8c5a80b35c3b66ac8c36084 100644
--- a/modules/mol/alg/src/local_dist_diff_test.hh
+++ b/modules/mol/alg/src/local_dist_diff_test.hh
@@ -20,19 +20,32 @@
#define OST_MOL_ALG_LOCAL_DIST_TEST_HH
#include <ost/mol/alg/module_config.hh>
+#include <ost/mol/entity_handle.hh>
#include <ost/seq/alignment_handle.hh>
#include <ost/mol/alg/distance_test_common.hh>
#include <ost/mol/alg/filter_clashes.hh>
namespace ost { namespace mol { namespace alg {
+struct StereoChemicalProps
+{
+ bool is_valid;
+ ost::mol::alg::StereoChemicalParams bond_table;
+ ost::mol::alg::StereoChemicalParams angle_table;
+ ost::mol::alg::ClashingDistances nonbonded_table;
+
+ StereoChemicalProps();
+ StereoChemicalProps(ost::mol::alg::StereoChemicalParams& init_bond_table,
+ ost::mol::alg::StereoChemicalParams& init_angle_table,
+ ost::mol::alg::ClashingDistances& init_nonbonded_table);
+};
+
struct lDDTSettings {
Real bond_tolerance;
Real angle_tolerance;
Real radius;
int sequence_separation;
String sel;
- String parameter_file_path;
bool structural_checks;
bool consistency_checks;
std::vector<Real> cutoffs;
@@ -44,12 +57,10 @@ struct lDDTSettings {
Real init_radius,
int init_sequence_separation,
String init_sel,
- String init_parameter_file_path,
bool init_structural_checks,
bool init_consistency_checks,
std::vector<Real>& init_cutoffs,
String init_label);
- void SetStereoChemicalParamsPath(const String& path="");
void PrintParameters();
std::string ToString();
};
@@ -76,10 +87,53 @@ struct lDDTLocalScore {
int init_total_dist);
String ToString(bool structural_checks) const;
+ String Repr() const;
static String GetHeader(bool structural_checks, int cutoffs_length);
};
+class lDDTScorer
+{
+ public:
+ lDDTSettings settings;
+ ost::mol::EntityHandle model;
+ std::vector<EntityHandle> references;
+ EntityView model_view;
+ std::vector<EntityView> references_view;
+ GlobalRDMap glob_dist_list;
+ StereoChemicalProps stereochemical_params;
+
+ lDDTScorer(std::vector<EntityHandle>& init_references,
+ ost::mol::EntityHandle& init_model,
+ lDDTSettings& init_settings);
+ lDDTScorer(std::vector<EntityHandle>& init_references,
+ ost::mol::EntityHandle& init_model,
+ lDDTSettings& init_settings,
+ StereoChemicalProps& init_stereochemical_params);
+ Real GetGlobalScore();
+ std::vector<lDDTLocalScore> GetLocalScores();
+ int GetNumConservedResidues(); // number of conserved distances in the model
+ int GetNumTotalResidues(); // the number of total distances in the reference structure
+ void PrintPerResidueStats();
+ bool IsValid();
+
+ private:
+ bool _score_calculated;
+ bool _score_valid;
+ bool _has_local_scores;
+ // number of conserved distances in the model and
+ // the number of total distances in the reference structure
+ int _num_con_res;
+ int _num_tot_res;
+ Real _global_score;
+ std::vector<lDDTLocalScore> _local_scores;
+ void _ComputelDDT();
+ void _GetLocallDDT();
+ void _CheckConsistency();
+ void _PrepareReferences();
+ void _PrepareGlobalRDMap();
+};
+
std::pair<int,int> DLLEXPORT_OST_MOL_ALG ComputeCoverage(const EntityView& v,const GlobalRDMap& glob_dist_list);
bool DLLEXPORT_OST_MOL_ALG IsResnumInGlobalRDMap(const ResNum& resnum, const GlobalRDMap& glob_dist_list);
@@ -205,7 +259,9 @@ 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,
- lDDTSettings& settings);
+ std::vector<Real>& cutoff_list,
+ int sequence_separation,
+ Real max_dist);
void DLLEXPORT_OST_MOL_ALG CheckStructure(EntityView& ent,
StereoChemicalParams& bond_table,
@@ -216,10 +272,12 @@ void DLLEXPORT_OST_MOL_ALG CheckStructure(EntityView& ent,
std::vector<lDDTLocalScore> DLLEXPORT_OST_MOL_ALG GetlDDTPerResidueStats(EntityHandle& model,
GlobalRDMap& glob_dist_list,
- lDDTSettings& settings);
+ bool structural_checks,
+ String label);
void DLLEXPORT_OST_MOL_ALG PrintlDDTPerResidueStats(std::vector<lDDTLocalScore>& scores,
- lDDTSettings& settings);
+ bool structural_checks,
+ int cutoffs_length);
}}}