PlotPSTHCategoriesGUI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots PSTHs and calculates basic neural measures for Chronic data. \n",
    "\n",
    "#### Run this notebook after creating the /PlaybackPkl files inside the site folder using PlotSpikeSortedCategoriesGUI.ipynb\n",
    "\n",
    "#### You should only have to modify the rootPath in cell 2 for this notebook to work.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dependencies \n",
    "import os\n",
    "import glob\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import pickle as pk\n",
    "from scipy.stats import t\n",
    "\n",
    "# Display and GUI\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# Data Processing Stuff\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "from IPython.display import display\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Paths\n",
    "Note that it is assumed that you have data directory accessible from the rootPath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "rootPath = '/Users/frederictheunissen/Code/songephys/'\n",
    "dataPath = 'data/birds/'\n",
    "\n",
    "birds = [\n",
    "    os.path.basename(bird)\n",
    "    for bird in sorted(glob.glob(os.path.join(rootPath+dataPath, \"*\")), key=os.path.getmtime, reverse=True)\n",
    "    if (os.path.isdir(bird))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input Widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Widget to pick bird\n",
    "bird_picker = widgets.Dropdown(\n",
    "    options=birds,\n",
    "    value=birds[0],\n",
    "    description=\"Bird\",\n",
    "    disabled=False,\n",
    ")\n",
    "\n",
    "site_picker = widgets.Dropdown(\n",
    "    options=['None'],\n",
    "    value='None',\n",
    "    description=\"Site\",\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "playbackPkl_picker = widgets.Dropdown(\n",
    "    options=['None'],\n",
    "    value='None',\n",
    "    description=\"ElectCluster\",\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "stimulus_picker = widgets.Dropdown(\n",
    "    options=['None'],\n",
    "    value='None',\n",
    "    description=\"Stimulus\",\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "save_button = widgets.Button(\n",
    "    description='Save PSTH',\n",
    "    disabled=False,\n",
    "    button_style='success', \n",
    "    tooltip='Click me to save figure', \n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<>:29: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "<>:101: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "<>:29: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "<>:101: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "/var/folders/zn/qkxxn59x2712k0l7j64f7ry40000gn/T/ipykernel_49333/1094220748.py:29: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "  if site is not 'None':\n",
      "/var/folders/zn/qkxxn59x2712k0l7j64f7ry40000gn/T/ipykernel_49333/1094220748.py:101: SyntaxWarning: \"is not\" with a literal. Did you mean \"!=\"?\n",
      "  if playbackPkl is not 'None':\n"
     ]
    }
   ],
   "source": [
    "def load_bird(change):\n",
    "    \n",
    "    bird = change['new']\n",
    "    if bird is not load_bird._loaded[0]:\n",
    "        sites = [ \n",
    "            os.path.basename(site)\n",
    "            for site in glob.glob(os.path.join(rootPath+dataPath+bird+'/sites/', \"*\"))\n",
    "            ]\n",
    "\n",
    "        load_bird._loaded = (bird, 'None', 'None', 'None')\n",
    "        site_picker.options=['None'] + sites\n",
    "        site_picker.value='None'\n",
    "        playbackPkl_picker.options=['None'] \n",
    "        playbackPkl_picker.value='None'\n",
    "        stimulus_picker.options=['None'] \n",
    "        stimulus_picker.value='None'\n",
    "        \n",
    "        \n",
    "    return load_bird._loaded\n",
    "\n",
    "\n",
    "def load_site(change):\n",
    "\n",
    "    global playPklPath    \n",
    "    \n",
    "    bird = bird_picker.value\n",
    "    site = change['new']\n",
    "    if (bird, site) is not (load_site._loaded[0],load_site._loaded[1]):\n",
    "        if site is not 'None':\n",
    "            # Set Paths\n",
    "            playPklPath = rootPath + dataPath + bird + '/sites/' + site + '/PlaybackPkl/'\n",
    "    \n",
    "            # Set up list of Playback files available for this site\n",
    "            playPkls = [ \n",
    "                os.path.basename(efile)\n",
    "                for efile in glob.glob(playPklPath + '*.pkl')\n",
    "                ]\n",
    "            playPkls.sort()\n",
    "            playbackPkl_picker.options=['None']+playPkls \n",
    "        else:\n",
    "            playPklPath = None\n",
    "            electrode_picker.options=['None']\n",
    "            \n",
    "        playbackPkl_picker.value='None'\n",
    "        stimulus_picker.options=['None'] \n",
    "        stimulus_picker.value='None'\n",
    "\n",
    "        load_site._loaded = (bird, site, 'None', 'None')\n",
    "    \n",
    "    return load_site._loaded\n",
    "\n",
    "def calc_zscore(stimName):\n",
    "    global dfRelTime\n",
    "    \n",
    "    rows = []\n",
    "    # This loop finds all rows that match but there should only be one.\n",
    "    for index, row in dfRelTime.iterrows():\n",
    "        if row['file'] == stimName:\n",
    "            rows.append(row)\n",
    "    if (len(rows) != 1):\n",
    "        print('Stimulus not found or too many')\n",
    "        return 0.0, 1.0\n",
    "    \n",
    "    # Plot the stimulus\n",
    "    row = rows[0]\n",
    "    \n",
    "    # Calculates response diff for each stim\n",
    "    spikeDiff = np.zeros(row['nTrials'])\n",
    "    for it in range(row['nTrials']):\n",
    "        spikeDiff[it] = np.sum((row['spikeTimes'][it] >= 0) & (row['spikeTimes'][it] < 0.5)) - np.sum((row['spikeTimes'][it] >= -0.5) & (row['spikeTimes'][it] < 0)) \n",
    "           \n",
    "    # Calculate z-score and pvalue\n",
    "    if (row['nTrials'] > 1) :\n",
    "        sdiffSD = np.std(spikeDiff, ddof=1)\n",
    "        if sdiffSD == 0:\n",
    "            spikeDiff[0] += 1\n",
    "            sdiffSD = np.std(spikeDiff, ddof=1)\n",
    "\n",
    "        zscore = np.mean(spikeDiff)/sdiffSD\n",
    "        if (zscore < 0.0):\n",
    "            pvalue = (t.cdf(zscore*np.sqrt(row['nTrials']), row['nTrials']))*2.0\n",
    "        else:\n",
    "            pvalue = (1.0 - t.cdf(zscore*np.sqrt(row['nTrials']), row['nTrials']))*2.0\n",
    "    else:\n",
    "        zscore = spikeDiff\n",
    "        pvalue = 1.0\n",
    "        \n",
    "    return zscore, pvalue\n",
    "\n",
    "    \n",
    "def load_playbackPkl(change):\n",
    "    global playPklPath\n",
    "    global unitInfo, dfRelTime\n",
    "    \n",
    "    bird = bird_picker.value\n",
    "    site = site_picker.value\n",
    "    playbackPkl = change['new']\n",
    "    \n",
    "    if (bird, site, playbackPkl) is not (load_playbackPkl._loaded[0],load_playbackPkl._loaded[1],load_playbackPkl._loaded[2]) :\n",
    "        \n",
    "        if playbackPkl is not 'None':            \n",
    "            # Load unitInfo and data frames \n",
    "            pklFile = playPklPath + playbackPkl\n",
    "            fileIn = open(pklFile, 'rb')\n",
    "            unitInfo = pk.load(fileIn)\n",
    "            dfAbsTime = pk.load(fileIn)\n",
    "            dfRelTime = pk.load(fileIn)\n",
    "            fileIn.close()\n",
    "            \n",
    "            stims = list(dfRelTime['file'])        \n",
    "            stimulus_picker.options=['None']+stims\n",
    "        else:\n",
    "            stimulus_picker.options=['None']\n",
    "        \n",
    "        stimulus_picker.value='None'       \n",
    "        load_playbackPkl._loaded = (bird, site, playbackPkl, 'None')\n",
    "        \n",
    "    return load_playbackPkl._loaded\n",
    "\n",
    "def load_stimulus(change):\n",
    "    global dfRelTime\n",
    "    \n",
    "    bird = bird_picker.value\n",
    "    site = site_picker.value\n",
    "    playbackPkl = playbackPkl_picker.value\n",
    "    stimulus = change['new']\n",
    "\n",
    "    load_stimulus._loaded = (bird, site, playbackPkl, stimulus)\n",
    "                        \n",
    "    return load_stimulus._loaded\n",
    "       \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize Input Widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_bird._loaded = ('None', 'None', 'None', 'None')\n",
    "load_site._loaded = ('None', 'None', 'None', 'None')\n",
    "load_playbackPkl._loaded = ('None', 'None', 'None', 'None')\n",
    "load_stimulus._loaded = ('None', 'None', 'None', 'None')\n",
    "\n",
    "load_bird({'type' : 'change', 'new': birds[0], 'old': 'None'} )\n",
    "bird_picker.observe(load_bird, 'value')\n",
    "site_picker.observe(load_site, 'value')\n",
    "playbackPkl_picker.observe(load_playbackPkl, 'value')\n",
    "stimulus_picker.observe(load_stimulus, 'value')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output Widget for Plot Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotPSTH_output = widgets.Output()\n",
    "plotMean_output = widgets.Output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "figPSTH = None\n",
    "def plot_PSTH(stimName):\n",
    "    \n",
    "    global dfRelTime, figPSTH\n",
    "    \n",
    "    \n",
    "    if (stimName == 'None'):\n",
    "        print('No stimulus chosen to plot')\n",
    "        return\n",
    "    \n",
    "    rows = []\n",
    "    # This loop finds all rows that match but there should only be one.\n",
    "    for index, row in dfRelTime.iterrows():\n",
    "        if row['file'] == stimName:\n",
    "            rows.append(row)\n",
    "    if (len(rows) != 1):\n",
    "        print('Stimulus not found or too many')\n",
    "        return\n",
    "    \n",
    "    # Plot the stimulus\n",
    "    row = rows[0]\n",
    "    if figPSTH is not None:\n",
    "        plt.close(figPSTH)\n",
    "        \n",
    "    figPSTH = plt.figure()\n",
    "    aSound = plt.axes([0.1, 0.85, 0.87, 0.15])\n",
    "    aSound.plot(row['tStim'], row['stimWav'])       \n",
    "        \n",
    "    # Plot Microphone and spikes for each trial \n",
    "    spikeDiff = np.zeros(row['nTrials'])\n",
    "    spikeTot = np.zeros(row['nTrials'])\n",
    "    for it in range(row['nTrials']):\n",
    "        aTrial = plt.axes([0.1, 0.85-0.6*(it+1)/row['nTrials'], 0.87, 0.6/row['nTrials']])\n",
    "        aTrial.plot(row['tMic'], row['micWav'][it], 'silver', lw=0.5)\n",
    "        aTrial.get_yaxis().set_ticks([])\n",
    "        \n",
    "        # Spike difference first 500 ms\n",
    "        spikeDiff[it] = np.sum((row['spikeTimes'][it] >= 0) & (row['spikeTimes'][it] < 0.5)) - np.sum((row['spikeTimes'][it] >= -0.5) & (row['spikeTimes'][it] < 0)) \n",
    "        \n",
    "        # Total number of spikes in first 100 ms\n",
    "        spikeTot[it] = np.sum((row['spikeTimes'][it] >= 0) & (row['spikeTimes'][it] < 1.0))\n",
    "        micMax = row['micWav'][it].max()\n",
    "        for tval in row['spikeTimes'][it]:\n",
    "            aTrial.plot([tval, tval], [0, micMax*0.8], 'k') \n",
    "            \n",
    "    # Calculate z-score and pvalue\n",
    "    if (row['nTrials'] > 1) :\n",
    "        sdiffSD = np.std(spikeDiff, ddof=1)\n",
    "        if sdiffSD == 0:\n",
    "            titleStr = 'Spike Diff is %f in all trials' % np.mean(spikeDiff)\n",
    "        else:\n",
    "            zscore = np.mean(spikeDiff)/sdiffSD\n",
    "            if (zscore < 0.0):\n",
    "                pvalue = (t.cdf(zscore*np.sqrt(row['nTrials']), row['nTrials']))*2.0\n",
    "            else:\n",
    "                pvalue = (1.0 - t.cdf(zscore*np.sqrt(row['nTrials']), row['nTrials']))*2.0\n",
    "            titleStr = 'z = %.3f p = %.4f'% (zscore, pvalue)\n",
    "    else:\n",
    "        titleStr = 'Spike Diff %f' % spikeDiff[0]\n",
    "    \n",
    "    plotPSTH_output.clear_output()\n",
    "    with plotPSTH_output:\n",
    "        print(titleStr)\n",
    "            \n",
    "    # Plot the average rate from KDE (smoothed PSTH)\n",
    "    aRate = plt.axes([0.1, 0.1, 0.87, 0.15])\n",
    "    aRate.plot(row['tKDE'], row['spikeKDE'])\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Spikes/s')\n",
    "    plt.show()\n",
    "    \n",
    "    # Check kde estimate\n",
    "    spikeIntegral = sum(row['spikeKDE'][(row['tKDE']>=0) & (row['tKDE']<1.0)])\n",
    "    spikeIntegral *= row['tKDE'][1]-row['tKDE'][0]\n",
    "\n",
    "    print('KDE rate estimate', spikeIntegral)\n",
    "    print('Actual rate', np.mean(spikeTot))\n",
    "    return\n",
    "\n",
    "def print_plot_PSTH(b):   \n",
    "    global playPklPath, figPSTH\n",
    "    \n",
    "    stimName = load_stimulus._loaded[3]\n",
    "    \n",
    "    if (stimName == 'None'):\n",
    "        return\n",
    "\n",
    "    clusterName = os.path.splitext(load_stimulus._loaded[2])[0]\n",
    "    stimName = os.path.splitext(load_stimulus._loaded[3])[0]\n",
    "    \n",
    "    figName = playPklPath + clusterName + '_' + stimName + '.pdf'        \n",
    "    figPSTH.savefig(figName)\n",
    "    \n",
    "\n",
    "def plot_meanRate(playbackPklName):  \n",
    "    global dfRelTime\n",
    "    kWidth = 0.03    # 30 ms width\n",
    "    \n",
    "    if (playbackPklName == 'None'):\n",
    "        print('No stimulus chosen to plot')\n",
    "        return\n",
    "    \n",
    "    allSpikes = []\n",
    "    # Loop through all stims and trials.\n",
    "    nTotal = 0\n",
    "    for index, row in dfRelTime.iterrows():\n",
    "        nTotal += row['nTrials']\n",
    "        for it in range(row['nTrials']):\n",
    "            allSpikes.append(list(row['spikeTimes'][it]))\n",
    "            \n",
    "    spikeDiff = np.zeros(nTotal)\n",
    "    itot = 0\n",
    "    for index, row in dfRelTime.iterrows():\n",
    "        for it in range(row['nTrials']):        \n",
    "            # Spike difference first 500 ms\n",
    "            spikeDiff[itot] = np.sum((row['spikeTimes'][it] >= 0) & (row['spikeTimes'][it] < 0.5)) - np.sum((row['spikeTimes'][it] >= -0.5) & (row['spikeTimes'][it] < 0)) \n",
    "            itot += 1\n",
    "\n",
    "    \n",
    "    if (len(allSpikes) == 0):\n",
    "        # When there is no data\n",
    "        tRate = np.linspace(0, 1.0, num=1001)\n",
    "        meanSpikeDens = np.zeros(tRate.shape)            \n",
    "    else:  \n",
    "        # KDE estimate of mean rate\n",
    "        allSpikes = np.concatenate(allSpikes).reshape(-1,1)\n",
    "        # The time axis for the mean rate - same for all rows - taking the last one\n",
    "        tRate = row['tKDE'].reshape(-1,1)\n",
    "        tInt = tRate[-1]-tRate[0]\n",
    "        sampRate = 1.0/(tRate[1]-tRate[0])\n",
    "        kdeMean = KernelDensity(kernel='gaussian', bandwidth=kWidth).fit(allSpikes)\n",
    "    \n",
    "        # Generate the kde rate.        \n",
    "        meanSpikeDens = np.exp(kdeMean.score_samples(tRate))\n",
    "        meanSpikeDens = meanSpikeDens*len(allSpikes)*sampRate/(sum(meanSpikeDens)*nTotal)\n",
    "\n",
    "    # Calculate z-score and pvalue\n",
    "    if (nTotal > 1) :\n",
    "        sdiffSD = np.std(spikeDiff, ddof=1)\n",
    "        if sdiffSD == 0:\n",
    "            titleStr = 'Spike Diff is %f in all trials' % np.mean(spikeDiff)\n",
    "        else:\n",
    "            zscore = np.mean(spikeDiff)/sdiffSD\n",
    "            if (zscore < 0.0):\n",
    "                pvalue = (t.cdf(zscore*np.sqrt(nTotal), nTotal))*2.0\n",
    "            else:\n",
    "                pvalue = (1.0 - t.cdf(zscore*np.sqrt(nTotal), nTotal))*2.0\n",
    "            titleStr = 'z = %.3f p = %.4f'% (zscore, pvalue)\n",
    "    elif nTotal == 1:\n",
    "        titleStr = 'Spike Diff %f' % spikeDiff[0]\n",
    "    else:\n",
    "        titleStr = 'No data'\n",
    "    \n",
    "    plotMean_output.clear_output()\n",
    "    with plotMean_output:\n",
    "        print(titleStr)\n",
    "\n",
    "    \n",
    "    plt.plot(tRate, meanSpikeDens)\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Spike Rate (spikes/s)')\n",
    "    plt.title('Playback Evoked Firing Rate')\n",
    "    plt.show()\n",
    "    return\n",
    "\n",
    "save_button.on_click(print_plot_PSTH)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output Widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Widget to plot rate with a given threhold\n",
    "PSTH_plot = widgets.interactive_output(\n",
    "    plot_PSTH, \n",
    "    {\n",
    "        \"stimName\" : stimulus_picker\n",
    "    }\n",
    ")\n",
    "\n",
    "mean_plot = widgets.interactive_output(\n",
    "    plot_meanRate, \n",
    "    {\n",
    "        \"playbackPklName\" : playbackPkl_picker\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1db684eedfd44061a3a79ee1b81ee94c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(VBox(children=(Dropdown(description='Bird', options=('ZF17M_ZF18F', 'ZF4F', 'ZF1…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "widgets.VBox([\n",
    "    widgets.HBox([widgets.VBox([\n",
    "        bird_picker,\n",
    "        site_picker,\n",
    "        playbackPkl_picker,\n",
    "        plotMean_output,\n",
    "        stimulus_picker,\n",
    "        plotPSTH_output\n",
    "    ]), mean_plot ]),\n",
    "    PSTH_plot,\n",
    "    save_button\n",
    "\n",
    "])"
   ]
  },
  {
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}