From 7ebe2b3d2347a2f76920a0d8bf68169a92eaff4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9dric=20Bouysset?= Date: Sun, 13 Nov 2022 17:55:01 +0100 Subject: [PATCH] Pi-stacking fix (#97) The Pi stacking interactions have been changed for a more accurate implementation: it now relies on the angle between the vector normal to a ring's plane and the vector between centroids (`normal_centroids_angles`) instead of the `shortest_distance` parameter. In addition, the `EdgeToFace` interaction makes sure that the intersection between both ring planes falls inside one of the rings. --- CHANGELOG.md | 5 +- prolif/fingerprint.py | 2 +- prolif/interactions.py | 161 +++++++++++++++++++++----- prolif/utils.py | 31 +++-- tests/notebooks/viz-pi-stacking.ipynb | 157 +++++++++++++++++++++++++ tests/test_interactions.py | 50 +++++++- 6 files changed, 362 insertions(+), 44 deletions(-) create mode 100644 tests/notebooks/viz-pi-stacking.ipynb diff --git a/CHANGELOG.md b/CHANGELOG.md index 1698b3d..d0ef46b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [1.1.0] - 2022-10-XX +## [1.1.0] - 2022-11-XX ### Added - `Fingerprint.run` now has a `converter_kwargs` parameter that can pass kwargs to the underlying RDKitConverter from MDAnalysis (Issue #57). @@ -21,7 +21,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 oxygen from esters, include some aromatic oxygen and nitrogen, - Anion: include resonance forms of carboxylic, sulfonic and phosphorus acids, - Cation: include amidine and guanidine, - - Metal ligand: exclude amides and some amines, + - Metal ligand: exclude amides and some amines. +- The Pi stacking interactions have been changed for a more accurate implementation. ### Fixed - Dead link in the quickstart notebook for the MDAnalysis quickstart (PR #75, @radifar). diff --git a/prolif/fingerprint.py b/prolif/fingerprint.py index 39b82be..d656674 100644 --- a/prolif/fingerprint.py +++ b/prolif/fingerprint.py @@ -442,7 +442,7 @@ def run(self, traj, lig, prot, residues=None, converter_kwargs=None, progress=Tr .. versionchanged:: 1.0.0 Added support for multiprocessing - .. versionadded:: 1.1.0 + .. versionchanged:: 1.1.0 Added support for passing kwargs to the RDKitConverter through the ``converter_kwargs`` parameter diff --git a/prolif/interactions.py b/prolif/interactions.py index 968ed70..383d102 100644 --- a/prolif/interactions.py +++ b/prolif/interactions.py @@ -134,6 +134,7 @@ class Hydrophobic(_Distance): distance : float Cutoff distance for the interaction + .. versionchanged:: 1.1.0 The initial SMARTS pattern was too broad. @@ -163,6 +164,7 @@ class _BaseHBond(Interaction): angles : tuple Min and max values for the ``[Donor]-[Hydrogen]...[Acceptor]`` angle + .. versionchanged:: 1.1.0 The initial SMARTS pattern was too broad. @@ -293,7 +295,7 @@ class _BaseIonic(_Distance): .. versionchanged:: 1.1.0 Handles resonance forms for common acids, amidine and guanidine. -""" + """ def __init__(self, cation="[+{1-},$([NX3&!$([NX3]-O)]-[C]=[NX3+])]", anion="[-{1-},$(O=[C,S,P]-[O-])]", @@ -329,6 +331,7 @@ class _BaseCationPi(Interaction): Min and max values for the angle between the vector normal to the ring plane and the vector going from the centroid to the cation + .. versionchanged:: 1.1.0 Handles resonance forms for amidine and guanidine as cations. @@ -385,29 +388,39 @@ def detect(self, ligand, residue): return super().detect(ligand, residue) -class PiStacking(Interaction): - """Pi-Stacking interaction between a ligand and a residue - +class _BasePiStacking(Interaction): + """Base class for Pi-Stacking interactions + Parameters ---------- centroid_distance : float Cutoff distance between each rings centroid - shortest_distance : float - Shortest distance allowed between the closest atoms of both rings plane_angles : tuple Min and max values for the angle between the ring planes + normal_centroids_angles : tuple + Min and max angles allowed between the vector normal to a ring's plane, + and the vector between the centroid of both rings. pi_ring : list List of SMARTS for aromatic rings + + + .. versionchanged:: 1.1.0 + The implementation now relies on the angle between the vector normal to a ring's + plane and the vector between centroids (``normal_centroids_angles``) instead of + the ``shortest_distance`` parameter. + """ def __init__(self, - centroid_distance=6.0, - shortest_distance=3.8, - plane_angles=(0, 90), + centroid_distance=5.5, + plane_angle=(0, 35), + normal_to_centroid_angle=(0, 30), pi_ring=("[a;r6]1:[a;r6]:[a;r6]:[a;r6]:[a;r6]:[a;r6]:1", "[a;r5]1:[a;r5]:[a;r5]:[a;r5]:[a;r5]:1")): self.pi_ring = [MolFromSmarts(s) for s in pi_ring] self.centroid_distance = centroid_distance - self.shortest_distance = shortest_distance**2 - self.plane_angles = tuple(radians(i) for i in plane_angles) + self.plane_angle = tuple(radians(i) for i in plane_angle) + self.normal_to_centroid_angle = tuple(radians(i) for i in normal_to_centroid_angle) + self.edge = False + self.edge_padding = 0 def detect(self, ligand, residue): for pi_rings in product(self.pi_ring, repeat=2): @@ -423,39 +436,124 @@ def detect(self, ligand, residue): cdist = lig_centroid.Distance(res_centroid) if cdist > self.centroid_distance: continue - squared_dist_matrix = np.add.outer( - (lig_pi_coords**2).sum(axis=-1), - (res_pi_coords**2).sum(axis=-1) - ) - 2*np.dot(lig_pi_coords, res_pi_coords.T) - shortest_dist = squared_dist_matrix.min().min() - if shortest_dist > self.shortest_distance: - continue # ligand lig_normal = get_ring_normal_vector(lig_centroid, lig_pi_coords) # residue res_normal = get_ring_normal_vector(res_centroid, res_pi_coords) - # angle between planes plane_angle = lig_normal.AngleTo(res_normal) - if angle_between_limits(plane_angle, *self.plane_angles, - ring=True): - return True, lig_match[0], res_match[0] + if not angle_between_limits( + plane_angle, *self.plane_angle, ring=True + ): + continue + c1c2 = lig_centroid.DirectionVector(res_centroid) + c2c1 = res_centroid.DirectionVector(lig_centroid) + n1c1c2 = lig_normal.AngleTo(c1c2) + n2c2c1 = res_normal.AngleTo(c2c1) + if angle_between_limits(n1c1c2, *self.normal_to_centroid_angle, ring=True): + tilted_plane_coords = res_pi_coords + tilted_normal = res_normal + plane_coords = lig_pi_coords + plane_normal = lig_normal + plane_centroid = lig_centroid + elif angle_between_limits(n2c2c1, *self.normal_to_centroid_angle, ring=True): + tilted_plane_coords = lig_pi_coords + tilted_normal = lig_normal + plane_coords = res_pi_coords + plane_normal = res_normal + plane_centroid = res_centroid + else: + continue + if self.edge: + # look for point of intersection between both ring planes + intersect_direction = plane_normal.CrossProduct(tilted_normal) + A = np.array([list(plane_normal), list(tilted_normal), list(intersect_direction)]) + if np.linalg.det(A) == 0: + continue + tilted_offset = tilted_normal.DotProduct(Geometry.Point3D(*tilted_plane_coords[0])) + plane_point = Geometry.Point3D(*plane_coords[0]) + plane_offset = plane_normal.DotProduct(plane_point) + d = np.array([[plane_offset], [tilted_offset], [0.]]) + intersect = np.linalg.solve(A, d).T + ring_radius = plane_centroid.Distance(plane_point) + if plane_centroid.Distance(Geometry.Point3D(*intersect[0])) > ring_radius + self.edge_padding: + continue + return True, lig_match[0], res_match[0] return False, None, None -class FaceToFace(PiStacking): +class FaceToFace(_BasePiStacking): """Face-to-face Pi-Stacking interaction between a ligand and a residue""" - def __init__(self, centroid_distance=4.5, plane_angles=(0, 40), **kwargs): - super().__init__(centroid_distance=centroid_distance, - plane_angles=plane_angles, **kwargs) + def __init__(self, + centroid_distance=5.5, + plane_angle=(0, 35), + normal_to_centroid_angle=(0, 30), + pi_ring=("[a;r6]1:[a;r6]:[a;r6]:[a;r6]:[a;r6]:[a;r6]:1", "[a;r5]1:[a;r5]:[a;r5]:[a;r5]:[a;r5]:1")): + super().__init__( + centroid_distance=centroid_distance, + plane_angle=plane_angle, + normal_to_centroid_angle=normal_to_centroid_angle, + pi_ring=pi_ring + ) -class EdgeToFace(PiStacking): - """Edge-to-face Pi-Stacking interaction between a ligand and a residue""" - def __init__(self, centroid_distance=6.0, plane_angles=(50, 90), **kwargs): - super().__init__(centroid_distance=centroid_distance, - plane_angles=plane_angles, **kwargs) +class EdgeToFace(_BasePiStacking): + """Edge-to-face Pi-Stacking interaction between a ligand and a residue + + .. versionchanged:: 1.1.0 + In addition to the changes made to the base pi-stacking interaction, this + implementation makes sure that the intersection between the perpendicular ring's + plane and the other's plane falls inside the ring. + + """ + def __init__(self, + centroid_distance=6.5, + plane_angle=(50, 90), + normal_to_centroid_angle=(0, 25), + edge_padding=0.3, + pi_ring=("[a;r6]1:[a;r6]:[a;r6]:[a;r6]:[a;r6]:[a;r6]:1", "[a;r5]1:[a;r5]:[a;r5]:[a;r5]:[a;r5]:1")): + super().__init__( + centroid_distance=centroid_distance, + plane_angle=plane_angle, + normal_to_centroid_angle=normal_to_centroid_angle, + pi_ring=pi_ring + ) + self.edge = True + self.edge_padding = edge_padding + + def detect(self, ligand, residue): + return super().detect(ligand, residue) + + +class PiStacking(Interaction): + """Pi-Stacking interaction between a ligand and a residue + + Parameters + ---------- + ftf_kwargs : dict + Parameters to pass to the underlying FaceToFace class + etf_kwargs : dict + Parameters to pass to the underlying EdgeToFace class + + + .. versionchanged:: 0.3.4 + `shortest_distance` has been replaced by `angle_normal_centroid` + + .. versionchanged:: 1.1.0 + The implementation now directly calls :class:`EdgeToFace` and :class:`FaceToFace` + instead of overwriting the default parameters with more generic ones. + + """ + def __init__(self, ftf_kwargs=None, etf_kwargs=None): + self.ftf = FaceToFace(**ftf_kwargs or {}) + self.etf = EdgeToFace(**etf_kwargs or {}) + + def detect(self, ligand, residue): + ftf = self.ftf.detect(ligand, residue) + if ftf[0] is True: + return ftf + return self.etf.detect(ligand, residue) class _BaseMetallic(_Distance): @@ -470,6 +568,7 @@ class _BaseMetallic(_Distance): distance : float Cutoff distance + .. versionchanged:: 1.1.0 The initial SMARTS pattern was too broad. diff --git a/prolif/utils.py b/prolif/utils.py index c0baaa3..11d0eda 100644 --- a/prolif/utils.py +++ b/prolif/utils.py @@ -54,15 +54,30 @@ def get_ring_normal_vector(centroid, coordinates): def angle_between_limits(angle, min_angle, max_angle, ring=False): - """Checks if an angle value is between min and max angles in radian. - If the angle to check involves a ring, include the angle that would be - obtained if we had used the other normal vector (same axis but opposite - direction) + """Checks if an angle value is between min and max angles in radian + + Parameters + ---------- + angle : float + Angle to check, in radians + min_angle : float + Lower bound angle, in radians + max_angle : float + Upper bound angle, in radians + ring : bool + Wether the angle being checked involves a ring or not + + Notes + ----- + When ``ring=True``, the angle is capped between 0 and 90, and so should be + the min and max angles. This is useful for angles involving a ring's plane + normal vector. """ - if ring and (angle > _90_deg_to_rad): - mirror_angle = _90_deg_to_rad - (angle % _90_deg_to_rad) - return (min_angle <= angle <= max_angle) or ( - min_angle <= mirror_angle <= max_angle) + if ring: + if angle >= pi: + angle %= _90_deg_to_rad + elif angle > _90_deg_to_rad: + angle = _90_deg_to_rad - (angle % _90_deg_to_rad) return (min_angle <= angle <= max_angle) diff --git a/tests/notebooks/viz-pi-stacking.ipynb b/tests/notebooks/viz-pi-stacking.ipynb new file mode 100644 index 0000000..ab81372 --- /dev/null +++ b/tests/notebooks/viz-pi-stacking.ipynb @@ -0,0 +1,157 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import nglview as nv\n", + "import MDAnalysis as mda\n", + "from MDAnalysis.transformations import translate, rotateby, center_in_box\n", + "import prolif as plf\n", + "from rdkit.Geometry import Point3D\n", + "from ipywidgets import interactive, HBox, Layout,VBox\n", + "from pprint import pprint" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "u1 = mda.Universe(plf.datafiles.datapath / \"benzene.mol2\")\n", + "elements = mda.topology.guessers.guess_types(u1.atoms.names)\n", + "u1.add_TopologyAttr(\"elements\", elements)\n", + "u1.segments.segids = np.array([\"U1\"], dtype=object)\n", + "u1.transfer_to_memory()\n", + "\n", + "def create(xyz=[0, 0, 0], rotation=[0,0,0]):\n", + " u2 = u1.copy()\n", + " u2.segments.segids = np.array([\"U2\"], dtype=object)\n", + " tr = translate(xyz)\n", + " rotx = rotateby(rotation[0], [1,0,0], ag=u2.atoms)\n", + " roty = rotateby(rotation[1], [0,1,0], ag=u2.atoms)\n", + " rotz = rotateby(rotation[2], [0,0,1], ag=u2.atoms)\n", + " u2.trajectory.add_transformations(tr, rotx, roty, rotz)\n", + " u2.transfer_to_memory()\n", + " u = mda.Merge(u1.atoms, u2.atoms)\n", + " return u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fp = plf.Fingerprint()\n", + "rad90 = np.pi/2\n", + "\n", + "def cap_angle(angle, cap=rad90):\n", + " if angle >= np.pi:\n", + " angle %= cap\n", + " elif angle > cap:\n", + " angle = (cap - (angle % cap))\n", + " return angle\n", + "\n", + "def measure(u):\n", + " ag1 = u.select_atoms(\"segid U1\")\n", + " ring1 = ag1.select_atoms(\"type C.2\").positions.astype(float)\n", + " c1 = plf.utils.get_centroid(ring1)\n", + " c1 = Point3D(*c1)\n", + " n1 = plf.utils.get_ring_normal_vector(c1, ring1)\n", + " \n", + " ag2 = u.select_atoms(\"segid U2\")\n", + " ring2 = ag2.select_atoms(\"type C.2\").positions.astype(float)\n", + " c2 = plf.utils.get_centroid(ring2)\n", + " c2 = Point3D(*c2)\n", + " n2 = plf.utils.get_ring_normal_vector(c2, ring2)\n", + " \n", + " planes_angle = n1.AngleTo(n2)\n", + " c1c2 = c1.DirectionVector(c2)\n", + " c2c1 = c2.DirectionVector(c1)\n", + " n1c1c2 = n1.AngleTo(c1c2)\n", + " n2c2c1 = n2.AngleTo(c2c1)\n", + " above = (plf.utils.angle_between_limits(n1c1c2, 0, np.radians(40), ring=True)\n", + " or\n", + " plf.utils.angle_between_limits(n2c2c1, 0, np.radians(40), ring=True))\n", + " \n", + " m1 = plf.Molecule.from_mda(ag1)\n", + " m2 = plf.Molecule.from_mda(ag2)\n", + " \n", + " print(f'''\n", + "centroid distance: {c1.Distance(c2):.3f}\n", + "planes: {np.degrees(cap_angle(planes_angle)):.3f}° above: {above}\n", + "n1c1c2: {np.degrees(cap_angle(n1c1c2)):.3f}° n2c2c1: {np.degrees(cap_angle(n2c2c1)):.3f}°\n", + "FTF: {fp.facetoface(m1, m2)} ETF: {fp.edgetoface(m1, m2)}''')\n", + " return c1, n1, c2, n2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "u = create(xyz=[0, 1.5, 4.5], rotation=[30, 0, 0])\n", + "v = nv.show_mdanalysis(u.atoms)\n", + "v.center(\"*\")\n", + "v._set_size(\"100%\", \"400px\")\n", + "v.camera = 'orthographic'\n", + "shapes = {}\n", + "\n", + "def view(dx=0, dy=1.5, dz=4.5, ax=30, ay=0, az=0):\n", + " new = create(xyz=[dx, dy, dz], rotation=[ax, ay, az])\n", + " u.atoms.positions = new.atoms.positions\n", + " v.set_coordinates({0: new.atoms.positions})\n", + " c1, n1, c2, n2 = measure(u)\n", + " try:\n", + " for comp in shapes.values():\n", + " comp.clear()\n", + " except:\n", + " pass\n", + " shapes[\"c1c2\"] = v.shape.add_cylinder(list(c1), list(c2), [1,0,0], .1)\n", + " shapes[\"n1\"] = v.shape.add_cylinder(list(c1), list(c1 + n1 + n1), [0,1,0], .1)\n", + " shapes[\"n2\"] = v.shape.add_cylinder(list(c2), list(c2 + n2 + n2), [0,0,1], .1)\n", + " \n", + "widget=interactive(view, \n", + " dx=(-7, 7, .5), dy=(-7, 7, .5), dz=(-7, 7, .5),\n", + " ax=(0, 180, 5), ay=(0, 180, 5), az=(0, 180, 5))\n", + "controls = HBox(widget.children[:-1], layout=Layout(flex_flow='row wrap'))\n", + "output = widget.children[-1]\n", + "display(VBox([controls, output, v]))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.8.13 ('prolif')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.13" + }, + "vscode": { + "interpreter": { + "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae" + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tests/test_interactions.py b/tests/test_interactions.py index 1e4a2bd..bde0bb1 100644 --- a/tests/test_interactions.py +++ b/tests/test_interactions.py @@ -1,6 +1,10 @@ import prolif import pytest +import numpy as np +import MDAnalysis as mda from MDAnalysis.topology.tables import vdwradii +from MDAnalysis.transformations import translate, rotateby +from rdkit import RDLogger from prolif.fingerprint import Fingerprint from prolif.interactions import (_INTERACTIONS, Interaction, VdWContact, get_mapindex) @@ -13,6 +17,11 @@ lg = RDLogger.logger() lg.setLevel(RDLogger.ERROR) +benzene = mda.Universe(prolif.datafiles.datapath / "benzene.mol2") +elements = mda.topology.guessers.guess_types(benzene.atoms.names) +benzene.add_TopologyAttr("elements", elements) +benzene.segments.segids = np.array(["U1"], dtype=object) +benzene.transfer_to_memory() interaction_instances = { name: cls() for name, cls in _INTERACTIONS.items() @@ -200,8 +209,10 @@ def test_vdwcontact_cache(self, lig_mol, prot_mol): ("_BaseCationPi.cation", "NC(=[NH2+])N", 3), ("_BaseCationPi.pi_ring", "c1ccccc1", 1), ("_BaseCationPi.pi_ring", "c1cocc1", 1), - ("PiStacking.pi_ring", "c1ccccc1", 1), - ("PiStacking.pi_ring", "c1cocc1", 1), + ("EdgeToFace.pi_ring", "c1ccccc1", 1), + ("EdgeToFace.pi_ring", "c1cocc1", 1), + ("FaceToFace.pi_ring", "c1ccccc1", 1), + ("FaceToFace.pi_ring", "c1cocc1", 1), ("_BaseMetallic.lig_pattern", "[Mg]", 1), ("_BaseMetallic.prot_pattern", "O", 1), ("_BaseMetallic.prot_pattern", "N", 1), @@ -221,3 +232,38 @@ def test_smarts_matches(self, interaction_qmol, smiles, expected): else: n_matches = len(mol.GetSubstructMatches(interaction_qmol)) assert n_matches == expected + + @pytest.mark.parametrize(["xyz", "rotation", "pi_type", "expected"], [ + ([0, 2.5, 4.5], [0, 0, 0], "facetoface", True), + ([0, 3, 4.5], [0, 0, 0], "facetoface", False), + ([0, 2, 4.5], [30, 0, 0], "facetoface", True), + ([0, 2, 4.5], [150, 0, 0], "facetoface", True), + ([0, 2, -4.5], [30, 0, 0], "facetoface", True), + ([0, 2, -4.5], [150, 0, 0], "facetoface", True), + ([1, 1.5, 3.5], [30, 15, 80], "facetoface", True), + ([0, 1.5, 4.5], [55, 0, 0], "edgetoface", True), + ([0, 1.5, 4.5], [90, 0, 0], "edgetoface", True), + ([0, 1.5, -4.5], [90, 0, 0], "edgetoface", True), + ([0, 6, -.5], [110, 0, 0], "edgetoface", True), + ([0, 4.5, -.5], [105, 0, 0], "edgetoface", True), + ([0, 1.5, 4.5], [115, 0, 0], "edgetoface", False), + ([0, 1.5, -4.5], [55, 0, 0], "edgetoface", False), + ]) + def test_pi_stacking(self, xyz, rotation, pi_type, expected, fingerprint): + r1, r2 = self.create_rings(xyz, rotation) + assert getattr(fingerprint, pi_type)(r1, r2) is expected + if expected is True: + other = "edgetoface" if pi_type == "facetoface" else "facetoface" + assert getattr(fingerprint, other)(r1, r2) is not expected + assert getattr(fingerprint, "pistacking")(r1, r2) is expected + + @staticmethod + def create_rings(xyz, rotation): + r2 = benzene.copy() + r2.segments.segids = np.array(["U2"], dtype=object) + tr = translate(xyz) + rotx = rotateby(rotation[0], [1,0,0], ag=r2.atoms) + roty = rotateby(rotation[1], [0,1,0], ag=r2.atoms) + rotz = rotateby(rotation[2], [0,0,1], ag=r2.atoms) + r2.trajectory.add_transformations(tr, rotx, roty, rotz) + return prolif.Molecule.from_mda(benzene), prolif.Molecule.from_mda(r2)