From c89a8ba6551df8c7cae5c7fca12ef14ee5ee34d6 Mon Sep 17 00:00:00 2001 From: Valerio Mariani <valerio.mariani@unibas.ch> Date: Thu, 11 Apr 2013 14:45:16 +0200 Subject: [PATCH] Moved logging logic out of functions and into the main lddt code Moved printed messages out of the FilterClashes and CheckStereoChemistry functions. Breaks those two functions --- modules/mol/alg/pymod/wrap_mol_alg.cc | 35 +++++++++-- modules/mol/alg/src/filter_clashes.cc | 30 ++++------ modules/mol/alg/src/filter_clashes.hh | 86 +++++++++++++++++++++++++-- modules/mol/alg/src/lddt.cc | 27 ++++++++- 4 files changed, 147 insertions(+), 31 deletions(-) diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc index b775b6396..7b403f02d 100644 --- a/modules/mol/alg/pymod/wrap_mol_alg.cc +++ b/modules/mol/alg/pymod/wrap_mol_alg.cc @@ -42,10 +42,10 @@ 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_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistDiffTest; -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; +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; +std::pair<mol::EntityView,mol::alg::StereoChemistryInfo> (*csc_a)(const mol::EntityView&, const mol::alg::StereoChemicalParams&, const mol::alg::StereoChemicalParams&, Real, Real, bool)=&mol::alg::CheckStereoChemistry; +std::pair<mol::EntityView,mol::alg::StereoChemistryInfo> (*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; @@ -170,5 +170,30 @@ BOOST_PYTHON_MODULE(_ost_mol_alg) def("ResidueNamesMatch",&mol::alg::ResidueNamesMatch, (arg("probe"), arg("reference"), arg("log_as_error")=false)); - + class_<mol::alg::BondLengthInfo> ("BondLengthInfo" ,init <>()) + .def(init<Real,Real,int>()) + .def("GetAvgLength",&mol::alg::BondLengthInfo::GetAvgLength) + .def("GetAvgZscore",&mol::alg::BondLengthInfo::GetAvgZscore) + .def("GetCount",&mol::alg::BondLengthInfo::GetCount) + ; + + class_<mol::alg::ClashingInfo> ("ClashingInfo" ,init <>()) + .def(init<int,Real>()) + .def("GetClashCount",&mol::alg::ClashingInfo::GetClashCount) + .def("GetAverageOffset",&mol::alg::ClashingInfo::GetAverageOffset) + ; + + class_<mol::alg::StereoChemistryInfo> ("StereoChemistryInfo" ,init <>()) + .def(init<Real,int,int,Real,int,int, const std::map<String,mol::alg::BondLengthInfo>&>()) + .def("GetAvgZscoreBonds",&mol::alg::StereoChemistryInfo::GetAvgZscoreBonds) + .def("GetBadBondCount",&mol::alg::StereoChemistryInfo::GetBadBondCount) + .def("GetBondCount",&mol::alg::StereoChemistryInfo::GetBondCount) + .def("GetAvgZscoreAngles",&mol::alg::StereoChemistryInfo::GetAvgZscoreAngles) + .def("GetBadAngleCount",&mol::alg::StereoChemistryInfo::GetBadAngleCount) + .def("GetAngleCount",&mol::alg::StereoChemistryInfo::GetAngleCount) + .def("GetAvgBondLengthInfo",&mol::alg::StereoChemistryInfo::GetAvgBondLengthInfo) + ; + + + } diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc index b7a22ba13..e544ce95b 100644 --- a/modules/mol/alg/src/filter_clashes.cc +++ b/modules/mol/alg/src/filter_clashes.cc @@ -329,7 +329,7 @@ ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_pro } -EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb) +std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb) { Real running_sum_zscore_bonds=0.0; Real running_sum_zscore_angles=0.0; @@ -487,10 +487,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam } Real avg_zscore_bonds = running_sum_zscore_bonds/static_cast<float>(bond_count); Real avg_zscore_angles = running_sum_zscore_angles/static_cast<float>(angle_count); - std::cout << "Average Z-Score for bond lengths: " << std::fixed << std::setprecision(5) << avg_zscore_bonds << std::endl; - std::cout << "Bonds outside of tolerance range: " << bad_bond_count << " out of " << bond_count << std::endl; - std::cout << "Bond\tAvg Length\tAvg zscore\tNum Bonds" << std::endl; - + std::map<String,BondLengthInfo> avg_bond_length_info; for (std::map<String,Real>::const_iterator bls_it=bond_length_sum.begin();bls_it!=bond_length_sum.end();++bls_it) { String key = (*bls_it).first; int counter=bond_counter_sum[key]; @@ -498,22 +495,23 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam Real sum_bond_zscore=bond_zscore_sum[key]; Real avg_length=sum_bond_length/static_cast<Real>(counter); Real avg_zscore=sum_bond_zscore/static_cast<Real>(counter); - std::cout << key << "\t" << std::fixed << std::setprecision(5) << std::left << std::setw(10) << avg_length << "\t" << std::left << std::setw(10) << avg_zscore << "\t" << counter << std::endl; + BondLengthInfo bond_info(avg_length,avg_zscore,counter); + avg_bond_length_info[key]=bond_info; } - std::cout << "Average Z-Score angle widths: " << std::fixed << std::setprecision(5) << avg_zscore_angles << std::endl; - std::cout << "Angles outside of tolerance range: " << bad_angle_count << " out of " << angle_count << std::endl; + StereoChemistryInfo info(avg_zscore_bonds, bad_bond_count, bond_count,avg_zscore_angles, + bad_angle_count, angle_count,avg_bond_length_info); filtered.AddAllInclusiveBonds(); - return filtered; + return std::make_pair<EntityView,StereoChemistryInfo>(filtered,info); } -EntityView CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb) +std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb) { return CheckStereoChemistry(ent.CreateFullView(), bond_table, angle_table, bond_tolerance, angle_tolerance, always_remove_bb); } -EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb) +std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb) { int distance_count = 0; int bad_distance_count = 0; @@ -550,7 +548,6 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis continue; } - // In theory, this should also trigger for disulfide bonds, but // since we don't detect disulfides correctly, we can't count on that // and we instead allow S-S distances down to 1.8. @@ -607,15 +604,14 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis if (bad_distance_count!=0) { average_offset = average_offset_sum / static_cast<Real>(bad_distance_count); } - std::cout << bad_distance_count << " 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) << average_offset << std::endl; + ClashingInfo info(bad_distance_count,average_offset); filtered.AddAllInclusiveBonds(); - return filtered; + return std::make_pair<EntityView,ClashingInfo>(filtered,info); } -EntityView FilterClashes(const EntityHandle& ent, - const ClashingDistances& min_distances, bool always_remove_bb) +std::pair<EntityView,ClashingInfo> FilterClashes(const EntityHandle& ent, + const ClashingDistances& min_distances, bool always_remove_bb) { return FilterClashes(ent.CreateFullView(), min_distances, always_remove_bb); } diff --git a/modules/mol/alg/src/filter_clashes.hh b/modules/mol/alg/src/filter_clashes.hh index e294cd489..d8e1c4c2e 100644 --- a/modules/mol/alg/src/filter_clashes.hh +++ b/modules/mol/alg/src/filter_clashes.hh @@ -24,6 +24,80 @@ namespace ost { namespace mol { namespace alg { +class BondLengthInfo +{ +public: + BondLengthInfo(): avg_length_(0),avg_zscore_(0),count_(0) {} + BondLengthInfo(Real avg_length,Real avg_zscore, int count): + avg_length_(avg_length), + avg_zscore_(avg_zscore), + count_(count){} + Real GetAvgLength() {return avg_length_;} + Real GetAvgZscore() {return avg_zscore_;} + int GetCount() {return count_;} + +private: + Real avg_length_; + Real avg_zscore_; + int count_; +}; + + +class ClashingInfo +{ + +public: + ClashingInfo(): clash_count_(0), average_offset_ (0) {} + ClashingInfo (int clash_count, Real average_offset): + clash_count_(clash_count), average_offset_ (average_offset) {} + int GetClashCount() {return clash_count_;} + Real GetAverageOffset() {return average_offset_;} + +private: + int clash_count_; + Real average_offset_; +}; + +class StereoChemistryInfo +{ +public: + StereoChemistryInfo(): + avg_zscore_bonds_(0), + bad_bond_count_(0), + bond_count_(0), + avg_zscore_angles_(0), + bad_angle_count_(0), + angle_count_(0), + avg_bond_length_info_(std::map<String,BondLengthInfo>()) {} + StereoChemistryInfo(Real avg_zscore_bonds, int bad_bond_count, int bond_count, + Real avg_zscore_angles, int bad_angle_count, int angle_count, + const std::map<String,BondLengthInfo>& avg_bond_length_info): + avg_zscore_bonds_(avg_zscore_bonds), + bad_bond_count_(bad_bond_count), + bond_count_(bond_count), + avg_zscore_angles_(avg_zscore_angles), + bad_angle_count_(bad_angle_count), + angle_count_(angle_count), + avg_bond_length_info_(avg_bond_length_info) {} + Real GetAvgZscoreBonds() {return avg_zscore_bonds_;} + int GetBadBondCount() {return bad_bond_count_;} + int GetBondCount() {return bond_count_;} + Real GetAvgZscoreAngles() {return avg_zscore_angles_;} + int GetBadAngleCount() {return bad_angle_count_;} + int GetAngleCount() {return angle_count_;} + std::map<String,BondLengthInfo> GetAvgBondLengthInfo() {return avg_bond_length_info_;} + +private: + Real avg_zscore_bonds_; + int bad_bond_count_; + int bond_count_; + Real avg_zscore_angles_; + int bad_angle_count_; + int angle_count_; + std::map<String,BondLengthInfo> avg_bond_length_info_; +}; + + /// \brief List of reference atom-atom distances to detect clashes between non-bonded atoms class DLLEXPORT_OST_MOL_ALG ClashingDistances { @@ -54,7 +128,7 @@ public: private: - std::map <String,std::pair<float,float> > min_distance_; + std::map <String,std::pair<Real,Real> > min_distance_; Real default_min_distance_; Real default_min_distance_tolerance_; bool valid_flag_; @@ -91,7 +165,7 @@ public: private: - std::map<std::pair<String,String>,std::pair<float,float> > params_; + std::map<std::pair<String,String>,std::pair<Real,Real> > params_; }; @@ -113,7 +187,7 @@ StereoChemicalParams DLLEXPORT_OST_MOL_ALG FillStereoChemicalParams(const String /// If a clash is detected in the backbone, all atoms in the residue are removed. This behavior is changed /// by the always_remove_bb flag: when the flag is set to true all atoms in the residue are removed even if /// a clash is just detected in the side-chain -EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityView& ent, +std::pair<EntityView,ClashingInfo> DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb=false); /// \brief Filters a structure based on detected clashes between non bonded atoms. Handle version @@ -123,7 +197,7 @@ EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityView& ent, /// If a clash is detected in the backbone, all atoms in the residue are removed. This behavior is changed /// by the always_remove_bb flag: when the flag is set to true all atoms in the residue are removed even if /// a clash is just detected in the side-chain -EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityHandle& ent, +std::pair<EntityView,ClashingInfo> DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityHandle& ent, const ClashingDistances& min_distances, bool always_remove_bb=false); /// \brief Filters a structure based on detected stereo-chemical violations. Entity version @@ -133,7 +207,7 @@ EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityHandle& ent, /// all atoms in the side chain are removed from the structure. If a violation is detected in the backbone, all /// atoms in the residue are removed. This behavior is changed by the always_remove_bb flag: when the flag is /// set to true all atoms in the residue are removed even if a violation is just detected in the side-chain -EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityView& ent, +std::pair<EntityView,StereoChemistryInfo> DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, @@ -147,7 +221,7 @@ EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityView& ent, /// all atoms in the side chain are removed from the structure. If a violation is detected in the backbone, all /// atoms in the residue are removed. This behavior is changed by the always_remove_bb flag: when the flag is /// set to true all atoms in the residue are removed even if a violation is just detected in the side-chain -EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityHandle& ent, +std::pair<EntityView,StereoChemistryInfo> DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc index 1787018b2..b98fa7d7c 100644 --- a/modules/mol/alg/src/lddt.cc +++ b/modules/mol/alg/src/lddt.cc @@ -373,21 +373,42 @@ int main (int argc, char **argv) std::cout << "Error reading the Clashing section of the stereo-chemical parameter file." << std::endl; exit(-1); } - // performs structural checks and filters the structure + // performs structural checks and filters the structure + StereoChemistryInfo stereo_chemistry_info; try { - v=alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance); + 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; 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 { - v=alg::FilterClashes(v,nonbonded_table); + 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; 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; -- GitLab