From b2947b24d374aa01bafe6e72d64f812e63739a55 Mon Sep 17 00:00:00 2001
From: Xavier Robin <xavier.robin@unibas.ch>
Date: Wed, 26 Jul 2023 13:17:02 +0200
Subject: [PATCH] feat: SCHWED-5954 report disconnected ligands.

Previous behavior was to raise an error if any ligand was disconnected.
However, now with the reporting it is more consistent to complete and
report why no score was computed.
---
 modules/mol/alg/pymod/ligand_scoring.py      | 37 +++++++++++++++-----
 modules/mol/alg/tests/test_ligand_scoring.py | 36 +++++++++++++++----
 2 files changed, 58 insertions(+), 15 deletions(-)

diff --git a/modules/mol/alg/pymod/ligand_scoring.py b/modules/mol/alg/pymod/ligand_scoring.py
index c06a99959..2a462b78d 100644
--- a/modules/mol/alg/pymod/ligand_scoring.py
+++ b/modules/mol/alg/pymod/ligand_scoring.py
@@ -652,6 +652,9 @@ class LigandScorer:
                             str(model_ligand), str(target_ligand)))
                         self._assignment_isomorphisms[target_i, model_i] = 0
                         continue
+                    except DisconnectedGraphError:
+                        # Disconnected graph is handled elsewhere
+                        continue
                     substructure_match = len(symmetries[0][0]) != len(
                         model_ligand.atoms)
 
@@ -1258,7 +1261,8 @@ class LigandScorer:
 
         Currently, the following reasons are reported:
 
-        * `no_ligand`: No ligand in the model.
+        * `no_ligand`: no ligand in the model.
+        * `disconnected`: ligand graph is disconnected.
         * `binding_site`: no residue in proximity of the target ligand.
         * `model_representation`: no representation of the reference binding
           site was found in the model
@@ -1291,7 +1295,8 @@ class LigandScorer:
 
         Currently, the following reasons are reported:
 
-        * `no_ligand`: No ligand in the target.
+        * `no_ligand`: no ligand in the target.
+        * `disconnected`: ligand graph is disconnected.
         * `binding_site`: no residue in proximity of the target ligand.
         * `model_representation`: no representation of the reference binding
           site was found in the model
@@ -1425,6 +1430,11 @@ class LigandScorer:
         if len(self.target_ligands) == 0:
             return ("no_ligand", "No ligand in the target")
 
+        # Is the ligand disconnected?
+        graph = _ResidueToGraph(ligand)
+        if not networkx.is_connected(graph):
+            return ("disconnected", "Ligand graph is disconnected")
+
         # Do we have isomorphisms with the target?
         for trg_lig_idx, assigned in enumerate(self._assignment_isomorphisms[:, ligand_idx]):
             if np.isnan(assigned):
@@ -1435,7 +1445,7 @@ class LigandScorer:
                         substructure_match=self.substructure_match,
                         by_atom_index=True,
                         return_symmetries=False)
-                except NoSymmetryError:
+                except (NoSymmetryError, DisconnectedGraphError):
                     assigned = 0.
                 else:
                     assigned = 1.
@@ -1481,6 +1491,11 @@ class LigandScorer:
         if len(self.model_ligands) == 0:
             return ("no_ligand", "No ligand in the model")
 
+        # Is the ligand disconnected?
+        graph = _ResidueToGraph(ligand)
+        if not networkx.is_connected(graph):
+            return ("disconnected", "Ligand graph is disconnected")
+
         # Is it because there was no valid binding site or no representation?
         if ligand in self._unassigned_target_ligands_reason:
             return self._unassigned_target_ligands_reason[ligand]
@@ -1495,7 +1510,7 @@ class LigandScorer:
                         substructure_match=self.substructure_match,
                         by_atom_index=True,
                         return_symmetries=False)
-                except NoSymmetryError:
+                except (NoSymmetryError, DisconnectedGraphError):
                     assigned = 0.
                 else:
                     assigned = 1.
