From 6ee5f25e044565a9283afdec6e8875520b88f26d Mon Sep 17 00:00:00 2001
From: Valerio Mariani <valerio.mariani@unibas.ch>
Date: Fri, 9 Nov 2012 13:40:03 +0100
Subject: [PATCH] Improved output of lddt
Conflicts:
modules/mol/alg/src/filter_clashes.cc
---
modules/mol/alg/src/filter_clashes.cc | 4 ++
modules/mol/alg/src/lddt.cc | 85 +++++++++++++++++++++------
2 files changed, 71 insertions(+), 18 deletions(-)
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index 396ebe39b..235464544 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -461,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) {
@@ -473,6 +474,7 @@ EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParam
filtered.AddAtom(atom);
}
}
+ res.SetBoolProp("stereo_chemical_violation_sidechain",true);
continue;
}
filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS);
@@ -576,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) {
@@ -588,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);
diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index c846af679..4d76ebd47 100644
--- a/modules/mol/alg/src/lddt.cc
+++ b/modules/mol/alg/src/lddt.cc
@@ -114,6 +114,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
@@ -244,8 +255,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();
@@ -261,6 +273,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);
@@ -297,8 +310,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);
@@ -363,8 +376,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) {
@@ -372,8 +383,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;
@@ -390,19 +399,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;
--
GitLab