diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 666cfe44d..72f0f50d7 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -927,6 +927,171 @@ def matchAtoms( ) +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property(element).num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped0 = set(mapping) + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def _matchAtoms( molecule0, molecule1, @@ -937,13 +1102,14 @@ def _matchAtoms( timeout=5 * _Units.Time.second, complete_rings_only=True, max_scoring_matches=1000, - roi=None, prune_perturbed_constraints=False, prune_crossing_constraints=False, prune_atom_types=False, property_map0={}, property_map1={}, mcs_kwargs={}, + *, + _check_mapping=True, ): import sys as _sys @@ -1035,7 +1201,10 @@ def _matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer. Keep + # the original, since the mapping check below re-enters this function, + # which expects a Time object. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1153,6 +1322,87 @@ def _matchAtoms( property_map1, ) + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. + if _check_mapping and not mcs_kwargs and not prematch and mappings: + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. + message = None + try: + best = mappings[0] + + # Attachment points where the MCS stopped on both sides. + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. + retry = _matchAtoms( + molecule0=molecule0, + molecule1=molecule1, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. The subset test is + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + message = ( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ @@ -1490,9 +1740,12 @@ def _roiMatch( ) mapping = None else: - mapping = matchAtoms( + mapping = _matchAtoms( res0_extracted, res1_extracted, + # The mapping check would report indices that are local to the + # extracted residue, not to molecule0 as its message claims. + _check_mapping=False, ) # Look up the absolute atom indices in the molecule if not using a custom ROI mapping. @@ -1711,11 +1964,14 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -1913,11 +2169,14 @@ def _flexAlign( # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. diff --git a/src/BioSimSpace/FreeEnergy/_atm.py b/src/BioSimSpace/FreeEnergy/_atm.py index b1f126f28..e7d237e51 100644 --- a/src/BioSimSpace/FreeEnergy/_atm.py +++ b/src/BioSimSpace/FreeEnergy/_atm.py @@ -593,7 +593,10 @@ def _makeSystemFromThree(protein, ligand_bound, ligand_free, displacement): BioSimSpace._SireWrappers.System The system for the ATM simulation. """ - from ..Align import matchAtoms as _matchAtoms + # The private form is used so that the mapping check can be switched + # off below. It is otherwise identical to 'matchAtoms', which is a + # dispatcher on 'roi' and is never given one here. + from ..Align._align import _matchAtoms from ..Align import rmsdAlign as _rmsdAlign from ..Types import Vector as _Vector @@ -683,7 +686,9 @@ def _findTranslationVector(system, displacement, protein, ligand): out_of_protein = displacement.value() * initial_normal_vector return out_of_protein - mapping = _matchAtoms(ligand_free, ligand_bound) + # ATM doesn't expose 'mcs_kwargs', nor a way to pass in a mapping, so + # the advice from the mapping check can't be acted on here. + mapping = _matchAtoms(ligand_free, ligand_bound, _check_mapping=False) ligand_free_aligned = _rmsdAlign(ligand_free, ligand_bound, mapping) prot_lig1 = (protein + ligand_bound).toSystem() diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 6ff60101f..db4675e34 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -703,6 +703,171 @@ def generateNetwork( return edges, scores +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property(element).num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped0 = set(mapping) + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def defaultMCSOptions(): """ Return the default options used for the RDKit maximum common substructure @@ -742,6 +907,8 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -948,7 +1115,9 @@ def matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer, keeping + # the original for any onward call that expects a Time. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1066,6 +1235,87 @@ def matchAtoms( property_map1, ) + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. + if _check_mapping and not mcs_kwargs and not prematch and mappings: + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. + message = None + try: + best = mappings[0] + + # Attachment points where the MCS stopped on both sides. + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. + retry = matchAtoms( + molecule0, + molecule1, + engine=engine, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + max_scoring_matches=max_scoring_matches, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. The subset test is + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + message = ( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS. if prune_perturbed_constraints: mappings = [ @@ -1177,6 +1427,9 @@ def rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_map molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -1334,6 +1587,9 @@ def flexAlign( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index ef7fc2804..dca7b31be 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest import sire as sr @@ -1397,3 +1398,248 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # 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") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. It should also never fire on + the ROI path, where the flagged indices would be local to the extracted + residue rather than to molecule0 as the message claims. + """ + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_check_suppressed_roi(protein_inputs): + # The ROI path maps each residue of interest separately, so any flagged + # indices would be local to that residue rather than to molecule0. + proteins, protein_mapping, roi = protein_inputs + p0 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_mut_peptide.pdb") + )[0] + p1 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_wt_peptide.pdb") + )[0] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + assert BSS.Align.matchAtoms(p0, p1, roi=roi) == protein_mapping + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 6c85a1ee0..6d265a166 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest from sire.legacy.Maths import Vector @@ -6,7 +7,12 @@ from sire.legacy.Mol import AtomIdx, Element, PartialMolecule import BioSimSpace.Sandpit.Exscientia as BSS -from tests.Sandpit.Exscientia.conftest import has_antechamber, has_openff, url +from tests.Sandpit.Exscientia.conftest import ( + has_antechamber, + has_openff, + has_tleap, + url, +) @pytest.fixture(scope="session") @@ -968,3 +974,230 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # 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") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Sandpit.Exscientia.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Sandpit.Exscientia.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Sandpit.Exscientia.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. + """ + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28)