Skip to content

Commit

Permalink
Fix pi stacking (#98)
Browse files Browse the repository at this point in the history 
Fix intersect calculation and adjust pi-stacking parameters
  • Loading branch information
@cbouy
cbouy authored Nov 17, 2022
1 parent 49e1f74 commit 1523f99
Show file tree
Hide file tree
Showing 3 changed files with 77 additions and 43 deletions.
76 changes: 48 additions & 28 deletions prolif/interactions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -420,7 +420,7 @@ def __init__(self,
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
self.ring_radius = 1.7

def detect(self, ligand, residue):
for pi_rings in product(self.pi_ring, repeat=2):
Expand Down Expand Up @@ -451,44 +451,64 @@ def detect(self, ligand, residue):
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:
if not (
angle_between_limits(
n1c1c2, *self.normal_to_centroid_angle, ring=True
) or angle_between_limits(
n2c2c1, *self.normal_to_centroid_angle, ring=True
)
):
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:
intersect = self._get_intersect_point(
lig_normal, lig_centroid, res_normal, res_centroid
)
if intersect is None:
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:
# check if intersection point falls ~within plane ring
intersect_dist = min(
lig_centroid.Distance(intersect),
res_centroid.Distance(intersect)
)
if intersect_dist > self.ring_radius:
continue
return True, lig_match[0], res_match[0]
return False, None, None

@staticmethod
def _get_intersect_point(
plane_normal, plane_centroid, tilted_normal, tilted_centroid,
):
# intersect line is orthogonal to both planes normal vectors
intersect_direction = plane_normal.CrossProduct(tilted_normal)
# setup system of linear equations to solve
A = np.array([list(plane_normal), list(tilted_normal), list(intersect_direction)])
if np.linalg.det(A) == 0:
return None
tilted_offset = tilted_normal.DotProduct(
Geometry.Point3D(*tilted_centroid)
)
plane_offset = plane_normal.DotProduct(
Geometry.Point3D(*plane_centroid)
)
d = np.array([[plane_offset], [tilted_offset], [0.]])
# point on intersect line
point = np.linalg.solve(A, d).T[0]
point = Geometry.Point3D(*point)
# find projection of centroid on intersect line using vector projection
vec = plane_centroid - point
intersect_direction.Normalize()
scalar_proj = intersect_direction.DotProduct(vec)
return point + intersect_direction * scalar_proj


class FaceToFace(_BasePiStacking):
"""Face-to-face Pi-Stacking interaction between a ligand and a residue"""
def __init__(self,
centroid_distance=5.5,
plane_angle=(0, 35),
normal_to_centroid_angle=(0, 30),
normal_to_centroid_angle=(0, 33),
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,
Expand All @@ -510,8 +530,8 @@ class EdgeToFace(_BasePiStacking):
def __init__(self,
centroid_distance=6.5,
plane_angle=(50, 90),
normal_to_centroid_angle=(0, 25),
edge_padding=0.3,
normal_to_centroid_angle=(0, 30),
ring_radius=1.5,
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,
Expand All @@ -520,7 +540,7 @@ def __init__(self,
pi_ring=pi_ring
)
self.edge = True
self.edge_padding = edge_padding
self.ring_radius = ring_radius

def detect(self, ligand, residue):
return super().detect(ligand, residue)
Expand Down
26 changes: 16 additions & 10 deletions tests/notebooks/viz-pi-stacking.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,10 @@
"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",
"from MDAnalysis.transformations import translate, rotateby\n",
"import prolif as plf\n",
"from rdkit.Geometry import Point3D\n",
"from ipywidgets import interactive, HBox, Layout,VBox\n",
"from pprint import pprint"
"from ipywidgets import interactive, HBox, Layout,VBox"
]
},
{
Expand Down Expand Up @@ -75,19 +74,18 @@
" 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",
" proj = plf.interactions.EdgeToFace._get_intersect_point(n1, c1, n2, c2)\n",
" pdist = min(c1.Distance(proj), c2.Distance(proj))\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",
"centroid distance: {c1.Distance(c2):.3f} pdist: {pdist:.3f}\n",
"planes: {np.degrees(cap_angle(planes_angle)):.3f}°\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"
" return c1, n1, c2, n2, proj"
]
},
{
Expand All @@ -109,7 +107,7 @@
" 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",
" c1, n1, c2, n2, proj = measure(u)\n",
" try:\n",
" for comp in shapes.values():\n",
" comp.clear()\n",
Expand All @@ -118,6 +116,7 @@
" 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",
" shapes[\"proj\"] = v.shape.add_sphere(list(proj), [.8,.2,.6], .3)\n",
" \n",
"widget=interactive(view, \n",
" dx=(-7, 7, .5), dy=(-7, 7, .5), dz=(-7, 7, .5),\n",
Expand All @@ -126,6 +125,13 @@
"output = widget.children[-1]\n",
"display(VBox([controls, output, v]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
18 changes: 13 additions & 5 deletions tests/test_interactions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
from rdkit import Chem, RDLogger

from . import mol2factory
from .test_base import ligand_mol
from .test_base import ligand_mol, protein_mol

# disable rdkit warnings
lg = RDLogger.logger()
Expand Down Expand Up @@ -234,20 +234,22 @@ def test_smarts_matches(self, interaction_qmol, smiles, expected):
assert n_matches == expected

@pytest.mark.parametrize(["xyz", "rotation", "pi_type", "expected"], [
([0, 2.5, 4.5], [0, 0, 0], "facetoface", True),
([0, 2.5, 4.0], [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),
([1, 2.5, 4.5], [30, 15, 65], "facetoface", True),
([0, 1.5, 4.5], [60, 0, 0], "edgetoface", True),
([0, 2, 5], [60, 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),
([0, 1.5, 4.5], [105, 0, 0], "edgetoface", False),
([0, 1.5, -4.5], [75, 0, 0], "edgetoface", False),
])
def test_pi_stacking(self, xyz, rotation, pi_type, expected, fingerprint):
r1, r2 = self.create_rings(xyz, rotation)
Expand All @@ -267,3 +269,9 @@ def create_rings(xyz, rotation):
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)

def test_edgetoface_phe331(self):
fp = Fingerprint()
lig, phe331 = ligand_mol[0], protein_mol["PHE331.B"]
assert fp.edgetoface(lig, phe331) is True
assert fp.pistacking(lig, phe331) is True

0 comments on commit 1523f99

Please sign in to comment.