Skip to content
Snippets Groups Projects
Commit 42761b45 authored by Valerio Mariani's avatar Valerio Mariani
Browse files

Added stere-chemical and clashing statistics to the output. Tweaked logging system

parent e59cf0f8
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
modules/mol/alg/src/filter_clashes.cc
+ 34
− 7
View file @ 42761b45
...@@ -306,6 +306,12 @@ ClashingDistances FillClashingDistances(std::vector<String>& stereo_chemical_pro ...@@ -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) 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") LOG_INFO("Checking stereo-chemistry")
EntityView filtered=ent.CreateEmptyView(); EntityView filtered=ent.CreateEmptyView();
ResidueViewList residues=ent.GetResidueList(); ResidueViewList residues=ent.GetResidueList();
...@@ -347,10 +353,10 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam ...@@ -347,10 +353,10 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
Real zscore = (blength - ref_length)/ref_stddev; Real zscore = (blength - ref_length)/ref_stddev;
if (blength < min_length || blength > max_length) { 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") LOG_INFO("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "FAIL")
bad_bond_count++;
remove_sc=true; remove_sc=true;
if (always_remove_bb==true) { if (always_remove_bb==true) {
remove_bb=true; remove_bb=true;
continue;
} }
String name=atom.GetName(); String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") { if (name=="CA" || name=="N" || name=="O" || name=="C") {
...@@ -358,7 +364,9 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam ...@@ -358,7 +364,9 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
} }
} else { } else {
LOG_VERBOSE("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "PASS") 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 ...@@ -399,10 +407,10 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
Real zscore = (awidth - ref_width)/ref_stddev; Real zscore = (awidth - ref_width)/ref_stddev;
if (awidth < min_width || awidth > max_width) { 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") LOG_INFO("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "FAIL")
bad_angle_count++;
remove_sc=true; remove_sc=true;
if (always_remove_bb==true) { if (always_remove_bb==true) {
remove_bb=true; remove_bb=true;
continue;
} }
String name=atom.GetName(); String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") { if (name=="CA" || name=="N" || name=="O" || name=="C") {
...@@ -411,6 +419,8 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam ...@@ -411,6 +419,8 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
} else { } else {
LOG_VERBOSE("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "PASS") 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 ...@@ -434,6 +444,12 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
} }
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS); 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; return filtered;
} }
...@@ -446,6 +462,8 @@ EntityView CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalPar ...@@ -446,6 +462,8 @@ EntityView CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalPar
EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb) 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") LOG_INFO("Filtering non-bonded clashes")
EntityView filtered=ent.CreateEmptyView(); EntityView filtered=ent.CreateEmptyView();
ResidueViewList residues=ent.GetResidueList(); ResidueViewList residues=ent.GetResidueList();
...@@ -461,6 +479,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis ...@@ -461,6 +479,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
for (AtomViewList::const_iterator for (AtomViewList::const_iterator
j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) { j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
AtomView atom=*j; AtomView atom=*j;
String ele1=atom.GetElement(); String ele1=atom.GetElement();
if (ele1=="H" || ele1=="D") { if (ele1=="H" || ele1=="D") {
continue; continue;
...@@ -469,6 +488,9 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis ...@@ -469,6 +488,9 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
for (AtomViewList::iterator for (AtomViewList::iterator
k=within.begin(), e3=within.end(); k!=e3; ++k) { k=within.begin(), e3=within.end(); k!=e3; ++k) {
AtomView atom2=*k; AtomView atom2=*k;
if (atom2==atom) { if (atom2==atom) {
continue; continue;
} }
...@@ -477,25 +499,28 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis ...@@ -477,25 +499,28 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
continue; continue;
} }
// In theory, this should also trigger for disulfide bonds, but // In theory, this should also trigger for disulfide bonds, but
// since we don't detect disulfides correctly, we can't count on that // 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()) { if (atom.GetHandle().FindBondToAtom(atom2.GetHandle()).IsValid()) {
continue; continue;
} }
Real d=geom::Length2(atom.GetPos()-atom2.GetPos()); Real d=geom::Length2(atom.GetPos()-atom2.GetPos());
std::pair <Real,Real> distance_tolerance=min_distances.GetClashingDistance(ele1, ele2); std::pair <Real,Real> distance_tolerance=min_distances.GetClashingDistance(ele1, ele2);
Real distance=distance_tolerance.first; Real distance=distance_tolerance.first;
Real tolerance=distance_tolerance.second; Real tolerance=distance_tolerance.second;
Real threshold=distance-tolerance; Real threshold=distance-tolerance;
if (d<threshold*threshold) { 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") 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; remove_sc=true;
if (always_remove_bb==true) { if (always_remove_bb==true) {
remove_bb=true; remove_bb=true;
continue;
} }
String name=atom.GetName(); String name=atom.GetName();
if (name=="CA" || name=="N" || name=="O" || name=="C") { if (name=="CA" || name=="N" || name=="O" || name=="C") {
...@@ -503,7 +528,8 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis ...@@ -503,7 +528,8 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
} }
} else { } 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") 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 ...@@ -525,6 +551,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis
} }
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS); 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; return filtered;
} }
......
modules/mol/alg/src/ldt.cc
+ 10
− 9
View file @ 42761b45
...@@ -106,7 +106,7 @@ int main (int argc, char **argv) ...@@ -106,7 +106,7 @@ int main (int argc, char **argv)
Real min_distance_tolerance = 0.0; Real min_distance_tolerance = 0.0;
Real bond_tolerance = 8.0; Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0; Real angle_tolerance = 8.0;
Real radius=8.0; Real radius=15.0;
IOProfile profile; IOProfile profile;
profile.bond_feasibility_check=false; profile.bond_feasibility_check=false;
...@@ -164,23 +164,22 @@ int main (int argc, char **argv) ...@@ -164,23 +164,22 @@ int main (int argc, char **argv)
exit(-1); exit(-1);
} }
int verbosity_level=0; int verbosity_level=0;
int ost_verbosity_level=2;
if (vm.count("verbosity")) { if (vm.count("verbosity")) {
verbosity_level=vm["verbosity"].as<int>(); verbosity_level=vm["verbosity"].as<int>();
if (verbosity_level==0) { if (verbosity_level==0) {
Logger::Instance().PushVerbosityLevel(0); ost_verbosity_level=2;
} else if (verbosity_level==1) { } else if (verbosity_level==1) {
Logger::Instance().PushVerbosityLevel(3); ost_verbosity_level=3;
} else if (verbosity_level==2) { } else if (verbosity_level==2) {
Logger::Instance().PushVerbosityLevel(4); ost_verbosity_level=4;
} else { } else {
std::cout << "Verbosity level " << verbosity_level << " is not available" << std::endl; std::cout << "Verbosity level " << verbosity_level << " is not available" << std::endl;
exit(-1); exit(-1);
} }
} }
if (verbosity_level>0) { Logger::Instance().PushVerbosityLevel(ost_verbosity_level);
LOG_INFO("LDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDist Difference Tolerance Status");
}
if (vm.count("bond_tolerance")) { if (vm.count("bond_tolerance")) {
bond_tolerance=vm["bond_tolerance"].as<Real>(); bond_tolerance=vm["bond_tolerance"].as<Real>();
...@@ -218,8 +217,8 @@ int main (int argc, char **argv) ...@@ -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("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("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("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) { for (size_t i=0; i<files.size(); ++i) {
EntityHandle model=load(files[i], profile); EntityHandle model=load(files[i], profile);
if (!model) { if (!model) {
...@@ -230,6 +229,7 @@ int main (int argc, char **argv) ...@@ -230,6 +229,7 @@ int main (int argc, char **argv)
} }
EntityView v=model.CreateFullView(); EntityView v=model.CreateFullView();
// The code in this following block is only used to make CASP9 models load correctly and normally commented out // 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 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)"); EntityView v1=model2.Select("not (rname==GLY and aname==CB)");
...@@ -249,6 +249,7 @@ int main (int argc, char **argv) ...@@ -249,6 +249,7 @@ int main (int argc, char **argv)
aitv.SetElement(firstletter); aitv.SetElement(firstletter);
} }
} }
v=v1;
std::cout << "File: " << files[i] << std::endl; std::cout << "File: " << files[i] << std::endl;
std::pair<int,int> cov = compute_coverage(v,glob_dist_list); 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; 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) ...@@ -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; std::cout << "Error reading the Clashing section of the stereo-chemical parameter file." << std::endl;
exit(-1); 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); 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; 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); v=alg::FilterClashes(v,nonbonded_table);
......
modules/mol/alg/src/local_dist_test.cc
+ 1
− 1
View file @ 42761b45
...@@ -132,7 +132,7 @@ std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list ...@@ -132,7 +132,7 @@ std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list
Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos()); Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
ReferenceDistance ref_dist=*ai; ReferenceDistance ref_dist=*ai;
if (within_tolerance(mdl_dist,ref_dist,tol).first) { 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() << " " << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
<< mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "PASS") << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "PASS")
overlap.first+=1; overlap.first+=1;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment