feat: expose resonant charges + test

docs/api.rst

@@ -113,6 +113,7 @@ Miscellaneous
     get_contact_residues
     find_atoms_by_pattern
     estimate_formal_charges
+    find_resonance_charges
     MatchWarning
     EvaluationWarning
     NoContactError

src/peppr/contacts.py

@@ -1,6 +1,7 @@
 __all__ = [
     "ContactMeasurement",
     "find_atoms_by_pattern",
+    "find_resonance_charges",
 ]
 
 from enum import IntEnum
@@ -283,13 +284,13 @@ def find_salt_bridges(
         When ``use_resonance=True``, Both oxygen atoms would be checked.
         """
         if use_resonance:
-            pos_mask, neg_mask, binding_site_conjugated_groups = (
-                _find_charged_atoms_in_resonance_structures(self._binding_site_mol)
+            pos_mask, neg_mask, binding_site_conjugated_groups = find_resonance_charges(
+                self._binding_site_mol
             )
             binding_site_pos_indices = np.where(pos_mask)[0]
             binding_site_neg_indices = np.where(neg_mask)[0]
-            pos_mask, neg_mask, ligand_conjugated_groups = (
-                _find_charged_atoms_in_resonance_structures(self._ligand_mol)
+            pos_mask, neg_mask, ligand_conjugated_groups = find_resonance_charges(
+                self._ligand_mol
             )
             ligand_pos_indices = np.where(pos_mask)[0]
             ligand_neg_indices = np.where(neg_mask)[0]
@@ -583,7 +584,7 @@ def _acceptable_angle(
     return abs(angle - ref_angle) <= tolerance
 
 
-def _find_charged_atoms_in_resonance_structures(
+def find_resonance_charges(
     mol: Chem.Mol,
 ) -> tuple[NDArray[np.bool_], NDArray[np.bool_], NDArray[np.int_]]:
     """

tests/test_contacts.py

@@ -5,6 +5,7 @@
 import biotite.structure.io.pdbx as pdbx
 import numpy as np
 import pytest
+from biotite.interface import rdkit as rdkit_interface
 import peppr
 from peppr.common import ACCEPTOR_PATTERN, DONOR_PATTERN, HBOND_DISTANCE_SCALING
 
@@ -301,3 +302,37 @@ def test_no_bonds(contact_method):
 
     contact_measurement = peppr.ContactMeasurement(receptor, ligand, 10.0)
     contact_method(contact_measurement)
+
+
+def test_find_resonance_charges():
+    """
+    Test finding charged atoms in resonance structures using a real ligand from a PDB file.
+    """
+    ligand = info.residue("ASP")
+    # standardize removes hydrogens - needed to estimate formal charges
+    ligand = peppr.standardize(ligand)
+
+    # set annotations for the benefit of finding charged atoms
+    ligand.set_annotation("charge", peppr.estimate_formal_charges(ligand, 7.4))
+    # create rdkit ligand object
+    ligand_mol = rdkit_interface.to_mol(ligand)
+    try:
+        peppr.sanitize(ligand_mol)
+    except Exception:
+        return np.nan
+    ligand_charged_atoms = np.where(ligand.charge != 0)[0]
+    assert np.equal(ligand_charged_atoms, [0, 7, 8]).all()
+
+    # get charged atoms and their resonance groups
+    pos_mask, neg_mask, ligand_conjugated_groups = peppr.find_resonance_charges(
+        ligand_mol
+    )
+    assert len(set(ligand_conjugated_groups)) < len(ligand_conjugated_groups), (
+        "Number of groups expected less than number of atoms as some are conjugated"
+    )
+    charged_atom_mask = pos_mask | neg_mask
+    ligand_charged_in_resonance_atoms = np.where(charged_atom_mask)[0]
+    assert set(ligand_charged_atoms) != set(ligand_charged_in_resonance_atoms), (
+        "Charged atoms do not match those found in resonance structures"
+    )
+    assert np.equal(ligand_charged_in_resonance_atoms, [0, 3, 6, 7, 8]).all()