Duncan - Data Explorer
Contents
Duncan - Data Explorer¶
Superficial Flexor: Postsynaptic Activity¶
Setup¶
Import and define functions
#@title {display-mode: "form" }
#@markdown Run this code cell to import packages and define functions
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy import ndimage
from scipy.optimize import curve_fit
from scipy.signal import hilbert,medfilt,resample, find_peaks, unit_impulse
import seaborn as sns
from datetime import datetime,timezone,timedelta
pal = sns.color_palette(n_colors=15)
pal = pal.as_hex()
import matplotlib.pyplot as plt
import random
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from pathlib import Path
from matplotlib.ticker import (AutoMinorLocator, MultipleLocator)
from ipywidgets import widgets, interact, interactive, interactive_output
%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")
print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
Mount Google Drive
#@title {display-mode: "form" }
#@markdown Run this cell to mount your Google Drive.
from google.colab import drive
drive.mount('/content/drive')
print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
Import data¶
Import data digitized with Nidaq USB6211 and recorded using Bonsai-rx as a .bin file
If you would like sample this Data Explorer, but do not have data, you can download the following examples (two channels digitized at 40000). Channel 0 is the signal measured from N3 and Channel 1 is the signal measured from the Superficial Flexor muscle.
#@title {display-mode: "form" }
#@markdown Specify the file path
#@markdown to your recorded data on Drive (find the filepath in the colab file manager:
filepath = "type/copy filepath here" #@param
# filepath = '/Users/kperks/Downloads/Duncan_final/Capstone_McKee0.bin'
#@markdown Specify the sampling rate and number of channels recorded.
sampling_rate = 40000 #@param
number_channels = 2 #@param
# downsample = False #@param
# newfs = 10000 #@param
#@markdown After you have filled out all form fields,
#@markdown run this code cell to load the data.
filepath = Path(filepath)
# No need to edit below this line
#################################
data = np.fromfile(Path(filepath), dtype = np.float64)
data = data.reshape(-1,number_channels)
data_dur = np.shape(data)[0]/sampling_rate
print('duration of recording was %0.2f seconds' %data_dur)
fs = sampling_rate
# if downsample:
## newfs = 10000 #downsample emg data
# chunksize = int(sampling_rate/newfs)
# data = data[0::chunksize,:]
# fs = int(np.shape(data)[0]/data_dur)
time = np.linspace(0,data_dur,np.shape(data)[0])
print('Data upload completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
Detect PSP peaks¶
For a more extensive RAW Data Explorer than the one provided in the above figure, use the DataExplorer.py application found in the howto section of the course website.
#@title {display-mode: "form" }
#@markdown Run this cell to create an interactive plot with a slider to scroll
#@markdown through the signal
#@markdown and set an appropriate event detection threshold and noise threshold.
slider_xrange = widgets.FloatRangeSlider(
min=0,
max=data_dur,
value=(0,data_dur),
step= 0.1,
readout=True,
continuous_update=False,
description='Time Range (s)',
style = {'description_width': '200px'})
slider_xrange.layout.width = '600px'
slider_yrange = widgets.FloatRangeSlider(
min=-2,
max=2,
value=[-0.9,0.1],
step=0.01,
continuous_update=False,
readout=True,
description='yrange',
style = {'description_width': '200px'})
slider_yrange.layout.width = '600px'
slider_text = widgets.Text(
value='0.04',
placeholder='0.04',
description='time to peak (s)',
style = {'description_width': '200px'},
disabled=False
)
select_channel = widgets.Select(
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)),
value=0,
#rows=10,
description='Channel',
style = {'description_width': '50px'},
disabled=False
)
select_channel.layout.width = '100px'
select_channel.layout.height = '50px'
# slider_threshold = widgets.FloatSlider(
# min=0,
# max=5,
# value=0.2,
# step=0.005,
# readout_format='.3f',
# continuous_update=False,
# readout=True,
# description='event detection threshold',
# style = {'description_width': '200px'})
# slider_threshold.layout.width = '600px'
radio_polarity = widgets.RadioButtons(
options=[1, -1],
value=1,
description='peaks polarity',
layout={'width': 'max-content'},
style = {'description_width': '100px'},
disabled=False
)
radio_polarity.layout.width = '300px'
# ui_peaks = widgets.HBox([select_channel, radio_polarity, slider_threshold])
# detect_type_radio = widgets.RadioButtons(
# options=['peak', 'level crossing positive', 'level crossing negative'],
# value='level crossing positive', # Defaults to 'level crossing'
# layout={'width': 'max-content'}, # If the items' names are long
# description='Type of event detection',
# style = {'description_width': '200px'},
# disabled=False
# )
iei_text = widgets.Text(
value='0.01',
placeholder='0.01',
description='minimum interval (s)',
style = {'description_width': '200px'},
disabled=False
)
noise_amp_text = widgets.Text(
value='0.06',
placeholder='0.06',
description='threshold amplitude (V)',
style = {'description_width': '200px'},
disabled=False
)
def update_plot(chan_,xrange,yrange,polarity,iei,noise_amp,onset):
fig, ax = plt.subplots(figsize=(10,5),num=1); #specify figure number so that it does not keep creating new ones
signal = data[:,chan_].reshape(-1)
# signal = signal-np.median(signal)
iei = float(iei)
noise_amp = float(noise_amp)
onset = float(onset)
prepeak_base = int(onset*fs)
if iei>1/fs:
d = iei*fs #minimum time allowed between distinct events
# if thresh_ >=0:
# r = find_peaks(signal,height=thresh_,distance=d)
# if thresh_ <0:
# r = find_peaks(-1*signal,height=-1*thresh_,distance=d)
# trial_times = r[0]/fs
if polarity == 1:
# r_thresh = find_peaks(signal,height=thresh_,distance=d)
r_prom = find_peaks(signal,distance=d,prominence=noise_amp)
r_ind = r_prom[0] #np.intersect1d(r_thresh[0],r_prom[0])
# lagging_peak_amp = r[1]['peak_heights']
peak_amp = signal[r_ind]
base_amp = signal[r_ind-prepeak_base]
# ax.hlines(thresh_, xrange[0],xrange[1],linestyle='--',color='green')
if polarity == -1:
# r = find_peaks(-1*lagging_signal,height=-1*thresh_,distance=d,prominence=float(noise_amp))
# r_thresh = find_peaks(-1*signal,height=thresh_,distance=d)
r_prom = find_peaks(-1*signal,distance=d,prominence=noise_amp)
r_ind = r_prom[0] #np.intersect1d(r_thresh[0],r_prom[0])
# lagging_peak_amp = -1*r[1]['peak_heights']
peak_amp = signal[r_ind]
base_amp = signal[r_ind-prepeak_base]
# ax.hlines(-thresh_, xrange[0],xrange[1],linestyle='--',color='green')
# peak_times = np.asarray([np.round(s/fs,2) for s in r_ind])
peak_times = np.asarray([s/fs for s in r_ind])
starti = int(xrange[0]*fs)+1
stopi = int(xrange[1]*fs)-1
ax.plot(time[starti:stopi], signal[starti:stopi], color='black')
inwin_inds = (peak_times>(xrange[0])) & (peak_times<(xrange[1]))
ax.scatter(peak_times[inwin_inds],peak_amp[inwin_inds], zorder=3,color='red',s=20)
ax.scatter(peak_times[inwin_inds]-onset,base_amp[inwin_inds], zorder=3,color='green',s=20)
ax.set_ylim(yrange[0],yrange[1])
ax.set_xlim(xrange[0],xrange[1])
ax.xaxis.set_minor_locator(AutoMinorLocator(5))
return peak_times, peak_amp-base_amp
# w_trials_ = interactive(update_plot, chan_=select_channel,
# xrange=slider_xrange,
# thresh_=slider_threshold, iei = iei_text, noise_amp = noise_amp_text);
# display(w_trials_)
w_peak_detect_all = interactive(update_plot, chan_=select_channel,
yrange=slider_yrange,xrange=slider_xrange,
polarity = radio_polarity, iei = iei_text,noise_amp=noise_amp_text,
onset=slider_text);
display(w_peak_detect_all)
# w_peak_detect_all = interactive(update_plot, {'chan_':select_channel,
# 'yrange':slider_yrange,
# 'xrange':slider_xrange,
# 'thresh_':slider_threshold,
# 'polarity':radio_polarity,
# 'iei':iei_text,
# 'noise_amp':noise_amp_text});
# display(ui_peaks,widgets.HBox([iei_text,noise_amp_text]),w_peak_detect_all,ui_range)
Define Bouts¶
To efficiently assess your data with this analysis, make sure to exclude any raw data that does not have a clean (low-noise) signal. For the simultaneously recorded pre- and post-synaptic signals, make sure to exclude raw data in which the post-synaptic electrode was not stably in the cell. The more data you are able to include, the better your spike sorting results will be.
#@title {display-mode: "form"}
#@markdown For this experiment, the entire file should be one long bout,
#@markdown but if there were regions that something got messed up or that you want to exclude, you can specify bouts with good data.
#@markdown Specify the list of bout ranges as follows: [[start of bout 0, end of bout 0],[start 1, end 1],...]] <br>
bouts_list = [[None,None]] #@param
#@markdown Then run this code cell to programatically define the list of bouts as 'bouts_list'.
Process data and save as .csv¶
Specify the “condition” for the data you just processed. The condition should be a string (ie. surrounded by quotes).
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.
#@title {display-mode: "form"}
condition = 'baseline' #@param
peak_times,peak_amp = w_peak_detect_all.result
mask = []
for b_ in bouts_list:
mask.append((peak_times>b_[0]) & (peak_times<b_[1]))
mask = sum(mask).astype(bool)
df = pd.DataFrame({'time':peak_times[mask],'amp':peak_amp[mask],'number':np.arange(0,sum(mask),1),'condition':condition})
df.to_csv('PSP_' + condition + '.csv')
Visualize the average raw signal time-locked to PSP peaks.¶
Average across unitary PSPs that are temporally isolated from each other (not temporally summated with another). Average is centered at the PSP peak. Only PSPs within the previously specified bouts are included.
#@title {display-mode:"form"}
#@markdown Specify which channel has the intracellular signal.
#@markdown Then run this cell to plot the average spike-triggered
#@markdown post-synaptic potential for each spike cluster you defined in Part I
intracellular_channel = 1 #@param
psp_duration = 0.04 #@param
slider_xrange = widgets.FloatRangeSlider(
min=-0.05,
max=0.10,
value=(-0.02,0.1),
step=0.005,
readout_format='.3f',
continuous_update=False,
readout=True,
description='xrange (s)'
)
slider_xrange.layout.width = '600px'
def update_plot(xrange):
# No need to edit below this line
#################################
windur = xrange[1]-xrange[0]
winsamps = int(windur * fs)
onset = int(xrange[0]*fs)
offset = int(xrange[1]*fs)
x = np.linspace(xrange[0],xrange[1],offset-onset)
hfig,ax = plt.subplots(figsize=(8,4))
ax.set_ylabel('volts recorded',fontsize=14)
ax.set_xlabel('seconds',fontsize=14)
# plt.xticks(fontsize=14)
# plt.yticks(fontsize=14)
trial_t = peak_times[mask]
trial_t_sorted = trial_t[((np.diff(trial_t,prepend=np.NaN)>psp_duration)&(np.diff(trial_t,append=np.NaN)>psp_duration))]
synwav = np.asarray([data[int(t*fs)+onset:int(t*fs)+offset,intracellular_channel] for t in trial_t_sorted
if (((int(t*fs)+onset)>0) & (int(t*fs)+offset<np.shape(data)[0]))])
wav_u = np.mean(synwav,0)
wav_std = np.std(synwav,0)
ax.plot(x,wav_u,linewidth = 3,color = 'black');
ax.fill_between(x, wav_u-wav_std, wav_u+wav_std, alpha = 0.25, color = 'purple')
# plt.legend(bbox_to_anchor=[1.5,1], fontsize=14);
ax.set_xlim(xrange[0],xrange[1])
w_psps_sorted_ = interactive(update_plot, xrange=slider_xrange);
display(w_psps_sorted_)
Compare across conditions¶
Provide a list of the csv files generated in the “Save processed data as .csv” 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.
#@title {display-mode: "form"}
#@markdown List the full path to each csv file created using the "adaptation analysis tool."
file_list = ['PSP_baseline.csv', 'PSP_dopamine.csv'] #@param
# file_list = ['adaptation_2_b_.csv','adaptation_5_b_.csv','adaptation_10_b_.csv']
#@markdown Then run this code cell to create a dataframe combining the data from all listed files.
df = pd.DataFrame({})
for f in file_list:
df_ = pd.read_csv(f)
df = pd.concat([df,df_])
A box plot (or box-and-whisker plot) shows the distribution of quantitative data in a way that facilitates comparisons between variables or across levels of a categorical variable. The box shows the quartiles of the dataset while the whiskers extend to show the rest of the distribution, except for points that are determined to be “outliers” using a method that is a function of the inter-quartile range.
#@title {display-mode: "form"}
#@markdown Run this code cell to create a bo xplot of condition by amplitude
x_param = 'condition' #@param
y_param = 'amp' #@param
hfig,ax = plt.subplots(figsize=(5,4))
sns.boxplot(data=df, x=x_param, y=y_param);
A violin plot plays a similar role as a box and whisker plot. It shows the distribution of quantitative data across several levels of one (or more) categorical variables such that those distributions can be compared. Unlike a box plot, in which all of the plot components correspond to actual datapoints, the violin plot features a kernel density estimation of the underlying distribution.
#@title {display-mode: "form"}
#@markdown Run this code cell to create a violin plot of condition by amplitude
x_param = 'condition' #@param
y_param = 'amp' #@param
hfig,ax = plt.subplots(figsize=(5,4))
sns.violinplot(data=df, x=x_param, y=y_param, bw=.15, scale="count");
Written by Dr. Krista Perks for courses taught at Wesleyan University.