diff --git a/modules/mol/alg/pymod/ligand_scoring_base.py b/modules/mol/alg/pymod/ligand_scoring_base.py
index 17a1ad00aeae5fa55eea1d7372b3376700440970..4cd4bb588bbe44dd0eb8473cd80c200a49f13fa8 100644
--- a/modules/mol/alg/pymod/ligand_scoring_base.py
+++ b/modules/mol/alg/pymod/ligand_scoring_base.py
@@ -6,7 +6,7 @@ from ost import LogWarning, LogScript, LogVerbose, LogDebug
from ost.mol.alg import chain_mapping
class LigandScorer:
- """ Scorer class to compute various small molecule ligand (non polymer) scores.
+ """ Scorer to compute various small molecule ligand (non polymer) scores.
.. note ::
Extra requirements:
@@ -16,7 +16,8 @@ class LigandScorer:
:class:`LigandScorer` is an abstract base class dealing with all the setup,
data storage, enumerating ligand symmetries and target/model ligand
- matching/assignment. But actual score computation is delegated to child classes.
+ matching/assignment. But actual score computation is delegated to child
+ classes.
At the moment, two such classes are available:
@@ -26,25 +27,10 @@ class LigandScorer:
* :class:`ost.mol.alg.ligand_scoring_scrmsd.SCRMSDScorer`
that computes a binding-site superposed, symmetry-corrected RMSD.
- By default, only exact matches between target and model ligands are
- considered. This is a problem when the target only contains a subset
- of the expected atoms (for instance if atoms are missing in an
- experimental structure, which often happens in the PDB). With
- `substructure_match=True`, complete model ligands can be scored against
- partial target ligands. One problem with this approach is that it is
- very easy to find good matches to small, irrelevant ligands like EDO, CO2
- or GOL. To counter that, the assignment algorithm considers the coverage,
- expressed as the fraction of atoms of the model ligand atoms covered in the
- target. Higher coverage matches are prioritized, but a match with a better
- score will be preferred if it falls within a window of `coverage_delta`
- (by default 0.2) of a worse-scoring match. As a result, for instance,
- with a delta of 0.2, a low-score match with coverage 0.96 would be
- preferred over a high-score match with coverage 0.70.
-
All versus all scores are available through the lazily computed
:attr:`score_matrix`. However, many things can go wrong... be it even
something as simple as two ligands not matching. Error states therefore
- encode scoring issues. Issues with a particular ligand are indicated by a
+ encode scoring issues. An Issue for a particular ligand is indicated by a
non-zero state in :attr:`model_ligand_states`/:attr:`target_ligand_states`.
This invalidates pairwise scores of such a ligand with all other ligands.
This and other issues in pairwise score computation are reported in
@@ -52,11 +38,29 @@ class LigandScorer:
Only if the respective location is 0, a valid pairwise score can be
expected. The states and their meaning can be explored with code::
- for state_code, (short_desc, desc) in scorer_obj.items():
+ for state_code, (short_desc, desc) in scorer_obj.state_decoding.items():
print(state_code)
print(short_desc)
print(desc)
+ A common use case is to derive a one-to-one mapping between ligands in
+ the model and the target for which :class:`LigandScorer` provides an
+ automated assignment procedure.
+ By default, only exact matches between target and model ligands are
+ considered. This is a problem when the target only contains a subset
+ of the expected atoms (for instance if atoms are missing in an
+ experimental structure, which often happens in the PDB). With
+ `substructure_match=True`, complete model ligands can be scored against
+ partial target ligands. One problem with this approach is that it is
+ very easy to find good matches to small, irrelevant ligands like EDO, CO2
+ or GOL. The assignment algorithm therefore considers the coverage,
+ expressed as the fraction of atoms of the model ligand atoms covered in the
+ target. Higher coverage matches are prioritized, but a match with a better
+ score will be preferred if it falls within a window of `coverage_delta`
+ (by default 0.2) of a worse-scoring match. As a result, for instance,
+ with a delta of 0.2, a low-score match with coverage 0.96 would be
+ preferred over a high-score match with coverage 0.70.
+
Assumptions:
:class:`LigandScorer` generally assumes that the
@@ -71,11 +75,11 @@ class LigandScorer:
scoring. :ref:`Molck <molck>` should be used with extra
care, as many of the options (such as `rm_non_std` or `map_nonstd_res`) can
cause ligands to be removed from the structure. If cleanup with Molck is
- needed, ligands should be kept aside and passed separately. Non-ligand residues
- should be valid compounds with atom names following the naming conventions
- of the component dictionary. Non-standard residues are acceptable, and if
- the model contains a standard residue at that position, only atoms with
- matching names will be considered.
+ needed, ligands should be kept aside and passed separately. Non-ligand
+ residues should be valid compounds with atom names following the naming
+ conventions of the component dictionary. Non-standard residues are
+ acceptable, and if the model contains a standard residue at that position,
+ only atoms with matching names will be considered.
Unlike most of OpenStructure, this class does not assume that the ligands
(either in the model or the target) are part of the PDB component
@@ -129,35 +133,47 @@ class LigandScorer:
# Setup scorer object and compute lDDT-PLI
model_ligands = [model_ligand.Select("ele != H")]
- ls = SCRMSDScorer(cleaned_model, cleaned_target, model_ligands)
+ sc = SCRMSDScorer(cleaned_model, cleaned_target, model_ligands)
+
+ # Perform assignment and read respective scores
+ for lig_pair in sc.assignment:
+ trg_lig = sc.target_ligands[lig_pair[0]]
+ mdl_lig = sc.model_ligands[lig_pair[1]]
+ score = sc.score_matrix[lig_pair[0], lig_pair[1]]
+ print(f"Score for {trg_lig} and {mdl_lig}: {score}")
:param model: Model structure - a deep copy is available as :attr:`model`.
No additional processing (ie. Molck), checks,
stereochemistry checks or sanitization is performed on the
input. Hydrogen atoms are kept.
:type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
- :param target: Target structure - a deep copy is available as :attr:`target`.
- No additional processing (ie. Molck), checks or sanitization
- is performed on the input. Hydrogen atoms are kept.
+ :param target: Target structure - a deep copy is available as
+ :attr:`target`. No additional processing (ie. Molck), checks
+ or sanitization is performed on the input. Hydrogen atoms are
+ kept.
:type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
:param model_ligands: Model ligands, as a list of
- :class:`~ost.mol.ResidueHandle` belonging to the model
- entity. Can be instantiated with either a :class:list of
- :class:`~ost.mol.ResidueHandle`/:class:`ost.mol.ResidueView`
- or of :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`.
- If `None`, ligands will be extracted based on the
- :attr:`~ost.mol.ResidueHandle.is_ligand` flag (this is
- normally set properly in entities loaded from mmCIF).
+ :class:`~ost.mol.ResidueHandle` belonging to the model
+ entity. Can be instantiated with either a :class:list
+ of :class:`~ost.mol.ResidueHandle`/
+ :class:`ost.mol.ResidueView` or of
+ :class:`ost.mol.EntityHandle`/
+ :class:`ost.mol.EntityView`.
+ If `None`, ligands will be extracted based on the
+ :attr:`~ost.mol.ResidueHandle.is_ligand` flag (this is
+ normally set properly in entities loaded from mmCIF).
:type model_ligands: :class:`list`
:param target_ligands: Target ligands, as a list of
- :class:`~ost.mol.ResidueHandle` belonging to the target
- entity. Can be instantiated either a :class:list of
- :class:`~ost.mol.ResidueHandle`/:class:`ost.mol.ResidueView`
- or of :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
- containing a single residue each.
- If `None`, ligands will be extracted based on the
- :attr:`~ost.mol.ResidueHandle.is_ligand` flag (this is
- normally set properly in entities loaded from mmCIF).
+ :class:`~ost.mol.ResidueHandle` belonging to the
+ target entity. Can be instantiated either a
+ :class:list of :class:`~ost.mol.ResidueHandle`/
+ :class:`ost.mol.ResidueView` or of
+ :class:`ost.mol.EntityHandle`/
+ :class:`ost.mol.EntityView` containing a single
+ residue each. If `None`, ligands will be extracted
+ based on the :attr:`~ost.mol.ResidueHandle.is_ligand`
+ flag (this is normally set properly in entities
+ loaded from mmCIF).
:type target_ligands: :class:`list`
:param resnum_alignments: Whether alignments between chemically equivalent
chains in *model* and *target* can be computed
@@ -174,14 +190,14 @@ class LigandScorer:
will be logged to the console with SCRIPT
level.
:type rename_ligand_chain: :class:`bool`
- :param substructure_match: Set this to True to allow incomplete (ie
+ :param substructure_match: Set this to True to allow incomplete (i.e.
partially resolved) target ligands.
:type substructure_match: :class:`bool`
:param coverage_delta: the coverage delta for partial ligand assignment.
:type coverage_delta: :class:`float`
:param max_symmetries: If more than that many isomorphisms exist for
- a target-ligand pair, it will be ignored and reported
- as unassigned.
+ a target-ligand pair, it will be ignored and reported
+ as unassigned.
:type max_symmetries: :class:`int`
"""
@@ -236,8 +252,8 @@ class LigandScorer:
self.__chain_mapper = None
# keep track of states
- # simple integers instead of enums - documentation of property describes
- # encoding
+ # simple integers instead of enums - encoding available in
+ # self.state_decoding
self._state_matrix = None
self._model_ligand_states = None
self._target_ligand_states = None
@@ -284,7 +300,8 @@ class LigandScorer:
Ligand pairs can be matched and a valid score can be expected if
respective location in this matrix is 0.
Target ligands are in rows, model ligands in columns. States are encoded
- as integers <= 9. Larger numbers encode errors for child classes.
+ as integers <= 9. Larger numbers encode errors for child classes.
+ Use something like ``self.state_decoding[3]`` to get a decscription.
:rtype: :class:`~numpy.ndarray`
"""
@@ -325,8 +342,8 @@ class LigandScorer:
Target ligands are in rows, model ligands in columns.
NaN values indicate that no value could be computed (i.e. different
- ligands). In other words: values are only valid if respective location
- :attr:`~state_matrix` is 0.
+ ligands). In other words: values are only valid if the respective
+ location in :attr:`~state_matrix` is 0.
:rtype: :class:`~numpy.ndarray`
"""
@@ -341,10 +358,10 @@ class LigandScorer:
Target ligands are in rows, model ligands in columns.
NaN values indicate that no value could be computed (i.e. different
- ligands). In other words: values are only valid if respective location
- :attr:`~state_matrix` is 0. If `substructure_match=False`, only full
- match isomorphisms are considered, and therefore only values of 1.0
- can be observed.
+ ligands). In other words: values are only valid if the respective
+ location in :attr:`~state_matrix` is 0. If `substructure_match=False`,
+ only full match isomorphisms are considered, and therefore only values
+ of 1.0 can be observed.
:rtype: :class:`~numpy.ndarray`
"""
@@ -362,7 +379,7 @@ class LigandScorer:
class to provide additional information for a scored ligand pair.
empty dictionaries indicate that the child class simply didn't return
anything or that no value could be computed (e.g. different ligands).
- In other words: values are only valid if respective location
+ In other words: values are only valid if respective location in the
:attr:`~state_matrix` is 0.
:rtype: :class:`~numpy.ndarray`
@@ -429,7 +446,7 @@ class LigandScorer:
""" Get a dictionary of score values, keyed by model ligand
Extract score with something like:
- `scorer.score[lig.GetChain().GetName()][lig.GetNumber()]`.
+ ``scorer.score[lig.GetChain().GetName()][lig.GetNumber()]``.
The returned scores are based on :attr:`~assignment`.
:rtype: :class:`dict`
@@ -451,7 +468,7 @@ class LigandScorer:
""" Get a dictionary of score details, keyed by model ligand
Extract dict with something like:
- `scorer.score[lig.GetChain().GetName()][lig.GetNumber()]`.
+ ``scorer.score[lig.GetChain().GetName()][lig.GetNumber()]``.
The returned info dicts are based on :attr:`~assignment`. The content is
documented in the respective child class.
@@ -767,23 +784,27 @@ class LigandScorer:
if rename_chain:
chain_ext = 2 # Extend the chain name by this
while True:
- new_chain_name = residue.chain.name + "_" + str(chain_ext)
+ new_chain_name = \
+ residue.chain.name + "_" + str(chain_ext)
new_chain = new_entity.FindChain(new_chain_name)
if new_chain.IsValid():
chain_ext += 1
continue
else:
- new_chain = new_editor.InsertChain(new_chain_name)
+ new_chain = \
+ new_editor.InsertChain(new_chain_name)
break
LogScript("Moved ligand residue %s to new chain %s" % (
residue.qualified_name, new_chain.name))
else:
- msg = "A residue number %s already exists in chain %s" % (
+ msg = \
+ "A residue number %s already exists in chain %s" % (
residue.number, residue.chain.name)
raise RuntimeError(msg)
# Add the residue with its original residue number
- new_res = new_editor.AppendResidue(new_chain, residue.name, residue.number)
+ new_res = new_editor.AppendResidue(new_chain, residue.name,
+ residue.number)
# Add atoms
for old_atom in residue.atoms:
new_editor.InsertAtom(new_res, old_atom.name, old_atom.pos,
@@ -802,28 +823,35 @@ class LigandScorer:
new_res = None
if res.entity.handle == old_entity.handle:
# Residue is part of the old_entity handle.
- # However, it may not be in the copied one, for instance it may have been a view
- # We try to grab it first, otherwise we copy it
+ # However, it may not be in the copied one, for instance it may
+ # have been a view We try to grab it first, otherwise we copy it
new_res = new_entity.FindResidue(res.chain.name, res.number)
if new_res and new_res.valid:
- LogVerbose("Ligand residue %s already in entity" % res.handle.qualified_name)
+ qname = res.handle.qualified_name
+ LogVerbose("Ligand residue %s already in entity" % qname)
else:
# Residue is not part of the entity, need to copy it first
new_res = _copy_residue(res, rename_chain)
- LogVerbose("Copied ligand residue %s" % res.handle.qualified_name)
+ qname = res.handle.qualified_name
+ LogVerbose("Copied ligand residue %s" % qname)
new_res.SetIsLigand(True)
return new_res
for ligand in ligands:
- if isinstance(ligand, mol.EntityHandle) or isinstance(ligand, mol.EntityView):
+ is_eh = isinstance(ligand, mol.EntityHandle)
+ is_ev = isinstance(ligand, mol.EntityView)
+ is_rh = isinstance(ligand, mol.ResidueHandle)
+ is_rv = isinstance(ligand, mol.ResidueView)
+ if is_eh or is_ev:
for residue in ligand.residues:
new_residue = _process_ligand_residue(residue, rename_chain)
extracted_ligands.append(new_residue)
- elif isinstance(ligand, mol.ResidueHandle) or isinstance(ligand, mol.ResidueView):
+ elif is_rh or is_rv:
new_residue = _process_ligand_residue(ligand, rename_chain)
extracted_ligands.append(new_residue)
else:
- raise RuntimeError("Ligands should be given as Entity or Residue")
+ raise RuntimeError("Ligands should be given as Entity or "
+ "Residue")
if new_editor is not None:
new_editor.UpdateICS()
@@ -846,22 +874,26 @@ class LigandScorer:
self._aux_matrix = np.empty(shape, dtype=dict)
# precompute ligand graphs
- target_graphs = [_ResidueToGraph(l, by_atom_index=True) for l in self.target_ligands]
- model_graphs = [_ResidueToGraph(l, by_atom_index=True) for l in self.model_ligands]
+ target_graphs = \
+ [_ResidueToGraph(l, by_atom_index=True) for l in self.target_ligands]
+ model_graphs = \
+ [_ResidueToGraph(l, by_atom_index=True) for l in self.model_ligands]
for g_idx, g in enumerate(target_graphs):
if not networkx.is_connected(g):
self._target_ligand_states[g_idx] = 4
self._state_matrix[g_idx,:] = 6
- LogVerbose("Disconnected graph observed for target ligand %s" % (
- str(self.target_ligands[g_idx])))
+ msg = "Disconnected graph observed for target ligand "
+ msg += str(self.target_ligands[g_idx])
+ LogVerbose(msg)
for g_idx, g in enumerate(model_graphs):
if not networkx.is_connected(g):
self._model_ligand_states[g_idx] = 4
self._state_matrix[:,g_idx] = 6
- LogVerbose("Disconnected graph observed for model ligand %s" % (
- str(self.model_ligands[g_idx])))
+ msg = "Disconnected graph observed for model ligand "
+ msg += str(self.model_ligands[g_idx])
+ LogVerbose(msg)
for target_id, target_ligand in enumerate(self.target_ligands):
@@ -1001,7 +1033,8 @@ class LigandScorer:
Returned score must be valid in this case (not None/NaN).
Child specific non-zero states must be >= 10.
"""
- raise NotImplementedError("_compute must be implemented by child class")
+ raise NotImplementedError("_compute must be implemented by child "
+ "class")
def _score_dir(self):
""" Return direction of score - defined by child class
@@ -1010,7 +1043,8 @@ class LigandScorer:
'+' for ascending scores, i.e. higher is better (lddt etc.)
'-' for descending scores, i.e. lower is better (rmsd etc.)
"""
- raise NotImplementedError("_score_dir must be implemented by child class")
+ raise NotImplementedError("_score_dir must be implemented by child "
+ "class")
def _ResidueToGraph(residue, by_atom_index=False):
@@ -1026,7 +1060,8 @@ def _ResidueToGraph(residue, by_atom_index=False):
:type by_atom_index: :class:`bool`
:rtype: :class:`~networkx.classes.graph.Graph`
- Nodes are labeled with the Atom's uppercase :attr:`~ost.mol.AtomHandle.element`.
+ Nodes are labeled with the Atom's uppercase
+ :attr:`~ost.mol.AtomHandle.element`.
"""
nxg = networkx.Graph()
@@ -1034,7 +1069,8 @@ def _ResidueToGraph(residue, by_atom_index=False):
nxg.add_node(atom.name, element=atom.element.upper())
# This will list all edges twice - once for every atom of the pair.
- # But as of NetworkX 3.0 adding the same edge twice has no effect, so we're good.
+ # But as of NetworkX 3.0 adding the same edge twice has no effect,
+ # so we're good.
nxg.add_edges_from([(
b.first.name,
b.second.name) for a in residue.atoms for b in a.GetBondList()])
@@ -1092,7 +1128,8 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
by_atom_index=by_atom_index)
if not networkx.is_connected(model_graph):
- raise DisconnectedGraphError("Disconnected graph for model ligand %s" % model_ligand)
+ msg = "Disconnected graph for model ligand %s" % model_ligand
+ raise DisconnectedGraphError(msg)
if target_graph is None:
@@ -1100,11 +1137,12 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
by_atom_index=by_atom_index)
if not networkx.is_connected(target_graph):
- raise DisconnectedGraphError("Disconnected graph for target ligand %s" % target_ligand)
+ msg = "Disconnected graph for target ligand %s" % target_ligand
+ raise DisconnectedGraphError(msg)
# Note the argument order (model, target) which differs from spyrmsd.
- # This is because a subgraph of model is isomorphic to target - but not the opposite
- # as we only consider partial ligands in the reference.
+ # This is because a subgraph of model is isomorphic to target - but not the
+ # opposite as we only consider partial ligands in the reference.
# Make sure to generate the symmetries correctly in the end
gm = networkx.algorithms.isomorphism.GraphMatcher(
model_graph, target_graph, node_match=lambda x, y:
@@ -1118,7 +1156,8 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
raise TooManySymmetriesError(
"Too many symmetries between %s and %s" % (
str(model_ligand), str(target_ligand)))
- symmetries.append((list(isomorphism.values()), list(isomorphism.keys())))
+ symmetries.append((list(isomorphism.values()),
+ list(isomorphism.keys())))
assert len(symmetries) > 0
LogDebug("Found %s isomorphic mappings (symmetries)" % len(symmetries))
elif gm.subgraph_is_isomorphic() and substructure_match:
@@ -1130,11 +1169,13 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
raise TooManySymmetriesError(
"Too many symmetries between %s and %s" % (
str(model_ligand), str(target_ligand)))
- symmetries.append((list(isomorphism.values()), list(isomorphism.keys())))
+ symmetries.append((list(isomorphism.values()),
+ list(isomorphism.keys())))
assert len(symmetries) > 0
# Assert that all the atoms in the target are part of the substructure
assert len(symmetries[0][0]) == len(target_ligand.atoms)
- LogDebug("Found %s subgraph isomorphisms (symmetries)" % len(symmetries))
+ n_sym = len(symmetries)
+ LogDebug("Found %s subgraph isomorphisms (symmetries)" % n_sym)
elif gm.subgraph_is_isomorphic():
LogDebug("Found subgraph isomorphisms (symmetries), but"
" ignoring because substructure_match=False")
diff --git a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
index b747921221949422f537b18328c156567823ea2c..f4ae9c4c6c27c4879119f2ade73ebe7f07e4fcbc 100644
--- a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
+++ b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
@@ -52,9 +52,9 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
* target_ligand: The actual target ligand for which the score was computed
* model_ligand: The actual model ligand for which the score was computed
* bs_ref_res: :class:`set` of residues with potentially non-zero
- contribution to score. That is every residue with at least one
- atom within *lddt_pli_radius* + max(*lddt_pli_thresholds*) of
- the ligand.
+ contribution to score. That is every residue with at least one
+ atom within *lddt_pli_radius* + max(*lddt_pli_thresholds*) of
+ the ligand.
* bs_mdl_res: Same for model
:param model: Passed to parent constructor - see :class:`LigandScorer`.
@@ -292,7 +292,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
if a is not None and b is not None:
flat_mapping[b] = a
- # for each mdl bs atom (as atom hash), the trg bs atoms (as index in scorer)
+ # for each mdl bs atom (as atom hash), the trg bs atoms (as index in
+ # scorer)
bs_atom_mapping = dict()
for mdl_cname, ref_cname in flat_mapping.items():
aln = cut_ref_mdl_alns[(ref_cname, mdl_cname)]
@@ -347,7 +348,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
# mask selects entries in trg_bs_indices that are not yet
# part of classic lDDT ref_indices for atom at trg_a_idx
# => added mdl contacts
- mask = np.isin(trg_bs_indices, ref_indices[ligand_start_idx + trg_a_idx],
+ mask = np.isin(trg_bs_indices,
+ ref_indices[ligand_start_idx + trg_a_idx],
assume_unique=True, invert=True)
added_indices = np.asarray([], dtype=np.int64)
added_distances = np.asarray([], dtype=np.float64)
@@ -355,17 +357,19 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
# compute ref distances on reference positions
added_indices = trg_bs_indices[mask]
tmp = scorer.positions.take(added_indices, axis=0)
- np.subtract(tmp, trg_ligand_pos[trg_a_idx][None, :], out=tmp)
+ np.subtract(tmp, trg_ligand_pos[trg_a_idx][None, :],
+ out=tmp)
np.square(tmp, out=tmp)
tmp = tmp.sum(axis=1)
- np.sqrt(tmp, out=tmp) # distances against all relevant atoms
+ np.sqrt(tmp, out=tmp)
added_distances = tmp
# extract the distances towards bs atoms that are symmetric
sym_mask = np.isin(added_indices, symmetric_atoms,
assume_unique=True)
- cache[(mdl_a_idx, trg_a_idx)] = (added_indices, added_distances,
+ cache[(mdl_a_idx, trg_a_idx)] = (added_indices,
+ added_distances,
added_indices[sym_mask],
added_distances[sym_mask])
@@ -401,7 +405,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
sym_idx_collector = [None] * scorer.n_atoms
sym_dist_collector = [None] * scorer.n_atoms
- for mapping, s_cache, p_cache in zip(chain_mappings, scoring_cache, penalty_cache):
+ for mapping, s_cache, p_cache in zip(chain_mappings, scoring_cache,
+ penalty_cache):
lddt_chain_mapping = dict()
lddt_alns = dict()
@@ -437,17 +442,17 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
for mdl_ch in mdl_chains:
if mdl_ch not in lddt_chain_mapping:
# check which chain in trg is closest
- chem_group_idx = None
+ chem_grp_idx = None
for i, m in enumerate(self._chem_mapping):
if mdl_ch in m:
- chem_group_idx = i
+ chem_grp_idx = i
break
- if chem_group_idx is None:
+ if chem_grp_idx is None:
raise RuntimeError("This should never happen... "
"ask Gabriel...")
closest_ch = None
closest_dist = None
- for trg_ch in self._chain_mapper.chem_groups[chem_group_idx]:
+ for trg_ch in self._chain_mapper.chem_groups[chem_grp_idx]:
if trg_ch not in lddt_chain_mapping.values():
if trg_ch not in already_mapped:
ch = self._chain_mapper.target.FindChain(trg_ch)
@@ -526,7 +531,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
conserved = np.sum(scorer._EvalAtoms(pos, ligand_at_indices,
self.lddt_pli_thresholds,
funky_ref_indices,
- funky_ref_distances), axis=0)
+ funky_ref_distances),
+ axis=0)
score = None
if N > 0:
score = np.mean(conserved/N)
@@ -692,7 +698,8 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
# the select statement also excludes the ligand in mdl_bs
# as it resides in a separate chain
mdl_cname = ch_tuple[1]
- mdl_bs_ch = mdl_bs.Select(f"cname={mol.QueryQuoteName(mdl_cname)}")
+ query = "cname=" + mol.QueryQuoteName(mdl_cname)
+ mdl_bs_ch = mdl_bs.Select(query)
for a in mdl_ligand_res.atoms:
close_atoms = \
mdl_bs_ch.FindWithin(a.GetPos(), self.lddt_pli_radius)
@@ -895,8 +902,10 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
self.__mappable_atoms = dict()
for (ref_cname, mdl_cname), aln in self._ref_mdl_alns.items():
self._mappable_atoms[(ref_cname, mdl_cname)] = set()
- ref_ch = self._chain_mapper.target.Select(f"cname={mol.QueryQuoteName(ref_cname)}")
- mdl_ch = self._chain_mapping_mdl.Select(f"cname={mol.QueryQuoteName(mdl_cname)}")
+ ref_query = f"cname={mol.QueryQuoteName(ref_cname)}"
+ mdl_query = f"cname={mol.QueryQuoteName(mdl_cname)}"
+ ref_ch = self._chain_mapper.target.Select(ref_query)
+ mdl_ch = self._chain_mapping_mdl.Select(mdl_query)
aln.AttachView(0, ref_ch)
aln.AttachView(1, mdl_ch)
for col in aln:
@@ -932,9 +941,9 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
if self.__ref_mdl_alns is None:
self.__ref_mdl_alns = \
chain_mapping._GetRefMdlAlns(self._chain_mapper.chem_groups,
- self._chain_mapper.chem_group_alignments,
- self._chem_mapping,
- self._chem_group_alns)
+ self._chain_mapper.chem_group_alignments,
+ self._chem_mapping,
+ self._chem_group_alns)
return self.__ref_mdl_alns
@property
diff --git a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
index 5006b9b30e036f2f8a4855a6945be11e72cef286..c76a23597af8c7d46b3cdcd7bf3772e3525cdecf 100644
--- a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
+++ b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
@@ -40,22 +40,21 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
* lddt_lp: lDDT of the binding site used for superposition
* bs_ref_res: :class:`list` of binding site residues in target
* bs_ref_res_mapped: :class:`list` of target binding site residues that
- are mapped to model
+ are mapped to model
* bs_mdl_res_mapped: :class:`list` of same length with respective model
- residues
+ residues
* bb_rmsd: Backbone RMSD (CA, C3' for nucleotides) for mapped residues
- given transform
+ given transform
* target_ligand: The actual target ligand for which the score was computed
* model_ligand: The actual model ligand for which the score was computed
* chain_mapping: :class:`dict` with a chain mapping of chains involved in
- binding site - key: trg chain name, value: mdl chain name
+ binding site - key: trg chain name, value: mdl chain name
* transform: :class:`geom.Mat4` to transform model binding site onto target
- binding site
+ binding site
* inconsistent_residues: :class:`list` of :class:`tuple` representing
- residues with inconsistent residue names upon
- mapping (which is given by bs_ref_res_mapped
- and bs_mdl_res_mapped). Tuples have two elements:
- 1) trg residue 2) mdl residue
+ residues with inconsistent residue names upon mapping (which is given by
+ bs_ref_res_mapped and bs_mdl_res_mapped). Tuples have two elements:
+ 1) trg residue 2) mdl residue
:param model: Passed to parent constructor - see :class:`LigandScorer`.
:type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
@@ -180,7 +179,8 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
rmsd = _SCRMSD_symmetries(symmetries, model_ligand,
target_ligand, transformation=r.transform)
- if best_rmsd_result["rmsd"] is None or rmsd < best_rmsd_result["rmsd"]:
+ if best_rmsd_result["rmsd"] is None or \
+ rmsd < best_rmsd_result["rmsd"]:
best_rmsd_result = {"rmsd": rmsd,
"lddt_lp": r.lDDT,
"bs_ref_res": r.substructure.residues,
@@ -191,7 +191,8 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
"model_ligand": model_ligand,
"chain_mapping": r.GetFlatChainMapping(),
"transform": r.transform,
- "inconsistent_residues": r.inconsistent_residues}
+ "inconsistent_residues":
+ r.inconsistent_residues}
target_ligand_state = 0
model_ligand_state = 0
@@ -203,13 +204,15 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
# try to identify error states
best_rmsd = np.nan
error_state = 20 # unknown error
- if self._get_target_binding_site(target_ligand).GetResidueCount() == 0:
+ N = self._get_target_binding_site(target_ligand).GetResidueCount()
+ if N == 0:
pair_state = 6 # binding_site
target_ligand_state = 10
elif len(representations) == 0:
pair_state = 11
- return (best_rmsd, pair_state, target_ligand_state, model_ligand_state, best_rmsd_result)
+ return (best_rmsd, pair_state, target_ligand_state, model_ligand_state,
+ best_rmsd_result)
def _score_dir(self):
""" Implements interface from parent
@@ -223,7 +226,8 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
# all possible binding sites, independent from actual model ligand
key = (target_ligand.handle.hash_code, 0)
else:
- key = (target_ligand.handle.hash_code, model_ligand.handle.hash_code)
+ key = (target_ligand.handle.hash_code,
+ model_ligand.handle.hash_code)
if key not in self._repr:
ref_bs = self._get_target_binding_site(target_ligand)
@@ -232,10 +236,12 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
ref_bs, self.model, inclusion_radius=self.lddt_lp_radius,
topn=self.binding_sites_topn)
else:
+ repr_in = self._get_get_repr_input(model_ligand)
+ radius = self.lddt_lp_radius
reprs = self._chain_mapper.GetRepr(ref_bs, self.model,
- inclusion_radius=self.lddt_lp_radius,
+ inclusion_radius=radius,
topn=self.binding_sites_topn,
- chem_mapping_result = self._get_get_repr_input(model_ligand))
+ chem_mapping_result=repr_in)
self._repr[key] = reprs
return self._repr[key]
@@ -245,35 +251,41 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
if target_ligand.handle.hash_code not in self._binding_sites:
# create view of reference binding site
- ref_residues_hashes = set() # helper to keep track of added residues
+ ref_residues_hashes = set() # helper to keep track of added residues
ignored_residue_hashes = {target_ligand.hash_code}
for ligand_at in target_ligand.atoms:
- close_atoms = self.target.FindWithin(ligand_at.GetPos(), self.bs_radius)
+ close_atoms = self.target.FindWithin(ligand_at.GetPos(),
+ self.bs_radius)
for close_at in close_atoms:
# Skip any residue not in the chain mapping target
ref_res = close_at.GetResidue()
h = ref_res.handle.GetHashCode()
if h not in ref_residues_hashes and \
h not in ignored_residue_hashes:
- if self._chain_mapper.target.ViewForHandle(ref_res).IsValid():
+ view = self._chain_mapper.target.ViewForHandle(ref_res)
+ if view.IsValid():
h = ref_res.handle.GetHashCode()
ref_residues_hashes.add(h)
elif ref_res.is_ligand:
- LogWarning("Ignoring ligand %s in binding site of %s" % (
- ref_res.qualified_name, target_ligand.qualified_name))
+ msg = f"Ignoring ligand {ref_res.qualified_name} "
+ msg += "in binding site of "
+ msg += str(target_ligand.qualified_name)
+ LogWarning(msg)
ignored_residue_hashes.add(h)
elif ref_res.chem_type == mol.ChemType.WATERS:
pass # That's ok, no need to warn
else:
- LogWarning("Ignoring residue %s in binding site of %s" % (
- ref_res.qualified_name, target_ligand.qualified_name))
+ msg = f"Ignoring residue {ref_res.qualified_name} "
+ msg += "in binding site of "
+ msg += str(target_ligand.qualified_name)
+ LogWarning(msg)
ignored_residue_hashes.add(h)
ref_bs = self.target.CreateEmptyView()
if ref_residues_hashes:
- # reason for doing that separately is to guarantee same ordering of
- # residues as in underlying entity. (Reorder by ResNum seems only
- # available on ChainHandles)
+ # reason for doing that separately is to guarantee same ordering
+ # of residues as in underlying entity. (Reorder by ResNum seems
+ # only available on ChainHandles)
for ch in self.target.chains:
for r in ch.residues:
if r.handle.GetHashCode() in ref_residues_hashes:
@@ -307,9 +319,9 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
if self.__ref_mdl_alns is None:
self.__ref_mdl_alns = \
chain_mapping._GetRefMdlAlns(self._chain_mapper.chem_groups,
- self._chain_mapper.chem_group_alignments,
- self._chem_mapping,
- self._chem_group_alns)
+ self._chain_mapper.chem_group_alignments,
+ self._chem_mapping,
+ self._chem_group_alns)
return self.__ref_mdl_alns
@property
@@ -383,17 +395,19 @@ def SCRMSD(model_ligand, target_ligand, transformation=geom.Mat4(),
generating at least that many isomorphisms and can take some time.
:type max_symmetries: :class:`int`
:rtype: :class:`float`
- :raises: :class:`NoSymmetryError` when no symmetry can be found,
- :class:`DisconnectedGraphError` when ligand graph is disconnected,
- :class:`TooManySymmetriesError` when more than `max_symmetries`
- isomorphisms are found.
+ :raises: :class:`ost.mol.alg.ligand_scoring_base.NoSymmetryError` when no
+ symmetry can be found,
+ :class:`ost.mol.alg.ligand_scoring_base.DisconnectedGraphError`
+ when ligand graph is disconnected,
+ :class:`ost.mol.alg.ligand_scoring_base.TooManySymmetriesError`
+ when more than *max_symmetries* isomorphisms are found.
"""
symmetries = ligand_scoring_base.ComputeSymmetries(model_ligand,
- target_ligand,
- substructure_match=substructure_match,
- by_atom_index=True,
- max_symmetries=max_symmetries)
+ target_ligand,
+ substructure_match=substructure_match,
+ by_atom_index=True,
+ max_symmetries=max_symmetries)
return _SCRMSD_symmetries(symmetries, model_ligand, target_ligand,
transformation)