@@ -1563,7 +1578,8 @@ def SCRMSD(model_ligand, target_ligand, transformation=geom.Mat4(),
                                ligand.
     :type substructure_match: :class:`bool`
     :rtype: :class:`float`
-    :raises: :class:`NoSymmetryError` when no symmetry can be found.
+    :raises: :class:`NoSymmetryError` when no symmetry can be found,
+             :class:`DisconnectedGraphError` when ligand graph is disconnected.
     """
 
     symmetries = _ComputeSymmetries(model_ligand, target_ligand,
@@ -1637,9 +1653,9 @@ def _ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
     target_graph = _ResidueToGraph(target_ligand, by_atom_index=by_atom_index)
 
     if not networkx.is_connected(model_graph):
-        raise RuntimeError("Disconnected graph for model ligand %s" % model_ligand)
+        raise DisconnectedGraphError("Disconnected graph for model ligand %s" % model_ligand)
     if not networkx.is_connected(target_graph):
-        raise RuntimeError("Disconnected graph for target ligand %s" % target_ligand)
+        raise DisconnectedGraphError("Disconnected graph for target ligand %s" % target_ligand)
 
     # 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
@@ -1683,5 +1699,10 @@ class NoSymmetryError(Exception):
     """
     pass
 
+class DisconnectedGraphError(Exception):
+    """Exception to be raised when the ligand graph is disconnected.
+    """
+    pass
+
 
-__all__ = ["LigandScorer", "SCRMSD", "NoSymmetryError"]
+__all__ = ["LigandScorer", "SCRMSD", "NoSymmetryError", "DisconnectedGraphError"]
diff --git a/modules/mol/alg/tests/test_ligand_scoring.py b/modules/mol/alg/tests/test_ligand_scoring.py
index 35de1ffb3..de7abda0c 100644
--- a/modules/mol/alg/tests/test_ligand_scoring.py
+++ b/modules/mol/alg/tests/test_ligand_scoring.py
@@ -467,6 +467,20 @@ class TestLigandScoring(unittest.TestCase):
         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 = [
             mdl.geometric_center,
@@ -494,7 +508,7 @@ class TestLigandScoring(unittest.TestCase):
         trg_zn = sc.target.FindResidue("H", 1)
         assert sc._find_unassigned_target_ligand_reason(trg_zn)[0] == "model_representation"
         assert sc.unassigned_target_ligands["H"][1][0] == "model_representation"
-        # CMO: not identical to anything in the model
+        # AFB: not identical to anything in the model
         trg_afb = sc.target.FindResidue("G", 1)
         assert sc._find_unassigned_target_ligand_reason(trg_afb)[0] == "identity"
         assert sc.unassigned_target_ligands["G"][1][0] == "identity"
@@ -502,6 +516,10 @@ class TestLigandScoring(unittest.TestCase):
         trg_fg3d1 = sc.target.FindResidue("F", 1)
         assert sc._find_unassigned_target_ligand_reason(trg_fg3d1)[0] == "stoichiometry"
         assert sc.unassigned_target_ligands["F"][1][0] == "stoichiometry"
+        # CMO: disconnected
+        trg_cmo1 = sc.target.FindResidue("L_CMO", 1)
+        sc._find_unassigned_target_ligand_reason(trg_cmo1)[0] == "disconnected"
+        assert sc.unassigned_target_ligands["L_CMO"][1][0] == "disconnected"
         # J.G3D1: assigned to L_2.G3D1 => error
         trg_jg3d1 = sc.target.FindResidue("J", 1)
         with self.assertRaises(RuntimeError):
@@ -537,6 +555,10 @@ class TestLigandScoring(unittest.TestCase):
         with self.assertRaises(RuntimeError):
             sc._find_unassigned_model_ligand_reason(mdl_mg_0)
         assert "L_MG_0" not in sc.unassigned_model_ligands
+        # CMO: disconnected
+        mdl_cmo1 = sc.model.FindResidue("L_CMO", 1)
+        sc._find_unassigned_model_ligand_reason(mdl_cmo1)[0] == "disconnected"
+        assert sc.unassigned_model_ligands["L_CMO"][1][0] == "disconnected"
         # Raises with an invalid ligand
         with self.assertRaises(ValueError):
             sc._find_unassigned_model_ligand_reason(sc.target_ligands[0])
@@ -552,7 +574,9 @@ class TestLigandScoring(unittest.TestCase):
             'L_OXY': {1: ('identity',
                           'Ligand was not found in the target (by isomorphism)')},
             'L_MG_2': {1: ('stoichiometry',
-                           'Ligand was assigned to an other target ligand (different stoichiometry)')}
+                           'Ligand was assigned to an other target ligand (different stoichiometry)')},
+            "L_CMO": {1: ('disconnected',
+                          'Ligand graph is disconnected')}
         }
         assert sc.unassigned_target_ligands == {
             'G': {1: ('identity',
@@ -564,7 +588,9 @@ class TestLigandScoring(unittest.TestCase):
             'K': {1: ('identity',
                       'Ligand was not found in the model (by isomorphism)')},
             'L_NA': {1: ('binding_site',
-                         'No residue in proximity of the target ligand')}
+                         'No residue in proximity of the target ligand')},
+            "L_CMO": {1: ('disconnected',
+                          'Ligand graph is disconnected')}
         }
         assert sc.lddt_pli["L_OXY"][1] is None
 
@@ -589,10 +615,6 @@ class TestLigandScoring(unittest.TestCase):
                               unassigned=True, rmsd_assignment=True)
 
 
-
-
-
-
 if __name__ == "__main__":
     from ost import testutils
     if testutils.DefaultCompoundLibIsSet():
-- 
GitLab