Calculate geodesic distance on the protein surface

protein_geodesic_distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "27560d78-04a2-4f68-a442-7e95f84e719e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Snippet for calculation of geodesic distance between \n",
    "# two points on the protein surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "af3b1f90-bfef-410f-9c2f-bb86f54f53f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymol\n",
    "import prody\n",
    "import potpourri3d as pp3d\n",
    "from sklearn.neighbors import NearestNeighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df77f7d4-675b-4419-87fb-24668ce397e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set filenames\n",
    "pdbfn = '/tmp/PDBDEV_00000015.cif'\n",
    "objfn = '/tmp/PDBDEV_00000015.obj'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "81e74c5a-9fde-4c6f-a131-b6941da43b5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " PyMOL not running, entering library mode (experimental)\n",
      " ExecutiveLoad-Detail: Detected mmCIF\n"
     ]
    }
   ],
   "source": [
    "# Clear pymol state\n",
    "pymol.cmd.reinitialize()\n",
    "# Load PDB or mmCIF file into pymol\n",
    "pymol.cmd.load(pdbfn)\n",
    "# hide everything\n",
    "pymol.cmd.hide('everything')\n",
    "# You might want to remove all-non protein\n",
    "# entities before following steps\n",
    "pymol.cmd.remove('hetatm')\n",
    "# generate molecular surface\n",
    "pymol.cmd.show('surface')\n",
    "# save surface as OBJ file\n",
    "pymol.cmd.save(objfn)\n",
    "# Clear pymol state to free the memory\n",
    "pymol.cmd.reinitialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "739ca70e-20b8-4a68-8ec9-31ac7f17a898",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read OBJ file as mesh\n",
    "V, F = pp3d.io.read_mesh(objfn)\n",
    "# Prepare distance solver from vertices\n",
    "solver = pp3d.PointCloudHeatSolver(V)\n",
    "# Construct KD-Tree for fast search of closes vertices\n",
    "NN = NearestNeighbors(n_neighbors=1).fit(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "538d266c-c06d-4f50-adc4-6286b23deb3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load PDB or mmCIF structure\n",
    "pdb = prody.parseMMCIF(pdbfn)\n",
    "# Select atom 1\n",
    "a1 = pdb.select('chain J and resnum 147 and name CA')\n",
    "# Select atom 2\n",
    "a2 = pdb.select('chain E and resnum 381 and name CA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "322cc79e-fe15-49ed-96cc-d9a583600764",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance between atom 1 an point 1 on the surface: 0.40 Å\n",
      "Distance between atom 2 an point 2 on the surface: 0.49 Å\n"
     ]
    }
   ],
   "source": [
    "# Find closest point on the surface for atom 1\n",
    "p1 = NN.kneighbors(a1.getCoords())[1][0][0]\n",
    "# Check the distance between points:\n",
    "ds1 = prody.calcDistance(a1, V[p1])[0]\n",
    "print(f'Distance between atom 1 an point 1 on the surface: {ds1:.2f} Å')\n",
    "      \n",
    "# Find closest point on the surface for atom 2\n",
    "p2 = NN.kneighbors(a2.getCoords())[1][0][0]\n",
    "# Check the distance between points:\n",
    "ds2 = prody.calcDistance(a2, V[p2])[0]\n",
    "print(f'Distance between atom 2 an point 2 on the surface: {ds2:.2f} Å')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9fa33072-fa65-40d3-84dd-c03932f9b0da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute geodesic distances from point 1 to every other point\n",
    "dists = solver.compute_distance(p1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0be186f6-decb-402c-be82-ca86262a7d84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Geodesic distance: 23.11 Å\n",
      "Euclidian distance (surface points): 16.09 Å\n",
      "Euclidian distance (atoms): 15.93 Å\n"
     ]
    }
   ],
   "source": [
    "# Find the distance to point 2\n",
    "dg = dists[p2]\n",
    "# Also find euclidian distances\n",
    "# Between surface points\n",
    "dep = prody.calcDistance(V[p1], V[p2])\n",
    "# Between atoms\n",
    "dea = prody.calcDistance(a1, a2)[0]\n",
    "# Print the distance\n",
    "print(f'Geodesic distance: {dg:.2f} Å')\n",
    "print(f'Euclidian distance (surface points): {dep:.2f} Å')\n",
    "print(f'Euclidian distance (atoms): {dea:.2f} Å')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a4d0481-d2b9-4f8d-ae35-ebcd75a836f8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}