02_tmaps.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import ListedColormap\n",
    "import joblib\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import rdMolDescriptors\n",
    "from mapchiral import mapchiral\n",
    "\n",
    "from tqdm import tqdm\n",
    "tqdm.pandas()\n",
    "\n",
    "import tmap as tm\n",
    "from faerun import Faerun\n",
    "\n",
    "plt.rcParams['font.sans-serif'] = \"Menlo\"\n",
    "plt.rcParams['font.family'] = \"sans-serif\"\n",
    "plt.rcParams.update({'font.size': 20})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Helper Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = ListedColormap(['#BCCFFF', '#648FFF', '#0042EF', '#D688AE', '#DC267F', '#A90052'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_tmap(query_name, savename):\n",
    "\n",
    "    # Extract the data\n",
    "    csv_files = [filename for filename in os.listdir('results/') if f'{query_name}' in filename and filename.endswith('.csv')]\n",
    "    csv_files = sorted(csv_files, key=lambda x: int(x.split('_')[3].split('.')[0]))\n",
    "\n",
    "    # Combine data from runs\n",
    "    df = pd.DataFrame()\n",
    "    for i, csv_file in enumerate(csv_files):\n",
    "        df_run = pd.read_csv(f'results/{csv_file}').head(1000)\n",
    "        df_run['run'] = int(csv_files[i].split('_')[3].split('.')[0]) + 1\n",
    "        df_run['fingerprint'] = csv_files[i].split('_')[2]\n",
    "        df = pd.concat([df, df_run]).reset_index(drop=True)\n",
    "    \n",
    "    # Clean data frame\n",
    "    df.sort_values(by=['fingerprint', 'generation'], inplace=True)\n",
    "    df['map4c'] = df['smiles'].apply(lambda x: mapchiral.encode(Chem.MolFromSmiles(x), max_radius=2, n_permutations=2048))\n",
    "    df['set'] = df['fingerprint'] + ': ' + df['run'].astype(str)\n",
    "\n",
    "    # Calculate TMAP layout \n",
    "    if not os.path.exists(f'tmaps/{savename}_tmap_layout.pkl'):\n",
    "        lf = tm.LSHForest(512, 32)\n",
    "        map4 = np.array(df['map4c'])\n",
    "        fps = []\n",
    "        for i in map4:\n",
    "            vec = tm.VectorUint(i)\n",
    "            fps.append(vec)\n",
    "        lf.batch_add(fps)\n",
    "        lf.index()\n",
    "        cfg = tm.LayoutConfiguration()\n",
    "        cfg.node_size = 1/32 \n",
    "        cfg.mmm_repeats = 2 \n",
    "        cfg.sl_extra_scaling_steps = 5 \n",
    "        cfg.k = 2000 \n",
    "        cfg.sl_scaling_type = tm.RelativeToAvgLength \n",
    "        x, y, s, t, _ = tm.layout_from_lsh_forest(lf, cfg)\n",
    "        tm_layout = {'x': list(x), 'y': list(y), 's': list(s), 't': list(t)}\n",
    "        joblib.dump(tm_layout, f'tmaps/{savename}_tmap_layout.pkl')\n",
    "    else:\n",
    "        tm_layout = joblib.load(f'tmaps/{savename}_tmap_layout.pkl')\n",
    "\n",
    "    # Define categories\n",
    "    set_labels, set_data = Faerun.create_categories(df['set'])\n",
    "\n",
    "    # Plot TMAP\n",
    "    f = Faerun(view=\"front\", coords=False, title=\"\", clear_color=\"#FFFFFF\", thumbnail_width=500)\n",
    "    f.add_scatter(f'{query_name}',{\"x\": tm.VectorFloat(tm_layout['x']), \"y\": tm.VectorFloat(tm_layout['y'] ), \"c\": [set_data], \"labels\": df['smiles'].values.tolist()}, shader=\"sphere\", point_scale=5, max_point_size=20, legend_labels=[set_labels], categorical=[True], colormap=[cmap], series_title=['Set'], has_legend=True)\n",
    "    f.add_tree(f'{query_name}_tree', {\"from\": tm.VectorUint(tm_layout['s']), \"to\": tm.VectorUint(tm_layout['t'])}, point_helper=f'{query_name}')\n",
    "    f.plot(f'tmaps/{savename}_tmap', template='smiles')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Epothilone A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/epothilone_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471bd82a310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/epothilone_tmap.html' target='_blank'>./tmaps/epothilone_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/epothilone_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('EpothiloneA', 'epothilone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cyclodextrin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/cyclodextrin_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471bd7ec110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/cyclodextrin_tmap.html' target='_blank'>./tmaps/cyclodextrin_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/cyclodextrin_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('aCyclodextrin', 'cyclodextrin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cyclosporin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/cyclosporin_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471bce0a310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/cyclosporin_tmap.html' target='_blank'>./tmaps/cyclosporin_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/cyclosporin_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('Cyclosporin', 'cyclosporin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonactin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/nonactin_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471b8b1b950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/nonactin_tmap.html' target='_blank'>./tmaps/nonactin_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/nonactin_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('Nonactin', 'nonactin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Onchidin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/onchidin_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471b8b2e250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/onchidin_tmap.html' target='_blank'>./tmaps/onchidin_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/onchidin_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('Onchidin', 'onchidin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Valinomycin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/valinomycin_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471e9552290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/valinomycin_tmap.html' target='_blank'>./tmaps/valinomycin_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/valinomycin_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_tmap('Valinomycin', 'valinomycin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### TMAP of random + self runs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import random runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_files = [filename for filename in os.listdir('results/') if f'PolymyxinB2_MAP4' in filename and 'GramicidinS' not in filename and 'PolymyxinB2_PolymyxinB2' not in filename and filename.endswith('.csv')]\n",
    "csv_files = sorted(csv_files, key=lambda x: ['1', '0', '2'].index(x.split('_')[3].split('.')[0]))\n",
    "\n",
    "df_rnd = pd.DataFrame()\n",
    "for i, csv_file in enumerate(csv_files):\n",
    "    df_run = pd.read_csv(f'results/{csv_file}')\n",
    "    df_run = df_run[df_run['dist'] < 0.5]\n",
    "    df_run['run'] = i \n",
    "    df_rnd = pd.concat([df_rnd, df_run])\n",
    "\n",
    "df_rnd.dropna(inplace=True)\n",
    "df_rnd.sort_values(by=['run', 'generation'], inplace=True)\n",
    "df_rnd['set'] = 'Random'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import self run (since they all generate the same compounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_self = pd.read_csv('results/20240523_PolymyxinB2_PolymyxinB2_MAP4_0.csv')\n",
    "df_self = df_self[df_self['dist'] < 0.5]\n",
    "df_self['run'] = 0\n",
    "df_self.dropna(inplace=True)\n",
    "df_self.sort_values(by=['run', 'generation'], inplace=True)\n",
    "df_self['set'] = 'Self'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.concat([df_rnd, df_self]).reset_index(drop=True)\n",
    "df['run'] = df['run'] + 1\n",
    "df['set_combined'] = df['set'] + ': ' + df['run'].astype(str)\n",
    "df.drop_duplicates(subset=['smiles'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine Levenshtein distance to Polymyxin B2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def levenshtein_distance(seq1, seq2):\n",
    "    blocks1 = seq1.split('-')\n",
    "    blocks2 = seq2.split('-')\n",
    "    len1 = len(blocks1)\n",
    "    len2 = len(blocks2)\n",
    "\n",
    "    dp = [[0] * (len2 + 1) for _ in range(len1 + 1)]\n",
    "    \n",
    "    for i in range(len1 + 1):\n",
    "        dp[i][0] = i\n",
    "    for j in range(len2 + 1):\n",
    "        dp[0][j] = j\n",
    "    \n",
    "    for i in range(1, len1 + 1):\n",
    "        for j in range(1, len2 + 1):\n",
    "            if blocks1[i - 1] == blocks2[j - 1]:\n",
    "                dp[i][j] = dp[i - 1][j - 1]\n",
    "            else:\n",
    "                dp[i][j] = min(dp[i - 1][j], dp[i][j - 1], dp[i - 1][j - 1]) + 1\n",
    "    \n",
    "    if \"c\" in seq1 and \"c\" in seq2:\n",
    "        return dp[len1][len2]\n",
    "    elif \"c\" in seq1 and \"c\" not in seq2:\n",
    "        return dp[len1][len2] - 1\n",
    "    elif \"c\" not in seq1 and \"c\" in seq2:\n",
    "        return dp[len1][len2] - 1\n",
    "    else:\n",
    "        return dp[len1][len2]\n",
    "    \n",
    "polymyxin_sequence = 'T011-BB045-BB012-BB045-b000-BB045-BB030-BB003-BB045-BB045-BB012'\n",
    "df['levenshtein'] = df['sequence'].apply(lambda x: levenshtein_distance(x, polymyxin_sequence))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate fingerprints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 12966/12966 [05:25<00:00, 39.83it/s]\n"
     ]
    }
   ],
   "source": [
    "df['map4c'] = df['smiles'].progress_apply(lambda x: mapchiral.encode(Chem.MolFromSmiles(x), max_radius=2, n_permutations=2048))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(f'tmaps/polymyxin_randself_tmap_layout.pkl'):\n",
    "    lf = tm.LSHForest(512, 32)\n",
    "    map4 = np.array(df['map4c'])\n",
    "    fps = []\n",
    "    for i in map4:\n",
    "        vec = tm.VectorUint(i)\n",
    "        fps.append(vec)\n",
    "    lf.batch_add(fps)\n",
    "    lf.index()\n",
    "    cfg = tm.LayoutConfiguration()\n",
    "    cfg.node_size = 1/32 \n",
    "    cfg.mmm_repeats = 2 \n",
    "    cfg.sl_extra_scaling_steps = 5 \n",
    "    cfg.k = 2000 \n",
    "    cfg.sl_scaling_type = tm.RelativeToAvgLength \n",
    "    x, y, s, t, _ = tm.layout_from_lsh_forest(lf, cfg)\n",
    "    tm_layout = {'x': list(x), 'y': list(y), 's': list(s), 't': list(t)}\n",
    "    joblib.dump(tm_layout, f'tmaps/polymyxin_randself_tmap_layout.pkl')\n",
    "else:\n",
    "    tm_layout = joblib.load(f'tmaps/polymyxin_randself_tmap_layout.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['highlights'] = ''\n",
    "df.loc[df['smiles'] == 'CC(C)CCCCC(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H]([C@@H](C)O)C(=O)N[C@@H]([C@H](C)O)C(=O)N[C@@H](CCN8)C(=O)N[C@@H](CC(C)C)C(=O)N[C@H](CC1=CC=CC=C1)C(=O)N[C@@H](CC(C)C)C(=O)N[C@@H](CC[NH3+])C(=O)N[C@H](CC[NH3+])C(=O)N[C@@H]([C@@H](C)O)C(=O)8', 'highlights'] = '7'\n",
    "df.loc[df['smiles'] == 'N9[C@@H](CC[NH3+])C(=O)N[C@@H](CC(C)C)C(=O)N[C@H](CC1=CC=CC=C1)C(=O)N[C@@H](CC1CCCCC1)C(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H]([C@H](C)O)C(=O)NCCCC(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H]([C@H](C)O)C(=O)N[C@@H](CC[NH3+])C(=O)9', 'highlights'] = '8'\n",
    "df.loc[df['smiles'] == 'CC(C)CCCCC(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H]([C@H](C)O)C(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H](CCN8)C(=O)N[C@@H](CC[NH3+])C(=O)N[C@H](CC1=CC=CC=C1)C(=O)N[C@@H](CC(C)C)C(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H](CC[NH3+])C(=O)N[C@@H]([C@H](C)O)C(=O)8', 'highlights'] = 'PMB2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"100%\"\n",
       "            height=\"500\"\n",
       "            src=\"./tmaps/polymyxin_randself_tmap.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7471bccda8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href='./tmaps/polymyxin_randself_tmap.html' target='_blank'>./tmaps/polymyxin_randself_tmap.html</a><br>"
      ],
      "text/plain": [
       "/home/markus/Developer/Code/10E60/tmaps/polymyxin_randself_tmap.html"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmap = ListedColormap(['#648FFF', '#DC267F', '#FFB000', '#CECECE'])\n",
    "set_labels, set_data = Faerun.create_categories(df['set_combined'])\n",
    "hihglight_labels, highlight_data = Faerun.create_categories(df['highlights'])\n",
    "\n",
    "# Plot TMAP\n",
    "f = Faerun(view=\"front\", coords=False, title=\"\", clear_color=\"#FFFFFF\", thumbnail_width=500)\n",
    "f.add_scatter(f'polymyxin_randself',\n",
    "              {\"x\": tm.VectorFloat(tm_layout['x']),\n",
    "                \"y\": tm.VectorFloat(tm_layout['y'] ), \n",
    "                \"c\": [set_data, \n",
    "                      df['dist'].values.tolist(),\n",
    "                      df['levenshtein'].values.tolist()],\n",
    "                \"labels\": df['smiles'].values.tolist()\n",
    "                }, \n",
    "            shader=\"sphere\", \n",
    "            point_scale=4, \n",
    "            max_point_size=20, \n",
    "            legend_labels=[set_labels, None, None], \n",
    "            categorical=[True, False, False], \n",
    "            colormap=[cmap, 'turbo', 'turbo'], \n",
    "            series_title=['Set', 'Jaccard Distance', 'Levenshtein Distance'], \n",
    "            has_legend=True)\n",
    "f.add_tree(f'polymyxin_randself_tree', {\"from\": tm.VectorUint(tm_layout['s']), \"to\": tm.VectorUint(tm_layout['t'])}, point_helper=f'polymyxin_randself')\n",
    "f.plot(f'tmaps/polymyxin_randself_tmap', template='smiles')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pdga2",
   "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}