Merge pull request #542 from OpenBioSim/fix_ring_break_merge

Author: lohedges

src/BioSimSpace/Align/_align.py

@@ -25,6 +25,7 @@
 __email__ = "lester.hedges@gmail.com"
 
 __all__ = [
+    "defaultMCSOptions",
     "generateNetwork",
     "matchAtoms",
     "viewMapping",
@@ -702,6 +703,29 @@ def generateNetwork(
     return edges, scores
 
 
+def defaultMCSOptions():
+    """
+    Return the default options used for the RDKit maximum common substructure
+    search. These can be overridden using the 'mcs_kwargs' argument of
+    :class:`matchAtoms <BioSimSpace.Align.matchAtoms>`.
+
+    Returns
+    -------
+
+    options : dict
+        The default RDKit MCS options.
+    """
+    return {
+        "atomCompare": _rdFMCS.AtomCompare.CompareAny,
+        "bondCompare": _rdFMCS.BondCompare.CompareAny,
+        "completeRingsOnly": True,
+        "ringMatchesRingOnly": True,
+        "matchChiralTag": False,
+        "matchValences": False,
+        "maximizeBonds": False,
+    }
+
+
 def matchAtoms(
     molecule0,
     molecule1,
@@ -719,6 +743,7 @@ def matchAtoms(
     prune_atom_types=False,
     property_map0={},
     property_map1={},
+    mcs_kwargs={},
 ):
     """
     Find mappings between atom indices in molecule0 to those in molecule1.
@@ -775,6 +800,11 @@ def matchAtoms(
         option is only relevant to MCS performed using RDKit and will be
         ignored when falling back on Sire.
 
+    mcs_kwargs : dict
+        A dictionary of keyword arguments used to override the defaults
+        passed to the RDKit MCS search. This option is only relevant to MCS
+        performed using RDKit and will be ignored when falling back on Sire.
+
     roi : list
         The region of interest to match.
         Consists of a list of ROI residue indices.
@@ -881,6 +911,7 @@ def matchAtoms(
             prune_atom_types=prune_atom_types,
             property_map0=property_map0,
             property_map1=property_map1,
+            mcs_kwargs=mcs_kwargs,
         )
     else:
         return _roiMatch(
@@ -912,6 +943,7 @@ def _matchAtoms(
     prune_atom_types=False,
     property_map0={},
     property_map1={},
+    mcs_kwargs={},
 ):
     import sys as _sys
 
@@ -975,6 +1007,9 @@ def _matchAtoms(
     if not isinstance(complete_rings_only, bool):
         raise TypeError("'complete_rings_only' must be of type 'bool'")
 
+    if not isinstance(mcs_kwargs, dict):
+        raise TypeError("'mcs_kwargs' must be of type 'dict'")
+
     if type(max_scoring_matches) is not int:
         raise TypeError("'max_scoring_matches' must be of type 'int'")
 
@@ -1012,24 +1047,21 @@ def _matchAtoms(
             _Convert.toRDKit(mol1, property_map=property_map1),
         ]
 
+        # Default MCS options, overridden by anything in 'mcs_kwargs'. The
+        # timeout is applied last so that it can't be overridden.
+        mcs_options = defaultMCSOptions()
+        mcs_options["completeRingsOnly"] = complete_rings_only
+        mcs_options.update(mcs_kwargs)
+        mcs_options["timeout"] = timeout
+
         # Generate the MCS match.
-        mcs = _rdFMCS.FindMCS(
-            mols,
-            atomCompare=_rdFMCS.AtomCompare.CompareAny,
-            bondCompare=_rdFMCS.BondCompare.CompareAny,
-            completeRingsOnly=complete_rings_only,
-            ringMatchesRingOnly=True,
-            matchChiralTag=False,
-            matchValences=False,
-            maximizeBonds=False,
-            timeout=timeout,
-        )
+        mcs = _rdFMCS.FindMCS(mols, **mcs_options)
 
         # Get the common substructure as a SMARTS string.
         mcs_smarts = _Chem.MolFromSmarts(mcs.smartsString)
 
-    except:
-        raise RuntimeError("RDKit MCS mapping failed!")
+    except Exception as e:
+        raise RuntimeError(f"RDKit MCS mapping failed: {e}")
 
     # Score the mappings and return them in sorted order (best to worst).
     mappings, scores = _score_rdkit_mappings(
@@ -1066,6 +1098,9 @@ def _matchAtoms(
                 "Using Sire MCS. Ignoring unsupported 'complete_rings_only' option!"
             )
 
+        if mcs_kwargs:
+            _warnings.warn("Using Sire MCS. Ignoring unsupported 'mcs_kwargs' options!")
+
         # Convert timeout to a Sire Unit.
         timeout = timeout * _SireUnits.second
 
@@ -2079,6 +2114,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    mcs_kwargs={},
     **kwargs,
 ):
     """
@@ -2130,6 +2166,11 @@ def merge(
         A dictionary that maps "properties" in molecule1 to their user
         defined values.
 
+    mcs_kwargs : dict
+        A dictionary of keyword arguments used to override the defaults
+        passed to the RDKit MCS search. This is only used when 'mapping'
+        is None, i.e. when a mapping is autogenerated.
+
     Returns
     -------
 
@@ -2208,6 +2249,7 @@ def merge(
             molecule1,
             property_map0=property_map0,
             property_map1=property_map1,
+            mcs_kwargs=mcs_kwargs,
         )
         molecule0 = rmsdAlign(molecule0, molecule1, mapping)
 

src/BioSimSpace/Align/_merge.py

@@ -1589,10 +1589,10 @@ def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size
     # Supplementary check for rings larger than max_path: find_paths may only
     # find the direct-bond path and miss the long way around the ring, giving
     # n=1 instead of n≥2. Sire's in_ring has no path-length limit and
-    # correctly identifies ring membership in macrocycles.
-    if (conn0.in_ring(idx0) and conn0.in_ring(idy0)) != (
-        conn1.in_ring(idx1) and conn1.in_ring(idy1)
-    ):
+    # correctly identifies ring membership in macrocycles. The two-atom
+    # overload asks whether the atoms share a ring, so a ring built entirely
+    # from dummy atoms, which breaks no bond between mapped atoms, is ignored.
+    if conn0.in_ring(idx0, idy0) != conn1.in_ring(idx1, idy1):
         return True, False
 
     # A direct bond was replaced by a ring path (or vice versa), leaving the

src/BioSimSpace/Sandpit/Exscientia/Align/_align.py

@@ -25,6 +25,7 @@
 __email__ = "lester.hedges@gmail.com"
 
 __all__ = [
+    "defaultMCSOptions",
     "generateNetwork",
     "matchAtoms",
     "viewMapping",
@@ -702,6 +703,29 @@ def generateNetwork(
     return edges, scores
 
 
+def defaultMCSOptions():
+    """
+    Return the default options used for the RDKit maximum common substructure
+    search. These can be overridden using the 'mcs_kwargs' argument of
+    :class:`matchAtoms <BioSimSpace.Align.matchAtoms>`.
+
+    Returns
+    -------
+
+    options : dict
+        The default RDKit MCS options.
+    """
+    return {
+        "atomCompare": _rdFMCS.AtomCompare.CompareAny,
+        "bondCompare": _rdFMCS.BondCompare.CompareAny,
+        "completeRingsOnly": True,
+        "ringMatchesRingOnly": True,
+        "matchChiralTag": False,
+        "matchValences": False,
+        "maximizeBonds": False,
+    }
+
+
 def matchAtoms(
     molecule0,
     molecule1,
@@ -717,6 +741,7 @@ def matchAtoms(
     max_scoring_matches=1000,
     property_map0={},
     property_map1={},
+    mcs_kwargs={},
 ):
     """
     Find mappings between atom indices in molecule0 to those in molecule1.
@@ -770,6 +795,11 @@ def matchAtoms(
         option is only relevant to MCS performed using RDKit and will be
         ignored when falling back on Sire.
 
+    mcs_kwargs : dict
+        A dictionary of keyword arguments used to override the defaults
+        passed to the RDKit MCS search. This option is only relevant to MCS
+        performed using RDKit and will be ignored when falling back on Sire.
+
     prune_perturbed_constraints : bool
         Whether to remove hydrogen atoms that are perturbed to heavy atoms
         from the mapping. This is True for AMBER by default and False for
@@ -930,24 +960,21 @@ def matchAtoms(
             _Convert.toRDKit(molecule1, property_map=property_map1),
         ]
 
+        # Default MCS options, overridden by anything in 'mcs_kwargs'. The
+        # timeout is applied last so that it can't be overridden.
+        mcs_options = defaultMCSOptions()
+        mcs_options["completeRingsOnly"] = complete_rings_only
+        mcs_options.update(mcs_kwargs)
+        mcs_options["timeout"] = timeout
+
         # Generate the MCS match.
-        mcs = _rdFMCS.FindMCS(
-            mols,
-            atomCompare=_rdFMCS.AtomCompare.CompareAny,
-            bondCompare=_rdFMCS.BondCompare.CompareAny,
-            completeRingsOnly=complete_rings_only,
-            ringMatchesRingOnly=True,
-            matchChiralTag=False,
-            matchValences=False,
-            maximizeBonds=False,
-            timeout=timeout,
-        )
+        mcs = _rdFMCS.FindMCS(mols, **mcs_options)
 
         # Get the common substructure as a SMARTS string.
         mcs_smarts = _Chem.MolFromSmarts(mcs.smartsString)
 
-    except:
-        raise RuntimeError("RDKit MCS mapping failed!")
+    except Exception as e:
+        raise RuntimeError(f"RDKit MCS mapping failed: {e}")
 
     # Score the mappings and return them in sorted order (best to worst).
     mappings, scores = _score_rdkit_mappings(
@@ -984,6 +1011,9 @@ def matchAtoms(
                 "Using Sire MCS. Ignoring unsupported 'complete_rings_only' option!"
             )
 
+        if mcs_kwargs:
+            _warnings.warn("Using Sire MCS. Ignoring unsupported 'mcs_kwargs' options!")
+
         # Convert timeout to a Sire Unit.
         timeout = timeout * _SireUnits.second
 
@@ -1370,6 +1400,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    mcs_kwargs={},
     **kwargs,
 ):
     """
@@ -1417,6 +1448,11 @@ def merge(
         A dictionary that maps "properties" in molecule1 to their user
         defined values.
 
+    mcs_kwargs : dict
+        A dictionary of keyword arguments used to override the defaults
+        passed to the RDKit MCS search. This is only used when 'mapping'
+        is None, i.e. when a mapping is autogenerated.
+
     Returns
     -------
 
@@ -1488,6 +1524,7 @@ def merge(
             molecule1,
             property_map0=property_map0,
             property_map1=property_map1,
+            mcs_kwargs=mcs_kwargs,
         )
         molecule0 = rmsdAlign(molecule0, molecule1, mapping)
 

src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py

@@ -1455,10 +1455,10 @@ def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size
     # Supplementary check for rings larger than max_path: find_paths may only
     # find the direct-bond path and miss the long way around the ring, giving
     # n=1 instead of n≥2. Sire's in_ring has no path-length limit and
-    # correctly identifies ring membership in macrocycles.
-    if (conn0.in_ring(idx0) and conn0.in_ring(idy0)) != (
-        conn1.in_ring(idx1) and conn1.in_ring(idy1)
-    ):
+    # correctly identifies ring membership in macrocycles. The two-atom
+    # overload asks whether the atoms share a ring, so a ring built entirely
+    # from dummy atoms, which breaks no bond between mapped atoms, is ignored.
+    if conn0.in_ring(idx0, idy0) != conn1.in_ring(idx1, idy1):
         return True, False
 
     # A direct bond was replaced by a ring path (or vice versa), leaving the