diff --git a/modules/seq/alg/src/variance_map.cc b/modules/seq/alg/src/variance_map.cc index 05d9eb04b1a814feeb8f961862975cdda725f0b3..38a27051e55f463e1f258b264deeade192ad5099 100644 --- a/modules/seq/alg/src/variance_map.cc +++ b/modules/seq/alg/src/variance_map.cc @@ -72,15 +72,16 @@ String GetJson(const T& data, uint num_rows, uint num_cols) { } -void FilllDDTValues(int len, const std::vector<Real>& dist_diff, - const std::vector<int>& below_15, +void FilllDDTValues(const std::vector<Real>& dist_diff, + const std::vector<int>& ref_distances, const std::vector<std::pair<int, int> >& idx_mapping, std::vector<Real>& local_lddt) { + int len = dist_diff.size(); std::vector<int> counts(len, 0); std::vector<int> total_counts(len, 0); - for(auto it = below_15.begin(); it != below_15.end(); ++it) { + for(auto it = ref_distances.begin(); it != ref_distances.end(); ++it) { Real diff = dist_diff[*it]; int fulfilled = 0; fulfilled += static_cast<int>(diff < static_cast<Real>(4.0)); @@ -106,7 +107,7 @@ void FilllDDTValues(int len, const std::vector<Real>& dist_diff, } -void FilllDDTValues(int len, const std::vector<Real>& dist_one, +void FilllDDTValues(const std::vector<Real>& dist_one, const std::vector<Real>& dist_two, const std::vector<int>& below_15_one, const std::vector<int>& below_15_two, @@ -114,13 +115,16 @@ void FilllDDTValues(int len, const std::vector<Real>& dist_one, std::vector<Real>& local_lddt_one, std::vector<Real>& local_lddt_two) { + // absolute dist differences are the same for both std::vector<Real> dist_diff(dist_one.size()); for(uint i = 0; i < dist_one.size(); ++i) { dist_diff[i] = std::abs(dist_one[i] - dist_two[i]); } - FilllDDTValues(len, dist_diff, below_15_two, idx_mapping, local_lddt_one); - FilllDDTValues(len, dist_diff, below_15_one, idx_mapping, local_lddt_two); + // estimate lDDT of structure 1 with reference distances (below 15) from + // structure 2 and vice versa + FilllDDTValues(dist_diff, below_15_two, idx_mapping, local_lddt_one); + FilllDDTValues(dist_diff, below_15_one, idx_mapping, local_lddt_two); } } // anon ns @@ -289,17 +293,20 @@ MeanlDDTPtr CreateMeanlDDTHA(const DistanceMapPtr dmap) { throw IntegrityError("empty distance map provided"); } - // the relevant pairwise distances are the ones that are off-diagonal from the - // pairwise distance matrices + // the relevant pairwise distances are the off-diagonal elements from the + // pairwise distance matrices -> the upper half without the diagonal elements uint n_off_diagonal = (len*len-len)/2; + // store that information in a linear layout for each structure std::vector<std::vector<Real> > dist_data(dmap->GetNumStructures()); + // keep track of distances <= 15 for lDDT calculation. Indices stored in here + // relate to the linear representations in dist_data std::vector<std::vector<int> > below_15(dmap->GetNumStructures()); for(uint s_idx = 0; s_idx < nstruc; ++s_idx) { - // the default value is simply a very large number - // this will trigger a distance difference above any relevant threshold + // the default value is simply a very large number this will trigger a + // distance difference above any relevant threshold if not set dist_data[s_idx].assign(n_off_diagonal, 10000000.0); } - + // keep track of which two residues are involved for each element in dist_data std::vector<std::pair<int,int> > off_diagonal_mapper(n_off_diagonal); uint off_diag_idx = 0; @@ -307,9 +314,8 @@ MeanlDDTPtr CreateMeanlDDTHA(const DistanceMapPtr dmap) { for (uint j_res = i_res + 1; j_res < len; ++j_res, ++off_diag_idx) { off_diagonal_mapper[off_diag_idx] = std::make_pair(i_res, j_res); const Distances& di = dmap->GetDistances(i_res, j_res); - const uint size = di.GetDataSize(); - for (uint k = 0; k < size; ++k) { - const std::pair<Real, int> ret = di.GetDataElement(k); + for (uint k = 0; k < di.GetDataSize(); ++k) { + const std::pair<Real, int>& ret = di.GetDataElement(k); dist_data[ret.second-1][off_diag_idx] = ret.first; if(ret.first <= 15.0) { below_15[ret.second-1].push_back(off_diag_idx); @@ -318,27 +324,26 @@ MeanlDDTPtr CreateMeanlDDTHA(const DistanceMapPtr dmap) { } } - std::vector<std::vector<Real> > avg_lddt_values(nstruc); - std::vector<std::vector<int> > avg_lddt_counts(nstruc); - + std::vector<std::vector<Real> > values(nstruc); + std::vector<std::vector<int> > counts(nstruc); for(uint i = 0; i < nstruc; ++i) { - avg_lddt_values[i].assign(len, 0); - avg_lddt_counts[i].assign(len, 0); + values[i].assign(len, 0); + counts[i].assign(len, 0); } std::vector<Real> lddt_i, lddt_j; for(uint i = 0; i < nstruc; ++i) { for(uint j = i+1; j < nstruc; ++j) { - FilllDDTValues(len, dist_data[i], dist_data[j], below_15[i], below_15[j], + FilllDDTValues(dist_data[i], dist_data[j], below_15[i], below_15[j], off_diagonal_mapper, lddt_i, lddt_j); for(uint k = 0; k < len; ++k) { if(!std::isnan(lddt_i[k])) { - avg_lddt_values[i][k] += lddt_i[k]; - avg_lddt_counts[i][k] += 1; + values[i][k] += lddt_i[k]; + counts[i][k] += 1; } if(!std::isnan(lddt_j[k])) { - avg_lddt_values[j][k] += lddt_j[k]; - avg_lddt_counts[j][k] += 1; + values[j][k] += lddt_j[k]; + counts[j][k] += 1; } } } @@ -347,9 +352,9 @@ MeanlDDTPtr CreateMeanlDDTHA(const DistanceMapPtr dmap) { MeanlDDTPtr return_ptr(new MeanlDDT(len, nstruc)); for(uint struc_idx = 0; struc_idx < nstruc; ++struc_idx) { for(uint res_idx = 0; res_idx < len; ++res_idx) { - if(avg_lddt_counts[struc_idx][res_idx] > 0) { + if(counts[struc_idx][res_idx] > 0) { return_ptr->Set(res_idx, struc_idx, - avg_lddt_values[struc_idx][res_idx]/avg_lddt_counts[struc_idx][res_idx]); + values[struc_idx][res_idx]/counts[struc_idx][res_idx]); } } }