diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index c4f2bfe8f708cfa38e249e1155679fbf91bfca36..a2d04c35d1f9ccc5650888a018e801392e5bd7f1 100644
--- a/modules/mol/alg/src/lddt.cc
+++ b/modules/mol/alg/src/lddt.cc
@@ -53,6 +53,9 @@ EntityHandle load(const String& file, const IOProfile& profile)
reader.Import(ent);
conop::Conopology& conop_inst=conop::Conopology::Instance();
conop_inst.ConnectAll(conop_inst.GetBuilder(), ent);
+ if (ent.GetChainList().size()!=1) {
+ std::cout << "WARNING: File " << file << "has more than one chain" << std::endl;
+ }
return ent;
}
std::cerr << "ERROR: '" << file << "' does not contain any ATOM records. "
@@ -78,6 +81,7 @@ void usage()
std::cerr << " -a <value> tolerance in stddevs for angles" << std::endl;
std::cerr << " -m <value> clashing distance for unknwon atom types" << std::endl;
std::cerr << " -r <value> distance inclusion radius" << std::endl;
+ std::cerr << " -i <value> sequence separation" << std::endl;
std::cerr << " -e print version" << std::endl;
}
@@ -108,6 +112,7 @@ int main (int argc, char **argv)
Real bond_tolerance = 8.0;
Real angle_tolerance = 8.0;
Real radius=15.0;
+ int sequence_separation = 0;
// creates the required loading profile
IOProfile profile;
@@ -128,6 +133,7 @@ int main (int argc, char **argv)
("bond_tolerance,b", po::value<Real>(), "tolerance in stddev for bonds")
("angle_tolerance,a", po::value<Real>(), "tolerance in stddev for angles")
("inclusion_radius,r", po::value<Real>(), "distance inclusion radius")
+ ("sequence_separation,i", po::value<int>(), "sequence separation")
("files", po::value< std::vector<String> >(), "input file(s)")
("reference",po::value<String>(),"reference(s)")
;
@@ -197,7 +203,10 @@ int main (int argc, char **argv)
if (vm.count("inclusion_radius")) {
radius=vm["inclusion_radius"].as<Real>();
}
-
+ if (vm.count("sequence_separation")) {
+ sequence_separation=vm["sequence_separation"].as<int>();
+ }
+
std::vector<Real> cutoffs;
cutoffs.push_back(0.5);
cutoffs.push_back(1.0);
@@ -231,7 +240,7 @@ int main (int argc, char **argv)
}
ref_list.push_back(ref.CreateFullView());
}
- glob_dist_list = CreateDistanceListFromMultipleReferences (ref_list,cutoffs,radius);
+ glob_dist_list = CreateDistanceListFromMultipleReferences (ref_list,cutoffs,sequence_separation,radius);
}
files.pop_back();
@@ -243,6 +252,7 @@ int main (int argc, char **argv)
std::cout << "Stereo-chemical and steric clash checks: Off " << std::endl;
}
std::cout << "Inclusion Radius: " << radius << std::endl;
+ std::cout << "Sequence separation: " << sequence_separation << std::endl;
if (structural_checks) {
std::cout << "Parameter filename: " << parameter_filename << std::endl;
std::cout << "Tolerance in stddevs for bonds: " << bond_tolerance << std::endl;
@@ -334,7 +344,7 @@ int main (int argc, char **argv)
// computes the lddt score
String label="localldt";
- std::pair<int,int> total_ov=alg::LocalDistDiffTest(v, glob_dist_list, cutoffs, label);
+ std::pair<int,int> total_ov=alg::LocalDistDiffTest(v, glob_dist_list, cutoffs, sequence_separation, 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;
diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc
index e69978fbc80e5d41d7f5a8d75bcd74bf1f37cd57..fdb1bc0a3e943c716bb82b21083465992f3b1b2f 100644
--- a/modules/mol/alg/src/local_dist_diff_test.cc
+++ b/modules/mol/alg/src/local_dist_diff_test.cc
@@ -79,7 +79,7 @@ std::pair<bool,Real> within_tolerance(Real mdl_dist, std::pair<Real,Real>& value
// helper function
std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_list, const ResNum& rnum,
- ChainView mdl_chain,
+ ChainView mdl_chain, int sequence_separation,
std::vector<Real>& tol_list, bool only_fixed,
bool swap,std::vector<std::pair<long int, long int> >& overlap_list, bool log )
{
@@ -94,7 +94,7 @@ std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_lis
AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
if (only_fixed) {
- if (first_atom.GetResNum()==second_atom.GetResNum()) {
+ if (std::abs(first_atom.GetResNum().GetNum()-second_atom.GetResNum().GetNum())<sequence_separation) {
continue;
}
if (swappable(second_atom.GetResidueName(), second_atom.GetAtomName())) {
@@ -102,7 +102,7 @@ std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_lis
}
}
if (!only_fixed) {
- if (first_atom.GetResNum()<=second_atom.GetResNum()) {
+ if (first_atom.GetResNum().GetNum()<=(second_atom.GetResNum().GetNum()+sequence_separation)) {
continue;
}
}
@@ -232,7 +232,7 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
// for each residue with multiple possible nomenclature conventions, checks which choice (switched or not)
// of atom nomenclature gives the highest lddt score then changes the naming convention of the input
// entity view accordingly
-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)
+void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, std::vector<Real> cutoff_list, int sequence_separation, std::vector<std::pair<long int, long int> > overlap_list)
{
ChainView mdl_chain=mdl.GetChainList()[0];
XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
@@ -249,11 +249,11 @@ void check_and_swap(const GlobalRDMap& glob_dist_list, const EntityView& mdl, st
if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || rname=="PHE" || rname=="LYS" || rname=="ARG")) {
continue;
}
- std::pair<Real, Real> ov1=calc_overlap1(i->second, rnum,mdl_chain,
+ std::pair<Real, Real> ov1=calc_overlap1(i->second, rnum,mdl_chain, sequence_separation,
cutoff_list, true,
false, overlap_list,false);
- std::pair<Real, Real> ov2=calc_overlap1(i->second, rnum, mdl_chain,
+ std::pair<Real, Real> ov2=calc_overlap1(i->second, rnum, mdl_chain, sequence_separation,
cutoff_list, true,
true, overlap_list,false);
@@ -472,7 +472,7 @@ GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist)
return dist_list;
}
-GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list, std::vector<Real>& cutoff_list, Real max_dist)
+GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list, std::vector<Real>& cutoff_list, int sequence_separation, Real max_dist)
{
int ref_counter=0;
ExistenceMap ex_map;
@@ -482,7 +482,7 @@ GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityVie
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<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,sequence_separation,overlap_list);
GlobalRDMap new_dist_list=CreateDistanceList(ref,max_dist);
merge_distance_lists(glob_dist_list,new_dist_list,ex_map,ref,ref_counter);
update_existence_map (ex_map,ref,ref_counter);
@@ -514,7 +514,7 @@ void PrintGlobalRDMap(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, int sequence_separation, const String& local_lddt_property_string)
{
if (!mdl.GetResidueCount()) {
LOG_WARNING("model structures doesn't contain any residues");
@@ -525,14 +525,14 @@ std::pair<long int,long int> LocalDistDiffTest(const EntityView& mdl, const Glob
return std::make_pair<long int,long int>(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,sequence_separation,overlap_list);
ChainView mdl_chain=mdl.GetChainList()[0];
overlap_list.clear();
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) {
- std::pair<long int, long int> ov1=calc_overlap1(i->second, rn, mdl_chain, cutoff_list,
+ std::pair<long int, long int> ov1=calc_overlap1(i->second, rn, mdl_chain, sequence_separation, cutoff_list,
false, false, overlap_list,true);
total_ov.first+=ov1.first;
total_ov.second+=ov1.second;
@@ -557,7 +557,8 @@ 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);
- std::pair<long int,long int> total_ov = LocalDistDiffTest(mdl, glob_dist_list, cutoffs, local_lddt_property_string);
+ int sequence_separation = 0;
+ std::pair<long int,long int> total_ov = LocalDistDiffTest(mdl, glob_dist_list, cutoffs, sequence_separation, local_lddt_property_string);
return static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
}
@@ -612,7 +613,7 @@ Real LocalDistDiffTest(const ost::seq::AlignmentHandle& aln,
return total_ov.first/(total_ov.second ? total_ov.second : 1);
}
-Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
+Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list, int sequence_separation)
{
std::vector<Real> cutoffs;
cutoffs.push_back(0.5);
@@ -620,10 +621,11 @@ Real LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
cutoffs.push_back(2.0);
cutoffs.push_back(4.0);
String label="locallddt";
- std::pair<long int,long int> total_ov=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, label);
+ std::pair<long int,long int> total_ov=alg::LocalDistDiffTest(v, global_dist_list, cutoffs, sequence_separation, label);
return static_cast<Real>(total_ov.first)/(static_cast<Real>(total_ov.second) ? static_cast<Real>(total_ov.second) : 1);
}
+// debugging code
/*
Real OldStyleLDDTHA(EntityView& v, const GlobalRDMap& global_dist_list)
{
diff --git a/modules/mol/alg/src/local_dist_diff_test.hh b/modules/mol/alg/src/local_dist_diff_test.hh
index 894d54605d4b4beccd1b5c18e54959b085790605..d60c4447b0c81b20ac38eea0bb71b2eb718f9146 100644
--- a/modules/mol/alg/src/local_dist_diff_test.hh
+++ b/modules/mol/alg/src/local_dist_diff_test.hh
@@ -95,7 +95,8 @@ typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
/// int properties named <string>_conserved and <string>_total.
std::pair<long int,long int> DLLEXPORT_OST_MOL_ALG LocalDistDiffTest(const EntityView& mdl,
const GlobalRDMap& dist_list,
- std::vector<Real> cutoff_list,
+ std::vector<Real> cutoff_list,
+ int sequence_separation = 0,
const String& local_ldt_property_string="");
/// \brief Calculates the Local Distance Difference Score for a given model with respect to a given target
@@ -131,7 +132,7 @@ Real DLLEXPORT_OST_MOL_ALG LocalDistDiffTest(const ost::seq::AlignmentHandle& al
///
/// Computes the Local Distance Difference High-Accuracy Test (with threshold 0.5,1,2 and 4 Angstrom)
/// Requires a list of distances to check and a model for which the score is computed
-Real DLLEXPORT_OST_MOL_ALG LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list);
+Real DLLEXPORT_OST_MOL_ALG LDDTHA(EntityView& v, const GlobalRDMap& global_dist_list, int sequence_separation=0);
/// \brief Creates a list of distances to check during a Local Difference Distance Test
///
@@ -144,7 +145,7 @@ GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist);
/// than the inclusion radius in all structures in which the two atoms are present, it will be included in the list
/// However, if the distance is longer than the inclusion radius in at least one of the structures, it
/// will not be considered a local interaction and will be exluded from the list
-GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list,std::vector<Real>& cutoff_list, Real max_dist);
+GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list,std::vector<Real>& cutoff_list, int sequence_separation, Real max_dist);
/// \brief Prints all distances in a global distance list to standard output
void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list);