diff --git a/modules/mol/alg/pymod/ligand_scoring_base.py b/modules/mol/alg/pymod/ligand_scoring_base.py
index c6e3d713fac3e8be389196cd4e88441312a6da2c..71223951e95d460b6802be630f5a989902121d27 100644
--- a/modules/mol/alg/pymod/ligand_scoring_base.py
+++ b/modules/mol/alg/pymod/ligand_scoring_base.py
@@ -222,6 +222,8 @@ class LigandScorer:
# simple integers instead of enums - documentation of property describes
# encoding
self._state_matrix = None
+ self._model_ligand_states = None
+ self._target_ligand_states = None
# score matrices
self._score_matrix = None
@@ -243,6 +245,7 @@ class LigandScorer:
iso = "subgraph isomorphism"
else:
iso = "full graph isomorphism"
+
self.state_decoding = \
{0: ("OK", "OK"),
1: ("identity", f"Ligands could not be matched (by {iso})"),
@@ -251,6 +254,10 @@ class LigandScorer:
3: ("no_iso", "No fully isomorphic match could be found - enabling "
"substructure_match might allow a match"),
4: ("disconnected", "Ligand graph is disconnected"),
+ 5: ("stoichiometry", "Ligand was already assigned to another ligand "
+ "(different stoichiometry)"),
+ 6: ("single_ligand_issue", "Cannot compute valid pairwise score as "
+ "either model or target ligand have non-zero state."),
9: ("unknown", "An unknown error occured in LigandScorer")}
@property
@@ -270,6 +277,32 @@ class LigandScorer:
self._compute_scores()
return self._state_matrix
+ @property
+ def model_ligand_states(self):
+ """ Encodes states of model ligands
+
+ Non-zero state in any of the model ligands invalidates the full
+ respective column in :attr:`~state_matrix`.
+
+ :rtype: :class:`~numpy.ndarray`
+ """
+ if self._model_ligand_states is None:
+ self._compute_scores()
+ return self._model_ligand_states
+
+ @property
+ def target_ligand_states(self):
+ """ Encodes states of target ligands
+
+ Non-zero state in any of the target ligands invalidates the full
+ respective row in :attr:`~state_matrix`.
+
+ :rtype: :class:`~numpy.ndarray`
+ """
+ if self._target_ligand_states is None:
+ self._compute_scores()
+ return self._target_ligand_states
+
@property
def score_matrix(self):
""" Get the matrix of scores.
@@ -427,6 +460,7 @@ class LigandScorer:
:rtype: :class:`list` of :class:`int`
"""
+ # compute on-the-fly, no need for caching
assigned = set([x[0] for x in self.assignment])
return [x for x in range(len(self.target_ligands)) if x not in assigned]
@@ -436,6 +470,7 @@ class LigandScorer:
:rtype: :class:`list` of :class:`int`
"""
+ # compute on-the-fly, no need for caching
assigned = set([x[1] for x in self.assignment])
return [x for x in range(len(self.model_ligands)) if x not in assigned]
@@ -448,7 +483,8 @@ class LigandScorer:
generated
:type trg_lig_idx: :class:`int`
"""
- return self._get_report(self.state_matrix[trg_lig_idx,:])
+ return self._get_report(self.target_ligand_states[trg_lig_idx],
+ self.state_matrix[trg_lig_idx,:])
def get_model_ligand_state_report(self, mdl_lig_idx):
""" Get summary of states observed with respect to all target ligands
@@ -459,20 +495,180 @@ class LigandScorer:
generated
:type mdl_lig_idx: :class:`int`
"""
- return self._get_report(self.state_matrix[:, mdl_lig_idx])
-
+ return self._get_report(self.model_ligand_states[mdl_lig_idx],
+ self.state_matrix[:, mdl_lig_idx])
- def _get_report(self, states):
+ def _get_report(self, ligand_state, pair_states):
""" Helper
"""
- report = list()
- for s in np.unique(states):
+ pair_report = list()
+ for s in np.unique(pair_states):
desc = self.state_decoding[s]
- report.append({"state": s,
- "short desc": desc[0],
- "desc": desc[1],
- "indices": np.flatnonzero(states == s).tolist()})
- return report
+ indices = np.flatnonzero(pair_states == s).tolist()
+ pair_report.append({"state": s,
+ "short desc": desc[0],
+ "desc": desc[1],
+ "indices": indices})
+
+ desc = self.state_decoding[ligand_state]
+ ligand_report = {"state": ligand_state,
+ "short desc": desc[0],
+ "desc": desc[1]}
+
+ return (ligand_report, pair_report)
+
+ def guess_target_ligand_unassigned_reason(self, trg_lig_idx):
+ """ Makes an educated guess why target ligand is not assigned
+
+ This either returns actual error states or custom states that are
+ derived from them.
+
+ :param trg_lig_idx: Index of target ligand
+ :type trg_lig_idx: :class:`int`
+ :returns: :class:`tuple` with two elements: 1) keyword 2) human readable
+ sentence describing the issue, (\"unknown\",\"unknown\") if
+ nothing obvious can be found.
+ :raises: :class:`RuntimeError` if specified target ligand is assigned
+ """
+ if trg_lig_idx not in self.unassigned_target_ligands:
+ raise RuntimeError("Specified target ligand is not unassigned")
+
+ # hardcoded tuple if there is simply nothing we can assign it to
+ if len(self.model_ligands) == 0:
+ return ("no_ligand", "No ligand in the model")
+
+ # if something with the ligand itself is wrong, we can be pretty sure
+ # thats why the ligand is unassigned
+ if self.target_ligand_states[trg_lig_idx] != 0:
+ return self.state_decoding[self.target_ligand_states[trg_lig_idx]]
+
+ # The next best guess comes from looking at pair states
+ tmp = np.unique(self.state_matrix[trg_lig_idx,:])
+
+ # In case of any 0, it could have been assigned so it's probably
+ # just not selected due to different stoichiometry - this is no
+ # defined state, we just return a hardcoded tuple in this case
+ if 0 in tmp:
+ return ("stoichiometry",
+ "Ligand was already assigned to an other "
+ "model ligand (different stoichiometry)")
+
+ # maybe its a symmetry issue?
+ if 2 in tmp:
+ return self.state_decoding[2]
+
+ # if the state is 6 (single_ligand_issue), there is an issue with its
+ # target counterpart.
+ if 6 in tmp:
+ mdl_idx = np.where(self.state_matrix[trg_lig_idx,:]==6)[0]
+ # we're reporting everything except disconnected error...
+ # don't ask...
+ for i in mdl_idx:
+ if self.model_ligand_states[i] == 0:
+ raise RuntimeError("This should never happen")
+ if self.model_ligand_states[i] != 4:
+ return self.state_decoding[self.model_ligand_states[i]]
+
+ # get rid of remaining single ligand issues (only disconnected error)
+ if 6 in tmp and len(tmp) > 1:
+ tmp = tmp[tmp!=6]
+
+ # prefer everything over identity state
+ if 1 in tmp and len(tmp) > 1:
+ tmp = tmp[tmp!=1]
+
+ # just return whatever is left
+ return self.state_decoding[tmp[0]]
+
+
+ def guess_model_ligand_unassigned_reason(self, mdl_lig_idx):
+ """ Makes an educated guess why model ligand is not assigned
+
+ This either returns actual error states or custom states that are
+ derived from them.
+
+ :param mdl_lig_idx: Index of model ligand
+ :type mdl_lig_idx: :class:`int`
+ :returns: :class:`tuple` with two elements: 1) keyword 2) human readable
+ sentence describing the issue, (\"unknown\",\"unknown\") if
+ nothing obvious can be found.
+ :raises: :class:`RuntimeError` if specified model ligand is assigned
+ """
+ if mdl_lig_idx not in self.unassigned_model_ligands:
+ raise RuntimeError("Specified model ligand is not unassigned")
+
+ # hardcoded tuple if there is simply nothing we can assign it to
+ if len(self.target_ligands) == 0:
+ return ("no_ligand", "No ligand in the target")
+
+ # if something with the ligand itself is wrong, we can be pretty sure
+ # thats why the ligand is unassigned
+ if self.model_ligand_states[mdl_lig_idx] != 0:
+ return self.state_decoding[self.model_ligand_states[mdl_lig_idx]]
+
+ # The next best guess comes from looking at pair states
+ tmp = np.unique(self.state_matrix[:,mdl_lig_idx])
+
+ # In case of any 0, it could have been assigned so it's probably
+ # just not selected due to different stoichiometry - this is no
+ # defined state, we just return a hardcoded tuple in this case
+ if 0 in tmp:
+ return ("stoichiometry",
+ "Ligand was already assigned to an other "
+ "model ligand (different stoichiometry)")
+
+ # maybe its a symmetry issue?
+ if 2 in tmp:
+ return self.state_decoding[2]
+
+ # if the state is 6 (single_ligand_issue), there is an issue with its
+ # target counterpart.
+ if 6 in tmp:
+ trg_idx = np.where(self.state_matrix[:,mdl_lig_idx]==6)[0]
+ # we're reporting everything except disconnected error...
+ # don't ask...
+ for i in trg_idx:
+ if self.target_ligand_states[i] == 0:
+ raise RuntimeError("This should never happen")
+ if self.target_ligand_states[i] != 4:
+ return self.state_decoding[self.target_ligand_states[i]]
+
+ # get rid of remaining single ligand issues (only disconnected error)
+ if 6 in tmp and len(tmp) > 1:
+ tmp = tmp[tmp!=6]
+
+ # prefer everything over identity state
+ if 1 in tmp and len(tmp) > 1:
+ tmp = tmp[tmp!=1]
+
+ # just return whatever is left
+ return self.state_decoding[tmp[0]]
+
+ @property
+ def unassigned_model_ligands_reasons(self):
+ return_dict = dict()
+ for i in self.unassigned_model_ligands:
+ lig = self.model_ligands[i]
+ cname = lig.GetChain().GetName()
+ rnum = lig.GetNumber()
+ if cname not in return_dict:
+ return_dict[cname] = dict()
+ return_dict[cname][rnum] = \
+ self.guess_model_ligand_unassigned_reason(i)[0]
+ return return_dict
+
+ @property
+ def unassigned_target_ligands_reasons(self):
+ return_dict = dict()
+ for i in self.unassigned_target_ligands:
+ lig = self.target_ligands[i]
+ cname = lig.GetChain().GetName()
+ rnum = lig.GetNumber()
+ if cname not in return_dict:
+ return_dict[cname] = dict()
+ return_dict[cname][rnum] = \
+ self.guess_target_ligand_unassigned_reason(i)[0]
+ return return_dict
@property
def _chain_mapper(self):
@@ -629,23 +825,58 @@ class LigandScorer:
shape = (len(self.target_ligands), len(self.model_ligands))
self._score_matrix = np.full(shape, np.nan, dtype=np.float32)
self._coverage_matrix = np.full(shape, np.nan, dtype=np.float32)
- self._state_matrix = np.full(shape, -1, dtype=np.int32)
+ self._state_matrix = np.full(shape, 0, dtype=np.int32)
+ self._model_ligand_states = np.zeros(len(self.model_ligands))
+ self._target_ligand_states = np.zeros(len(self.target_ligands))
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]
+
+ 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])))
+
+ 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])))
+
+
for target_id, target_ligand in enumerate(self.target_ligands):
LogVerbose("Analyzing target ligand %s" % target_ligand)
+
+ if self._target_ligand_states[target_id] == 4:
+ # Disconnected graph - already updated states and reported
+ # to LogVerbose
+ continue
+
for model_id, model_ligand in enumerate(self.model_ligands):
LogVerbose("Compare to model ligand %s" % model_ligand)
#########################################################
# Compute symmetries for given target/model ligand pair #
#########################################################
+
+ if self._model_ligand_states[model_id] == 4:
+ # Disconnected graph - already updated states and reported
+ # to LogVerbose
+ continue
+
try:
symmetries = ComputeSymmetries(
model_ligand, target_ligand,
substructure_match=self.substructure_match,
by_atom_index=True,
- max_symmetries=self.max_symmetries)
+ max_symmetries=self.max_symmetries,
+ model_graph = model_graphs[model_id],
+ target_graph = target_graphs[target_id])
LogVerbose("Ligands %s and %s symmetry match" % (
str(model_ligand), str(target_ligand)))
except NoSymmetryError:
@@ -667,38 +898,57 @@ class LigandScorer:
self._state_matrix[target_id, model_id] = 3
continue
except DisconnectedGraphError:
+ # this should never happen as we guard against
+ # DisconnectedGraphError when precomputing the graph
LogVerbose("Disconnected graph observed for %s and %s" % (
str(model_ligand), str(target_ligand)))
- self._state_matrix[target_id, model_id] = 4
+ # just set both ligand states to 4
+ self._model_ligand_states[model_id] = 4
+ self._model_ligand_states[target_id] = 4
+ self._state_matrix[target_id, model_id] = 6
continue
#####################################################
# Compute score by calling the child class _compute #
#####################################################
- score, state, aux = self._compute(symmetries, target_ligand,
- model_ligand)
+ score, pair_state, target_ligand_state, model_ligand_state, aux\
+ = self._compute(symmetries, target_ligand, model_ligand)
############
# Finalize #
############
- if state != 0:
- # non-zero error states up to 9 are reserved for base class
- if state <= 9:
- raise RuntimeError("Child returned reserved err. state")
-
- # Ensure that returned state is associated with a
- # description. This is a requirement when subclassing
- # LigandScorer => state_decoding dict from base class must
- # be modified in subclass constructor
- if state not in self.state_decoding:
- raise RuntimeError(f"Subclass returned state "
- f"\"{state}\" for which no "
- f"description is available. Point "
- f"the developer of the used scorer "
- f"to this error message.")
-
- self._state_matrix[target_id, model_id] = state
- if state == 0:
+
+ # Ensure that returned states are associated with a
+ # description. This is a requirement when subclassing
+ # LigandScorer => state_decoding dict from base class must
+ # be modified in subclass constructor
+ if pair_state not in self.state_decoding or \
+ target_ligand_state not in self.state_decoding or \
+ model_ligand_state not in self.state_decoding:
+ raise RuntimeError(f"Subclass returned state "
+ f"\"{state}\" for which no "
+ f"description is available. Point "
+ f"the developer of the used scorer "
+ f"to this error message.")
+
+ # if any of the ligand states is non-zero, this must be marked
+ # by the scorer subclass by specifying a specific pair state
+ if target_ligand_state != 0 and pair_state != 6:
+ raise RuntimeError("Observed non-zero target-ligand "
+ "state which must trigger certain "
+ "pair state. Point the developer of "
+ "the used scorer to this error message")
+
+ if model_ligand_state != 0 and pair_state != 6:
+ raise RuntimeError("Observed non-zero model-ligand "
+ "state which must trigger certain "
+ "pair state. Point the developer of "
+ "the used scorer to this error message")
+
+ self._state_matrix[target_id, model_id] = pair_state
+ self._target_ligand_states[target_id] = target_ligand_state
+ self._model_ligand_states[model_id] = model_ligand_state
+ if pair_state == 0:
if score is None or np.isnan(score):
raise RuntimeError("LigandScorer returned invalid "
"score despite 0 error state")
@@ -784,7 +1034,8 @@ def _ResidueToGraph(residue, by_atom_index=False):
def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
by_atom_index=False, return_symmetries=True,
- max_symmetries=1e6):
+ max_symmetries=1e6, model_graph = None,
+ target_graph = None):
"""Return a list of symmetries (isomorphisms) of the model onto the target
residues.
@@ -820,11 +1071,18 @@ def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
"""
# Get the Graphs of the ligands
- model_graph = _ResidueToGraph(model_ligand, by_atom_index=by_atom_index)
- target_graph = _ResidueToGraph(target_ligand, by_atom_index=by_atom_index)
+ if model_graph is None:
+ model_graph = _ResidueToGraph(model_ligand,
+ by_atom_index=by_atom_index)
if not networkx.is_connected(model_graph):
raise DisconnectedGraphError("Disconnected graph for model ligand %s" % model_ligand)
+
+
+ if target_graph is None:
+ target_graph = _ResidueToGraph(target_ligand,
+ by_atom_index=by_atom_index)
+
if not networkx.is_connected(target_graph):
raise DisconnectedGraphError("Disconnected graph for target ligand %s" % target_ligand)
diff --git a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
index 7ddd8f5eea7676a2ca4706b86f473f51764cfd63..5566253611d860f3adccf2cf3b05cb5ac0f3b9ba 100644
--- a/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
+++ b/modules/mol/alg/pymod/ligand_scoring_lddtpli.py
@@ -153,18 +153,23 @@ class LDDTPLIScorer(ligand_scoring_base.LigandScorer):
target_ligand,
model_ligand)
- state = 0
+ pair_state = 0
score = result["lddt_pli"]
if score is None or np.isnan(score):
if result["lddt_pli_n_contacts"] == 0:
# it's a space ship!
- state = 10
+ pair_state = 10
else:
# unknwon error state
- state = 20
+ pair_state = 20
- return (score, state, result)
+ # the ligands get a zero-state...
+ target_ligand_state = 0
+ model_ligand_state = 0
+
+ return (score, pair_state, target_ligand_state, model_ligand_state,
+ result)
def _score_dir(self):
""" Implements interface from parent
diff --git a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
index e300301e368ad4144ef44f91db0bb9678dce9ab4..19c51a5da134fe9ea03a8ae0f769f0fda73f5b03 100644
--- a/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
+++ b/modules/mol/alg/pymod/ligand_scoring_scrmsd.py
@@ -193,19 +193,23 @@ class SCRMSDScorer(ligand_scoring_base.LigandScorer):
"transform": r.transform,
"inconsistent_residues": r.inconsistent_residues}
+ target_ligand_state = 0
+ model_ligand_state = 0
+ pair_state = 0
+
if best_rmsd_result["rmsd"] is not None:
best_rmsd = best_rmsd_result["rmsd"]
- error_state = 0
else:
# try to identify error states
best_rmsd = np.nan
error_state = 20 # unknown error
if self._get_target_binding_site(target_ligand).GetResidueCount() == 0:
- error_state = 10 # binding_site
+ pair_state = 6 # binding_site
+ target_ligand_state = 10
elif len(representations) == 0:
- error_state = 11 # model_representation
+ pair_state = 11
- return (best_rmsd, error_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
diff --git a/modules/mol/alg/tests/test_ligand_scoring.py b/modules/mol/alg/tests/test_ligand_scoring.py
index 1bd7971c9237fa3264b38a9c03fd3bbd6b38cb84..55d1cb2665638b7c71dadc210cca9abd6913f97d 100644
--- a/modules/mol/alg/tests/test_ligand_scoring.py
+++ b/modules/mol/alg/tests/test_ligand_scoring.py
@@ -507,10 +507,11 @@ class TestLigandScoring(unittest.TestCase):
mdl_ed.UpdateICS()
trg_ed.UpdateICS()
+ print("1")
sc = LigandScorer(mdl, trg, None, None, unassigned=True,
full_bs_search=True,
add_mdl_contacts=False)
-
+ print("done")
# Check unassigned targets
# NA: not in contact with target
trg_na = sc.target.FindResidue("L_NA", 1)
diff --git a/modules/mol/alg/tests/test_ligand_scoring_fancy.py b/modules/mol/alg/tests/test_ligand_scoring_fancy.py
index 15000f5188772a96f7d8abd89c6a944481621f6e..934ef783407872ce57f5e4785505367a3d062f11 100644
--- a/modules/mol/alg/tests/test_ligand_scoring_fancy.py
+++ b/modules/mol/alg/tests/test_ligand_scoring_fancy.py
@@ -629,13 +629,18 @@ class TestLigandScoringFancy(unittest.TestCase):
self.assertEqual(sc.unassigned_target_ligands, [0])
self.assertEqual(sc.unassigned_model_ligands, [0])
- trg_report = sc.get_target_ligand_state_report(0)
- mdl_report = sc.get_model_ligand_state_report(0)
+ trg_report, trg_pair_report = sc.get_target_ligand_state_report(0)
+ mdl_report, mdl_pair_report = sc.get_model_ligand_state_report(0)
- self.assertEqual(len(trg_report), 1)
- self.assertEqual(len(mdl_report), 1)
- self.assertEqual(trg_report[0]["short desc"], "symmetries")
- self.assertEqual(mdl_report[0]["short desc"], "symmetries")
+ # the individual ligands are OK
+ self.assertEqual(trg_report["short desc"], "OK")
+ self.assertEqual(mdl_report["short desc"], "OK")
+
+ # but there are too many symmetries
+ self.assertEqual(len(trg_pair_report), 1)
+ self.assertEqual(len(mdl_pair_report), 1)
+ self.assertEqual(trg_pair_report[0]["short desc"], "symmetries")
+ self.assertEqual(mdl_pair_report[0]["short desc"], "symmetries")
def test_no_binding_site(self):
"""
@@ -665,18 +670,24 @@ class TestLigandScoringFancy(unittest.TestCase):
self.assertTrue(np.isnan(sc.score_matrix[0, 3]))
- trg_report = sc.get_target_ligand_state_report(0)
+ trg_report, trg_pair_report = sc.get_target_ligand_state_report(0)
+
+ exp_lig_report = {'state': 10.0,
+ 'short desc': 'binding_site',
+ 'desc': 'No residues were in proximity of the target ligand.'}
- expected_report = [{'state': 1, 'short desc': 'identity',
+ exp_pair_report = [{'state': 1, 'short desc': 'identity',
'desc': 'Ligands could not be matched (by subgraph isomorphism)',
'indices': [0, 1, 2, 4, 5, 6]},
- {'state': 10, 'short desc': 'binding_site',
- 'desc': 'No residues were in proximity of the target ligand.',
+ {'state': 6, 'short desc': 'single_ligand_issue',
+ 'desc': 'Cannot compute valid pairwise score as either model or target ligand have non-zero state.',
'indices': [3]}]
- # order of report is fix
- self.assertDictEqual(trg_report[0], expected_report[0])
- self.assertDictEqual(trg_report[1], expected_report[1])
+ # order of report is fix
+ self.assertDictEqual(trg_report, exp_lig_report)
+ self.assertDictEqual(trg_pair_report[0], exp_pair_report[0])
+ self.assertDictEqual(trg_pair_report[1], exp_pair_report[1])
+
def test_no_lddt_pli_contact(self):
"""
@@ -717,12 +728,12 @@ class TestLigandScoringFancy(unittest.TestCase):
'indices': [0, 1]}]
# both target ligands are expected to have the same report => no_contact
- report_1 = sc.get_target_ligand_state_report(0)
- report_2 = sc.get_target_ligand_state_report(1)
- self.assertEqual(len(report_1), 1)
- self.assertEqual(len(report_2), 1)
- self.assertDictEqual(report_1[0], expected_report[0])
- self.assertDictEqual(report_2[0], expected_report[0])
+ report_1, report_pair_1 = sc.get_target_ligand_state_report(0)
+ report_2, report_pair_2 = sc.get_target_ligand_state_report(1)
+ self.assertEqual(len(report_pair_1), 1)
+ self.assertEqual(len(report_pair_2), 1)
+ self.assertDictEqual(report_pair_1[0], expected_report[0])
+ self.assertDictEqual(report_pair_2[0], expected_report[0])
# However RMSD should be assigned
@@ -742,14 +753,190 @@ class TestLigandScoringFancy(unittest.TestCase):
'indices': [0, 1]}]
# both target ligands are expected to have the same report => no_contact
- report_1 = sc.get_target_ligand_state_report(0)
- report_2 = sc.get_target_ligand_state_report(1)
- self.assertEqual(len(report_1), 1)
- self.assertEqual(len(report_2), 1)
- self.assertDictEqual(report_1[0], expected_report[0])
- self.assertDictEqual(report_2[0], expected_report[0])
-
+ report_1, report_pair_1 = sc.get_target_ligand_state_report(0)
+ report_2, report_pair_2 = sc.get_target_ligand_state_report(1)
+ self.assertEqual(len(report_pair_1), 1)
+ self.assertEqual(len(report_pair_2), 1)
+ self.assertDictEqual(report_pair_1[0], expected_report[0])
+ self.assertDictEqual(report_pair_2[0], expected_report[0])
+
+ def test_unassigned_reasons(self):
+ """Test reasons for being unassigned."""
+ trg = _LoadMMCIF("1r8q.cif.gz")
+ mdl = _LoadMMCIF("P84080_model_02.cif.gz")
+ def _AppendResidueWithBonds(ed, chain, old_res):
+ new_res = ed.AppendResidue(chain, old_res.name)
+ for old_atom in old_res.atoms:
+ ed.InsertAtom(new_res, old_atom.name, old_atom.pos, old_atom.element,
+ old_atom.occupancy, old_atom.b_factor, old_atom.is_hetatom)
+ for old_bond in old_atom.bonds:
+ new_first = new_res.FindAtom(old_bond.first.name)
+ new_second = new_res.FindAtom(old_bond.second.name)
+ ed.Connect(new_first, new_second)
+ return new_res
+
+ # Add interesting ligands to model and target
+ mdl_ed = mdl.EditXCS()
+ trg_ed = trg.EditXCS()
+
+ # Add ZN: representation in the model (chain missing in model)
+ new_chain = mdl_ed.InsertChain("L_ZN")
+ mdl_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = _AppendResidueWithBonds(mdl_ed, new_chain, trg.Select("rname=ZN").residues[0].handle)
+ new_res.is_ligand = True
+
+ # Add NA: not in contact with target
+ new_chain = trg_ed.InsertChain("L_NA")
+ trg_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = trg_ed.AppendResidue(new_chain, "NA")
+ new_atom = trg_ed.InsertAtom(new_res, "NA", geom.Vec3(100, 100, 100), "NA")
+ new_res.is_ligand = True
+ new_chain = mdl_ed.InsertChain("L_NA")
+ mdl_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = mdl_ed.AppendResidue(new_chain, "NA")
+ new_atom = mdl_ed.InsertAtom(new_res, "NA", geom.Vec3(100, 100, 100), "NA")
+ new_res.is_ligand = True
+
+ # Add OXY: no symmetry/ not identical -
+ new_chain = mdl_ed.InsertChain("L_OXY")
+ mdl_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = mdl_ed.AppendResidue(new_chain, "OXY")
+ new_atom1 = mdl_ed.InsertAtom(new_res, "O1", geom.Vec3(0, 0, 0), "O")
+ new_atom2 = mdl_ed.InsertAtom(new_res, "O2", geom.Vec3(1, 1, 1), "O")
+ mdl_ed.Connect(new_atom1, new_atom2)
+ new_res.is_ligand = True
+
+ # Add CMO: disconnected
+ new_chain = mdl_ed.InsertChain("L_CMO")
+ mdl_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = mdl_ed.AppendResidue(new_chain, "CMO")
+ new_atom1 = mdl_ed.InsertAtom(new_res, "O", geom.Vec3(0, 0, 0), "O")
+ new_atom2 = mdl_ed.InsertAtom(new_res, "C", geom.Vec3(1, 1, 1), "O")
+ new_res.is_ligand = True
+ new_chain = trg_ed.InsertChain("L_CMO")
+ trg_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = trg_ed.AppendResidue(new_chain, "CMO")
+ new_atom1 = trg_ed.InsertAtom(new_res, "O", geom.Vec3(0, 0, 0), "O")
+ new_atom2 = trg_ed.InsertAtom(new_res, "C", geom.Vec3(1, 1, 1), "O")
+ new_res.is_ligand = True
+
+ # Add 3 MG in model: assignment/stoichiometry
+ mg_pos = [
+ geom.Vec3(3.871, 12.343, 44.485),
+ geom.Vec3(3.871, 12.343, 44.485) + 1,
+ geom.Vec3(3.871, 12.343, 44.485) + 100
+ ]
+ for i in range(3):
+ new_chain = mdl_ed.InsertChain("L_MG_%d" % i)
+ mdl_ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
+ new_res = mdl_ed.AppendResidue(new_chain, "MG")
+ new_atom = mdl_ed.InsertAtom(new_res, "MG", mg_pos[i], "MG")
+ new_res.is_ligand = True
+
+ mdl_ed.UpdateICS()
+ trg_ed.UpdateICS()
+
+ sc = ligand_scoring_lddtpli.LDDTPLIScorer(mdl, trg, None, None)
+
+ # Check unassigned targets
+ # NA: not in contact with target
+ trg_na = sc.target.FindResidue("L_NA", 1)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["L_NA"][1], "no_contact")
+ # AFB: not identical to anything in the model
+ trg_afb = sc.target.FindResidue("G", 1)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["G"][1], "identity")
+ # F.G3D1: J.G3D1 assigned instead
+ trg_fg3d1 = sc.target.FindResidue("F", 1)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["F"][1], "stoichiometry")
+ # CMO: disconnected
+ trg_cmo1 = sc.target.FindResidue("L_CMO", 1)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["L_CMO"][1], "disconnected")
+ # J.G3D1: assigned to L_2.G3D1 => check if it is assigned
+ self.assertTrue(5 not in sc.unassigned_target_ligands)
+ self.assertNotIn("J", sc.unassigned_target_ligands_reasons)
+
+ # Check unassigned models
+ # OXY: not identical to anything in the model
+ mdl_oxy = sc.model.FindResidue("L_OXY", 1)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_OXY"][1], "identity")
+ self.assertTrue("L_OXY" not in sc.score)
+ # NA: not in contact with target
+ mdl_na = sc.model.FindResidue("L_NA", 1)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_NA"][1], "no_contact")
+ self.assertTrue("L_NA" not in sc.score)
+
+ # MG in L_MG_2 has stupid coordinates and is not assigned
+ mdl_mg_2 = sc.model.FindResidue("L_MG_2", 1)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_MG_2"][1], "stoichiometry")
+ self.assertTrue("L_MG_2" not in sc.score)
+
+ self.assertNotIn("L_MG_0", sc.unassigned_model_ligands_reasons)
+ # CMO: disconnected
+ mdl_cmo1 = sc.model.FindResidue("L_CMO", 1)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_CMO"][1], "disconnected")
+
+ # Should work with rmsd_assignment too
+ sc = ligand_scoring_scrmsd.SCRMSDScorer(mdl, trg, None, None,
+ full_bs_search=True)
+
+
+
+ self.assertDictEqual(sc.unassigned_model_ligands_reasons, {
+ 'L_ZN': {1: 'model_representation'},
+ 'L_NA': {1: 'binding_site'},
+ 'L_OXY': {1: 'identity'},
+ 'L_MG_2': {1: 'stoichiometry'},
+ "L_CMO": {1: 'disconnected'}
+ })
+ self.assertDictEqual(sc.unassigned_target_ligands_reasons, {
+ 'G': {1: 'identity'},
+ 'H': {1: 'model_representation'},
+ 'J': {1: 'stoichiometry'},
+ 'K': {1: 'identity'},
+ 'L_NA': {1: 'binding_site'},
+ "L_CMO": {1: 'disconnected'}
+ })
+ self.assertTrue("L_OXY" not in sc.score)
+
+ # With missing ligands
+ sc = ligand_scoring_lddtpli.LDDTPLIScorer(mdl.Select("cname=A"), trg, None, None)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["E"][1], 'no_ligand')
+
+ sc = ligand_scoring_lddtpli.LDDTPLIScorer(mdl, trg.Select("cname=A"), None, None)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_2"][1], 'no_ligand')
+
+ sc = ligand_scoring_scrmsd.SCRMSDScorer(mdl.Select("cname=A"), trg, None, None)
+ self.assertEqual(sc.unassigned_target_ligands_reasons["E"][1], 'no_ligand')
+
+ sc = ligand_scoring_scrmsd.SCRMSDScorer(mdl, trg.Select("cname=A"), None, None)
+ self.assertEqual(sc.unassigned_model_ligands_reasons["L_2"][1], 'no_ligand')
+
+ # However not everything must be missing
+ with self.assertRaises(ValueError):
+ sc = LigandScorer(mdl.Select("cname=A"), trg.Select("cname=A"), None, None)
+
+ # Test with partial bs search (full_bs_search=True)
+ # Here we expect L_MG_2 to be unassigned because of stoichiometry
+ # rather than model_representation, as it no longer matters so far from
+ # the binding site.
+ sc = ligand_scoring_scrmsd.SCRMSDScorer(mdl, trg, None, None,
+ full_bs_search=True)
+ self.assertEqual(sc.unassigned_model_ligands_reasons, {
+ 'L_ZN': {1: 'model_representation'},
+ 'L_NA': {1: 'binding_site'},
+ 'L_OXY': {1: 'identity'},
+ 'L_MG_2': {1: 'stoichiometry'},
+ "L_CMO": {1: 'disconnected'}
+ })
+ self.assertEqual(sc.unassigned_target_ligands_reasons, {
+ 'G': {1: 'identity'},
+ 'H': {1: 'model_representation'},
+ 'J': {1: 'stoichiometry'},
+ 'K': {1: 'identity'},
+ 'L_NA': {1: 'binding_site'},
+ "L_CMO": {1: 'disconnected'}
+ })
if __name__ == "__main__":
from ost import testutils