diff --git a/docs/api.rst b/docs/api.rst index 3973d25..cc03375 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -113,6 +113,7 @@ Miscellaneous get_contact_residues find_atoms_by_pattern estimate_formal_charges + find_resonance_charges MatchWarning EvaluationWarning NoContactError diff --git a/src/peppr/contacts.py b/src/peppr/contacts.py index 308a987..ce103cb 100644 --- a/src/peppr/contacts.py +++ b/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_]]: """ diff --git a/tests/test_contacts.py b/tests/test_contacts.py index 6a67644..3e2459e 100644 --- a/tests/test_contacts.py +++ b/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()