diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc index 27ffb844853f7d68821d3356ec6895ef4389b4f1..712672c216fd28e749dc81a538e1cceb7cd7c353 100644 --- a/modules/mol/alg/src/filter_clashes.cc +++ b/modules/mol/alg/src/filter_clashes.cc @@ -19,6 +19,7 @@ #include <ost/log.hh> #include <ost/mol/mol.hh> #include <sstream> +#include <iomanip> #include <math.h> #include "filter_clashes.hh" #include <ost/units.hh> @@ -460,6 +461,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam if (remove_bb) { LOG_INFO("ACTION: removing whole residue " << res); + res.SetBoolProp("stereo_chemical_violation_backbone",true); continue; } if (remove_sc) { @@ -472,13 +474,14 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam filtered.AddAtom(atom); } } + res.SetBoolProp("stereo_chemical_violation_sidechain",true); continue; } 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("Average Z-Score for bond lengths: " << std::fixed << std::setprecision(5) << avg_zscore_bonds); LOG_SCRIPT("Bonds outside of tolerance range: " << bad_bond_count << " out of " << bond_count); LOG_SCRIPT("Bond\tAvg Length\tAvg zscore\tNum Bonds") @@ -489,9 +492,9 @@ 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); - LOG_SCRIPT(key << "\t" << avg_length << "\t" << avg_zscore << "\t" << counter); + LOG_SCRIPT(key << "\t" << std::fixed << std::setprecision(5) << std::left << std::setw(10) << avg_length << "\t" << std::left << std::setw(10) << avg_zscore << "\t" << counter); } - LOG_SCRIPT("Average Z-Score angle widths: " << avg_zscore_angles); + LOG_SCRIPT("Average Z-Score angle widths: " << std::fixed << std::setprecision(5) << avg_zscore_angles); LOG_SCRIPT("Angles outside of tolerance range: " << bad_angle_count << " out of " << angle_count); return filtered; } @@ -575,6 +578,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis if (remove_bb) { LOG_VERBOSE("ACTION: removing whole residue " << res); + res.SetBoolProp("steric_clash",true); continue; } if (remove_sc) { @@ -587,6 +591,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis filtered.AddAtom(atom); } } + res.SetBoolProp("steric_clash",true); continue; } filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS); @@ -596,7 +601,7 @@ EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_dis average_offset = average_offset_sum / static_cast<Real>(bad_distance_count); } LOG_SCRIPT(bad_distance_count << " non-bonded short-range distances shorter than tolerance distance"); - LOG_SCRIPT("Distances shorter than tolerance are on average shorter by: " << average_offset); + LOG_SCRIPT("Distances shorter than tolerance are on average shorter by: " << std::fixed << std::setprecision(5) << average_offset); return filtered; } diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc index a7615007b548256067a4c2c093c7f248fc1b7012..23069def36cedb5560cbd63d21ca0822dfaf8e12 100644 --- a/modules/mol/alg/src/lddt.cc +++ b/modules/mol/alg/src/lddt.cc @@ -111,6 +111,17 @@ std::pair<int,int> compute_coverage (const EntityView& v,const GlobalRDMap& glob return std::make_pair<int,int>(first,second); } +bool is_resnum_in_globalrdmap(const ResNum& resnum, const GlobalRDMap& glob_dist_list) +{ + for (GlobalRDMap::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) { + ResNum rn = i->first; + if (rn==resnum) { + return true; + } + } + return false; +} + int main (int argc, char **argv) { // sets some default values for parameters @@ -241,8 +252,9 @@ int main (int argc, char **argv) if (!ref) { exit(-1); } - ref_list.push_back(ref.CreateFullView()); - glob_dist_list = CreateDistanceList(ref.CreateFullView(),radius); + EntityView refview=ref.GetChainList()[0].Select("peptide=true"); + ref_list.push_back(refview); + glob_dist_list = CreateDistanceList(refview,radius); } else { std::cout << "Multi-reference mode: On" << std::endl; for (std::vector<StringRef>::const_iterator ref_file_split_sr_it = ref_file_split_sr.begin(); @@ -258,6 +270,7 @@ int main (int argc, char **argv) exit(-1); } } + EntityView refview=ref.GetChainList()[0].Select("peptide=true"); ref_list.push_back(ref.CreateFullView()); } glob_dist_list = CreateDistanceListFromMultipleReferences (ref_list,cutoffs,sequence_separation,radius); @@ -293,8 +306,8 @@ int main (int argc, char **argv) } continue; } - EntityView v=model.CreateFullView(); - + EntityView v=model.GetChainList()[0].Select("peptide=true"); + EntityView outv=model.GetChainList()[0].Select("peptide=true"); for (std::vector<EntityView>::const_iterator ref_list_it = ref_list.begin(); ref_list_it != ref_list.end(); ++ref_list_it) { bool cons_check = ResidueNamesMatch(v,*ref_list_it); @@ -355,8 +368,6 @@ int main (int argc, char **argv) std::cout << e.what() << std::endl; exit(-1); } - 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; try { v=alg::FilterClashes(v,nonbonded_table); } catch (std::exception& e) { @@ -364,8 +375,6 @@ int main (int argc, char **argv) std::cout << e.what() << std::endl; exit(-1); } - cov = compute_coverage(v,glob_dist_list); - std::cout << "Coverage after clashing checks: " << (float(cov.first)/float(cov.second)) << " (" << cov.first << " out of " << cov.second << " residues)" << std::endl; } if (cov.first==0) { std::cout << "Global LDDT score: 0.0" << std::endl; @@ -382,19 +391,59 @@ int main (int argc, char **argv) // prints the residue-by-residue statistics std::cout << "Local LDDT Score:" << std::endl; - std::cout << "Chain\tResName\tResNum\tScore\t(Conserved/Total, over " << cutoffs.size() << " thresholds)" << std::endl; - for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){ - ResidueView ritv = *rit; - Real lddt_local = 0; - int conserved_dist = 0; - int total_dist = 0; - if (ritv.HasProp(label)) { + if (structural_checks) { + std::cout << "Chain\tResName\tResNum\tAsses.\tQ.Prob.\tScore\t(Conserved/Total, over " << cutoffs.size() << " thresholds)" << std::endl; + } else { + std::cout << "Chain\tResName\tResNum\tAsses.\tScore\t(Conserved/Total, over " << cutoffs.size() << " thresholds)" << std::endl; + } + for (ResidueViewIter rit=outv.ResiduesBegin();rit!=outv.ResiduesEnd();++rit){ + ResidueView ritv=*rit; + ResNum rnum = ritv.GetNumber(); + bool assessed = false; + String assessed_string="No"; + bool quality_problems = false; + String quality_problems_string="No"; + Real lddt_local = -1; + String lddt_local_string="-"; + int conserved_dist = -1; + int total_dist = -1; + String dist_string = "-"; + if (is_resnum_in_globalrdmap(rnum,glob_dist_list)) { + assessed = true; + assessed_string="Yes"; + } + if (ritv.HasProp("stereo_chemical_violation_sidechain") || ritv.HasProp("steric_clash_sidechain")) { + quality_problems = true; + quality_problems_string="Yes"; + } + if (ritv.HasProp("stereo_chemical_violation_backbone") || ritv.HasProp("steric_clash_backbone")) { + quality_problems = true; + quality_problems_string="Yes+"; + } + + if (assessed==true) { + if (ritv.HasProp(label)) { lddt_local=ritv.GetFloatProp(label); + std::stringstream stkeylddt; + stkeylddt << std::fixed << std::setprecision(4) << lddt_local; + lddt_local_string=stkeylddt.str(); conserved_dist=ritv.GetIntProp(label+"_conserved"); - total_dist=ritv.GetIntProp(label+"_total"); + total_dist=ritv.GetIntProp(label+"_total"); + std::stringstream stkeydist; + stkeydist << "("<< conserved_dist << "/" << total_dist << ")"; + dist_string=stkeydist.str(); + } else { + lddt_local = 0; + lddt_local_string="0.0000"; + conserved_dist = 0; + total_dist = 0; + dist_string="(0/0)"; + } } - if (lddt_local!=0) { - std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << std::setprecision(4) << lddt_local << "\t" << "("<< conserved_dist << "/" << total_dist << ")" <<std::endl; + if (structural_checks) { + std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << quality_problems_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl; + } else { + std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl; } } std::cout << std::endl;