{ "cells": [ { "cell_type": "markdown", "id": "32647951-7dcf-4553-b42a-5d8508d95cf9", "metadata": {}, "source": [ "\"Open   " ] }, { "cell_type": "markdown", "id": "8f701b03-52cb-4215-979a-1ccaa3839703", "metadata": {}, "source": [ "# Data Explorer\n", "\n", "\n", "\n", "If you would like to explore the analysis for this lab, but do not have data, you can download examples for the following experiments using the linked shared files:" ] }, { "cell_type": "markdown", "id": "f07451ff-cf56-4763-b4ac-bd099ba1f145", "metadata": {}, "source": [ "\n", "# Table of Contents\n", "\n", "- [Introduction](#intro)\n", "- [Setup](#setup)\n", "- [Analysis tool 1](#one)\n", "- [Analysis tool 2](#two)\n" ] }, { "cell_type": "markdown", "id": "488d476b-9387-43d7-b4bf-0528d9f4950a", "metadata": {}, "source": [ "\n", "# Setup\n", "\n", "[toc](#toc)" ] }, { "cell_type": "markdown", "id": "b94f344a-8c2a-4080-b972-8e7a30c91e66", "metadata": {}, "source": [ "Import and define functions" ] }, { "cell_type": "code", "execution_count": null, "id": "be70dacf-b862-4afa-8fc0-dce8c363fdf9", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown Run this code cell to import packages and define functions \n", "from pathlib import Path\n", "import random\n", "import numpy as np\n", "import pandas as pd\n", "import plotly.express as px\n", "import plotly.graph_objects as go\n", "from plotly.subplots import make_subplots\n", "from scipy import ndimage, optimize, signal\n", "from scipy.optimize import curve_fit\n", "from scipy.signal import hilbert,medfilt,resample, find_peaks,butter\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "from sklearn.decomposition import PCA\n", "from sklearn.cluster import KMeans\n", "from datetime import datetime,timezone,timedelta\n", "pal = sns.color_palette(n_colors=15)\n", "pal = pal.as_hex()\n", "from numpy import NaN\n", "\n", "from matplotlib.ticker import (AutoMinorLocator, MultipleLocator)\n", "from ipywidgets import widgets, interact, interactive,interactive_output\n", "%config InlineBackend.figure_format = 'retina'\n", "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n", "\n", "def monoExp(x, m, t, b):\n", " return m * np.exp(-x / t) + b\n", "\n", "\n", "print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))" ] }, { "cell_type": "markdown", "id": "f91a09e1-c305-4809-8302-c1fad4140efa", "metadata": {}, "source": [ "Mount Google Drive" ] }, { "cell_type": "code", "execution_count": null, "id": "2b3a4500-470b-44e3-8572-a0d55f3bc1cb", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown Run this cell to mount your Google Drive.\n", "\n", "from google.colab import drive\n", "drive.mount('/content/drive')\n", "\n", "print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))" ] }, { "cell_type": "markdown", "id": "a85d4e3e-b775-4dc2-8b8f-89878c97e750", "metadata": {}, "source": [ "\n", "\n", "# Analysis tool 1: adaptation rate\n", "\n", "Use this analysis tool for experiment files from steps \\#2-4 in the Lab Manual (the effect of stimulus rate on the ERG response). \n", "\n", "First, you will tweak the parameters of the analysis tool so that the peak ERG amplitude on each trial is automatically calculated for the entire file at once. \n", " > One of the parameters you need for this processing step (\"time to peak (s)\") comes from using the Dash Data Explorer to determine the rise time of the ERG to a 10msec pulse. \n", "\n", "Then, you will programatically (ie. automatically) write the processing results to a csv file for each stimulus rate condition. \n", "\n", "Finally, you will programatically load the csv files from all conditions to create an overlaid scatter plot of the results.\n", "\n" ] }, { "cell_type": "markdown", "id": "59f3dc95-c2fe-4221-bcd5-39335118d74e", "metadata": {}, "source": [ "## 1. Import raw data and measure ERG amplitude\n", "\n", "Raw data should be in a format digitized with *Nidaq USB6211* and recorded using *Bonsai-rx* as a *.bin* file" ] }, { "cell_type": "code", "execution_count": null, "id": "8e71acc5-b8df-4e0f-a74c-2f27d0f5be59", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "\n", "#@markdown Specify the file path \n", "#@markdown to your recorded data in the colab runtime (find the filepath in the colab file manager):\n", "\n", "filepath = \"full filepath goes here\" #@param \n", "# filepath = '/Volumes/Untitled/BIOL247/data/crayfish-erg/KP_20221104/rate-5hz_0.bin'\n", "\n", "#@markdown Specify the sampling rate and number of channels recorded.\n", "\n", "sampling_rate = None #@param\n", "number_channels = None #@param\n", "\n", "# sampling_rate = 30000 #@param\n", "# number_channels = 2 #@param\n", "\n", "# downsample = False #@param\n", "# newfs = 10000 #@param\n", "\n", "#@markdown After you have filled out all form fields, \n", "#@markdown run this code cell to load the data. \n", "\n", "filepath = Path(filepath)\n", "\n", "# No need to edit below this line\n", "#################################\n", "data = np.fromfile(Path(filepath), dtype = np.float64)\n", "data = data.reshape(-1,number_channels)\n", "data_dur = np.shape(data)[0]/sampling_rate\n", "print('duration of recording was %0.2f seconds' %data_dur)\n", "\n", "fs = sampling_rate\n", "# if downsample:\n", "# # newfs = 10000 #downsample emg data\n", "# chunksize = int(sampling_rate/newfs)\n", "# data = data[0::chunksize,:]\n", "# fs = int(np.shape(data)[0]/data_dur)\n", "\n", "time = np.linspace(0,data_dur,np.shape(data)[0])\n", "\n", "print('Data upload completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))\n", "\n", "#@markdown Run this cell to create an interactive plot with a slider to scroll \n", "#@markdown through the signal\n", "#@markdown and set an appropriate event detection threshold and noise threshold. \n", "\n", "slider_xrange = widgets.FloatRangeSlider(\n", " min=0,\n", " max=data_dur,\n", " value=(0,data_dur),\n", " step= 0.5,\n", " readout=True,\n", " continuous_update=False,\n", " description='Time Range (s)',\n", " style = {'description_width': '200px'})\n", "slider_xrange.layout.width = '600px'\n", "\n", "slider_yrange = widgets.FloatRangeSlider(\n", " min=-5,\n", " max=5,\n", " value=[-2,1],\n", " step=0.05,\n", " continuous_update=False,\n", " readout=True,\n", " description='yrange',\n", " style = {'description_width': '200px'})\n", "slider_yrange.layout.width = '600px'\n", "\n", "slider_text = widgets.Text(\n", " value='0.04',\n", " placeholder='0.04',\n", " description='time to peak (s)',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "select_channel = widgets.Select(\n", " options=np.arange(np.shape(data)[1]), # start with a single trial on a single bout... it will update when runs ; old: np.arange(len(trial_times)),\n", " value=0,\n", " #rows=10,\n", " description='Channel',\n", " style = {'description_width': '50px'},\n", " disabled=False\n", ")\n", "select_channel.layout.width = '100px'\n", "select_channel.layout.height = '50px'\n", "\n", "# slider_threshold = widgets.FloatSlider(\n", "# min=0,\n", "# max=5,\n", "# value=0.2,\n", "# step=0.005,\n", "# readout_format='.3f',\n", "# continuous_update=False,\n", "# readout=True,\n", "# description='event detection threshold',\n", "# style = {'description_width': '200px'})\n", "# slider_threshold.layout.width = '600px'\n", "\n", "radio_polarity = widgets.RadioButtons(\n", " options=[1, -1],\n", " value=-1,\n", " description='peaks polarity',\n", " layout={'width': 'max-content'},\n", " style = {'description_width': '100px'},\n", " disabled=False\n", ")\n", "radio_polarity.layout.width = '300px'\n", "\n", "# ui_peaks = widgets.HBox([select_channel, radio_polarity, slider_threshold])\n", "\n", "# detect_type_radio = widgets.RadioButtons(\n", "# options=['peak', 'level crossing positive', 'level crossing negative'],\n", "# value='level crossing positive', # Defaults to 'level crossing'\n", "# layout={'width': 'max-content'}, # If the items' names are long\n", "# description='Type of event detection',\n", "# style = {'description_width': '200px'},\n", "# disabled=False\n", "# )\n", "\n", "iei_text = widgets.Text(\n", " value='0.01',\n", " placeholder='0.01',\n", " description='minimum interval (s)',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "\n", "noise_amp_text = widgets.Text(\n", " value='0.06',\n", " placeholder='0.06',\n", " description='threshold amplitude (V)',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "\n", "def update_plot(chan_,xrange,yrange,polarity,iei,noise_amp,onset):\n", " fig, ax = plt.subplots(figsize=(10,5),num=1); #specify figure number so that it does not keep creating new ones\n", " \n", " signal = data[:,chan_].reshape(-1)\n", " # signal = signal-np.median(signal)\n", "\n", " \n", " iei = float(iei)\n", " noise_amp = float(noise_amp)\n", " onset = float(onset)\n", " prepeak_base = int(onset*fs)\n", " \n", " if iei>1/fs:\n", " d = iei*fs #minimum time allowed between distinct events\n", "\n", " # if thresh_ >=0:\n", " # r = find_peaks(signal,height=thresh_,distance=d)\n", " # if thresh_ <0:\n", " # r = find_peaks(-1*signal,height=-1*thresh_,distance=d)\n", " # trial_times = r[0]/fs\n", "\n", " # r_thresh = find_peaks(signal,height=thresh_,distance=d)\n", " r_prom = find_peaks(polarity * signal,distance=d,prominence=noise_amp)\n", " r_ind = r_prom[0] #np.intersect1d(r_thresh[0],r_prom[0])\n", " # lagging_peak_amp = r[1]['peak_heights']\n", " peak_amp = signal[r_ind]\n", " base_amp = signal[r_ind-prepeak_base]\n", " # ax.hlines(thresh_, xrange[0],xrange[1],linestyle='--',color='green')\n", "\n", " peak_times = np.asarray([np.round(s/fs,2) for s in r_ind])\n", "\n", " starti = int(xrange[0]*fs)+1\n", " stopi = int(xrange[1]*fs)-1\n", " ax.plot(time[starti:stopi], signal[starti:stopi], color='black')\n", "\n", "\n", " inwin_inds = (peak_times>(xrange[0])) & (peak_times<(xrange[1]))\n", " ax.scatter(peak_times[inwin_inds],peak_amp[inwin_inds], zorder=3,color='red',s=20)\n", " ax.scatter(peak_times[inwin_inds]-onset,base_amp[inwin_inds], zorder=3,color='green',s=20)\n", "\n", " ax.set_ylim(yrange[0],yrange[1])\n", " ax.set_xlim(xrange[0],xrange[1])\n", "\n", " ax.xaxis.set_minor_locator(AutoMinorLocator(5))\n", "\n", "\n", " return peak_times-np.min(peak_times), peak_amp-base_amp\n", "\n", "# w_trials_ = interactive(update_plot, chan_=select_channel, \n", "# xrange=slider_xrange, \n", "# thresh_=slider_threshold, iei = iei_text, noise_amp = noise_amp_text);\n", "# display(w_trials_)\n", "w_peak_detect_all = interactive(update_plot, chan_=select_channel, \n", " yrange=slider_yrange,xrange=slider_xrange,\n", " polarity = radio_polarity, iei = iei_text,noise_amp=noise_amp_text,\n", " onset=slider_text);\n", "display(w_peak_detect_all)\n", "\n", "# w_peak_detect_all = interactive(update_plot, {'chan_':select_channel,\n", "# 'yrange':slider_yrange, \n", "# 'xrange':slider_xrange,\n", "# 'thresh_':slider_threshold,\n", "# 'polarity':radio_polarity,\n", "# 'iei':iei_text,\n", "# 'noise_amp':noise_amp_text});\n", "# display(ui_peaks,widgets.HBox([iei_text,noise_amp_text]),w_peak_detect_all,ui_range)" ] }, { "cell_type": "markdown", "id": "3fd41f4c-1b4e-420c-9118-037970743e42", "metadata": {}, "source": [ "## 2. Save processed data as .csv\n", "\n", "Specify the stimulus rate for the data you just processed. The rate should be a string (ie. surrounded by quotes).\n", "> It will take a minute or so to process and actually save, but the .csv file will be saved to the local \"content\" folder of Google Collaboratory. You can download them individually to your own computer (or move them to your Google Drive) to use again later without re-processing the data. " ] }, { "cell_type": "code", "execution_count": null, "id": "5674bf5d-b89c-4d6e-a2d6-98aac49333ea", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "rate = 'None' #@param\n", "\n", "peak_times,peak_amp = w_peak_detect_all.result\n", "\n", "df = pd.DataFrame({'time':peak_times,'amp':-peak_amp,'number':np.arange(0,len(peak_amp),1),'rate':rate})\n", "\n", "df.to_csv('adaptation_' + rate + '.csv')" ] }, { "cell_type": "markdown", "id": "2972d387-fde3-4c49-9ca1-0accfcf6afce", "metadata": { "tags": [] }, "source": [ "## 3. Compare across conditions \n", "\n", "Provide a list of the csv files generated in the last step. The code will iterate through each data file and plot the results. You will also see the result of a mono-exponential fit to the data." ] }, { "cell_type": "code", "execution_count": null, "id": "915bf086-0617-42b5-9fc6-a64438a6f4fa", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown list the full path to each csv file created using the \"adaptation analysis tool.\"\n", "file_list = ['file1_.csv', 'file1_.csv', 'file1_.csv'] #@param\n", "# file_list = ['adaptation_2_b_.csv','adaptation_5_b_.csv','adaptation_10_b_.csv']\n", "\n", "x_axis_limit_number = None #@param\n", "x_axis_limit_time = None #@param\n", "\n", "df = pd.DataFrame({})\n", "for f in file_list:\n", " df_ = pd.read_csv(f)\n", " df = pd.concat([df,df_])\n", "\n", "hfig,ax = plt.subplots(figsize=(10,4))\n", "g=sns.scatterplot(x='number',y='amp',data=df, hue='rate');\n", "for rate in df['rate'].unique():\n", " xs = df[df['rate']==rate]['number']\n", " ys = df[df['rate']==rate]['amp']\n", " # p0 = [10, 2, np.min(ys)] # start with values near those we expect\n", " params_, cv = curve_fit(monoExp, xs, ys,maxfev = 5000) #monoexponential\n", " m, t, c = params_\n", " y_fit = monoExp(xs, m, t, c)#+np.min(ys)\n", "\n", " # ax.scatter(xs,-ys)\n", " plt.plot(xs,y_fit,label=f'm = {m:.2f}, tau = {t:.2f}, c = {c:.2f}',linestyle='--')\n", " plt.legend()\n", "g.legend(loc='center left', bbox_to_anchor=(1, 0.5), ncol=1);\n", "ax.set_xlim(-1,x_axis_limit_number)\n", "\n", "hfig,ax = plt.subplots(figsize=(10,4))\n", "g=sns.scatterplot(x='time',y='amp',data=df, hue='rate');\n", "for rate in df['rate'].unique():\n", " xs = df[df['rate']==rate]['time']\n", " ys = df[df['rate']==rate]['amp']\n", " # p0 = [10, 2, np.min(ys)] # start with values near those we expect\n", " params_, cv = curve_fit(monoExp, xs, ys,maxfev = 5000) #monoexponential\n", " m, t, c = params_\n", " y_fit = monoExp(xs, m, t, c)#+np.min(ys)\n", "\n", " # ax.scatter(xs,-ys)\n", " plt.plot(xs,y_fit,label=f'm = {m:.2f}, tau = {t:.2f}, c = {c:.2f}',linestyle='--')\n", " plt.legend()\n", "g.legend(loc='center left', bbox_to_anchor=(1, 0.5), ncol=1);\n", "ax.set_xlim(-1,x_axis_limit_time)" ] }, { "cell_type": "markdown", "id": "4dc1abe5-679c-493d-ad28-5840d8965523", "metadata": {}, "source": [ "# Analysis tool 2: Individual trial analysis\n", "\n", "Use this analysis tool for experiment files from steps \\#5 and 6 in the Lab Manual." ] }, { "cell_type": "markdown", "id": "8b30f077-df2c-47fb-8021-7f2370064fa3", "metadata": { "tags": [] }, "source": [ "\n", "\n", "## 1. Import data \n", "\n", "Import data digitized with *Nidaq USB6211* and recorded using *Bonsai-rx* as a *.bin* file" ] }, { "cell_type": "code", "execution_count": null, "id": "ac805ac8-8b1b-4ad5-9cd7-134e2d6ad519", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\" }\n", "\n", "#@markdown Specify the file path \n", "#@markdown to your recorded data in the colab runtime (find the filepath in the colab file manager):\n", "\n", "filepath = \"full filepath goes here\" #@param \n", "# filepath = '/Volumes/Untitled/BIOL247/data/crayfish-erg/KP_20221104/paired-pulse_0.bin'\n", "\n", "#@markdown Specify the sampling rate and number of channels recorded.\n", "\n", "sampling_rate = None #@param\n", "number_channels = None #@param\n", "\n", "# sampling_rate = 30000 #@param\n", "# number_channels = 2 #@param\n", "\n", "# downsample = False #@param\n", "# newfs = 10000 #@param\n", "\n", "#@markdown After you have filled out all form fields, \n", "#@markdown run this code cell to load the data. \n", "\n", "filepath = Path(filepath)\n", "\n", "# No need to edit below this line\n", "#################################\n", "data = np.fromfile(Path(filepath), dtype = np.float64)\n", "data = data.reshape(-1,number_channels)\n", "data_dur = np.shape(data)[0]/sampling_rate\n", "print('duration of recording was %0.2f seconds' %data_dur)\n", "\n", "fs = sampling_rate\n", "# if downsample:\n", "# # newfs = 10000 #downsample emg data\n", "# chunksize = int(sampling_rate/newfs)\n", "# data = data[0::chunksize,:]\n", "# fs = int(np.shape(data)[0]/data_dur)\n", "\n", "time = np.linspace(0,data_dur,np.shape(data)[0])\n", "\n", "print('Data upload completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))" ] }, { "cell_type": "markdown", "id": "71afe800-9977-450b-b2b3-821dbedced22", "metadata": {}, "source": [ "## Plot raw data\n", "\n", "This visualization tool is not necessary to go on, but provides a sanity check that you have loaded the data correctly before going on to process and analyze the results. " ] }, { "cell_type": "code", "execution_count": null, "id": "a35f0827-e477-41b8-af14-9d8e7171dff5", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown Run this code cell to plot the imported data.
\n", "#@markdown Use the range slider to scroll through the data in time.\n", "#@markdown Use the channel slider to choose which channel to plot\n", "#@markdown Be patient with the range refresh... the more data you are plotting the slower it will be. \n", "\n", "slider_xrange = widgets.FloatRangeSlider(\n", " min=0,\n", " max=data_dur,\n", " value=(0,1),\n", " step= 1,\n", " readout=True,\n", " continuous_update=False,\n", " description='Time Range (s)')\n", "slider_xrange.layout.width = '600px'\n", "\n", "slider_chan = widgets.IntSlider(\n", " min=0,\n", " max=number_channels-1,\n", " value=0,\n", " step= 1,\n", " continuous_update=False,\n", " description='channel')\n", "slider_chan.layout.width = '300px'\n", "\n", "# a function that will modify the xaxis range\n", "def update_plot(x,chan):\n", " fig, ax = plt.subplots(figsize=(10,5),num=1); #specify figure number so that it does not keep creating new ones\n", " starti = int(x[0]*fs)\n", " stopi = int(x[1]*fs)\n", " ax.plot(time[starti:stopi], data[starti:stopi,chan])\n", "\n", "w = interact(update_plot, x=slider_xrange, chan=slider_chan);" ] }, { "cell_type": "markdown", "id": "3d7706ec-acec-4728-b127-0ac07432dc07", "metadata": {}, "source": [ "For a more extensive ***RAW*** Data Explorer than the one provided in the above figure, use the [DataExplorer.py](https://raw.githubusercontent.com/neurologic/Neurophysiology-Lab/main/howto/Data-Explorer.py) application found in the [howto section](https://neurologic.github.io/Neurophysiology-Lab/howto/Dash-Data-Explorer.html) of the course website." ] }, { "cell_type": "markdown", "id": "301914d7-595e-42ae-aba0-87e77d5cbe98", "metadata": {}, "source": [ "## 2. Define bout and stimulus times\n", "\n", "The time between stimulus onset and the evoked neural event (in this case the LFP), the duration of each stimulus pulse, and the time between different stimulus pulses are critical parameters of the data on each trial. \n", "\n", "Our first task in processing and analyzing data from this experiment is to figure out the stimulus onset times. You can then segment the data in to separate bouts if the raw recording was not one continuous successful protocol. " ] }, { "cell_type": "markdown", "id": "0a7b5528-f6fd-4a8a-9cd5-bceb29af6324", "metadata": {}, "source": [ "### Define stimulus times\n", "\n", "In this experiment, stimulus times should be defined using the stimulus monitor because you will be analyzing the retinal LFP signal evoked by each stimulus pulse. " ] }, { "cell_type": "code", "execution_count": null, "id": "cdaef195-8789-40d4-a65a-95acfc8a3588", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown Run this cell to create an interactive plot with a slider to scroll \n", "#@markdown through the signal\n", "#@markdown and set an appropriate event detection threshold \n", "#@markdown (you can do so based on level crossing or peaks). \n", "\n", "slider_xrange = widgets.FloatRangeSlider(\n", " min=0,\n", " max=data_dur,\n", " value=(0,1),\n", " step= 0.05,\n", " readout=True,\n", " continuous_update=False,\n", " description='Time Range (s)',\n", " style = {'description_width': '200px'})\n", "slider_xrange.layout.width = '600px'\n", "\n", "# slider_yrange = widgets.FloatRangeSlider(\n", "# min=np.min(stim)-0.5,\n", "# max=np.max(stim)+0.5,\n", "# value=[np.min(stim),np.max(stim)],\n", "# step=0.05,\n", "# continuous_update=False,\n", "# readout=True,\n", "# description='yrange',\n", "# style = {'description_width': '200px'})\n", "# slider_yrange.layout.width = '600px'\n", "\n", "select_channel = widgets.Select(\n", " options=np.arange(np.shape(data)[1]), # start with a single trial on a single bout... it will update when runs ; old: np.arange(len(trial_times)),\n", " value=0,\n", " #rows=10,\n", " description='Channel used to detect events',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "\n", "slider_threshold = widgets.FloatSlider(\n", " min=-2,\n", " max=2,\n", " value=0.2,\n", " step=0.001,\n", " readout_format='.3f',\n", " continuous_update=False,\n", " readout=True,\n", " description='event detection threshold',\n", " style = {'description_width': '200px'})\n", "slider_threshold.layout.width = '600px'\n", "\n", "detect_type_radio = widgets.RadioButtons(\n", " options=['peak', 'level crossing'],\n", " value='peak', # Defaults to 'level crossing'\n", " layout={'width': 'max-content'}, # If the items' names are long\n", " description='Type of event detection',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "\n", "radio_polarity = widgets.RadioButtons(\n", " options=[1, -1],\n", " value=-1,\n", " description='peaks polarity',\n", " disabled=False,\n", " style = {'description_width': '200px'}\n", ")\n", "\n", "iei_text = widgets.Text(\n", " value='0.005',\n", " placeholder='0.005',\n", " description='min IEI (seconds)',\n", " style = {'description_width': '200px'},\n", " disabled=False\n", ")\n", "\n", "def update_plot(chan_,xrange,thresh_,detect_type,polarity,iei):\n", " fig, ax = plt.subplots(figsize=(10,5),num=1); #specify figure number so that it does not keep creating new ones\n", " \n", " signal = data[:,chan_]\n", " signal = signal-np.median(signal)\n", " \n", " iei = float(iei)\n", " \n", " if iei>0.001:\n", " d = iei*fs #minimum time allowed between distinct events\n", " \n", " if detect_type == 'peak':\n", " r = find_peaks(signal*polarity,height=thresh_,distance=d)\n", " trial_times = r[0]/fs\n", " # print(r[1])\n", " ax.scatter(trial_times,[r[1]['peak_heights']*polarity],marker='^',s=300,color='purple',zorder=3)\n", " \n", " if detect_type == 'level crossing':\n", " # get the changes in bool value for a bool of signal greater than threshold\n", " # if polarity == 1:\n", " threshold_crossings = np.diff(signal*polarity > thresh_, prepend=False)\n", " # get indices where threshold crossings are true\n", " tcross = np.argwhere(threshold_crossings)[:,0]\n", " # get a mask for only positive level crossings\n", " mask_ = [signal[t]-signal[t-1] > 0 for t in tcross]\n", " # if polarity == -1:\n", " # threshold_crossings = np.diff(signal*polarity < thresh_*polarity, prepend=False)\n", " # # get indices where threshold crossings are true\n", " # tcross = np.argwhere(threshold_crossings)[:,0]\n", " # # get a mask for only positive level crossings\n", " # mask_ = [signal[t]-signal[t-1] > 0 for t in tcross]\n", " \n", " # trial times are positive level crossings\n", " trial_times = tcross[mask_]/fs\n", " ax.scatter(trial_times,[thresh_*polarity]*len(trial_times),marker='^',s=300,color='purple',zorder=3)\n", "\n", " starti = int(xrange[0]*fs)+1\n", " stopi = int(xrange[1]*fs)-1\n", " ax.plot(time[starti:stopi], signal[starti:stopi], color='black')\n", " \n", " # ax.plot(tmp,color='black')\n", " ax.hlines(thresh_*polarity, time[starti],time[stopi],linestyle='--',color='green')\n", " \n", " # ax.set_ylim(yrange[0],yrange[1])\n", " ax.set_xlim(xrange[0],xrange[1])\n", " \n", "\n", " ax.xaxis.set_minor_locator(AutoMinorLocator(5))\n", "\n", " \n", " return trial_times\n", "\n", "w_trials_ = interactive(update_plot, chan_=select_channel, \n", " xrange=slider_xrange, \n", " thresh_=slider_threshold, detect_type = detect_type_radio, \n", " polarity=radio_polarity, iei = iei_text);\n", "display(w_trials_)" ] }, { "cell_type": "code", "execution_count": null, "id": "942d64fd-5da8-46ed-9dc5-d3dbdc6a21da", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown Run this cell to finalize the list of stimulus times \n", "#@markdown after settling on a channel and threshold in the interactive plot.
\n", "#@markdown This stores the event times in an array called 'event_times'.
\n", "#@markdown Additionally, the stimulus pulse duration on each trial will be calculated \n", "#@markdown from the raw data so that you can access that information in the interactive trial-by-trial data viewer.\n", "\n", "trial_times, peakthresh_, stimchan = w_trials_.result\n", "\n", "# get stimulus pulse durations for square-wave stimuli\n", "# assumes that each pulse has an onset and offset.... so it finds offset-onset\n", "signal = data[:,stimchan]\n", "# get the changes in bool value for a bool of signal greater than threshold\n", "threshold_crossings = np.diff(signal > peakthresh_, prepend=False)\n", "# get indices where threshold crossings are true\n", "tcross = np.argwhere(threshold_crossings)[:,0]\n", "\n", "stimdur = (np.diff(tcross)*1000)[::2]\n", "stimdur = [np.round(d/fs,2) for d in stimdur]" ] }, { "cell_type": "markdown", "id": "39e97ab3-7300-4821-adb7-1b9416b10a15", "metadata": {}, "source": [ "### Define Bouts" ] }, { "cell_type": "code", "execution_count": null, "id": "55fbe04b-ca90-4dd1-b465-829c49ff7335", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode: \"form\"}\n", "\n", "#@markdown For this experiment, the entire file should be one long bout, \n", "#@markdown but if there were regions that something got messed up or that you want to exclude, you can specify bouts with good data.\n", "#@markdown Specify the list of bout ranges as follows: [[start of bout 0, end of bout 0],[start 1, end 1],...]]
\n", "\n", "bouts_list = [[None,None]] #@param\n", "\n", "#@markdown Then run this code cell to programatically define the list of bouts as 'bouts_list'." ] }, { "cell_type": "markdown", "id": "90d583aa-317f-4543-a0e0-b3bc7a3228d1", "metadata": {}, "source": [ "## 3. Measure the Raw Data time-locked to individual stimuli\n", "\n", "Obtain information from the raw signal time-locked to each stimulus that you need to answer Response prompts 3 and 4." ] }, { "cell_type": "code", "execution_count": null, "id": "d07423b4-8547-4147-b176-45d1fbbd2f72", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode:\"form\"}\n", "\n", "#@markdown Run this code cell to create an interactive plot to \n", "#@markdown examine the raw signal time-locked to individual stimuli (event_times).\n", "#@markdown You can overlay multple channels by selecting more than one.\n", "#@markdown You can overlay multiple stimulus times by selecting more than one. \n", "#@markdown (To select more than one item from an option menu, press the control/command key \n", "#@markdown while mouse clicking or shift while using up/down arrows)\n", "\n", "slider_xrange = widgets.FloatRangeSlider(\n", " min=-0.5,\n", " max=2,\n", " value=(-0.1,1),\n", " step=0.001,\n", " continuous_update=False,\n", " readout=True,\n", " readout_format='.4f',\n", " description='xrange (s)'\n", ")\n", "slider_xrange.layout.width = '600px'\n", "\n", "slider_yrange = widgets.FloatRangeSlider(\n", " min=-5,\n", " max=5, # normal range for earthworm experiments\n", " value=(-2,2),\n", " step=0.01,\n", " continuous_update=False,\n", " readout=True,\n", " description='yrange'\n", ")\n", "slider_yrange.layout.width = '600px'\n", "\n", "ui_range = widgets.VBox([slider_xrange, slider_yrange])\n", "\n", "# trials in bout 0 to start...\n", "trials_t = trial_times[(trial_times>bouts_list[0][0]) & (trial_timesbouts_list[bout_][0]) & (trial_times0) & ((int(fs*t_)+win_1)) 0:\n", " iti = 1000*(trials_t[trial_list[0]] - trials_t[trial_list[0]-1])\n", " trial_readout.value=f'time since last event is (ms): {iti:.2f}'\n", " trial_abs_readout.value=f'time of this event is (s): {trials_t[trial_list[0]]}'\n", " \n", " for lt_ in lagging_peak_times:\n", " ax.scatter(lagging_peak_times,lagging_peak_amp,color='black',s=50,zorder=3)\n", " \n", "\n", " ax.set_ylim(yrange[0],yrange[1]);\n", " ax.set_xlabel('seconds')\n", " # ax.vlines(0,yrange[0],yrange[1],color='black')\n", "\n", " \n", "# # Change major ticks to show every 20.\n", " # ax_pwm.xaxis.set_major_locator(MultipleLocator(5))\n", " # ax_pwm.yaxis.set_major_locator(MultipleLocator(5))\n", "\n", " # # Change minor ticks to show every 5. (20/4 = 5)\n", " # ax_mro.yaxis.set_minor_locator(AutoMinorLocator(10))\n", " ax.xaxis.set_minor_locator(AutoMinorLocator(10))\n", " # ax_pwm.yaxis.set_minor_locator(AutoMinorLocator(5))\n", "\n", "# # Turn grid on for both major and minor ticks and style minor slightly\n", "# # # differently.\n", " ax.grid(which='major', color='gray', linestyle='-')\n", " ax.grid(which='minor', color='gray', linestyle=':')\n", "# ax_pwm.grid(which='major', color='gray', linestyle='-')\n", "# ax_pwm.grid(which='minor', color='gray', linestyle=':')\n", "\n", "\n", "# display(w)\n", "w = interactive_output(update_plot, {'chan_list':select_channels,\n", " 'trial_list':select_trials, \n", " 'bout_':select_bouts,\n", " 'yrange':slider_yrange, \n", " 'xrange':slider_xrange,\n", " 'lagging_chan_':detect_chan_radio,\n", " 'thresh_':slider_threshold,\n", " 'polarity':radio_polarity,\n", " 'noise_amp':noise_amp_text});\n", "display(trial_abs_readout,stimdur_readout,trial_readout,lagging_time_readout,lagging_amp_readout,\n", " ui_trials,ui_peaks,noise_amp_text,w,ui_range)\n" ] }, { "cell_type": "markdown", "id": "948f937f-720a-4c50-8fff-f3dc72ef8f01", "metadata": {}, "source": [ "## 4. Plot processed data\n", "\n", "Make a \"scatter plot\" of processed data from a csv file. \n" ] }, { "cell_type": "code", "execution_count": null, "id": "cdc62101-c92d-4fb8-99bc-b1d894cd2e5a", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#@title {display-mode:\"form\"}\n", "\n", "#@markdown Specify the filepath to a csv file\n", "filepath = 'filepath to csv' #@param\n", "# filepath = '/Users/kperks/Documents/Teaching/Neurophysiology-Lab/modules/crayfish-erg/duration_.csv'\n", "\n", "#@markdown Specify the single header name of the column you want for the x-axis data. \n", "x_column = ['column header for x axis'] #@param\n", "# x_column = ['peakamp']\n", "\n", "# #@markdown Specify categorical bins using np.arange(start,stop,step) if the x_column is a continuous variable. \n", "# #@markdown Use empty brackets if not needed.\n", "# categorical_bins = np.arange(0,25,0.5) #@param\n", "\n", "#@markdown Specify the header name(s) of the column(s)s you want for your y points. \n", "#@markdown If more than one header is specified (separated by commas), each will be plotted overlaid in a different color for a scatter plot\n", "y_column = ['column header(s) for y axis 1'] #@param\n", "# y_column = ['steady state','peak']\n", "# y_column = ['m1','m2']\n", "# y_column = ['peakt']\n", "\n", "# #@markdown Specify the plot type ('scatter' or 'violin'). Note that for a 'violin' plot, only the 'x_column' data would be used.\n", "# plot_type = 'plot type' #@param\n", "\n", "# @markdown Specify the plot type ('scatter' or 'point'). \n", "# @markdown Note that for a 'point' plot, the x-axis is categorical and the y-axis data has error bars.\n", "# plot_type = 'point'\n", "plot_type = 'scatter' \n", "\n", "df = pd.read_csv(filepath)\n", "\n", "hfig,ax = plt.subplots(figsize=(8,5))\n", "\n", "if plot_type == 'scatter':\n", " df_melted = df[y_column+x_column].melt(x_column[0],var_name='headers')\n", " g = sns.scatterplot(data=df_melted,x=x_column[0],y='value',hue='headers',s=100,alpha=0.75, ax = ax);\n", " \n", " \n", "if plot_type == 'point':\n", " df_melted = df[y_column+x_column].melt(x_column[0],var_name='headers')\n", " \n", "# if len(categorical_bins)>0:\n", "# df_melted[x_column[0]] = pd.cut(df_melted[x_column[0]],bins=categorical_bins,labels=categorical_bins[1:])\n", " \n", " g = sns.pointplot(data=df_melted,x=x_column[0],y='value',hue='headers',alpha=0.75);\n", "\n", " \n", "g.legend(loc='center left', bbox_to_anchor=(1.1, 0.5), ncol=1);\n", "# ax.set_xlim(x_lim[0],x_lim[1])\n" ] }, { "cell_type": "markdown", "id": "d451b557-6777-4492-955e-a41130ad39b9", "metadata": {}, "source": [ "
\n", "Written by Dr. Krista Perks for courses taught at Wesleyan University." ] }, { "cell_type": "code", "execution_count": null, "id": "5c707588-d5f6-4142-8289-da7ba6a97dbc", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "57a8a68c-055e-4af0-ba4d-67d0d67148f1", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "a350712e-e145-4475-9588-e3f19469b5ce", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "86fd2b4a-b890-465f-bf33-c2057738789e", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "b05de192-2baf-4d46-8848-6db1e3b50bff", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "e6519277-fc86-44fd-a75b-f74cad0efcf1", "metadata": {}, "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" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }