interactive-multiparameter-screening-plot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## jGCaMP8 interactive multiparameter plot\n",
    "\n",
    "Created with [plotly](https://plotly.com/)\n",
    "\n",
    "### Instructions\n",
    " Click on `Voila` button on the top banner to run and display the interactive plot\n",
    "\n",
    "### Operation\n",
    "\n",
    "All features can be plotted on the X axis, Y axis, or incorporated into the colormap. You can zoom in, pan, and scale axes using the pop-up menu on the top right of the plot. To return to the default view, click the \"Reset axes\" button on the pop-up menu. Click on any construct to show all construct features in a table at the bottom of the page. Controls, jGCaMP8 series, and XCaMP series constructs are highlighted in red.\n",
    "\n",
    "**All features are normalized to in-plate GCaMP6s controls**. \n",
    "\n",
    "For example, the table below (generated by clicking on any construct) should be interpreted as \"The DF/F (1 AP) of construct 500.656 is 4.96-fold higher than GCaMP6s. The half-rise time (1 AP) of construct 500.656 is 0.28-fold of GCaMP6s (i.e. 500.656 is 3.6x faster). \n",
    "\n",
    "| | 500.656|\n",
    "| ----------- | ----------- |\n",
    "| **DF/F (1 AP)** | 4.96 |\n",
    "|**Half-rise time (1 AP)** | 0.28 |\n",
    "\n",
    "\n",
    "\n",
    "### Widget controls\n",
    "\n",
    "* **X axis / Y axis / color**: Set what to plot on each axis\n",
    "\n",
    "* **X scale, Y scale**: linear or logarithmic axes\n",
    "\n",
    "* **Show all construct names**: Turn on to show all construct names. _Note: the construct names may not show up immediately. You may need to pan or zoom once to have them appear.__\n",
    "\n",
    "\n",
    "### Contact\n",
    "ilya kolb ([email](kolbi@hhmi.org))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "multi-parameter screening plots\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import plotly.express as px\n",
    "import pandas as pd\n",
    "import plotly.io as pio\n",
    "import plotly.graph_objects as go\n",
    "from ipywidgets import widgets\n",
    "from utils import condition_df, filter_construct_range\n",
    "\n",
    "pio.renderers.default='browser'\n",
    "\n",
    "plottableVars = ['DF/F (1 AP)', 'DF/F (3 AP)', 'DF/F (10 AP)', 'DF/F (160 AP)',\n",
    "       'Half-rise time (1 AP)', 'Half-rise time (3 AP)', 'Half-rise time (10 AP)', 'Half-rise time (160 AP)',\n",
    "       'Time to peak (1 AP)', 'Time to peak (3 AP)', 'Time to peak (10 AP)',\n",
    "       'Time to peak (160 AP)', 'Half-decay time (1 AP)', 'Half-decay time (3 AP)',\n",
    "       'Half-decay time (10 AP)', 'Half-decay time (160 AP)', \n",
    "       'SNR (1 AP)', 'SNR (3 AP)', 'SNR (10 AP)', 'SNR (160 AP)',\n",
    "        'dprime (1 AP)', 'dprime (3 AP)', 'dprime (10 AP)', 'dprime (160 AP)',\n",
    "       'Norm. F0']\n",
    "\n",
    "# data conditioning (filter out everything that does not pass these criteria)\n",
    "min_dff = 0\n",
    "max_rise = 4\n",
    "min_rise = 0.1\n",
    "min_decay = 0.01\n",
    "timetopeak_max = 3\n",
    "\n",
    "csv_dir = r'./data/data_all_GCaMP6scontrol.csv'\n",
    "id_seq_dir = r'./data/jG8-id-and-seq.csv'\n",
    "failed_wells_dir = r'./data/failed_wells.csv'\n",
    "\n",
    "mapping = {'GCaMP6s': '10.641' , 'jGCaMP8f': '500.456', 'jGCaMP8m': '500.686', 'jGCaMP8s': '500.688','jGCaMP8.712': '500.712',\n",
    "           'GCaMP6f': '10.693', 'jGCaMP7f': '10.921', 'jGCaMP7s': '10.1473', 'jGCaMP7c': '10.1513', 'jGCaMP7b': '10.1561',\n",
    "           'XCaMP-Gf': '538.1', 'XCaMP-G': '538.2', 'XCaMP-Gf0': '538.3'}\n",
    "mapping_swapped = dict([(value, key) for key, value in mapping.items()])\n",
    "data = pd.read_csv(csv_dir, na_values = 'NaN', sep='\\t') # data from screen\n",
    "\n",
    "#load and process id-sequence map table\n",
    "data_id_seq = pd.read_csv(id_seq_dir)\n",
    "data_id_seq['ID'] = data_id_seq['ID'].str.replace('dot', '.')\n",
    "data_id_seq = data_id_seq.set_index('ID', verify_integrity=True, drop=True)\n",
    "\n",
    "# load and process id-longname map table\n",
    "\n",
    "\n",
    "# load list of failed wells\n",
    "data_failed_wells = pd.read_csv(failed_wells_dir)\n",
    "\n",
    "'''\n",
    "add constructs that failed segmentation to the list, i.e. constructs in data_failed_wells that are NOT in the full dataset\n",
    "remove TE (dummy wells)\n",
    "concatenate to main dataset\n",
    "'''\n",
    "failed_seg_constructs = [d for d in data_failed_wells['Construct'] if not any(data['construct'].str.fullmatch(d))]\n",
    "failed_seg_constructs = [f for f in failed_seg_constructs if ('TE' not in f)]\n",
    "data = pd.concat([data, pd.DataFrame(data={'construct':failed_seg_constructs})])\n",
    "\n",
    "data = data.set_index('construct', drop=True, verify_integrity=True)\n",
    "\n",
    "# condition data table\n",
    "data = condition_df(data)\n",
    "\n",
    "# drop 500.723-1000\n",
    "data = data[[not filter_construct_range(i, '500', 722, 1000) for i in data.index]]\n",
    "\n",
    "data.dropna(axis = 0, how = 'any') # remove NaNs\n",
    "data = data.drop(columns=['d_fmax_f0', 'd_fmax_f0_p', 'es50', 'first_assay_date', 'last_assay_date', 'variant_type'\\\n",
    "                         ]) # remove extra columns\n",
    "\n",
    "# add aa sequence information\n",
    "data = data.join(data_id_seq) \n",
    "data['DNA sequence'] = data['DNA sequence'].replace(np.nan, 'Not available')\n",
    "data['AA sequence'] = data['AA sequence'].replace(np.nan, 'Not available')\n",
    "data['mutations'] = data['mutations'].replace(np.nan, 'Not available')\n",
    "\n",
    "# pretty column names\n",
    "data = data.rename(columns={'construct':'Construct', 'replicate_number': 'Number of wells', 'sequence': 'DNA sequence', 'aa_sequence': 'AA sequence'})\n",
    "\n",
    "# add mutation name (long name)\n",
    "data['Clone name'] = data['Clone name'].replace(np.nan, 'Not available')\n",
    "\n",
    "## condition data to eliminate metrics of sensors that didn't pass QC\n",
    "data_filt = data.copy()\n",
    "data_filt.loc[(data['Half-rise time (1 AP)'] >= max_rise) | \n",
    "         (data['Half-rise time (3 AP)'] >= max_rise) | \n",
    "         (data_filt['DF/F (1 AP)'] <= min_dff) | \n",
    "         (data_filt['DF/F (3 AP)'] <= min_dff) | \n",
    "         (data_filt['DF/F (10 AP)'] <= min_dff) | \n",
    "         (data_filt['DF/F (160 AP)'] <= min_dff) | \n",
    "         (data_filt['Half-rise time (1 AP)'] <= min_rise) |\n",
    "         (data_filt['Half-rise time (3 AP)'] <= min_rise) |\n",
    "         (data_filt['Half-rise time (10 AP)'] <= min_rise) |\n",
    "         (data_filt['Half-rise time (160 AP)'] <= min_rise) |\n",
    "         (data_filt['Half-decay time (1 AP)'] <= min_decay) |\n",
    "         (data_filt['Half-decay time (3 AP)'] <= min_decay) |\n",
    "         (data_filt['Half-decay time (10 AP)'] <= min_decay) |\n",
    "         (data_filt['Half-decay time (160 AP)'] <= min_decay) | \n",
    "         (data_filt['Time to peak (1 AP)'] >= timetopeak_max) | \n",
    "         (data_filt['Time to peak (3 AP)'] >= timetopeak_max) |\n",
    "         (data_filt['Time to peak (10 AP)'] >= timetopeak_max) | \n",
    "         (data_filt['Time to peak (160 AP)'] >= timetopeak_max), plottableVars + list(data.columns[data.columns.str.endswith('_p')].values)]= np.nan\n",
    "\n",
    "\n",
    "# for public use, filter by DF/F (1 AP), add colloquial names\n",
    "data_filt_public = data_filt.sort_values(by='DF/F (1 AP)', ascending=False).copy()\n",
    "data_filt_public = data_filt_public.rename(index=mapping_swapped)\n",
    "\n",
    "# add failed segmentation label\n",
    "data_filt_public.loc[data_filt_public.index.isin(failed_seg_constructs), 'DF/F (1 AP)'] = 'Failed segmentation'\n",
    "\n",
    "# add No detectable response label\n",
    "data_filt_public.loc[data_filt_public['DF/F (1 AP)'].isna(), 'DF/F (1 AP)'] = 'No detectable response'\n",
    "\n",
    "# save filtered data\n",
    "print('Saving data_filt pickle... ', end='')\n",
    "data_filt_public.to_pickle('outputs/data_filt_public.pkl')\n",
    "print('Done')\n",
    "\n",
    "# save filtered data csv\n",
    "print('Saving data_filt csv... ', end='')\n",
    "data_filt_public.to_csv(r'outputs/data_filt_public.csv')\n",
    "print('Done')\n",
    "\n",
    "# hightlight_txt_array for plotting highlight names in mapping. '' for non-highlighted, name for highlighted\n",
    "hightlight_txt_array = [mapping_swapped.get(c_id) if (c_id in mapping_swapped.keys()) else '' for c_id in data_filt.index ]\n",
    "highlight_TF_array = np.logical_not(np.array(hightlight_txt_array) == '')\n",
    "\n",
    "print('Total constructs: ' + str(len(data)))\n",
    "\n",
    "failure_reasons = ['Failed segmentation', 'No detectable response']\n",
    "\n",
    "for f in failure_reasons:\n",
    "    print(f + ': ' + str((data_filt_public['DF/F (1 AP)'] == f).sum()))\n",
    "\n",
    "g = go.FigureWidget({\n",
    "    'data': [{'customdata': data_filt.index,\n",
    "              'hovertemplate': '<b>%{customdata}</b><br>DF/F (1 AP)=%{x:.3f}<br>Half-rise time (1 AP)=%{y:.3f}<br>Half-decay time (1 AP)=%{marker.color:.3f}', \n",
    "              # ('%{x}<br>Half-rise (1FP)=%{y}<br' ... '{customdata[0]}<extra></extra>'),\n",
    "              'legendgroup': '',\n",
    "              'marker': {'color': data_filt['Half-decay time (1 AP)'],\n",
    "                         'coloraxis': 'coloraxis',\n",
    "                         'size': 10, # data_filt['Decay (1FP)'],\n",
    "                         'sizemode': 'area',\n",
    "                         'sizeref': 0.02,\n",
    "                         'symbol': 'circle',\n",
    "                         'opacity': 0.4,\n",
    "                         'line' : {\n",
    "                             'color': 'red',\n",
    "                             'width': 2 * highlight_TF_array,\n",
    "                            }\n",
    "                        },\n",
    "              'text': hightlight_txt_array,\n",
    "              'textfont': {'color': 'red'},\n",
    "              'textposition': 'top center',\n",
    "              'mode': 'text+markers',\n",
    "              'name': '',\n",
    "              'orientation': 'v',\n",
    "              'showlegend': False,\n",
    "              'type': 'scatter',\n",
    "              'x': data_filt['DF/F (1 AP)'],\n",
    "              'xaxis': 'x',\n",
    "              'y': data_filt['Half-rise time (1 AP)'],\n",
    "              'yaxis': 'y'}],\n",
    "    'layout': {'coloraxis': {'colorbar': {'title': {'text': 'Half-decay time (1 AP)'}},\n",
    "                             'colorscale': [[0.0, '#0d0887'], [0.1111111111111111,\n",
    "                                            '#46039f'], [0.2222222222222222,\n",
    "                                            '#7201a8'], [0.3333333333333333,\n",
    "                                            '#9c179e'], [0.4444444444444444,\n",
    "                                            '#bd3786'], [0.5555555555555556,\n",
    "                                            '#d8576b'], [0.6666666666666666,\n",
    "                                            '#ed7953'], [0.7777777777777778,\n",
    "                                            '#fb9f3a'], [0.8888888888888888,\n",
    "                                            '#fdca26'], [1.0, '#f0f921']]},\n",
    "               'legend': {'itemsizing': 'constant', 'tracegroupgap': 0},\n",
    "               'margin': {'t': 60},\n",
    "               'height': 500,\n",
    "               'width' : 700,\n",
    "               'xaxis': {'anchor': 'y', 'domain': [0.0, 1.0], 'title': {'text': 'DF/F (1 AP)'}},\n",
    "               'yaxis': {'anchor': 'x', 'domain': [0.0, 1.0], 'title': {'text': 'Half-rise time (1 AP)'}}}\n",
    "})\n",
    "\n",
    "\n",
    "# GUI elements\n",
    "x_dropdown = widgets.Dropdown(\n",
    "    options=plottableVars,\n",
    "    value='DF/F (1 AP)',\n",
    "    description='X axis:',\n",
    ")\n",
    "y_dropdown = widgets.Dropdown(\n",
    "    options=plottableVars,\n",
    "    value='Half-rise time (1 AP)',\n",
    "    description='Y axis:',\n",
    ")\n",
    "color_dropdown = widgets.Dropdown(\n",
    "    options=plottableVars,\n",
    "    value='Half-decay time (1 AP)',\n",
    "    description='color:',\n",
    ")\n",
    "xscale_radio = widgets.RadioButtons(\n",
    "    options=['linear', 'log'],\n",
    "    description='X scale:',\n",
    "    disabled=False\n",
    ")\n",
    "yscale_radio = widgets.RadioButtons(\n",
    "    options=['linear', 'log'],\n",
    "    description='Y scale:',\n",
    "    disabled=False\n",
    ")\n",
    "show_names_chkbx = widgets.Checkbox(\n",
    "    value=False,\n",
    "    description='Show all construct names',\n",
    "    disabled=False,\n",
    "    indent=False\n",
    ")\n",
    "\n",
    "# construct info table\n",
    "construct_table = widgets.Output(layout={'border': '1px solid black', 'width':'30%'})\n",
    "construct_info = widgets.Output(layout={'border': '1px solid black', 'width':'40%'})\n",
    "\n",
    "def response(change):\n",
    "    with g.batch_update():\n",
    "        x_val = data_filt[x_dropdown.value]\n",
    "        y_val = data_filt[y_dropdown.value]\n",
    "        \n",
    "        g.data[0]['x'] = x_val# np.log10(x_val) if is_log_xaxis else x_val\n",
    "        \n",
    "        g.data[0]['y'] = y_val# np.log10(y_val) if is_log_yaxis else y_val\n",
    "        g.data[0].marker.color = data_filt[color_dropdown.value]\n",
    "        \n",
    "        g.layout.xaxis.title.text = x_dropdown.value\n",
    "        g.layout.yaxis.title.text = y_dropdown.value\n",
    "\n",
    "        g.layout.coloraxis.colorbar.title.text = color_dropdown.value\n",
    "        g.layout.xaxis.type = xscale_radio.value\n",
    "        g.layout.yaxis.type = yscale_radio.value\n",
    "        \n",
    "        # update construct text\n",
    "        if show_names_chkbx.value:\n",
    "            # show all\n",
    "            hightlight_txt_array = data_filt.index\n",
    "            \n",
    "        else:\n",
    "            # show only highlights\n",
    "            hightlight_txt_array = [mapping_swapped.get(c_id) if (c_id in mapping_swapped.keys()) else '' for c_id in data_filt.index ]\n",
    "        g.data[0].text = hightlight_txt_array\n",
    "        \n",
    "        # update hover text\n",
    "        g.data[0].hovertemplate = '<b>%{customdata}</b><br>' + x_dropdown.value + '=%{x:.3f}<br>' + y_dropdown.value + '=%{y:.3f}<br>' + color_dropdown.value + '=%{marker.color:.3f}'\n",
    "    \n",
    "# click behavior (https://plotly.com/python/click-events/)\n",
    "def update_point(trace, points, selector):\n",
    "    construct_table.clear_output()\n",
    "    construct_info.clear_output()\n",
    "    construct_table.append_display_data(pd.DataFrame(data_filt.iloc[points.point_inds[0]][plottableVars]))\n",
    "    construct_info.append_display_data('Name: ' + data_filt.iloc[points.point_inds[0]]['Clone name'])\n",
    "    construct_info.append_display_data('Sequence: ' + data_filt.iloc[points.point_inds[0]]['AA sequence'])\n",
    "    \n",
    "    \n",
    "g.data[0].on_click(update_point)\n",
    "\n",
    "x_dropdown.observe(response, names=\"value\")\n",
    "y_dropdown.observe(response, names=\"value\")\n",
    "color_dropdown.observe(response, names=\"value\")\n",
    "xscale_radio.observe(response, names=\"value\")\n",
    "yscale_radio.observe(response, names=\"value\")\n",
    "show_names_chkbx.observe(response, names='value')\n",
    "\n",
    "scale_wdgets = widgets.HBox([xscale_radio, yscale_radio])\n",
    "dropdown_wdgts = widgets.HBox([x_dropdown, y_dropdown, color_dropdown])\n",
    "out_widgts = widgets.HBox([construct_table, construct_info])\n",
    "\n",
    "v = widgets.VBox([dropdown_wdgts, \n",
    "                  scale_wdgets, \n",
    "                  show_names_chkbx,\n",
    "                  g,\n",
    "                 out_widgts])\n",
    "v"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "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": 2
}