analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('tableau-colorblind10')\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptm_torsion = {\n",
    "    \"SEP\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"OG\"],\n",
    "        [\"CA\", \"CB\", \"OG\", \"P\"],\n",
    "        [\"CB\", \"OG\", \"P\", \"O1P\"]],\n",
    "    \"TPO\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"OG1\"],\n",
    "        [\"CA\", \"CB\", \"OG1\", \"P\"],\n",
    "        [\"CB\", \"OG1\", \"P\", \"O1P\"]],\n",
    "    \"PTR\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD1\"],\n",
    "        [\"CE1\", \"CZ\", \"OH\", \"P\"],\n",
    "        [\"CZ\", \"OH\", \"P\", \"O1P\"]],\n",
    "    \"M3L\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD\"],\n",
    "        [\"CB\", \"CG\", \"CD\", \"CE\"],\n",
    "        [\"CG\", \"CD\", \"CE\", \"NZ\"],\n",
    "        [\"CD\", \"CE\", \"NZ\", \"CM1\"]\n",
    "        ],\n",
    "    \"ALY\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD\"],\n",
    "        [\"CB\", \"CG\", \"CD\", \"CE\"],\n",
    "        [\"CG\", \"CD\", \"CE\", \"NZ\"],\n",
    "        [\"CD\", \"CE\", \"NZ\", \"CH\"],\n",
    "        ],\n",
    "}\n",
    "\n",
    "nchis = {k: len(v) for k, v in ptm_torsion.items()}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polar_histogram(data, nbins=36):\n",
    "    if len(data) == 0:\n",
    "        bins = np.linspace(-np.pi, np.pi, nbins, endpoint=False)\n",
    "        return np.zeros(nbins - 1), np.zeros(nbins - 1), (bins[:-1] + bins[1:]) / 2\n",
    "    hist, bins = np.histogram(data, bins=np.linspace(-np.pi, np.pi, nbins+1))\n",
    "    sqrt_hist = np.sqrt(hist)/np.max(np.sqrt(hist))\n",
    "    return sqrt_hist, hist/np.sum(hist) * 100, (bins[:-1] + bins[1:]) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import i0\n",
    "\n",
    "def vonmises_kde(data, kappa, n_bins=36):\n",
    "    x = np.linspace(-np.pi, np.pi, n_bins)\n",
    "    # integrate vonmises kernels\n",
    "    kde = np.exp(kappa * np.cos(x[:, None] - data[None, :])).sum(1) / (2*np.pi*i0(kappa))\n",
    "    kde /= np.trapz(kde, x=x)\n",
    "    return x, kde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PTM data analysis; SEP as an example for visualizing histograms across backbone regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read ptm dataframes\n",
    "ptm_data = []\n",
    "for ptm in ptm_torsion:\n",
    "    ptm_data.append(pd.read_csv(f\"../data/{ptm}_df.csv\"))\n",
    "ptm_data = pd.concat(ptm_data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resn = \"SEP\"\n",
    "data = ptm_data[ptm_data.residue==resn]\n",
    "# Normalize the data to [0, 1] for color mapping\n",
    "norm = colors.Normalize(0, 100)\n",
    "cmap = cm.winter_r\n",
    "\n",
    "fig, ax = plt.subplots(4, 4, figsize=(10, 10), subplot_kw={'projection': 'polar'})\n",
    "for j in range(4):\n",
    "    if j == 0: \n",
    "        # backbone independent\n",
    "        angles = data\n",
    "        xlabel = \"bb. indp.\"\n",
    "    elif j == 1:\n",
    "        # helix\n",
    "        angles = data[(-105 <= data.phi).values & (data.phi < -45).values & \\\n",
    "                            (-60 <= data.psi).values & (data.psi < 30).values]\n",
    "        xlabel = r\"$\\alpha$-helix\"\n",
    "    elif j == 2:\n",
    "        # beta\n",
    "        angles = data[(data.phi.values < -105) & (data.psi.values > 90)]\n",
    "        xlabel = r\"$\\beta$-sheet\"\n",
    "    else:\n",
    "        # PII\n",
    "        angles = data[(-105 <= data.phi).values  & (data.phi < -45).values & \\\n",
    "                        (data.psi.values > 105)]\n",
    "        xlabel = \"polyproline\"\n",
    "    for i in range(4):\n",
    "        ax[i, j].set_theta_zero_location('E')\n",
    "        ax[i, j].yaxis.grid(False)\n",
    "        ax[i, j].set_yticks([])\n",
    "        if i < 3:\n",
    "            ax[i, j].set_xlabel(None)\n",
    "        if j > 0:\n",
    "            ax[i, j].set_ylabel(None)\n",
    "\n",
    "        sqrth, h, b = polar_histogram(angles[f\"chi{i+1}\"].apply(np.radians).values, nbins=36)\n",
    "        if i in [1, 2] and resn in [\"PTR\", \"TYR\"]:\n",
    "            sqrth, h, b = polar_histogram(angles[f\"chi{i+1}\"].apply(np.radians).values, nbins=24)\n",
    "            mask = np.logical_and(np.pi/4 < b, b <= 3*np.pi/4)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(15), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(-3*np.pi/4 < b, b <= -np.pi/4)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(15), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(-np.pi/4 < b, b <= np.pi/4)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(15), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_or(-3*np.pi/4 >= b, b > 3*np.pi/4)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(15), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "\n",
    "            ax[i, j].set_xticks([np.pi/4, 3*np.pi/4, 5*np.pi/4, 7*np.pi/4], [r\"$45^\\circ$\", r\"$135^\\circ$\", r\"$-135^\\circ$\", r\"$-45^\\circ$\"])  \n",
    "            continue\n",
    "        elif i == 4 and resn == \"ALY\":\n",
    "            sqrth, h, b = polar_histogram(angles[f\"chi{i+1}\"].apply(np.radians).values, nbins=36)\n",
    "            mask = np.logical_and(np.pi/6 < b, b <= np.pi/2)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(np.pi/2 < b, b <= 5*np.pi/6)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(-np.pi/6 < b, b <= np.pi/6)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(-np.pi/2 < b, b <= -np.pi/6)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_and(-5*np.pi/6 < b, b <= -np.pi/2)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "            mask = np.logical_or(-5*np.pi/6 >= b, b > 5*np.pi/6)\n",
    "            ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "\n",
    "            ax[i, j].set_xticks([np.pi/6, np.pi/2, 5*np.pi/6, 7*np.pi/6, 3*np.pi/2, 11*np.pi/6], \n",
    "                                [r\"$30^\\circ$\", r\"$90^\\circ$\", r\"$150^\\circ$\", r\"$-150^\\circ$\", r\"$-90^\\circ$\", r\"$-30^\\circ$\"])  \n",
    "            continue\n",
    "        \n",
    "        mask = np.logical_and(0 < b, b <= 2*np.pi/3)\n",
    "        ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "        mask = np.logical_and(-2*np.pi/3 < b, b <= 0)\n",
    "        ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "        mask = np.logical_or(2*np.pi/3 < b, b <= -2*np.pi/3)\n",
    "        ax[i, j].bar(b[mask], sqrth[mask], width=np.radians(10), alpha=0.7, color=cmap(norm(h[mask].sum())))\n",
    "        ax[i, j].set_xticks([0, 2*np.pi/3, 4*np.pi/3], [r\"$0^\\circ$\", r\"$120^\\circ$\", r\"$-120^\\circ$\"])  \n",
    "\n",
    "# Add colorbar\n",
    "#sm = cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "#sm.set_array([])\n",
    "#cbar = plt.colorbar(sm)\n",
    "#cbar.ax.set_ylabel('Rotamer population (%)')  \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize chi1 distributions on ramachandran plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 4, figsize=(14, 3))\n",
    "norm = colors.Normalize(-180, 180)\n",
    "cmap = cm.twilight_shifted\n",
    "ax[0].scatter(ptm_data[ptm_data.residue==\"SEP\"].phi, ptm_data[ptm_data.residue==\"SEP\"].psi, \n",
    "                c=ptm_data[ptm_data.residue==\"SEP\"].chi1, cmap=cmap, alpha=0.7)\n",
    "ax[1].scatter(ptm_data[ptm_data.residue==\"TPO\"].phi, ptm_data[ptm_data.residue==\"TPO\"].psi, \n",
    "                c=ptm_data[ptm_data.residue==\"TPO\"].chi1, cmap=cmap, alpha=0.7)\n",
    "ax[2].scatter(ptm_data[ptm_data.residue==\"PTR\"].phi, ptm_data[ptm_data.residue==\"PTR\"].psi, \n",
    "                c=ptm_data[ptm_data.residue==\"PTR\"].chi1, cmap=cmap, alpha=0.7)\n",
    "ax[3].scatter(ptm_data[ptm_data.residue==\"M3L\"].phi, ptm_data[ptm_data.residue==\"M3L\"].psi, \n",
    "                c=ptm_data[ptm_data.residue==\"M3L\"].chi1, cmap=cmap, alpha=0.7)\n",
    "\n",
    "for i in range(4):\n",
    "    ax[i].set_xlim(-180, 180)\n",
    "    ax[i].set_ylim(-180, 180)\n",
    "    # Highlight grid lines at x = [-120, 0, 120] and y = [-120, 0, 120] with magenta\n",
    "    for val in [30, 90, 150, -30, -90, -150]:\n",
    "        ax[i].axvline(x=val, color='grey', linestyle='--', linewidth=0.6)\n",
    "        ax[i].axhline(y=val, color='grey', linestyle='--', linewidth=0.6)\n",
    "    if i > 0:\n",
    "        ax[i].set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Helper function for returning backbone angles of residues forming local hydrogen bonds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.PDB import calc_dihedral, PDBParser\n",
    "\n",
    "parser = PDBParser()\n",
    "def phosphate_hbond(pdb, ptm):\n",
    "    struct = parser.get_structure(\"t\", pdb)\n",
    "    backbone = []\n",
    "    for chain in struct[0]:\n",
    "        for res in chain:\n",
    "            if res.get_resname() == ptm:\n",
    "                try:\n",
    "                    dist = min(res[\"N\"] - res[f\"O{i+1}P\"] for i in range(3))\n",
    "                    phi = calc_dihedral(chain[res.id[1] - 1][\"C\"].get_vector(), \n",
    "                            res[\"N\"].get_vector(), res[\"CA\"].get_vector(), res[\"C\"].get_vector())\n",
    "                    psi = calc_dihedral(res[\"N\"].get_vector(), \n",
    "                            res[\"CA\"].get_vector(), res[\"C\"].get_vector(), chain[res.id[1] + 1][\"N\"].get_vector())\n",
    "                    chi1 = calc_dihedral(*(res[atom].get_vector() for atom in ptm_torsion[ptm][0]))\n",
    "                    chi2 = calc_dihedral(*(res[atom].get_vector() for atom in ptm_torsion[ptm][1]))\n",
    "                    if dist <= 3.5 or min(chain[res.id[1] + 1][\"N\"] - res[f\"O{i+1}P\"] for i in range(3)) <= 3.5:\n",
    "                        backbone.append((np.degrees(phi), np.degrees(psi), np.degrees(chi1), np.degrees(chi2)))\n",
    "                except KeyError as e:\n",
    "                    continue\n",
    "                    #raise KeyError(pdb, chain.id, res.id)\n",
    "    return backbone\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mda",
   "language": "python",
   "name": "mda"
  },
  "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}