better detecting of cross-residue bonds both by distance and by deposited bonds

Author: stefdoerr

moleculekit/rdkittools.py

@@ -316,7 +316,7 @@ def _try_strip_unmatched_terminals(
 _BONDTYPE_ORDER = {"1": 1, "2": 2, "3": 3, "4": 4, "5": 5, "6": 6, "ar": 1, "un": 1, "mc": 1}
 
 
-def _detect_cross_residue_bonds(mol, selidx):
+def _detect_interresidue_bonds(mol, selidx):
     """Return ``(cross_bonds, sel_start, sel_end)`` for covalent bonds linking
     the selected residue to the rest of ``mol``.
 
@@ -342,55 +342,39 @@ def _detect_cross_residue_bonds(mol, selidx):
             else:
                 cross_bonds.append((b - sel_start, a, bt))
 
-    # PDB inputs often arrive without explicit peptide / phosphodiester bonds
-    # (CONECT records cover only HET groups). Without them, the boundary
-    # atoms look free-standing and AddHs over-protonates them. Mirror the
-    # proximity check used by ``_has_peptide_neighbour`` in
-    # nonstandard_residues so a caller doesn't have to run ``mol.guessBonds``
-    # first.
-    sel_names = mol.name[selidx]
-    sel_elems = mol.element[selidx]
-
-    def _add_cross_by_proximity(local_idx, other_mask_full, threshold):
-        if local_idx in {li for li, _, _ in cross_bonds}:
-            return
-        own_pos = mol.coords[selidx[local_idx], :, mol.frame]
-        other_mask = other_mask_full.copy()
-        other_mask[selidx] = False
-        candidates = np.where(other_mask)[0]
-        if not len(candidates):
-            return
-        d = np.linalg.norm(
-            mol.coords[candidates, :, mol.frame] - own_pos, axis=1
-        )
-        within = d < threshold
-        if not within.any():
-            return
-        partner = int(candidates[np.argmin(np.where(within, d, np.inf))])
-        cross_bonds.append((local_idx, partner, "1"))
-
-    # Peptide N-C bonds (~1.33 A, threshold 1.6 A)
-    if {"N", "CA", "C"}.issubset(sel_names):
-        for own_side, other_name in (("N", "C"), ("C", "N")):
-            hits = np.where(sel_names == own_side)[0]
-            if len(hits):
-                _add_cross_by_proximity(
-                    int(hits[0]), mol.name == other_name, threshold=1.6
-                )
-
-    # Nucleic acid phosphodiester P-O3' bonds (~1.6 A, threshold 1.8 A).
-    # Two directions: (1) own P to external O3' of previous residue,
-    # (2) own O3' to external P of next residue. The O3' check also runs for
-    # 5'-terminal residues that lack their own P but still bond to next.
-    if "P" in sel_elems:
-        for own_p_idx in np.where(sel_elems == "P")[0]:
-            _add_cross_by_proximity(
-                int(own_p_idx), mol.element == "O", threshold=1.8
-            )
-    for own_o3_idx in np.where(sel_names == "O3'")[0]:
-        _add_cross_by_proximity(
-            int(own_o3_idx), mol.element == "P", threshold=1.8
+    # PDB inputs often arrive without explicit peptide / phosphodiester /
+    # isopeptide bonds (CONECT records cover only HET groups). Without them the
+    # boundary atoms look free-standing and AddHs over-protonates them (or leaves
+    # an isopeptide carbon over-valent). Recover them from geometry against the
+    # file-adjacent residues using the shared ``geometric_interresidue_links`` primitive, so
+    # the definition and thresholds match autoSegment / detect / systemPrepare.
+    # This covers the standard peptide / phosphodiester backbone AND the
+    # non-standard side-chain isopeptide (e.g. microcystin's ACB.CG -> ARG.N) in
+    # one pass.
+    from moleculekit.tools.nonstandard_residues import geometric_interresidue_links
+
+    def _residue_atoms_of(rep):
+        mask = (
+            (mol.resid == mol.resid[rep])
+            & (mol.chain == mol.chain[rep])
+            & (mol.insertion == mol.insertion[rep])
+            & (mol.segid == mol.segid[rep])
         )
+        return np.where(mask)[0]
+
+    neighbors = []
+    if sel_start > 0:
+        neighbors.append(_residue_atoms_of(sel_start - 1))
+    if sel_end + 1 < mol.numAtoms:
+        neighbors.append(_residue_atoms_of(sel_end + 1))
+
+    have_cross = {li for li, _, _ in cross_bonds}
+    for neighbor in neighbors:
+        for ia, ib, _kind in geometric_interresidue_links(mol, selidx, neighbor):
+            local_idx = int(ia) - sel_start
+            if local_idx not in have_cross:
+                cross_bonds.append((local_idx, int(ib), "1"))
+                have_cross.add(local_idx)
 
     return cross_bonds, sel_start, sel_end
 
@@ -645,7 +629,7 @@ def template_residue_from_smiles(
             "The selection contains gaps in the atom indexes. Please select a single molecule residue only."
         )
 
-    cross_bonds, sel_start, sel_end = _detect_cross_residue_bonds(mol, selidx)
+    cross_bonds, sel_start, sel_end = _detect_interresidue_bonds(mol, selidx)
 
     residue = mol.copy(sel=selidx)
     if guessBonds:
@@ -912,7 +896,7 @@ def template_residue_from_molecule(
         )
         return
 
-    cross_bonds, sel_start, sel_end = _detect_cross_residue_bonds(mol, selidx)
+    cross_bonds, sel_start, sel_end = _detect_interresidue_bonds(mol, selidx)
 
     residue = mol.copy(sel=selidx)
     if guessBonds:

moleculekit/tools/autosegment.py

@@ -102,24 +102,42 @@ def _polymer_linked(
     ``nucleic_fallback_cutoff``. If no usable atom pair exists, returns False
     (treated as a break).
     """
+    from moleculekit.tools.nonstandard_residues import (
+        geometric_interresidue_links,
+        _CANONICAL_RESNAMES,
+    )
+
     if category == "protein":
-        c = _residue_atom_coord(mol, prev_idx, "C")
-        n = _residue_atom_coord(mol, curr_idx, "N")
-        if c is not None and n is not None:
-            if float(np.linalg.norm(c - n)) <= protein_cutoff:
+        # Geometric backbone continuation via the shared primitive: a peptide
+        # bond always links; a side-chain isopeptide links only at a NON-canonical
+        # junction (geometric inference is less certain than a deposited bond, so
+        # two canonical residues are not merged on a coincidental side-chain
+        # contact - e.g. microcystin's ACB.CG->ARG.N links, a Gln-Lys crosslink
+        # between two standard residues does not).
+        noncanon = (
+            str(mol.resname[prev_idx[0]]) not in _CANONICAL_RESNAMES
+            or str(mol.resname[curr_idx[0]]) not in _CANONICAL_RESNAMES
+        )
+        for _, _, k in geometric_interresidue_links(
+            mol, prev_idx, curr_idx, amide_dist=protein_cutoff
+        ):
+            if k == "peptide" or (k == "isopeptide" and noncanon):
                 return True
-        else:
+        # CA-CA fallback ONLY when the backbone C/N needed for the precise check
+        # is missing (an incomplete residue) - never for present-but-far C/N,
+        # which is a real chain gap.
+        if (
+            _residue_atom_coord(mol, prev_idx, "C") is None
+            or _residue_atom_coord(mol, curr_idx, "N") is None
+        ):
             ca0 = _residue_atom_coord(mol, prev_idx, "CA")
             ca1 = _residue_atom_coord(mol, curr_idx, "CA")
             if ca0 is not None and ca1 is not None:
                 if float(np.linalg.norm(ca0 - ca1)) <= ca_fallback_cutoff:
                     return True
-        # Distance alone misses NON-standard backbone links - the backbone N of
-        # ``curr`` acylated by a side-chain carbonyl of ``prev`` (e.g. the
-        # gamma-glutamyl isopeptide FGA.CD->DAM.N in microcystin), or prev's
-        # backbone C bonded into curr. Honor such an explicit bond so the
-        # connected polymer stays one segment.
-        return _has_interresidue_backbone_bond(
+        # Last resort: honor an explicit deposited backbone bond the geometry
+        # above did not catch (a stretched or non-amide CONECT link).
+        return _has_deposited_backbone_bond(
             mol, prev_idx, curr_idx, prev_names=("C",), curr_names=("N",)
         )
 
@@ -134,23 +152,23 @@ def _polymer_linked(
             if c3 is not None and p is not None:
                 if float(np.linalg.norm(c3 - p)) <= nucleic_fallback_cutoff:
                     return True
-        return _has_interresidue_backbone_bond(
+        return _has_deposited_backbone_bond(
             mol, prev_idx, curr_idx,
             prev_names=("O3'", "O3*", "C3'", "C3*"), curr_names=("P",),
         )
 
     return False
 
 
-def _has_interresidue_backbone_bond(mol, prev_idx, curr_idx, prev_names, curr_names):
-    """Return True if an explicit bond joins the two residues at the backbone:
-    a bond from ``curr``'s backbone atom (any of ``curr_names``) to any atom of
-    ``prev``, or from ``prev``'s backbone atom (any of ``prev_names``) to any
-    atom of ``curr``. This catches non-standard backbone linkages (e.g. an
-    isopeptide acylating curr's N) that a pure distance check misses, while
-    NOT merging side-chain-only crosslinks (e.g. a disulfide), which touch
-    neither backbone atom. Only consulted as a fallback when the distance
-    checks fail, so the per-pair bond scan runs at most at chain breaks.
+def _has_deposited_backbone_bond(mol, prev_idx, curr_idx, prev_names, curr_names):
+    """Return True if an EXPLICIT (deposited) bond joins the two residues at the
+    backbone: a bond from ``curr``'s backbone atom (any of ``curr_names``) to any
+    atom of ``prev``, or from ``prev``'s backbone atom (any of ``prev_names``) to
+    any atom of ``curr``. This honors a deposited CONECT link the geometric
+    :func:`geometric_interresidue_links` check could miss (a stretched or non-amide bond),
+    while NOT merging side-chain-only crosslinks (e.g. a disulfide), which touch
+    neither backbone atom. Only consulted as a fallback when the geometric checks
+    in :func:`_polymer_linked` fail, so the bond scan runs at most at chain breaks.
     """
     if mol.bonds is None or len(mol.bonds) == 0:
         return False