diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index b775b6396338c2ee45357bcfd56d259e6704b002..7b403f02d410525986db487beec6397c678c4fc3 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 b7a22ba13b0e30852669e132cf77a089d0874428..e544ce95b724bd0d01372bd18696b7d1ee1d0e3f 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 e294cd4895ad9790a3b001a73873f8b7b4812cc8..d8e1c4c2ea7978536e7e4ad8073bf0e08e1b87a3 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 1787018b2c6b61df01827811b4aa355d0e37e5da..b98fa7d7c14217edb0cabcc47d7b3b695f52c9c9 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;