capillary_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Capillary Electrophoresis Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook reads in, formats and displays capillary electrophoresis trace files from an Agilent Fragment Analyzer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### File Prep and Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable, AxesGrid\n",
    "from IPython.display import display, clear_output\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import Layout, Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# make sure to reconfigure these file paths for your system\n",
    "base_path = r'SET INPUT FOLDER'\n",
    "base_output_loc = 'SET OUTPUT LOCATION'\n",
    "final_data = os.path.join(base_path, r'SET TO CAPILLARY ELECTROPHORESIS DATA FOLDER')\n",
    "\n",
    "print('Data available:', os.path.isdir(final_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Labelling - all labels are specific to particular wells/samples\n",
    "base_data = final_data\n",
    "prepend = '2021 10 11 14H 11M'  # results from paper\n",
    "available_samples = {'A1': 'Ladder',\n",
    "                     'A2': 'Single BRCA1\\n Payload',\n",
    "#                      'A3': 'Dual BRCA1\\n Payloads ' +   r\"$\\bf{(DB)}$\",\n",
    "                     'A3': 'Dual BRCA1\\n Payloads (DB)',\n",
    "                     'A4': 'DB +\\n Target (5 \\u03BCM)',\n",
    "                     'A5': 'DB +\\n Target (0.436 \\u03BCM)',\n",
    "                     'A6': 'DB +\\n Target (50 nM)',\n",
    "                     'A7': 'DB +\\n Mistoe (50 nM)',\n",
    "                     'A8': 'DB +\\n Misbody (50 nM)',\n",
    "                     'A9': 'DB +\\n Mismany (50 nM)',\n",
    "                     'A10':'DB +\\n 2-MT (50 nM)',\n",
    "                     'B1': 'Single miR-141 RNA\\n Payload',\n",
    "                     'B2': 'Dual miR-141 RNA\\n Payloads (DR)',\n",
    "                     'B3': 'DR +\\n Target (5 \\u03BCM)',\n",
    "                     'B4': 'DR +\\n Target (0.5 \\u03BCM)',\n",
    "                     'B5': 'DR +\\n Target (50 nM)',\n",
    "                     'C1': 'Single Aldosterone\\n Payload',\n",
    "                     'C2': 'Dual Aldosterone\\n Payloads (DA)',\n",
    "                     'C3': 'DA +\\n Target (1.5mM)',\n",
    "                     'C4': 'DA +\\n Target (150\\u03BCM)',\n",
    "                     'C5': 'DA +\\n Target (0.5\\u03BCM)',\n",
    "                     'C6': 'DA +\\n DNA Target (5\\u03BCM)',\n",
    "                     'D1': 'Dual HIV Block\\n (Pure)',\n",
    "                     'D2': 'Dual Aldosterone (DA)\\n + Dual BRCA1 (DB)\\n Payloads',\n",
    "                     'D3': 'DB + DA\\n+ 5uM BRCA target',\n",
    "                     'D4': 'DB + DA\\n+ 1.5mM Aldosterone',\n",
    "                     'D5': 'Quad BRCA1\\n Payloads (QB)',\n",
    "                     'D6': 'DB + QB',\n",
    "                     'D7': 'DB + QB\\n+ 5\\u03BCM BRCA1 target',\n",
    "                     'E1': 'Dual Thrombin\\n Payloads',\n",
    "                     'E2': 'Dual Thrombin\\n Payloads\\n+ Target (5\\u03BCM)',\n",
    "                     'F1': 'Impure\\n U4U2Thrco',\n",
    "                     'F2': 'Impure\\n QB',\n",
    "                     'F3': 'Impure\\n HR',\n",
    "                     'F4': 'Impure\\n U3',\n",
    "                     'H12': 'Capillary\\n Ladder'}\n",
    "\n",
    "presets = {\n",
    "    'BRCA tests':['A2', 'A3', 'A4', 'A5', 'A6', 'A7', 'A8', 'A9', 'A10'],\n",
    "    'RNA tests':['B1', 'B2', 'B3', 'B4', 'B5'],\n",
    "    'Aldos tests': ['C1', 'C2', 'C3', 'C4', 'C5', 'C6'],\n",
    "    'Thrombin tests': ['E1', 'E2'],\n",
    "    'Multi-Block Tests': ['D2','D3','D4','D5','D6','D7'],\n",
    "    'Others':['D1', 'F1', 'F2', 'F3', 'F4']\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data loading and formatting \n",
    "files = []\n",
    "dfs = []\n",
    "\n",
    "for key, label in available_samples.items():\n",
    "    file = '%s %s  Samp%s.csv' % (prepend, key, key)  # all files have this naming format\n",
    "    file_loc = os.path.join(base_data, file)\n",
    "    files.append(file_loc)\n",
    "    dfs.append(pd.read_csv(file_loc, header=0).add_suffix('_%s' % key))  \n",
    "    \n",
    "concatenated_df = pd.concat(dfs,axis=1).abs()\n",
    "intensity_df = concatenated_df[[c for c in concatenated_df.columns if 'Intensity' in c]]\n",
    "s_array = intensity_df.to_numpy()  # this is the final array containing all the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# main interface for selecting lanes and plotting results\n",
    "%matplotlib inline\n",
    "fig = None\n",
    "plt.rcParams.update({'font.sans-serif':'Helvetica'})  # consistent figure formatting\n",
    "mark_over_limits = True  # set to true to mark which bands have intensities above the selected limit\n",
    "plot_elution_time = False  # set to true to print elution time on y-axis\n",
    "\n",
    "\n",
    "def plot_preset(btn):  # plot one of the pre-selected groups of data\n",
    "    selection = presets[sel_preset.value[0]]\n",
    "    data_sel = intensity_df[[c for c in intensity_df.columns if c.split('_')[-1] in selection]].to_numpy()\n",
    "    plot_capillary(data_sel, selection)\n",
    "\n",
    "def plot_samples(btn):  # plot specific selected samples\n",
    "    selection = sel_mul.value\n",
    "    data_sel = intensity_df[[c for c in intensity_df.columns if c.split('_')[-1] in selection]].to_numpy()\n",
    "    plot_capillary(data_sel, selection)\n",
    "\n",
    "def plot_capillary(data, selection):  # vertical graph seems to be easier to read\n",
    "    if plot_type.value == 'vertical':\n",
    "        plot_capillary_vertical(data, selection)\n",
    "    else:\n",
    "        plot_capillary_horizontal(data, selection)\n",
    "\n",
    "def plot_capillary_vertical(cap_data, data_keys, normalization=False, level_shift=1800, max_level=12000):\n",
    "    # level_shift determines how much of the data to chop from the beginning of the x-axis, max_level indicates max intensity of colorbar\n",
    "    \n",
    "    global fig, out\n",
    "    \n",
    "    cmap = plt.cm.Blues\n",
    "    cmap.set_bad((1, 0, 0, 1))  # colour pixels outside of the specified range red for easy viewing\n",
    "    if mark_over_limits:\n",
    "        cap_data[cap_data > max_level] = np.nan  # marking illegal pixels for easy viewing\n",
    "    \n",
    "    cap_data = cap_data[level_shift:,:]  # shift data as first section is always blank\n",
    "    \n",
    "    title_font = 23  # specific graph formatting\n",
    "    label_font = 20\n",
    "    tick_font = 16\n",
    "    \n",
    "    with out:  # using output context to allow for clearing and re-drawing figure if necessary\n",
    "        if fig is not None:\n",
    "            clear_output()\n",
    "            \n",
    "        # Individual Normalization\n",
    "        if normalization:\n",
    "            for i in range(cap_data.shape[1]):\n",
    "                cap_data[:, i] = (cap_data[:, i]-cap_data[:,i].min())/(cap_data[:,i].max()-cap_data[:,i].min())\n",
    "\n",
    "        fig, ax = plt.subplots(sharey=True, ncols=cap_data.shape[1], figsize=(1.5*len(data_keys), 8))  # main plot shares y axis throughout\n",
    "\n",
    "        if cap_data.shape[1] == 1:\n",
    "            ax = [ax]\n",
    "            \n",
    "        dx = 80/300  # manual offset to center labels\n",
    "        # dx = 150/300  # manual offset to center labels\n",
    "\n",
    "        offset = mpl.transforms.ScaledTranslation(dx, 0, fig.dpi_scale_trans)  # graph label offset\n",
    "\n",
    "        for ind in range(cap_data.shape[1]):  # preparing heatmaps and visual elements\n",
    "            if data_keys[ind] == 'DA +\\n Target (1.5mM)':\n",
    "                print(cap_data[:,ind])\n",
    "            heatmap = ax[ind].imshow(np.flip(cap_data[:,ind][np.newaxis,:]).T, cmap=cmap, aspect=\"auto\", vmin=0, vmax = max_level)\n",
    "            # setting vmin=0,vmax=1.0 not required if data is normalized\n",
    "            ax[ind].tick_params(axis='x',which='both', bottom=False, labelleft=False, right=False, left=False, top=False)\n",
    "            ax[ind].xaxis.set_label_position('top')\n",
    "            \n",
    "            label = ax[ind].set_xlabel(available_samples[data_keys[ind]], fontsize=label_font, rotation=60, labelpad=10, ha='center')\n",
    "            \n",
    "            ax[ind].get_yaxis().set_visible(False)\n",
    "            label.set_transform(label.get_transform() + offset)\n",
    "            ax[ind].set_xticks([])  # blanking out x ticks\n",
    "            \n",
    "        if len(ax) == 1:\n",
    "            c_axis = ax[0]\n",
    "        else:\n",
    "            c_axis = ax.ravel().tolist()\n",
    "\n",
    "        cbar = fig.colorbar(heatmap, ax=c_axis, location='bottom', pad=0.05)\n",
    "\n",
    "        cbar.ax.set_xlabel('Signal Intensity (a.u.)', fontsize=title_font)\n",
    "        cbar.ax.tick_params(labelsize=tick_font)\n",
    "#         cbar.ax.get_yaxis().labelpad = 50\n",
    "        \n",
    "        ax[0].get_yaxis().set_visible(True)\n",
    "        \n",
    "        if plot_elution_time:\n",
    "            ax[0].yaxis.set_ticks([0, cap_data.shape[0]])\n",
    "            ax[0].set_yticklabels(['%.1f minutes' % ((cap_data.shape[0]+level_shift)/60), '30 minutes'],fontsize=tick_font, rotation=90, va='center')\n",
    "            ax[0].tick_params(axis='y', direction='out', length=12, width=1, which='both')\n",
    "            ax[0].set_ylabel('Elution Time',fontsize=title_font, labelpad=-15)\n",
    "        else:\n",
    "            ax[0].set_yticks([])\n",
    "            ax[0].set_ylabel('Elution Time',fontsize=title_font, labelpad=15)\n",
    "            \n",
    "        plt.show()\n",
    "    \n",
    "def plot_capillary_horizontal(cap_data, data_keys, normalization=False, level_shift=1800, max_level=12000):\n",
    "    global fig, out\n",
    "    \n",
    "    cmap = plt.cm.Blues\n",
    "    cmap.set_bad((1, 0, 0, 1))  # colour pixels outside of the specified range red for easy viewing\n",
    "    cap_data[cap_data > max_level] = np.nan  # marking illegal pixels for easy viewing\n",
    "    cap_data = cap_data[level_shift:,:]  # shift data as first section is always blank\n",
    "    title_font = 23\n",
    "    label_font = 20\n",
    "    tick_font = 16\n",
    "    \n",
    "    with out:\n",
    "        if fig is not None:\n",
    "            clear_output()\n",
    "            \n",
    "        # Individual Normalization\n",
    "        if normalization:\n",
    "            for i in range(cap_data.shape[1]):\n",
    "                cap_data[:, i] = (cap_data[:, i]-cap_data[:,i].min())/(cap_data[:,i].max()-cap_data[:,i].min())\n",
    "\n",
    "        fig, ax = plt.subplots(sharex=True, nrows=cap_data.shape[1], figsize=(12,len(data_keys)))\n",
    "\n",
    "        if cap_data.shape[1] == 1:\n",
    "            ax = [ax]\n",
    "    \n",
    "        for ind in range(cap_data.shape[1]):\n",
    "            heatmap = ax[ind].imshow(np.flip(cap_data[:,ind][np.newaxis,:]), cmap=cmap, aspect=\"auto\", vmin=0, vmax = max_level)\n",
    "            # setting vmin=0,vmax=1.0 not required if data is normalized\n",
    "            ax[ind].tick_params(left=False, labelleft=False, right=False)\n",
    "            ax[ind].set_ylabel(available_samples[data_keys[ind]], fontsize=label_font, rotation=0, va=\"center\", ha='center', labelpad=90)\n",
    "            ax[ind].get_xaxis().set_visible(False)\n",
    "\n",
    "        if len(ax) == 1:\n",
    "            c_axis = ax[0]\n",
    "        else:\n",
    "            c_axis = ax.ravel().tolist()\n",
    "        cbar = fig.colorbar(heatmap, ax=c_axis)\n",
    "\n",
    "        cbar.ax.set_ylabel('Signal Intensity\\n (Arbitrary Units)', fontsize=title_font, rotation=270)\n",
    "        cbar.ax.tick_params(labelsize=tick_font)\n",
    "        cbar.ax.get_yaxis().labelpad = 50\n",
    "        \n",
    "        ax[-1].xaxis.set_ticks([0, cap_data.shape[0]])\n",
    "        ax[-1].set_xticklabels(['%.1f minutes' % ((cap_data.shape[0]+level_shift)/60), '30 minutes'],fontsize=tick_font)\n",
    "        ax[-1].get_xaxis().set_visible(True)\n",
    "        ax[-1].set_xlabel('Elution Time',fontsize=title_font, labelpad=-15)\n",
    "        ax[-1].tick_params(axis='x', direction='out', length=12, width=1, which='both')\n",
    "        \n",
    "        plt.show()\n",
    "        \n",
    "\n",
    "def save_plot(btn):\n",
    "    fig.savefig(save_loc.value, bbox_inches='tight', dpi=300)\n",
    "    \n",
    "# GUI LAYOUT AND FORMATTING \n",
    "layout = widgets.Layout(width='auto') #set width and height\n",
    "sel_mul = widgets.SelectMultiple(\n",
    "                                    options=list(available_samples.keys()),\n",
    "                                    value=['A1'],\n",
    "                                    rows=10,\n",
    "                                    description='Available Samples',\n",
    "                                    layout=Layout(width=\"20%\",  justify_content=\"center\"),\n",
    "                                    disabled=False,\n",
    "                                    style={'description_width': 'initial'}\n",
    "                                 )\n",
    "sel_preset = widgets.SelectMultiple(\n",
    "                                    options=list(presets.keys()),\n",
    "                                    value=['BRCA tests'],\n",
    "                                    description='Preset Combinations',\n",
    "                                    layout=Layout(width=\"30%\"),\n",
    "                                    disabled=False,\n",
    "                                    style={'description_width': 'initial'}\n",
    "                                 )\n",
    "save_loc = widgets.Text(\n",
    "#                             value = os.path.join(os.getcwd(), 'cap_output.pdf'),\n",
    "                            value = os.path.join(base_output_loc, 'cap_output.png'),\n",
    "                            placeholder='Enter full filename here',\n",
    "                            description='',\n",
    "                            disabled=False,\n",
    "                            style={'description_width': 'initial'}\n",
    "                        ) \n",
    "\n",
    "plot_type = widgets.Dropdown(\n",
    "                                options=['vertical', 'horizontal'],\n",
    "                                value='vertical',\n",
    "                                description='',\n",
    "                                disabled=False,\n",
    "                                style={'description_width': 'initial'}\n",
    "                            )\n",
    "\n",
    "save_stack = widgets.VBox([widgets.Label('Output Filename:', layout=widgets.Layout(display='flex',justify_content='center')),\n",
    "                           save_loc,\n",
    "                           widgets.Label('Plot Type:', layout=widgets.Layout(display='flex',justify_content='center')),\n",
    "                           plot_type])\n",
    "\n",
    "sample_btn = widgets.Button(button_style='success',description = 'View Selected Samples', layout=layout)\n",
    "sample_btn.on_click(plot_samples)\n",
    "\n",
    "preset_btn = widgets.Button(button_style='success',description = 'View Selected Pre-set', layout=layout)\n",
    "preset_btn.on_click(plot_preset)\n",
    "\n",
    "save_btn = widgets.Button(button_style='info',description = 'Save Figure', layout=layout)\n",
    "save_btn.on_click(save_plot)\n",
    "btns = widgets.VBox([sample_btn, preset_btn, save_btn])\n",
    "\n",
    "display(widgets.HBox([sel_mul,sel_preset, btns, save_stack]))\n",
    "\n",
    "out = Output()  # this setup allows me to clear the figure without having to clear the widgets too\n",
    "display(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}