Improved output of lddt

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;
 }
 

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;