diff --git a/modules/mol/alg/pymod/ligand_scoring.py b/modules/mol/alg/pymod/ligand_scoring.py
index 25bea3916daebb66386f9e61b6e9e10253db3354..594492045052277e72df761e82e7c5f9d431196c 100644
--- a/modules/mol/alg/pymod/ligand_scoring.py
+++ b/modules/mol/alg/pymod/ligand_scoring.py
@@ -856,16 +856,18 @@ class LigandScorer:
if mdl_a_hash_code in bs_atom_mapping:
trg_bs_indices.append(bs_atom_mapping[mdl_a_hash_code])
elif mdl_a_hash_code not in mdl_lig_hashes:
- at_key = (a.GetResidue().GetNumber(), a.name)
- cname = a.GetChain().name
- cname_key = (flat_mapping[cname], cname)
- if at_key in self.mappable_atoms[cname_key]:
- # Its a contact in the model but not part of trg_bs.
- # It can still be mapped using the global
- # mdl_ch/ref_ch alignment
- # dist in ref > self.lddt_pli_radius + max_thresh
- # => guaranteed to be non-fulfilled contact
- n_penalty += 1
+ if a.GetChain().GetName() in flat_mapping:
+ # Its in a mapped chain
+ at_key = (a.GetResidue().GetNumber(), a.name)
+ cname = a.GetChain().name
+ cname_key = (flat_mapping[cname], cname)
+ if at_key in self.mappable_atoms[cname_key]:
+ # Its a contact in the model but not part of
+ # trg_bs. It can still be mapped using the
+ # global mdl_ch/ref_ch alignment
+ # d in ref > self.lddt_pli_radius + max_thresh
+ # => guaranteed to be non-fulfilled contact
+ n_penalty += 1
n_penalties.append(n_penalty)
@@ -964,20 +966,20 @@ class LigandScorer:
if mdl_ch not in lddt_chain_mapping:
# check which chain in trg is closest
chem_group_idx = None
- for i, m in self.chem_mapping:
+ for i, m in enumerate(self.chem_mapping):
if mdl_ch in m:
chem_group_idx = i
break
if chem_group_idx is None:
raise RuntimeError("This should never happen... "
"ask Gabriel...")
- mdl_ch = self.chain_mapping_mdl.FindChain(mdl_ch)
- mdl_center = mdl_ch.geometric_center
+ mdl_ch_view = self.chain_mapping_mdl.FindChain(mdl_ch)
+ mdl_center = mdl_ch_view.geometric_center
closest_ch = None
closest_dist = None
- for trg_ch in self.chem_groups[chem_group_idx]:
- if trg_ch not in lddt_mapping.values():
- c = self.target.FindChain(trg_ch).geometric_center
+ for trg_ch in self.chain_mapper.chem_groups[chem_group_idx]:
+ if trg_ch not in lddt_chain_mapping.values():
+ c = self.chain_mapper.target.FindChain(trg_ch).geometric_center
d = geom.Distance(mdl_center, c)
if closest_dist is None or d < closest_dist:
closest_dist = d
diff --git a/modules/mol/alg/tests/test_ligand_scoring.py b/modules/mol/alg/tests/test_ligand_scoring.py
index e0e87fe392f5df97be7dfc1dda1d3600a8ed1c3d..c60c7cba15d3955e28f8ace0a1fb23f5098a0f86 100644
--- a/modules/mol/alg/tests/test_ligand_scoring.py
+++ b/modules/mol/alg/tests/test_ligand_scoring.py
@@ -371,7 +371,7 @@ class TestLigandScoring(unittest.TestCase):
# lddt_pli
self.assertAlmostEqual(sc.lddt_pli["J"][mol.ResNum(1)], 0.9127105666156202, 5)
- self.assertAlmostEqual(sc.lddt_pli["F"][mol.ResNum(1)], 0.915929203539823, 6)
+ self.assertAlmostEqual(sc.lddt_pli["F"][mol.ResNum(1)], 0.915929203539823, 5)
# lddt_pli_details
self.assertEqual(sc.lddt_pli_details["J"][mol.ResNum(1)]["lddt_pli_n_contacts"], 653)
self.assertEqual(len(sc.lddt_pli_details["J"][mol.ResNum(1)]["bs_ref_res"]), 15)
@@ -382,6 +382,21 @@ class TestLigandScoring(unittest.TestCase):
self.assertEqual(sc.lddt_pli_details["F"][mol.ResNum(1)]["target_ligand"].qualified_name, 'K.G3D1')
self.assertEqual(sc.lddt_pli_details["F"][mol.ResNum(1)]["model_ligand"].qualified_name, 'F.G3D1')
+ # lddt_pli with added mdl contacts
+ sc = LigandScorer(trg, trg_4c0a, None, None, check_resnames=False,
+ add_mdl_contacts=True)
+ self.assertAlmostEqual(sc.lddt_pli["J"][mol.ResNum(1)], 0.8988340192043895, 5)
+ self.assertAlmostEqual(sc.lddt_pli["F"][mol.ResNum(1)], 0.9039735099337749, 5)
+ # lddt_pli_details
+ self.assertEqual(sc.lddt_pli_details["J"][mol.ResNum(1)]["lddt_pli_n_contacts"], 729)
+ self.assertEqual(len(sc.lddt_pli_details["J"][mol.ResNum(1)]["bs_ref_res"]), 63)
+ self.assertEqual(sc.lddt_pli_details["J"][mol.ResNum(1)]["target_ligand"].qualified_name, 'I.G3D1')
+ self.assertEqual(sc.lddt_pli_details["J"][mol.ResNum(1)]["model_ligand"].qualified_name, 'J.G3D1')
+ self.assertEqual(sc.lddt_pli_details["F"][mol.ResNum(1)]["lddt_pli_n_contacts"], 755)
+ self.assertEqual(len(sc.lddt_pli_details["F"][mol.ResNum(1)]["bs_ref_res"]), 62)
+ self.assertEqual(sc.lddt_pli_details["F"][mol.ResNum(1)]["target_ligand"].qualified_name, 'K.G3D1')
+ self.assertEqual(sc.lddt_pli_details["F"][mol.ResNum(1)]["model_ligand"].qualified_name, 'F.G3D1')
+
def test_rmsd_assignment(self):
"""Test that the RMSD-based assignment works.
@@ -775,6 +790,86 @@ class TestLigandScoring(unittest.TestCase):
self.assertAlmostEqual(sc2.rmsd['00001_FE'][1], 2.1703, 4)
self.assertAlmostEqual(sc2.rmsd['00001_FE_2'][1], 2.1703, 4)
+ def test_added_mdl_contacts(self):
+
+ # binding site for ligand in chain G consists of chains A and B
+ prot = _LoadMMCIF("1r8q.cif.gz")
+
+ # model has the full binding site
+ mdl = mol.CreateEntityFromView(prot.Select("cname=A,B,G"), True)
+
+ # chain C has same sequence as chain A but is not in contact
+ # with ligand in chain G
+ # target has thus incomplete binding site only from chain B
+ trg = mol.CreateEntityFromView(prot.Select("cname=B,C,G"), True)
+
+ # if added model contacts are not considered, the incomplete binding
+ # site only from chain B is perfectly reproduced by model which also has
+ # chain B
+ sc = LigandScorer(mdl, trg, add_mdl_contacts=False)
+ self.assertAlmostEqual(sc.lddt_pli_matrix[0,0], 1.0, 5)
+
+ # if added model contacts are considered, contributions from chain B are
+ # perfectly reproduced but all contacts of ligand towards chain A are
+ # added as penalty
+ sc = LigandScorer(mdl, trg, add_mdl_contacts=True)
+
+ lig = prot.Select("cname=G")
+ A_count = 0
+ B_count = 0
+ for a in lig.atoms:
+ close_atoms = mdl.FindWithin(a.GetPos(), sc.lddt_pli_radius)
+ for ca in close_atoms:
+ cname = ca.GetChain().GetName()
+ if cname == "G":
+ pass # its a ligand atom...
+ elif cname == "A":
+ A_count += 1
+ elif cname == "B":
+ B_count += 1
+
+ self.assertAlmostEqual(sc.lddt_pli_matrix[0,0],
+ B_count/(A_count + B_count), 5)
+
+ # Same as before but additionally we remove residue TRP.66
+ # from chain C in the target to test mapping magic...
+ # Chain C is NOT in contact with the ligand but we only
+ # add contacts from chain A as penalty that are mappable
+ # to the closest chain with same sequence. That would be
+ # chain C
+ query = "cname=B,G or (cname=C and rnum!=66)"
+ trg = mol.CreateEntityFromView(prot.Select(query), True)
+ sc = LigandScorer(mdl, trg, add_mdl_contacts=True)
+
+ TRP66_count = 0
+ for a in lig.atoms:
+ close_atoms = mdl.FindWithin(a.GetPos(), sc.lddt_pli_radius)
+ for ca in close_atoms:
+ cname = ca.GetChain().GetName()
+ if cname == "A" and ca.GetResidue().GetNumber().GetNum() == 66:
+ TRP66_count += 1
+
+ self.assertEqual(TRP66_count, 134)
+
+ # remove TRP66_count from original penalty
+ self.assertAlmostEqual(sc.lddt_pli_matrix[0,0],
+ B_count/(A_count + B_count - TRP66_count), 5)
+
+ # Move a random atom in the model from chain B towards the ligand center
+ # chain B is also present in the target and interacts with the ligand,
+ # but that atom would be far away and thus adds to the penalty. Since
+ # the ligand is small enough, the number of added contacts should be
+ # exactly the number of ligand atoms.
+ mdl_ed = mdl.EditXCS()
+ at = mdl.FindResidue("B", mol.ResNum(8)).FindAtom("NZ")
+ mdl_ed.SetAtomPos(at, lig.geometric_center)
+ sc = LigandScorer(mdl, trg, add_mdl_contacts=True)
+
+ # compared to the last assertAlmostEqual, we add the number of ligand
+ # atoms as additional penalties
+ self.assertAlmostEqual(sc.lddt_pli_matrix[0,0],
+ B_count/(A_count + B_count - TRP66_count + \
+ lig.GetAtomCount()), 5)
if __name__ == "__main__":
from ost import testutils