From 42761b45f5de985f9ccacf8897ec345ec15a5e5c Mon Sep 17 00:00:00 2001
From: Valerio Mariani <valerio.mariani@unibas.ch>
Date: Fri, 25 May 2012 12:56:59 +0200
Subject: [PATCH] Added stere-chemical and clashing statistics to the output.
Tweaked logging system
---
modules/mol/alg/src/filter_clashes.cc | 41 +++++++++++++++++++++-----
modules/mol/alg/src/ldt.cc | 19 ++++++------
modules/mol/alg/src/local_dist_test.cc | 2 +-
3 files changed, 45 insertions(+), 17 deletions(-)
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index 4b81ed10b..1b163d176 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -306,6 +306,12 @@ 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)
{
+ Real running_sum_zscore_bonds=0.0;
+ Real running_sum_zscore_angles=0.0;
+ int bond_count = 0;
+ int bad_bond_count = 0;
+ int angle_count = 1;
+ int bad_angle_count = 0;
LOG_INFO("Checking stereo-chemistry")
EntityView filtered=ent.CreateEmptyView();
ResidueViewList residues=ent.GetResidueList();
@@ -347,10 +353,10 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
Real zscore = (blength - ref_length)/ref_stddev;
if (blength < min_length || blength > max_length) {
LOG_INFO("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "FAIL")
+ bad_bond_count++;
remove_sc=true;
if (always_remove_bb==true) {
remove_bb=true;
- continue;
}
String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") {
@@ -358,7 +364,9 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
}
} else {
LOG_VERBOSE("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "PASS")
- }
+ }
+ bond_count++;
+ running_sum_zscore_bonds+=zscore;
}
}
@@ -399,10 +407,10 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
Real zscore = (awidth - ref_width)/ref_stddev;
if (awidth < min_width || awidth > max_width) {
LOG_INFO("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "FAIL")
+ bad_angle_count++;
remove_sc=true;
if (always_remove_bb==true) {
remove_bb=true;
- continue;
}
String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") {
@@ -411,6 +419,8 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
} else {
LOG_VERBOSE("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "PASS")
}
+ angle_count++;
+ running_sum_zscore_angles+=zscore;
}
}
}
@@ -434,6 +444,12 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
}
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS);
}
+ 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);
+ LOG_SCRIPT("Average Z-Score for bond lengths: " << avg_zscore_bonds);
+ LOG_SCRIPT("Bonds outside of tolerance range: " << bad_bond_count << " out of " << bond_count);
+ LOG_SCRIPT("Average Z-Score angle widths: " << avg_zscore_angles);
+ LOG_SCRIPT("Angles outside of tolerance range: " << bad_angle_count << " out of " << angle_count);
return filtered;
}
@@ -446,6 +462,8 @@ EntityView CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalPar
EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb)
{
+ int distance_count = 0;
+ int bad_distance_count = 0;
LOG_INFO("Filtering non-bonded clashes")
EntityView filtered=ent.CreateEmptyView();
ResidueViewList residues=ent.GetResidueList();
@@ -461,6 +479,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
for (AtomViewList::const_iterator
j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
AtomView atom=*j;
+
String ele1=atom.GetElement();
if (ele1=="H" || ele1=="D") {
continue;
@@ -469,6 +488,9 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
for (AtomViewList::iterator
k=within.begin(), e3=within.end(); k!=e3; ++k) {
AtomView atom2=*k;
+
+
+
if (atom2==atom) {
continue;
}
@@ -477,25 +499,28 @@ 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.
+ // and we instead allow S-S distances down to 1.8.
if (atom.GetHandle().FindBondToAtom(atom2.GetHandle()).IsValid()) {
continue;
}
+
+
Real d=geom::Length2(atom.GetPos()-atom2.GetPos());
std::pair <Real,Real> distance_tolerance=min_distances.GetClashingDistance(ele1, ele2);
Real distance=distance_tolerance.first;
Real tolerance=distance_tolerance.second;
Real threshold=distance-tolerance;
+
if (d<threshold*threshold) {
LOG_INFO(atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetName() << " " << atom2.GetResidue().GetNumber() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "FAIL")
-
+ bad_distance_count++;
remove_sc=true;
if (always_remove_bb==true) {
remove_bb=true;
- continue;
}
String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") {
@@ -503,7 +528,8 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
}
} else {
LOG_VERBOSE("CLASH:" << " " << atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetNumber() << " " << atom2.GetResidue().GetName() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "PASS")
- }
+ }
+ distance_count++;
}
}
@@ -525,6 +551,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
}
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS);
}
+ LOG_SCRIPT(bad_distance_count << " out of " << distance_count << " non-bonded short-range distances checked shorter than tolerance distance");
return filtered;
}
diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc
index 5dbab61cf..48dfc9a6c 100644
--- a/modules/mol/alg/src/ldt.cc
+++ b/modules/mol/alg/src/ldt.cc
@@ -106,7 +106,7 @@ int main (int argc, char **argv)
Real min_distance_tolerance = 0.0;
Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0;
- Real radius=8.0;
+ Real radius=15.0;
IOProfile profile;
profile.bond_feasibility_check=false;
@@ -164,23 +164,22 @@ int main (int argc, char **argv)
exit(-1);
}
int verbosity_level=0;
+ int ost_verbosity_level=2;
if (vm.count("verbosity")) {
verbosity_level=vm["verbosity"].as<int>();
if (verbosity_level==0) {
- Logger::Instance().PushVerbosityLevel(0);
+ ost_verbosity_level=2;
} else if (verbosity_level==1) {
- Logger::Instance().PushVerbosityLevel(3);
+ ost_verbosity_level=3;
} else if (verbosity_level==2) {
- Logger::Instance().PushVerbosityLevel(4);
+ ost_verbosity_level=4;
} else {
std::cout << "Verbosity level " << verbosity_level << " is not available" << std::endl;
exit(-1);
}
}
- if (verbosity_level>0) {
- LOG_INFO("LDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDist Difference Tolerance Status");
- }
+ Logger::Instance().PushVerbosityLevel(ost_verbosity_level);
if (vm.count("bond_tolerance")) {
bond_tolerance=vm["bond_tolerance"].as<Real>();
@@ -218,8 +217,8 @@ int main (int argc, char **argv)
LOG_INFO("BOND INFO FORMAT: Chain Residue ResNum Bond Min Max Observed Z-score Status");
LOG_INFO("ANGLE INFO FORMAT: Chain Residue ResNum Angle Min Max Observed Z-score Status");
LOG_INFO("CLASH INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 Observed Difference Status");
- LOG_INFO("LDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDist Difference Tolerance Status");
}
+ LOG_INFO("LDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDist Difference Tolerance Status");
for (size_t i=0; i<files.size(); ++i) {
EntityHandle model=load(files[i], profile);
if (!model) {
@@ -230,6 +229,7 @@ int main (int argc, char **argv)
}
EntityView v=model.CreateFullView();
+
// The code in this following block is only used to make CASP9 models load correctly and normally commented out
EntityView model2=model.Select("aname!=CEN,NV,OT1,OT,CAY,CY,OXT,1OCT,NT,OT2,2OCT,OVL1,OC1,O1,OC2,O2,OVU1");
EntityView v1=model2.Select("not (rname==GLY and aname==CB)");
@@ -249,6 +249,7 @@ int main (int argc, char **argv)
aitv.SetElement(firstletter);
}
}
+ v=v1;
std::cout << "File: " << files[i] << std::endl;
std::pair<int,int> cov = compute_coverage(v,glob_dist_list);
std::cout << "Coverage: " << (float(cov.first)/float(cov.second)) << " (" << cov.first << " out of " << cov.second << " residues)" << std::endl;
@@ -284,7 +285,7 @@ int main (int argc, char **argv)
std::cout << "Error reading the Clashing section of the stereo-chemical parameter file." << std::endl;
exit(-1);
}
- EntityView v=alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance);
+ v=alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance);
cov = compute_coverage(v,glob_dist_list);
std::cout << "Coverage after stereo-chemical checks: " << (float(cov.first)/float(cov.second)) << " (" << cov.first << " out of " << cov.second << " residues)" << std::endl;
v=alg::FilterClashes(v,nonbonded_table);
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
index 83d842226..430dd31be 100644
--- a/modules/mol/alg/src/local_dist_test.cc
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -132,7 +132,7 @@ std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list
Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
ReferenceDistance ref_dist=*ai;
if (within_tolerance(mdl_dist,ref_dist,tol).first) {
- LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
+ LOG_INFO("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
<< " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
<< mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "PASS")
overlap.first+=1;
--
GitLab