diff --git a/modules/mol/alg/pymod/ligand_scoring.py b/modules/mol/alg/pymod/ligand_scoring.py
index 84d1291ac67f5e6e3b91d770f0db7ddbf0d325d9..db2e1d307c6cffc8be4de9be3a29310659a1c4d3 100644
--- a/modules/mol/alg/pymod/ligand_scoring.py
+++ b/modules/mol/alg/pymod/ligand_scoring.py
@@ -4,6 +4,7 @@ import numpy as np
import networkx
from ost import mol
+from ost import geom
from ost import LogError, LogWarning, LogScript, LogInfo, LogVerbose, LogDebug
from ost.mol.alg import chain_mapping
@@ -303,4 +304,142 @@ def ResidueToGraph(residue):
return nxg
-__all__ = ["LigandScorer", "ResidueToGraph"]
+def LigandRMSD(model_ligand, target_ligand, transformation=geom.Mat4(),
+ substructure_match=False):
+ """Calculate symmetry-corrected RMSD.
+
+ Binding site superposition must be computed separately and passed as
+ `transformation`.
+
+ :param model_ligand: The model ligand
+ :type model_ligand: :class:`ost.mol.ResidueHandle` or
+ :class:`ost.mol.ResidueView`
+ :param target_ligand: The target ligand
+ :type target_ligand: :class:`ost.mol.ResidueHandle` or
+ :class:`ost.mol.ResidueView`
+ :param transformation: Optional transformation to apply on each atom
+ position of model_ligand.
+ :type transformation: :class:`ost.geom.Mat4`
+ :param substructure_match: Set this to True to allow partial target
+ ligand.
+ :type substructure_match: :class:`bool`
+ :rtype: :class:`float`
+ :raises: :class:`NoSymmetryError` when no symmetry can be found.
+ """
+
+ symmetries = _ComputeSymmetries(model_ligand, target_ligand,
+ substructure_match=substructure_match,
+ by_atom_index=True)
+
+ best_rmsd = float('inf')
+ for i, (trg_sym, mdl_sym) in enumerate(symmetries):
+ # Prepare Entities for RMSD
+ trg_lig_rmsd_ent = mol.CreateEntity()
+ trg_lig_rmsd_editor = trg_lig_rmsd_ent.EditXCS()
+ trg_lig_rmsd_chain = trg_lig_rmsd_editor.InsertChain("_")
+ trg_lig_rmsd_res = trg_lig_rmsd_editor.AppendResidue(trg_lig_rmsd_chain, "LIG")
+
+ mdl_lig_rmsd_ent = mol.CreateEntity()
+ mdl_lig_rmsd_editor = mdl_lig_rmsd_ent.EditXCS()
+ mdl_lig_rmsd_chain = mdl_lig_rmsd_editor.InsertChain("_")
+ mdl_lig_rmsd_res = mdl_lig_rmsd_editor.AppendResidue(mdl_lig_rmsd_chain, "LIG")
+
+ for mdl_anum, trg_anum in zip(mdl_sym, trg_sym):
+ # Rename model atoms according to symmetry
+ trg_atom = target_ligand.atoms[trg_anum]
+ mdl_atom = model_ligand.atoms[mdl_anum]
+ # Add atoms in the correct order to the RMSD entities
+ trg_lig_rmsd_editor.InsertAtom(trg_lig_rmsd_res, trg_atom.name, trg_atom.pos)
+ mdl_lig_rmsd_editor.InsertAtom(mdl_lig_rmsd_res, mdl_atom.name, mdl_atom.pos)
+
+ trg_lig_rmsd_editor.UpdateICS()
+ mdl_lig_rmsd_editor.UpdateICS()
+
+ rmsd = mol.alg.CalculateRMSD(mdl_lig_rmsd_ent.CreateFullView(),
+ trg_lig_rmsd_ent.CreateFullView(),
+ transformation)
+ if rmsd < best_rmsd:
+ best_rmsd = rmsd
+
+ return best_rmsd
+
+
+def _ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
+ by_atom_index=False):
+ """Return a list of symmetries (isomorphisms) of the model onto the target
+ residues.
+
+ :param model_ligand: The model ligand
+ :type model_ligand: :class:`ost.mol.ResidueHandle` or
+ :class:`ost.mol.ResidueView`
+ :param target_ligand: The target ligand
+ :type target_ligand: :class:`ost.mol.ResidueHandle` or
+ :class:`ost.mol.ResidueView`
+ :param substructure_match: Set this to True to allow partial ligands
+ in the reference.
+ :type substructure_match: :class:`bool`
+ :param by_atom_index: Set this parameter to True if you need the symmetries
+ to refer to atom index (within the residue).
+ Otherwise, if False, the symmetries refer to atom
+ names.
+ :raises: :class:`NoSymmetryError` when no symmetry can be found.
+
+ """
+
+ # Get the Graphs of the ligands
+ model_graph = ResidueToGraph(model_ligand)
+ target_graph = ResidueToGraph(target_ligand)
+
+ if by_atom_index:
+ networkx.relabel_nodes(model_graph,
+ {a: b for a, b in zip(
+ [a.name for a in model_ligand.atoms],
+ range(len(model_ligand.atoms)))},
+ False)
+ networkx.relabel_nodes(target_graph,
+ {a: b for a, b in zip(
+ [a.name for a in target_ligand.atoms],
+ range(len(target_ligand.atoms)))},
+ False)
+
+ # 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.
+ # Make sure to generate the symmetries correctly in the end
+ GM = networkx.algorithms.isomorphism.GraphMatcher(
+ model_graph, target_graph, node_match=lambda x, y:
+ x["element"] == y["element"])
+ if GM.is_isomorphic():
+ symmetries = [
+ (list(isomorphism.values()), list(isomorphism.keys()))
+ for isomorphism in GM.isomorphisms_iter()]
+ assert len(symmetries) > 0
+ LogDebug("Found %s isomorphic mappings (symmetries)" % len(symmetries))
+ elif GM.subgraph_is_isomorphic() and substructure_match:
+ symmetries = [(list(isomorphism.values()), list(isomorphism.keys())) for isomorphism in
+ GM.subgraph_isomorphisms_iter()]
+ 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))
+ elif GM.subgraph_is_isomorphic():
+ LogDebug("Found subgraph isomorphisms (symmetries), but"
+ " ignoring because substructure_match=False")
+ raise NoSymmetryError("No symmetry between %s and %s" % (
+ str(model_ligand), str(target_ligand)))
+ else:
+ LogDebug("Found no isomorphic mappings (symmetries)")
+ raise NoSymmetryError("No symmetry between %s and %s" % (
+ str(model_ligand), str(target_ligand)))
+
+ return symmetries
+
+
+class NoSymmetryError(Exception):
+ """Exception to be raised when no symmetry can be found
+
+ Those are cases we might want to capture for default behavior.
+ """
+ pass
+
+__all__ = ["LigandScorer", "ResidueToGraph", "LigandRMSD", "NoSymmetryError"]
diff --git a/modules/mol/alg/tests/test_ligand_scoring.py b/modules/mol/alg/tests/test_ligand_scoring.py
index baf8b2bf6244dc23ea038f9aebc2093aeaad7736..98eda28319a212684dc94efe4e787a919e50f157 100644
--- a/modules/mol/alg/tests/test_ligand_scoring.py
+++ b/modules/mol/alg/tests/test_ligand_scoring.py
@@ -1,9 +1,10 @@
import unittest, os, sys
-from ost import io, mol
+from ost import io, mol, geom
# check if we can import: fails if numpy or scipy not available
try:
from ost.mol.alg.ligand_scoring import *
+ import ost.mol.alg.ligand_scoring
except ImportError:
print("Failed to import ligand_scoring.py. Happens when numpy, scipy or "
"networkx is missing. Ignoring test_ligand_scoring.py tests.")
@@ -122,6 +123,86 @@ class TestLigandScoring(unittest.TestCase):
graph = ResidueToGraph(mdl_lig.residues[0])
assert len(graph.edges) == 34
assert len(graph.nodes) == 32
+ # Check an arbitrary node
+ assert [a for a in graph.adj["14"].keys()] == ["13", "29"]
+
+ def test__ComputeSymmetries(self):
+ """Test that _ComputeSymmetries works.
+ """
+ trg, trg_seqres = io.LoadMMCIF(os.path.join('testfiles', "1r8q.cif.gz"), seqres=True)
+ mdl, mdl_seqres = io.LoadMMCIF(os.path.join('testfiles', "P84080_model_02.cif.gz"), seqres=True)
+
+ trg_mg1 = trg.FindResidue("E", 1)
+ trg_g3d1 = trg.FindResidue("F", 1)
+ trg_afb1 = trg.FindResidue("G", 1)
+ trg_g3d2 = trg.FindResidue("J", 1)
+ mdl_g3d = mdl.FindResidue("L_2", 1)
+
+ sym = ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d, trg_g3d1)
+ assert len(sym) == 72
+
+ sym = ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d, trg_g3d1, by_atom_index=True)
+ assert len(sym) == 72
+
+ # Test that it works with views and only consider atoms in the view
+ # Skip PA, PB and O[1-3]A and O[1-3]B in target and model
+ # We assume atom index are fixed and won't change
+ trg_g3d1_sub = trg_g3d1.Select("aindex>6019").residues[0]
+ mdl_g3d_sub = mdl_g3d.Select("aindex>1447").residues[0]
+
+ sym = ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d_sub, trg_g3d1_sub)
+ assert len(sym) == 6
+
+ sym = ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d_sub, trg_g3d1_sub, by_atom_index=True)
+ assert len(sym) == 6
+
+ # Substructure matches
+ self.skipTest("Substructure matches don't work yet")
+ sym = ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d, trg_g3d1_sub, substructure_match=True)
+ assert len(sym) == 6
+
+ # Missing atoms only allowed in target, not in model
+ with self.assertRaises(NoSymmetryError):
+ ost.mol.alg.ligand_scoring._ComputeSymmetries(mdl_g3d_sub, trg_g3d1, substructure_match=True)
+
+ def test_LigandRMSD(self):
+ """Test that LigandRMSD works.
+ """
+ trg, trg_seqres = io.LoadMMCIF(os.path.join('testfiles', "1r8q.cif.gz"), seqres=True)
+ mdl, mdl_seqres = io.LoadMMCIF(os.path.join('testfiles', "P84080_model_02.cif.gz"), seqres=True)
+
+ trg_mg1 = trg.FindResidue("E", 1)
+ trg_g3d1 = trg.FindResidue("F", 1)
+ trg_afb1 = trg.FindResidue("G", 1)
+ trg_g3d2 = trg.FindResidue("J", 1)
+ mdl_g3d = mdl.FindResidue("L_2", 1)
+
+ rmsd = LigandRMSD(mdl_g3d, trg_g3d1)
+ self.assertAlmostEqual(rmsd, 2.21341e-06, 10)
+ rmsd = LigandRMSD(mdl_g3d, trg_g3d2)
+ self.assertAlmostEqual(rmsd, 61.21325, 4)
+
+ # Ensure we raise a NoSymmetryError if the ligand is wrong
+ with self.assertRaises(NoSymmetryError):
+ LigandRMSD(mdl_g3d, trg_mg1)
+ with self.assertRaises(NoSymmetryError):
+ LigandRMSD(mdl_g3d, trg_afb1)
+
+ # Assert that transform works
+ trans = geom.Mat4(-0.999256, 0.00788487, -0.0377333, -15.4397,
+ 0.0380652, 0.0473315, -0.998154, 29.9477,
+ -0.00608426, -0.998848, -0.0475963, 28.8251,
+ 0, 0, 0, 1)
+ rmsd = LigandRMSD(mdl_g3d, trg_g3d2, transformation=trans)
+ self.assertAlmostEqual(rmsd, 0.293972, 5)
+
+ # Assert that substructure matches work
+ self.skipTest("Substructure matches don't work yet")
+ trg_g3d1_sub = trg_g3d1.Select("aindex>6019").residues[0] # Skip PA, PB and O[1-3]A and O[1-3]B.
+ mdl_g3d_sub = mdl_g3d.Select("aindex>1447").residues[0] # Skip PA, PB and O[1-3]A and O[1-3]B.
+ with self.assertRaises(NoSymmetryError):
+ LigandRMSD(mdl_g3d, trg_g3d1_sub) # no full match
+ rmsd = LigandRMSD(mdl_g3d, trg_g3d1_sub, substructure_match=True)
if __name__ == "__main__":