feat: SCHWED-5783 implement symmetry-corrected RMSD

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"]

modules/mol/alg/tests/test_ligand_scoring.py

 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__":