diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 2602bb5cc..666cfe44d 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/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 `. + + 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) diff --git a/src/BioSimSpace/Align/_merge.py b/src/BioSimSpace/Align/_merge.py index 24b50649e..cac08e0e2 100644 --- a/src/BioSimSpace/Align/_merge.py +++ b/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 diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 2f7170fcb..6ff60101f 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/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 `. + + 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) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py index c4102611f..f14002dca 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py +++ b/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 diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index be01565e5..ef7fc2804 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -674,10 +674,10 @@ def test_roi_flex_align(protein_inputs): def test_empty_custom_roi_mapping(): # mut contains a proline mutation at position 15 wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] # use the custom_roi_map to specify that residue 15 in the WT protein should be @@ -691,15 +691,16 @@ def test_empty_custom_roi_mapping(): for atom_idx in roi_res_idx: assert atom_idx not in mapping.keys() + @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") def test_custom_roi_ring_break_merge(): # wt contains a leucine at position 15 # mut contains a proline at position 15 wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() @@ -743,13 +744,14 @@ def test_custom_roi_ring_break_merge(): assert n_bonds_created == 1 assert n_bonds_annihilated == 0 + @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") def test_custom_roi_map_invalid_outside_roi(): wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() @@ -761,7 +763,6 @@ def test_custom_roi_map_invalid_outside_roi(): molecule0=wt, molecule1=mut, roi=[15], - custom_roi_map={ 0: 0, 1: 1, @@ -1313,9 +1314,9 @@ def test_ring_breaking_cross_bond_cleanup(): mol_info.atom_idx(p.atom3()).value(), } for a, b in changing: - assert not ( - a in atoms and b in atoms - ), f"improper{suffix} spans absent bond ({a},{b})" + assert not (a in atoms and b in atoms), ( + f"improper{suffix} spans absent bond ({a},{b})" + ) # Check that the ring-breaking and ring-making bond properties are set. def _read_pairs(prop_name): @@ -1326,9 +1327,73 @@ def _read_pairs(prop_name): stored_breaking = _read_pairs("ring_breaking_bonds") stored_making = _read_pairs("ring_making_bonds") - assert ( - stored_breaking == ring_breaking - ), f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}" - assert ( - stored_making == ring_making - ), f"ring_making_bonds property mismatch: {stored_making} != {ring_making}" + assert stored_breaking == ring_breaking, ( + f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}" + ) + assert stored_making == ring_making, ( + f"ring_making_bonds property mismatch: {stored_making} != {ring_making}" + ) + + +@pytest.fixture(scope="session") +def ejm31(): + return BSS.IO.readMolecules( + [f"{url}/lig_ejm31.prm7.bz2", f"{url}/lig_ejm31.rst7.bz2"] + ).getMolecules()[0] + + +@pytest.fixture(scope="session") +def jmc28(): + return BSS.IO.readMolecules( + [f"{url}/lig_jmc28.prm7.bz2", f"{url}/lig_jmc28.rst7.bz2"] + ).getMolecules()[0] + + +def test_default_mcs_options(): + # The MCS defaults should be discoverable, and ring matching is on. + options = BSS.Align.defaultMCSOptions() + assert options["ringMatchesRingOnly"] is True + assert options["completeRingsOnly"] is True + + # The returned dictionary is a copy, so mutating it has no side effects. + options["ringMatchesRingOnly"] = False + assert BSS.Align.defaultMCSOptions()["ringMatchesRingOnly"] is True + + +def test_mcs_kwargs_ring_matches_ring_only(ejm31, jmc28): + # Perturbing a methyl to a 2-methylcyclopropyl. Atom 19 is the methyl + # carbon in ejm31 and the ring carbon bonded to the carbonyl in jmc28. + + # By default an acyclic atom can't map onto a ring atom, so the whole + # substituent is unmapped. + mapping = BSS.Align.matchAtoms(ejm31, jmc28) + assert 19 not in mapping + + # Allowing the match maps the two carbons onto each other, along with one + # of the methyl hydrogens. + mapping = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False} + ) + assert mapping[19] == 19 + assert len(mapping) == 30 + + # Only two hydrogens are removed and the ring is grown from dummy atoms. + assert sorted(set(range(32)) - set(mapping)) == [27, 28] + + +def test_mcs_kwargs_merge(ejm31, jmc28): + # The options are used when merge autogenerates a mapping. + merged = BSS.Align.merge(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + sire_mol = merged._sire_object + + # A ring grown entirely from dummy atoms breaks no bond between mapped + # atoms, so the merge doesn't require 'allow_ring_breaking' and the end + # states have the same number of bonds. + assert sire_mol.num_atoms() == 41 + assert len(sire_mol.property("bond0").potentials()) == len( + sire_mol.property("bond1").potentials() + ) + + # No ring is broken or made, so neither property is set. + assert not sire_mol.has_property("ring_breaking_bonds") + assert not sire_mol.has_property("ring_making_bonds") diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 18f6c86f0..6c85a1ee0 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -904,3 +904,67 @@ def test_ring_opening_and_size_change(ligands, mapping): BSS.Align.merge( m0, m1, mapping, allow_ring_breaking=True, allow_ring_size_change=True ) + + +@pytest.fixture(scope="session") +def ejm31(): + return BSS.IO.readMolecules( + [f"{url}/lig_ejm31.prm7.bz2", f"{url}/lig_ejm31.rst7.bz2"] + ).getMolecules()[0] + + +@pytest.fixture(scope="session") +def jmc28(): + return BSS.IO.readMolecules( + [f"{url}/lig_jmc28.prm7.bz2", f"{url}/lig_jmc28.rst7.bz2"] + ).getMolecules()[0] + + +def test_default_mcs_options(): + # The MCS defaults should be discoverable, and ring matching is on. + options = BSS.Align.defaultMCSOptions() + assert options["ringMatchesRingOnly"] is True + assert options["completeRingsOnly"] is True + + # The returned dictionary is a copy, so mutating it has no side effects. + options["ringMatchesRingOnly"] = False + assert BSS.Align.defaultMCSOptions()["ringMatchesRingOnly"] is True + + +def test_mcs_kwargs_ring_matches_ring_only(ejm31, jmc28): + # Perturbing a methyl to a 2-methylcyclopropyl. Atom 19 is the methyl + # carbon in ejm31 and the ring carbon bonded to the carbonyl in jmc28. + + # By default an acyclic atom can't map onto a ring atom, so the whole + # substituent is unmapped. + mapping = BSS.Align.matchAtoms(ejm31, jmc28) + assert 19 not in mapping + + # Allowing the match maps the two carbons onto each other, along with one + # of the methyl hydrogens. + mapping = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False} + ) + assert mapping[19] == 19 + assert len(mapping) == 30 + + # Only two hydrogens are removed and the ring is grown from dummy atoms. + assert sorted(set(range(32)) - set(mapping)) == [27, 28] + + +def test_mcs_kwargs_merge(ejm31, jmc28): + # The options are used when merge autogenerates a mapping. + merged = BSS.Align.merge(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + sire_mol = merged._sire_object + + # A ring grown entirely from dummy atoms breaks no bond between mapped + # atoms, so the merge doesn't require 'allow_ring_breaking' and the end + # states have the same number of bonds. + assert sire_mol.num_atoms() == 41 + assert len(sire_mol.property("bond0").potentials()) == len( + sire_mol.property("bond1").potentials() + ) + + # No ring is broken or made, so neither property is set. + assert not sire_mol.has_property("ring_breaking_bonds") + assert not sire_mol.has_property("ring_making_bonds")