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

Lddt now reports exact count of conserved and toatal distances 

Even residue by residue!
parent aacdb58b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
modules/mol/alg/pymod/wrap_mol_alg.cc
+ 1
− 1
View file @ 72cb7cad
...@@ -36,7 +36,7 @@ void export_entity_to_density(); ...@@ -36,7 +36,7 @@ void export_entity_to_density();
namespace { namespace {
Real (*lddt_a)(const mol::EntityView&, const mol::alg::GlobalRDMap& , std::vector<Real>, const String&)=&mol::alg::LocalDistDiffTest; std::pair<long int,long int> (*lddt_a)(const mol::EntityView&, const mol::alg::GlobalRDMap& , std::vector<Real>, const String&)=&mol::alg::LocalDistDiffTest;
Real (*lddt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistDiffTest; Real (*lddt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistDiffTest;
Real (*lddt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistDiffTest; Real (*lddt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistDiffTest;
mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes; mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes;
......
modules/mol/alg/src/lddt.cc
+ 20
− 6
View file @ 72cb7cad
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
// along with this library; if not, write to the Free Software Foundation, Inc., // along with this library; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
#include <iomanip>
#if defined (_MSC_VER) #if defined (_MSC_VER)
#define BOOST_ALL_DYN_LINK 1 #define BOOST_ALL_DYN_LINK 1
#endif #endif
...@@ -101,7 +102,7 @@ int main (int argc, char **argv) ...@@ -101,7 +102,7 @@ int main (int argc, char **argv)
String version = "Beta - 2012-06-13"; String version = "Beta - 2012-06-13";
Real bond_tolerance = 8.0; Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0; Real angle_tolerance = 8.0;
Real radius=15.0; Real radius=15.0;
IOProfile profile; IOProfile profile;
profile.bond_feasibility_check=false; profile.bond_feasibility_check=false;
...@@ -296,19 +297,32 @@ int main (int argc, char **argv) ...@@ -296,19 +297,32 @@ int main (int argc, char **argv)
std::cout << "Global LDDT score: 0.0" << std::endl; std::cout << "Global LDDT score: 0.0" << std::endl;
return 0; return 0;
} }
Real lddt=LDDTHA(v, glob_dist_list);
std::vector<Real> cutoffs;
cutoffs.push_back(0.5);
cutoffs.push_back(1.0);
cutoffs.push_back(2.0);
cutoffs.push_back(4.0);
String label="localldt";
std::pair<int,int> total_ov=alg::LocalDistDiffTest(v, glob_dist_list, cutoffs, label);
Real lddt = static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
std::cout << "Global LDDT score: " << lddt << std::endl; std::cout << "Global LDDT score: " << lddt << std::endl;
std::cout << "(" << std::fixed << total_ov.first << " conserved distances in the model out of " << total_ov.second << " checked)" << std::endl;
std::cout << "Local LDDT Score:" << std::endl; std::cout << "Local LDDT Score:" << std::endl;
std::cout << "Chain\tResName\tResNum\tScore" << std::endl; std::cout << "Chain\tResName\tResNum\tScore\t(Conserved/Total)" << std::endl;
for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){ for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
ResidueView ritv = *rit; ResidueView ritv = *rit;
Real lddt_local = 0; Real lddt_local = 0;
if (ritv.HasProp("locallddt")) { int conserved_dist = 0;
lddt_local=ritv.GetFloatProp("locallddt"); int total_dist = 0;
if (ritv.HasProp(label)) {
lddt_local=ritv.GetFloatProp(label);
conserved_dist=ritv.GetIntProp(label+"_conserved");
total_dist=ritv.GetIntProp(label+"_total");
} }
if (lddt_local!=0) { if (lddt_local!=0) {
std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << lddt_local << std::endl; std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << lddt_local << "\t" << "("<< conserved_dist << "/" << total_dist << ")" <<std::endl;
} }
} }
std::cout << std::endl; std::cout << std::endl;
......
modules/mol/alg/src/local_dist_diff_test.cc
+ 25
− 20
View file @ 72cb7cad
...@@ -78,7 +78,7 @@ std::pair<bool,Real> within_tolerance(Real mdl_dist, std::pair<Real,Real>& value ...@@ -78,7 +78,7 @@ std::pair<bool,Real> within_tolerance(Real mdl_dist, std::pair<Real,Real>& value
std::pair<Real, Real> calc_overlap1(const ResidueRDMap& res_distance_list, const ResNum& rnum, std::pair<Real, Real> calc_overlap1(const ResidueRDMap& res_distance_list, const ResNum& rnum,
ChainView mdl_chain, ChainView mdl_chain,
std::vector<Real>& tol_list, bool only_fixed, std::vector<Real>& tol_list, bool only_fixed,
bool swap,std::vector<std::pair<Real, Real> >& overlap_list, bool log ) bool swap,std::vector<std::pair<long int, long int> >& overlap_list, bool log )
{ {
std::pair<Real, Real> overlap(0.0, 0.0); std::pair<Real, Real> overlap(0.0, 0.0);
ResidueView mdl_res=mdl_chain.FindResidue(rnum); ResidueView mdl_res=mdl_chain.FindResidue(rnum);
...@@ -226,7 +226,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq, ...@@ -226,7 +226,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
return overlap; return overlap;
} }
void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, std::vector<Real> cutoff_list, std::vector<std::pair<Real, Real> > overlap_list) void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, std::vector<Real> cutoff_list, std::vector<std::pair<long int, long int> > overlap_list)
{ {
ChainView mdl_chain=mdl.GetChainList()[0]; ChainView mdl_chain=mdl.GetChainList()[0];
XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT); XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
...@@ -471,7 +471,7 @@ GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityVie ...@@ -471,7 +471,7 @@ GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityVie
ref_counter++; ref_counter++;
for (std::vector<EntityView>::const_iterator ref_list_it=ref_list.begin()+1;ref_list_it!=ref_list.end();++ref_list_it) { for (std::vector<EntityView>::const_iterator ref_list_it=ref_list.begin()+1;ref_list_it!=ref_list.end();++ref_list_it) {
EntityView ref = *ref_list_it; EntityView ref = *ref_list_it;
std::vector<std::pair<Real, Real> > overlap_list(ref.GetResidueCount(), std::pair<Real, Real>(0.0, 0.0)); std::vector<std::pair<long int, long int> > overlap_list(ref.GetResidueCount(), std::pair<long int, long int>(0, 0));
check_and_swap(glob_dist_list,ref,cutoff_list,overlap_list); check_and_swap(glob_dist_list,ref,cutoff_list,overlap_list);
GlobalRDMap new_dist_list=CreateDistanceList(ref,max_dist); GlobalRDMap new_dist_list=CreateDistanceList(ref,max_dist);
merge_distance_lists(glob_dist_list,new_dist_list,ex_map,ref,ref_counter); merge_distance_lists(glob_dist_list,new_dist_list,ex_map,ref,ref_counter);
...@@ -503,22 +503,22 @@ void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list){ ...@@ -503,22 +503,22 @@ void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list){
Real LocalDistDiffTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list, std::pair<long int,long int> LocalDistDiffTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list,
std::vector<Real> cutoff_list, const String& local_lddt_property_string) std::vector<Real> cutoff_list, const String& local_lddt_property_string)
{ {
if (!mdl.GetResidueCount()) { if (!mdl.GetResidueCount()) {
LOG_WARNING("model structures doesn't contain any residues"); LOG_WARNING("model structures doesn't contain any residues");
return 0.0; return std::make_pair<long int,long int>(0,0);
} }
if (glob_dist_list.size()==0) { if (glob_dist_list.size()==0) {
LOG_WARNING("global reference list is empty"); LOG_WARNING("global reference list is empty");
return 0.0; return std::make_pair<long int,long int>(0,0);
} }
std::vector<std::pair<Real, Real> > overlap_list(mdl.GetResidueCount(), std::pair<Real, Real>(0.0, 0.0)); std::vector<std::pair<long int, long int> > overlap_list(mdl.GetResidueCount(), std::pair<long int, long int>(0, 0));
check_and_swap(glob_dist_list,mdl,cutoff_list,overlap_list); check_and_swap(glob_dist_list,mdl,cutoff_list,overlap_list);
ChainView mdl_chain=mdl.GetChainList()[0]; ChainView mdl_chain=mdl.GetChainList()[0];
overlap_list.clear(); overlap_list.clear();
std::pair<Real, Real> total_ov(0.0, 0.0); std::pair<long int, long int> total_ov(0, 0);
for (GlobalRDMap::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) { for (GlobalRDMap::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
ResNum rn = i->first; ResNum rn = i->first;
if (i->second.size()!=0) { if (i->second.size()!=0) {
...@@ -531,13 +531,15 @@ Real LocalDistDiffTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list, ...@@ -531,13 +531,15 @@ Real LocalDistDiffTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list,
ResidueView mdlr=mdl_chain.FindResidue(rn); ResidueView mdlr=mdl_chain.FindResidue(rn);
if (mdlr.IsValid()) { if (mdlr.IsValid()) {
int mdl_res_index =mdlr.GetIndex(); int mdl_res_index =mdlr.GetIndex();
Real local_lddt=overlap_list[mdl_res_index].first/(overlap_list[mdl_res_index].second ? overlap_list[mdl_res_index].second : 1); Real local_lddt=static_cast<Real>(overlap_list[mdl_res_index].first)/(static_cast<Real>(overlap_list[mdl_res_index].second) ? static_cast<Real>(overlap_list[mdl_res_index].second) : 1);
mdlr.SetFloatProp(local_lddt_property_string, local_lddt); mdlr.SetFloatProp(local_lddt_property_string, local_lddt);
mdlr.SetIntProp(local_lddt_property_string+"_conserved", overlap_list[mdl_res_index].first);
mdlr.SetIntProp(local_lddt_property_string+"_total", overlap_list[mdl_res_index].second);
} }
} }
} }
overlap_list.clear(); overlap_list.clear();
return total_ov.first/(total_ov.second ? total_ov.second : 1); return std::make_pair<long int,long int>(total_ov.first,total_ov.second);
} }
Real LocalDistDiffTest(const EntityView& mdl, const EntityView& target, Real cutoff, Real max_dist, const String& local_lddt_property_string) Real LocalDistDiffTest(const EntityView& mdl, const EntityView& target, Real cutoff, Real max_dist, const String& local_lddt_property_string)
...@@ -545,11 +547,11 @@ Real LocalDistDiffTest(const EntityView& mdl, const EntityView& target, Real cut ...@@ -545,11 +547,11 @@ Real LocalDistDiffTest(const EntityView& mdl, const EntityView& target, Real cut
std::vector<Real> cutoffs; std::vector<Real> cutoffs;
cutoffs.push_back(cutoff); cutoffs.push_back(cutoff);
GlobalRDMap glob_dist_list = CreateDistanceList(target,max_dist); GlobalRDMap glob_dist_list = CreateDistanceList(target,max_dist);
return LocalDistDiffTest(mdl, glob_dist_list, cutoffs, local_lddt_property_string); std::pair<long int,long int> total_ov = LocalDistDiffTest(mdl, glob_dist_list, cutoffs, local_lddt_property_string);
return static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
} }
Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln, Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln,
Real cutoff, Real max_dist, int ref_index, int mdl_index) Real cutoff, Real max_dist, int ref_index, int mdl_index)
{ {
...@@ -597,8 +599,7 @@ Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln, ...@@ -597,8 +599,7 @@ Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln,
total_ov.first+=ov1.first; total_ov.first+=ov1.first;
total_ov.second+=ov1.second; total_ov.second+=ov1.second;
} }
return total_ov.first/(total_ov.second ? total_ov.second : 1); return total_ov.first/(total_ov.second ? total_ov.second : 1);
return 0.0;
} }
Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list) Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
...@@ -609,8 +610,8 @@ Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list) ...@@ -609,8 +610,8 @@ Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
cutoffs.push_back(2.0); cutoffs.push_back(2.0);
cutoffs.push_back(4.0); cutoffs.push_back(4.0);
String label="locallddt"; String label="locallddt";
Real lddt=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, label); std::pair<long int,long int> total_ov=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, label);
return lddt; return static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
} }
Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list) Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
...@@ -618,16 +619,20 @@ Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list) ...@@ -618,16 +619,20 @@ Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
Real lddt =0; Real lddt =0;
std::vector<Real> cutoffs05; std::vector<Real> cutoffs05;
cutoffs05.push_back(0.5); cutoffs05.push_back(0.5);
Real lddt05=alg::LocalDistDiffTest(v, global_dist_list, cutoffs05, "locallddt0.5"); std::pair<Real,Real> lddt05o=alg::LocalDistDiffTest(v, global_dist_list, cutoffs05, "locallddt0.5");
std::vector<Real> cutoffs1; std::vector<Real> cutoffs1;
cutoffs1.push_back(1.0); cutoffs1.push_back(1.0);
Real lddt1=alg::LocalDistDiffTest(v, global_dist_list, cutoffs1, "locallddt1"); std::pair<Real,Real> lddt1o=alg::LocalDistDiffTest(v, global_dist_list, cutoffs1, "locallddt1");
std::vector<Real> cutoffs2; std::vector<Real> cutoffs2;
cutoffs2.push_back(2.0); cutoffs2.push_back(2.0);
Real lddt2=alg::LocalDistDiffTest(v, global_dist_list, cutoffs2, "locallddt2"); std::pair<Real,Real> lddt2o=alg::LocalDistDiffTest(v, global_dist_list, cutoffs2, "locallddt2");
std::vector<Real> cutoffs4; std::vector<Real> cutoffs4;
cutoffs4.push_back(4.0); cutoffs4.push_back(4.0);
Real lddt4=alg::LocalDistDiffTest(v, global_dist_list, cutoffs4, "locallddt4"); std::pair<Real,Real> lddt4o=alg::LocalDistDiffTest(v, global_dist_list, cutoffs4, "locallddt4");
Real lddt05 = lddt05o.first/(lddt05o.second ? lddt05o.second : 1);
Real lddt1 = lddt1o.first/(lddt1o.second ? lddt1o.second : 1);
Real lddt2 = lddt2o.first/(lddt2o.second ? lddt2o.second : 1);
Real lddt4 = lddt4o.first/(lddt4o.second ? lddt4o.second : 1);
lddt = (lddt05+lddt1+lddt2+lddt4)/4.0; lddt = (lddt05+lddt1+lddt2+lddt4)/4.0;
for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){ for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
ResidueView ritv = *rit; ResidueView ritv = *rit;
......
modules/mol/alg/src/local_dist_diff_test.hh
+ 1
− 1
View file @ 72cb7cad
...@@ -55,7 +55,7 @@ typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair< ...@@ -55,7 +55,7 @@ typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<
typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap; typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
typedef std::map<UniqueAtomIdentifier,int> ExistenceMap; typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
Real DLLEXPORT_OST_MOL_ALG LocalDistDiffTest(const EntityView& mdl, std::pair<long int,long int> DLLEXPORT_OST_MOL_ALG LocalDistDiffTest(const EntityView& mdl,
const GlobalRDMap& dist_list, const GlobalRDMap& dist_list,
std::vector<Real> cutoff_list, std::vector<Real> cutoff_list,
const String& local_ldt_property_string=""); const String& local_ldt_property_string="");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment