diff --git a/modules/io/src/mol/CMakeLists.txt b/modules/io/src/mol/CMakeLists.txt index e0792e93ee64ec6038451cad5a283a5c5c12144c..c88d99db8690af8efafc44a896f1cd017d03ff19 100644 --- a/modules/io/src/mol/CMakeLists.txt +++ b/modules/io/src/mol/CMakeLists.txt @@ -20,6 +20,7 @@ star_parser.cc mmcif_reader.cc mmcif_info.cc pdb_str.cc +stereochemical_params_reader.cc PARENT_SCOPE ) @@ -47,5 +48,6 @@ surface_io_handler.hh load_surface.hh surface_io_msms_handler.hh pdb_str.hh +stereochemical_params_reader.hh PARENT_SCOPE ) diff --git a/modules/io/src/mol/stereochemical_params_reader.cc b/modules/io/src/mol/stereochemical_params_reader.cc new file mode 100644 index 0000000000000000000000000000000000000000..e934177436ea4ca49f977144f2b7b3dd4081d01a --- /dev/null +++ b/modules/io/src/mol/stereochemical_params_reader.cc @@ -0,0 +1,64 @@ +//------------------------------------------------------------------------------ +// This file is part of the OpenStructure project <www.openstructure.org> +// +// Copyright (C) 2008-2011 by the OpenStructure authors +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation; either version 3.0 of the License, or (at your option) +// any later version. +// This library is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this library; if not, write to the Free Software Foundation, Inc., +// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA +//------------------------------------------------------------------------------ +#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_propsd.txt"; +} + +StereoChemicalParamsReader::StereoChemicalParamsReader(const String& init_filename): filename(init_filename) {} + +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); +} +}} // ns + diff --git a/modules/io/src/mol/stereochemical_params_reader.hh b/modules/io/src/mol/stereochemical_params_reader.hh new file mode 100644 index 0000000000000000000000000000000000000000..f8990f55268de150d6f17a28160f53e695781bfa --- /dev/null +++ b/modules/io/src/mol/stereochemical_params_reader.hh @@ -0,0 +1,40 @@ +//------------------------------------------------------------------------------ +// This file is part of the OpenStructure project <www.openstructure.org> +// +// Copyright (C) 2008-2011 by the OpenStructure authors +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation; either version 3.0 of the License, or (at your option) +// any later version. +// This library is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this library; if not, write to the Free Software Foundation, Inc., +// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA +//------------------------------------------------------------------------------ +#ifndef OST_IO_STEREOCHEMICAL_PARAMS_READER_H +#define OST_IO_STEREOCHEMICAL_PARAMS_READER_H + +#include <ost/io/module_config.hh> +#include <ost/mol/alg/filter_clashes.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); +}; + +}} // ns + +#endif \ 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 690a10aff11a238a45663186ac6b98c3ace29db5..455d1e7dbd8114aeb7ef972cc50691de2e9a0e8c 100644 --- a/modules/mol/alg/src/filter_clashes.cc +++ b/modules/mol/alg/src/filter_clashes.cc @@ -201,7 +201,7 @@ bool StereoChemicalParams::IsEmpty() const return false; } -StereoChemicalParams FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file) +StereoChemicalParams FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file, bool check) { StereoChemicalParams table; bool found=false; @@ -276,11 +276,18 @@ StereoChemicalParams FillStereoChemicalParams(const String& header, std::vector< if (found==false) { std::cout << "Could not find the relevant section in the stereo-chemical parameter file" << std::endl; return StereoChemicalParams(); - }; + }; + if (check) { + if (table.IsEmpty()) { + std::stringstream serr; + serr << "Error reading the " << header << " section of the stereo-chemical parameter file."; + throw ost::Error(serr.str()); + } + } return table; }; -ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_props_file) +ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_props_file, bool check) { ClashingDistances table; bool found=false; @@ -340,7 +347,14 @@ ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_pro if (found==false) { std::cout << "Could not find the relevant section in the stereo-chemical parameter file" << std::endl; return ClashingDistances(); - } + } + if (check) { + if (table.IsEmpty()) { + std::stringstream serr; + serr << "Error reading the Clashing section of the stereo-chemical parameter file."; + throw ost::Error(serr.str()); + } + } return table; } diff --git a/modules/mol/alg/src/filter_clashes.hh b/modules/mol/alg/src/filter_clashes.hh index f53f123dc09922937382a5bb67502b0add7603c2..fe6d12c75a7303ec637e194d341831c1200ac4b0 100644 --- a/modules/mol/alg/src/filter_clashes.hh +++ b/modules/mol/alg/src/filter_clashes.hh @@ -244,13 +244,13 @@ private: /// \brief Fills a list of reference clashing distances from the content of a parameter file /// /// Requires a list of strings holding the contents of a parameter file, one line per string -ClashingDistances DLLEXPORT_OST_MOL_ALG FillClashingDistances(std::vector<String>& stereo_chemical_props_file); +ClashingDistances DLLEXPORT_OST_MOL_ALG FillClashingDistances(std::vector<String>& stereo_chemical_props_file, bool check=false); /// \brief Fills a list of stereo-chemical statistics from the content of a parameter file /// /// Requires a list of strings holding the contents of a parameter file, one line per string /// The header can be 'Bonds' to read bond statistics or 'Angles' to read angle statistics -StereoChemicalParams DLLEXPORT_OST_MOL_ALG FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file); +StereoChemicalParams DLLEXPORT_OST_MOL_ALG FillStereoChemicalParams(const String& header, std::vector<String>& stereo_chemical_props_file, bool check=false); /// \brief Filters a structure based on detected clashes between non bonded atoms. Entity version /// diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc index dcde44f0f8111a01094e26da90ca84d0deb54731..406647f479ad8324d9f3486e300edd72235ad38a 100644 --- a/modules/mol/alg/src/lddt.cc +++ b/modules/mol/alg/src/lddt.cc @@ -32,6 +32,7 @@ #include <ost/mol/alg/filter_clashes.hh> #include <ost/io/mol/pdb_reader.hh> #include <ost/io/io_exception.hh> +#include <ost/io/stereochemical_params_reader.hh> #include <ost/conop/conop.hh> #include <ost/conop/compound_lib.hh> #include <ost/string_ref.hh> @@ -390,74 +391,27 @@ int main (int argc, char **argv) } if (structural_checks) { - // reads in parameter files - boost::filesystem::path loc(parameter_filename); - boost::filesystem::ifstream infile(loc); - if (!infile) { - std::cout << "Could not find " << parameter_filename << std::endl; - exit(-1); - } - std::vector<String> stereo_chemical_props; - String line; - while (std::getline(infile, line)) - { - std::stringstream line_stream(line); - stereo_chemical_props.push_back(line); - } - 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); - } - 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); - } - - ClashingDistances nonbonded_table = FillClashingDistances(stereo_chemical_props); - - if (nonbonded_table.IsEmpty()) { - std::cout << "Error reading the Clashing section of the stereo-chemical parameter file." << std::endl; - exit(-1); - } - // performs structural checks and filters the structure - StereoChemistryInfo stereo_chemistry_info; - try { - std::pair<EntityView,StereoChemistryInfo> csc_result = alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance); - v = csc_result.first; - stereo_chemistry_info = csc_result.second; - } catch (std::exception& e) { - std::cout << "An error occurred during the structure quality checks, stage 1:" << std::endl; + StereoChemicalParamsReader stereochemical_params(parameter_filename); + try { + stereochemical_params.Read(true); + } catch (std::exception& e) { std::cout << e.what() << std::endl; exit(-1); } - std::cout << "Average Z-Score for bond lengths: " << std::fixed << std::setprecision(5) << stereo_chemistry_info.GetAvgZscoreBonds() << std::endl; - std::cout << "Bonds outside of tolerance range: " << stereo_chemistry_info.GetBadBondCount() << " out of " << stereo_chemistry_info.GetBondCount() << std::endl; - std::cout << "Bond\tAvg Length\tAvg zscore\tNum Bonds" << std::endl; - std::map<String,BondLengthInfo> avg_bond_length_info = stereo_chemistry_info.GetAvgBondLengthInfo(); - for (std::map<String,BondLengthInfo>::const_iterator abli_it=avg_bond_length_info.begin();abli_it!=avg_bond_length_info.end();++abli_it) { - String key = (*abli_it).first; - BondLengthInfo bond_length_info = (*abli_it).second; - std::cout << key << "\t" << std::fixed << std::setprecision(5) << std::left << std::setw(10) << - bond_length_info.GetAvgLength() << "\t" << std::left << std::setw(10) << bond_length_info.GetAvgZscore() << "\t" << bond_length_info.GetCount() << std::endl; - } - std::cout << "Average Z-Score angle widths: " << std::fixed << std::setprecision(5) << stereo_chemistry_info.GetAvgZscoreAngles() << std::endl; - std::cout << "Angles outside of tolerance range: " << stereo_chemistry_info.GetBadAngleCount() << " out of " << stereo_chemistry_info.GetAngleCount() << std::endl; - ClashingInfo clash_info; + try { - std::pair<EntityView,ClashingInfo> fc_result = alg::FilterClashes(v,nonbonded_table); - v = fc_result.first; - clash_info = fc_result.second; - } catch (std::exception& e) { - std::cout << "An error occurred during the structure quality checks, stage 2:" << std::endl; + CheckStructure(v, + stereochemical_params.bond_table, + stereochemical_params.angle_table, + stereochemical_params.nonbonded_table, + bond_tolerance, + angle_tolerance); + } catch (std::exception& e) { std::cout << e.what() << std::endl; exit(-1); } - std::cout << clash_info.GetClashCount() << " non-bonded short-range distances shorter than tolerance distance" << std::endl; - std::cout << "Distances shorter than tolerance are on average shorter by: " << std::fixed << std::setprecision(5) << clash_info.GetAverageOffset() << std::endl; - } + if (cov.first==0) { std::cout << "Global LDDT score: 0.0" << std::endl; 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 b6e0ecc046cca187cedc3473657a7de62e7c0a31..3ad22db00bcea5b2d078e004e10fe909912d62b6 100644 --- a/modules/mol/alg/src/local_dist_diff_test.cc +++ b/modules/mol/alg/src/local_dist_diff_test.cc @@ -1,3 +1,4 @@ +#include <iomanip> #include <ost/log.hh> #include <ost/mol/mol.hh> #include "local_dist_diff_test.hh" @@ -566,6 +567,49 @@ GlobalRDMap PreparelDDTGlobalRDMap(const std::vector<EntityView>& ref_list, return glob_dist_list; } +void CheckStructure(EntityView& ent, + StereoChemicalParams& bond_table, + StereoChemicalParams& angle_table, + ClashingDistances& nonbonded_table, + Real bond_tolerance, + Real angle_tolerance){ + // performs structural checks and filters the structure + StereoChemistryInfo stereo_chemistry_info; + try { + std::pair<EntityView,StereoChemistryInfo> csc_result = alg::CheckStereoChemistry(ent,bond_table,angle_table,bond_tolerance,angle_tolerance); + ent = csc_result.first; + stereo_chemistry_info = csc_result.second; + } catch (std::exception& e) { + std::cout << "An error occurred during the structure quality checks, stage 1:" << std::endl; + std::cout << e.what() << std::endl; + exit(-1); + } + std::cout << "Average Z-Score for bond lengths: " << std::fixed << std::setprecision(5) << stereo_chemistry_info.GetAvgZscoreBonds() << std::endl; + std::cout << "Bonds outside of tolerance range: " << stereo_chemistry_info.GetBadBondCount() << " out of " << stereo_chemistry_info.GetBondCount() << std::endl; + std::cout << "Bond\tAvg Length\tAvg zscore\tNum Bonds" << std::endl; + std::map<String,BondLengthInfo> avg_bond_length_info = stereo_chemistry_info.GetAvgBondLengthInfo(); + for (std::map<String,BondLengthInfo>::const_iterator abli_it=avg_bond_length_info.begin();abli_it!=avg_bond_length_info.end();++abli_it) { + String key = (*abli_it).first; + BondLengthInfo bond_length_info = (*abli_it).second; + std::cout << key << "\t" << std::fixed << std::setprecision(5) << std::left << std::setw(10) << + bond_length_info.GetAvgLength() << "\t" << std::left << std::setw(10) << bond_length_info.GetAvgZscore() << "\t" << bond_length_info.GetCount() << std::endl; + } + std::cout << "Average Z-Score angle widths: " << std::fixed << std::setprecision(5) << stereo_chemistry_info.GetAvgZscoreAngles() << std::endl; + std::cout << "Angles outside of tolerance range: " << stereo_chemistry_info.GetBadAngleCount() << " out of " << stereo_chemistry_info.GetAngleCount() << std::endl; + ClashingInfo clash_info; + try { + std::pair<EntityView,ClashingInfo> fc_result = alg::FilterClashes(ent,nonbonded_table); + ent = fc_result.first; + clash_info = fc_result.second; + } catch (std::exception& e) { + std::stringstream serr; + serr << "An error occurred during the structure quality checks, stage 2: " << e.what(); + throw ost::Error(serr.str()); + } + std::cout << clash_info.GetClashCount() << " non-bonded short-range distances shorter than tolerance distance" << std::endl; + std::cout << "Distances shorter than tolerance are on average shorter by: " << std::fixed << std::setprecision(5) << clash_info.GetAverageOffset() << std::endl; +} + // debugging code /* Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_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 9eb5829fe6d9790b0030beab56067666469fbca7..05b9fd887c1631ef9e3c70123f6fb8667ead6757 100644 --- a/modules/mol/alg/src/local_dist_diff_test.hh +++ b/modules/mol/alg/src/local_dist_diff_test.hh @@ -22,6 +22,7 @@ #include <ost/mol/alg/module_config.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 { @@ -144,6 +145,13 @@ GlobalRDMap DLLEXPORT_OST_MOL_ALG PreparelDDTGlobalRDMap( int sequence_separation, Real max_dist); +void DLLEXPORT_OST_MOL_ALG CheckStructure(EntityView& ent, + StereoChemicalParams& bond_table, + StereoChemicalParams& angle_table, + ClashingDistances& nonbonded_table, + Real bond_tolerance, + Real angle_tolerance); + }}} #endif