Skip to content

Improve ring-break merge detection and expose RDKit MCS kwargs #542

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 3 commits into from
Jul 29, 2026
Merged

Improve ring-break merge detection and expose RDKit MCS kwargs #542

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b8760dd
Use two-atom in_ring so unrelated rings don't flag a ring break.
lohedges Jul 29, 2026
646fbef
Add mcs_kwargs to allow the user to configure the RDKit MCS.
lohedges Jul 29, 2026
ab5eb3b
Add unit tests for mcs_kwargs and updated ring-break check.
lohedges Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 55 additions & 13 deletions src/BioSimSpace/Align/_align.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
__email__ = "lester.hedges@gmail.com"

__all__ = [
"defaultMCSOptions",
"generateNetwork",
"matchAtoms",
"viewMapping",
Expand Down Expand Up @@ -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,
Expand All @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -912,6 +943,7 @@ def _matchAtoms(
prune_atom_types=False,
property_map0={},
property_map1={},
mcs_kwargs={},
):
import sys as _sys

Expand Down Expand Up @@ -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'")

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -2079,6 +2114,7 @@ def merge(
roi=None,
property_map0={},
property_map1={},
mcs_kwargs={},
**kwargs,
):
"""
Expand Down Expand Up @@ -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
-------

Expand Down Expand Up @@ -2208,6 +2249,7 @@ def merge(
molecule1,
property_map0=property_map0,
property_map1=property_map1,
mcs_kwargs=mcs_kwargs,
)
molecule0 = rmsdAlign(molecule0, molecule1, mapping)

Expand Down
8 changes: 4 additions & 4 deletions src/BioSimSpace/Align/_merge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
63 changes: 50 additions & 13 deletions src/BioSimSpace/Sandpit/Exscientia/Align/_align.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
__email__ = "lester.hedges@gmail.com"

__all__ = [
"defaultMCSOptions",
"generateNetwork",
"matchAtoms",
"viewMapping",
Expand Down Expand Up @@ -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,
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -1370,6 +1400,7 @@ def merge(
roi=None,
property_map0={},
property_map1={},
mcs_kwargs={},
**kwargs,
):
"""
Expand Down Expand Up @@ -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
-------

Expand Down Expand Up @@ -1488,6 +1524,7 @@ def merge(
molecule1,
property_map0=property_map0,
property_map1=property_map1,
mcs_kwargs=mcs_kwargs,
)
molecule0 = rmsdAlign(molecule0, molecule1, mapping)

Expand Down
8 changes: 4 additions & 4 deletions src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading