added support for templating a residue through an SDF file

Author: stefdoerr

moleculekit/molecule.py

@@ -3113,17 +3113,24 @@ def templateResidueFromMolecule(
         _logger=True,
     ):
         """Assign bonds, bond orders, formal charges and (optionally) hydrogens
-        to a residue from a reference Molecule template, matched by atom name.
+        to a residue from a reference Molecule template.
 
         Like :meth:`templateResidueFromSmiles`, but the template is a reference
         Molecule (e.g. loaded from a CIF) that already carries bonds, bond
-        orders and formal charges. The selected residue's heavy atoms are
-        mapped onto the reference by NAME and the reference's bond orders and
-        formal charges are transferred verbatim, so they must already be
-        correct in the reference (this function does not re-derive them). The
-        reference is used only as a template and is never appended. The
-        molecule is mutated in place. The residue and reference must share the
-        same set of heavy-atom names.
+        orders and formal charges. When the reference's heavy-atom names are
+        unique and equal to the residue's, the atoms are mapped by NAME and the
+        reference's bond orders and formal charges are transferred verbatim
+        (ideal for CIF references: unambiguous under molecular symmetry). When
+        the names do not match (for example a reference read from an SDF file,
+        whose atom names are only element symbols), the reference is converted
+        to a SMILES and matched by element and connectivity instead; for a
+        symmetric residue this can place a charge or double bond on an
+        equivalent atom differently while giving the same molecule. Either way
+        the reference is used only as a template and is never appended, and the
+        molecule is mutated in place. The reference must describe the same
+        residue: it may carry extra terminal atoms (such as a free amino acid's
+        OXT or a covalent leaving group), which are stripped, but a reference
+        missing heavy atoms of the residue raises.
 
         Parameters
         ----------
@@ -3132,9 +3139,9 @@ def templateResidueFromMolecule(
             of the residue(s) to template. May span multiple copies.
         refmol : Molecule
             The reference template. Its bonds, bond orders and formal charges
-            are copied onto the residue verbatim, so they must already be
-            correct; heavy-atom names must be unique and must match the
-            residue's.
+            must already be correct. When its heavy-atom names are unique and
+            match the residue's they are used directly; otherwise it is matched
+            by SMILES, so unnamed references (e.g. from SDF) are supported.
         addHs : bool
             If True, add hydrogens after bond orders are transferred.
         onlyOnAtoms : str or np.ndarray
@@ -3145,7 +3152,7 @@ def templateResidueFromMolecule(
         Examples
         --------
         >>> mol = Molecule("complex.pdb")  # doctest: +SKIP
-        >>> mol.templateResidueFromMolecule("resname LIG", Molecule("LIG.cif"), addHs=True, guessBonds=True)
+        >>> mol.templateResidueFromMolecule("resname LIG", Molecule("LIG.cif"), addHs=True, guessBonds=True)  # doctest: +SKIP
         """
         from moleculekit.rdkittools import template_residue_from_molecule
 

moleculekit/rdkittools.py

@@ -746,6 +746,23 @@ def template_residue_from_smiles(
                 "mol.guessBonds() before templating, or rely on LINK/CONECT records)."
             )
 
+    # Every residue heavy atom must be covered by the match. The check above
+    # handles a template that is a superset of the residue (extra template
+    # atoms are stripped); this handles the opposite case, an incomplete
+    # template that would leave a residue heavy atom with its guessed bonds.
+    # Erroring here keeps a reference missing heavy atoms (routed in from
+    # template_residue_from_molecule) from producing structurally wrong output.
+    res_heavy_idx = [a.GetIdx() for a in rmol.GetAtoms() if a.GetSymbol() != "H"]
+    unmatched_res = sorted(set(res_heavy_idx) - set(at1))
+    if unmatched_res:
+        unmatched_names = [str(residue.name[i]) for i in unmatched_res]
+        raise RuntimeError(
+            f"The SMILES template '{smiles}' does not cover all heavy atoms of "
+            f"the residue; unmatched residue atoms: {unmatched_names}. The "
+            "template must describe the same residue (it may carry extra "
+            "terminal atoms, but must not be missing any)."
+        )
+
     _apply_template_mapping(
         mol,
         selidx,
@@ -859,6 +876,42 @@ def template_residue_from_molecule(
             "The selection contains gaps in the atom indexes. Please select a single molecule residue only."
         )
 
+    # Decide the matching strategy for this single residue. Verbatim name
+    # matching is only viable when both sides have unique heavy-atom names and
+    # the name sets are equal (e.g. a CIF reference); it is unambiguous under
+    # molecular symmetry. Otherwise (e.g. an SDF reference, whose atom names are
+    # only element symbols) convert the reference to a SMILES and match by
+    # element and connectivity through the SMILES path.
+    res_heavy_names = mol.name[selidx][mol.element[selidx] != "H"]
+    ref_heavy_names = refmol.name[refmol.element != "H"]
+    name_match_viable = (
+        len(set(res_heavy_names)) == len(res_heavy_names)
+        and len(set(ref_heavy_names)) == len(ref_heavy_names)
+        and set(res_heavy_names) == set(ref_heavy_names)
+    )
+    if not name_match_viable:
+        ref_rdkit = refmol.toRDKitMol(
+            sanitize=True, kekulize=False, assignStereo=False, _logger=False
+        )
+        smiles = Chem.MolToSmiles(ref_rdkit)
+        if _logger:
+            logger.info(
+                "Reference atom names do not match the residue (e.g. an SDF "
+                "reference); matching by SMILES instead. Reference SMILES: "
+                f"{smiles!r}"
+            )
+        template_residue_from_smiles(
+            mol,
+            selidx,
+            smiles,
+            sanitizeSmiles=True,
+            addHs=addHs,
+            onlyOnAtoms=onlyOnAtoms,
+            guessBonds=guessBonds,
+            _logger=_logger,
+        )
+        return
+
     cross_bonds, sel_start, sel_end = _detect_cross_residue_bonds(mol, selidx)
 
     residue = mol.copy(sel=selidx)

tests/test_rdkittools.py

@@ -106,6 +106,50 @@ def _cmp(mol, ref_mol):
         )
 
 
+def test_templateResidueFromMolecule_sdf_reference(tmp_path):
+    """A reference Molecule without matching atom names (e.g. read from an SDF
+    file, whose atom names are only element symbols) is matched by converting it
+    to SMILES and using MCS, giving the same result as the name-matched path."""
+    import numpy as np
+
+    testdir = os.path.join(curr_dir, "test_molecule", "test_templating")
+    ref_named = Molecule(os.path.join(testdir, "BEN_pH7.4.cif"))
+
+    # Round-trip the reference through SDF. SDF keeps bonds, bond orders and
+    # formal charges but assigns element-only atom names, so name matching is no
+    # longer viable and the SMILES route is taken.
+    sdf_path = os.path.join(tmp_path, "ref.sdf")
+    ref_named.write(sdf_path)
+    ref_sdf = Molecule(sdf_path)
+    # The SDF reference's names are not unique (element symbols repeat), so the
+    # name path cannot apply and the SMILES route must be exercised.
+    assert len(np.unique(ref_sdf.name)) < ref_sdf.numAtoms
+
+    mol = Molecule(os.path.join(testdir, "BEN.pdb"))
+    mol.templateResidueFromMolecule(
+        "resname BEN", ref_sdf, addHs=True, guessBonds=True
+    )
+
+    # The SMILES / MCS path is symmetry-tolerant: benzamidine's two amidine
+    # nitrogens are interchangeable, so which one carries the +1 and the double
+    # bond can differ from the name-matched reference by a resonance swap. It is
+    # the same molecule, so compare by canonical SMILES (invariant under that
+    # automorphism) plus atom count, composition and net charge.
+    from rdkit import Chem
+
+    def _canonical(m):
+        rd = m.toRDKitMol(
+            sanitize=True, kekulize=False, assignStereo=False, _logger=False
+        )
+        return Chem.MolToSmiles(Chem.RemoveHs(rd))
+
+    ben = mol.copy(sel="resname BEN")
+    assert ben.numAtoms == ref_named.numAtoms
+    assert sorted(ben.element) == sorted(ref_named.element)
+    assert int(ben.formalcharge.sum()) == int(ref_named.formalcharge.sum())
+    assert _canonical(ben) == _canonical(ref_named)
+
+
 def test_templateResidueFromSmiles_multiresidue():
     """A selection that spans multiple residues with the same resname
     (e.g. ``resname BEN`` when there are several BEN copies in the
@@ -199,6 +243,16 @@ def test_templateResidueFromSmiles_incorrect_smiles():
     assert "2" in ben.bondtype
 
 
+def test_templateResidueFromSmiles_incomplete_template_errors():
+    """A SMILES that does not cover every residue heavy atom must raise rather
+    than silently leaving atoms untemplated. Benzene matches only BEN's ring,
+    leaving the amidine C and two N atoms unmatched."""
+    testdir = os.path.join(curr_dir, "test_molecule", "test_templating")
+    mol = Molecule(os.path.join(testdir, "BEN.pdb"))
+    with pytest.raises(RuntimeError):
+        mol.templateResidueFromSmiles("resname BEN", "c1ccccc1", guessBonds=True)
+
+
 @pytest.mark.parametrize(
     "smiles",
     (