diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index f00fea65b459a1e4634ec4cee561dbe7b43a7ea8..deb570db8bd52444686c1a783d739b3f28509feb 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -36,7 +36,7 @@ void export_entity_to_density();
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_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;
diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index 4009cd3a86eb104a5e309a74f443b7465c9e29ab..9c2c1bf6121e6137e330ba03de6267753086fba3 100644
--- a/modules/mol/alg/src/lddt.cc
+++ b/modules/mol/alg/src/lddt.cc
@@ -16,6 +16,7 @@
// along with this library; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
+#include <iomanip>
#if defined (_MSC_VER)
#define BOOST_ALL_DYN_LINK 1
#endif
@@ -101,7 +102,7 @@ int main (int argc, char **argv)
String version = "Beta - 2012-06-13";
Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0;
- Real radius=15.0;
+ Real radius=15.0;
IOProfile profile;
profile.bond_feasibility_check=false;
@@ -296,19 +297,32 @@ int main (int argc, char **argv)
std::cout << "Global LDDT score: 0.0" << std::endl;
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 << "(" << 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 << "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){
ResidueView ritv = *rit;
Real lddt_local = 0;
- if (ritv.HasProp("locallddt")) {
- lddt_local=ritv.GetFloatProp("locallddt");
+ int conserved_dist = 0;
+ 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) {
- 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;
diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc
index 37b38d89a7d4a0d55ceac1609d1fc847964d3462..0511e71ef0fb397f135e3f18f86afab6c7653f6d 100644
--- a/modules/mol/alg/src/local_dist_diff_test.cc
+++ b/modules/mol/alg/src/local_dist_diff_test.cc
@@ -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,
ChainView mdl_chain,
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);
ResidueView mdl_res=mdl_chain.FindResidue(rnum);
@@ -226,7 +226,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
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];
XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
@@ -471,7 +471,7 @@ GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityVie
ref_counter++;
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;
- 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);
GlobalRDMap new_dist_list=CreateDistanceList(ref,max_dist);
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){
-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)
{
if (!mdl.GetResidueCount()) {
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) {
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);
ChainView mdl_chain=mdl.GetChainList()[0];
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) {
ResNum rn = i->first;
if (i->second.size()!=0) {
@@ -531,13 +531,15 @@ Real LocalDistDiffTest(const EntityView& mdl, const GlobalRDMap& glob_dist_list,
ResidueView mdlr=mdl_chain.FindResidue(rn);
if (mdlr.IsValid()) {
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.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();
- 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)
@@ -545,11 +547,11 @@ Real LocalDistDiffTest(const EntityView& mdl, const EntityView& target, Real cut
std::vector<Real> cutoffs;
cutoffs.push_back(cutoff);
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 cutoff, Real max_dist, int ref_index, int mdl_index)
{
@@ -597,8 +599,7 @@ Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln,
total_ov.first+=ov1.first;
total_ov.second+=ov1.second;
}
- return total_ov.first/(total_ov.second ? total_ov.second : 1);
- return 0.0;
+ return total_ov.first/(total_ov.second ? total_ov.second : 1);
}
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(4.0);
String label="locallddt";
- Real lddt=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, label);
- return lddt;
+ std::pair<long int,long int> total_ov=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, label);
+ 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)
@@ -618,16 +619,20 @@ Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
Real lddt =0;
std::vector<Real> cutoffs05;
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;
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;
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;
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;
for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
ResidueView ritv = *rit;
diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh
index 6e1c28a8088b1fade9918ce2c8d46653e9e6531f..7ed7e46586a6a0becd36c60682176f8a7323f440 100644
--- a/modules/mol/alg/src/local_dist_diff_test.hh
+++ b/modules/mol/alg/src/local_dist_diff_test.hh
@@ -55,7 +55,7 @@ typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<
typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
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,
std::vector<Real> cutoff_list,
const String& local_ldt_property_string="");