diff --git a/modules/io/pymod/export_omf_io.cc b/modules/io/pymod/export_omf_io.cc
index fb98237b0ffbb2e5c6e5cd6d00ceceeefde080dd..43355a4a908f97ae3f567daba57f5f0e739c34c5 100644
--- a/modules/io/pymod/export_omf_io.cc
+++ b/modules/io/pymod/export_omf_io.cc
@@ -77,6 +77,7 @@ void export_omf_io() {
.def("GetChainNames", &wrap_get_chain_names)
.def("GetPositions", &OMF::GetPositions, return_value_policy<reference_existing_object>(),(arg("cname")))
.def("GetBFactors", &OMF::GetBFactors, return_value_policy<reference_existing_object>(),(arg("cname")))
+ .def("GetAvgBFactors", &OMF::GetAvgBFactors,(arg("cname")))
.def("GetSequence", &OMF::GetSequence, (arg("cname")))
;
}
diff --git a/modules/io/src/mol/omf.cc b/modules/io/src/mol/omf.cc
index 973595161ed6295631bf6e26593d9984b81b7d1e..91fdca413267c42bf41582e4422a37149bdbb59c 100644
--- a/modules/io/src/mol/omf.cc
+++ b/modules/io/src/mol/omf.cc
@@ -3833,6 +3833,31 @@ const std::vector<Real>& OMF::GetBFactors(const String& cname) const {
return it->second->bfactors;
}
+std::vector<Real> OMF::GetAvgBFactors(const String& cname) const {
+ auto it = chain_data_.find(cname);
+ if(it == chain_data_.end()) {
+ throw ost::Error("Provided chain name not in OMF structure");
+ }
+ const std::vector<Real>& bfactors = it->second->bfactors;
+ const std::vector<int>& res_def_indices = it->second->res_def_indices;
+ std::vector<Real> avg_bfactors;
+ avg_bfactors.reserve(it->second->res_def_indices.size());
+ int current_atom_idx = 0;
+ for(auto i = res_def_indices.begin(); i != res_def_indices.end(); ++i) {
+ int size = residue_definitions_[*i].anames.size();
+ Real summed_bfac = 0.0;
+ for(int j = 0; j < size; ++j) {
+ summed_bfac += bfactors[current_atom_idx];
+ ++current_atom_idx;
+ }
+ if(size > 0) {
+ summed_bfac /= size;
+ }
+ avg_bfactors.push_back(summed_bfac);
+ }
+ return avg_bfactors;
+}
+
String OMF::GetSequence(const String& cname) const {
auto it = chain_data_.find(cname);
if(it == chain_data_.end()) {
diff --git a/modules/io/src/mol/omf.hh b/modules/io/src/mol/omf.hh
index f5bb8b9966b0b5c74d9d114f745f212c16b52799..9503b1f541a928b1a3f3730a39fffc60bbab36c4 100644
--- a/modules/io/src/mol/omf.hh
+++ b/modules/io/src/mol/omf.hh
@@ -199,6 +199,8 @@ public:
const std::vector<Real>& GetBFactors(const String& cname) const;
+ std::vector<Real> GetAvgBFactors(const String& cname) const;
+
String GetSequence(const String& cname) const;
private:
diff --git a/modules/mol/alg/pymod/chain_mapping.py b/modules/mol/alg/pymod/chain_mapping.py
index c8bfa7c7789ca7440cdb84c55ce5ea6bf03709cb..01b0554a022328fb2b825192935a902b0af92acc 100644
--- a/modules/mol/alg/pymod/chain_mapping.py
+++ b/modules/mol/alg/pymod/chain_mapping.py
@@ -26,12 +26,14 @@ class MappingResult:
Constructor is directly called within the functions, no need to construct
such objects yourself.
"""
- def __init__(self, target, model, chem_groups, mapping, alns):
+ def __init__(self, target, model, chem_groups, mapping, alns,
+ opt_score=None):
self._target = target
self._model = model
self._chem_groups = chem_groups
self._mapping = mapping
self._alns = alns
+ self._opt_score = opt_score
@property
def target(self):
@@ -89,6 +91,18 @@ class MappingResult:
"""
return self._alns
+ @property
+ def opt_score(self):
+ """ Placeholder property without any guarantee of being set
+
+ Different scores get optimized in the various chain mapping algorithms.
+ Some of them may set their final optimal score in that property.
+ Consult the documentation of the respective chain mapping algorithm
+ for more information. Won't be in the return dict of
+ :func:`JSONSummary`.
+ """
+ return self._opt_score
+
def GetFlatMapping(self, mdl_as_key=False):
""" Returns flat mapping as :class:`dict` for all mapable chains
@@ -760,6 +774,9 @@ class ChainMapper:
scoring ones are extend by *block_seed_size* chains and the best
scoring one is exhaustively extended.
+ Sets :attr:`MappingResult.opt_score` in case of no trivial one-to-one
+ mapping.
+
:param model: Model to map
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param inclusion_radius: Inclusion radius for lDDT
@@ -827,11 +844,13 @@ class ChainMapper:
alns)
mapping = None
+ opt_lddt = None
if strategy == "naive":
- mapping = _lDDTNaive(self.target, mdl, inclusion_radius, thresholds,
- self.chem_groups, chem_mapping, ref_mdl_alns,
- self.n_max_naive)
+ mapping, opt_lddt = _lDDTNaive(self.target, mdl, inclusion_radius,
+ thresholds, self.chem_groups,
+ chem_mapping, ref_mdl_alns,
+ self.n_max_naive)
else:
# its one of the greedy strategies - setup greedy searcher
the_greed = _lDDTGreedySearcher(self.target, mdl, self.chem_groups,
@@ -846,6 +865,8 @@ class ChainMapper:
elif strategy == "greedy_block":
mapping = _lDDTGreedyBlock(the_greed, block_seed_size,
block_blocks_per_chem_group)
+ # cached => lDDT computation is fast here
+ opt_lddt = the_greed.lDDT(self.chem_groups, mapping)
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, mapping):
@@ -857,7 +878,7 @@ class ChainMapper:
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, mapping,
- alns)
+ alns, opt_score = opt_lddt)
def GetQSScoreMapping(self, model, contact_d = 12.0, strategy = "naive",
@@ -894,6 +915,9 @@ class ChainMapper:
scoring ones are extend by *block_seed_size* chains and the block with
with best QS score is exhaustively extended.
+ Sets :attr:`MappingResult.opt_score` in case of no trivial one-to-one
+ mapping.
+
:param model: Model to map
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param contact_d: Max distance between two residues to be considered as
@@ -935,18 +959,21 @@ class ChainMapper:
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, one_to_one,
alns)
+ mapping = None
+ opt_qsscore = None
if strategy == "naive":
- mapping = _QSScoreNaive(self.target, mdl, self.chem_groups,
- chem_mapping, ref_mdl_alns, contact_d,
- self.n_max_naive)
+ mapping, opt_qsscore = _QSScoreNaive(self.target, mdl,
+ self.chem_groups,
+ chem_mapping, ref_mdl_alns,
+ contact_d, self.n_max_naive)
else:
# its one of the greedy strategies - setup greedy searcher
-
- the_greed = _QSScoreGreedySearcher(self.target, mdl, self.chem_groups,
- chem_mapping, ref_mdl_alns,
- contact_d = contact_d,
- steep_opt_rate=steep_opt_rate)
+ the_greed = _QSScoreGreedySearcher(self.target, mdl,
+ self.chem_groups,
+ chem_mapping, ref_mdl_alns,
+ contact_d = contact_d,
+ steep_opt_rate=steep_opt_rate)
if strategy == "greedy_fast":
mapping = _QSScoreGreedyFast(the_greed)
elif strategy == "greedy_full":
@@ -954,6 +981,9 @@ class ChainMapper:
elif strategy == "greedy_block":
mapping = _QSScoreGreedyBlock(the_greed, block_seed_size,
block_blocks_per_chem_group)
+ # cached => QSScore computation is fast here
+ opt_qsscore = the_greed.Score(mapping, check=False)
+
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, mapping):
@@ -965,7 +995,7 @@ class ChainMapper:
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, mapping,
- alns)
+ alns, opt_score = opt_qsscore)
def GetRigidMapping(self, model, strategy = "greedy_single_gdtts",
single_chain_gdtts_thresh=0.4, subsampling=None,
@@ -2225,7 +2255,7 @@ def _lDDTNaive(trg, mdl, inclusion_radius, thresholds, chem_groups,
best_mapping = mapping
best_lddt = lDDT
- return best_mapping
+ return (best_mapping, best_lddt)
def _lDDTGreedyFast(the_greed):
@@ -2648,7 +2678,7 @@ def _QSScoreNaive(trg, mdl, chem_groups, chem_mapping, ref_mdl_alns, contact_d,
if score_result.QS_global > best_score:
best_mapping = mapping
best_score = score_result.QS_global
- return best_mapping
+ return (best_mapping, best_score)
def _QSScoreGreedyFast(the_greed):