diff --git a/actions/stimulus-spike-lfp-response-other-drive/attributes.yaml b/actions/stimulus-spike-lfp-response-other-drive/attributes.yaml new file mode 100644 index 000000000..d9f8669f7 --- /dev/null +++ b/actions/stimulus-spike-lfp-response-other-drive/attributes.yaml @@ -0,0 +1,5 @@ +registered: '2020-12-18T15:14:54' +data: + results: results.csv + notebook: 20_stimulus-spike-lfp-response.ipynb + html: 20_stimulus-spike-lfp-response.html diff --git a/actions/stimulus-spike-lfp-response-other-drive/data/10-calculate-stimulus-spike-lfp-response-other-drive.html b/actions/stimulus-spike-lfp-response-other-drive/data/10-calculate-stimulus-spike-lfp-response-other-drive.html new file mode 100644 index 000000000..e027a21b9 --- /dev/null +++ b/actions/stimulus-spike-lfp-response-other-drive/data/10-calculate-stimulus-spike-lfp-response-other-drive.html @@ -0,0 +1,14186 @@ + + +
+ +%load_ext autoreload
+%autoreload 2
+
import matplotlib.pyplot as plt
+%matplotlib inline
+import spatial_maps as sp
+import septum_mec.analysis.data_processing as dp
+import septum_mec.analysis.registration
+import expipe
+import os
+import pathlib
+import scipy
+import scipy.signal
+import numpy as np
+import exdir
+import pandas as pd
+import optogenetics as og
+import quantities as pq
+import shutil
+from distutils.dir_util import copy_tree
+import elephant as el
+import neo
+from scipy.signal import find_peaks
+from scipy.interpolate import interp1d
+from matplotlib import mlab
+
+from tqdm import tqdm_notebook as tqdm
+from tqdm._tqdm_notebook import tqdm_notebook
+tqdm_notebook.pandas()
+
data_loader = dp.Data()
+actions = data_loader.actions
+project = data_loader.project
+
output = pathlib.Path('output/stimulus-spike-lfp-response-other-drive')
+(output / 'data').mkdir(parents=True, exist_ok=True)
+
identify_neurons = actions['identify-neurons']
+# sessions = pd.read_csv(identify_neurons.data_path('sessions'))
+units = pd.read_csv(identify_neurons.data_path('units'))
+
def get_lim(action_id):
+ stim_times = data_loader.stim_times(action_id)
+ if stim_times is None:
+ return [0, np.inf]
+ stim_times = np.array(stim_times)
+ return [stim_times.min(), stim_times.max()]
+
+def get_mask(lfp, lim):
+ return (lfp.times >= lim[0]) & (lfp.times <= lim[1])
+
+def zscore(a):
+ return (a - a.mean(0)) / a.std(0)
+
+def compute_stim_freq(action_id):
+ stim_times = data_loader.stim_times(action_id)
+ if stim_times is None:
+ return np.nan
+ stim_times = np.array(stim_times)
+ return 1 / np.mean(np.diff(stim_times))
+
def signaltonoise(a, axis=0, ddof=0):
+ a = np.asanyarray(a)
+ m = a.mean(axis)
+ sd = a.std(axis=axis, ddof=ddof)
+ return np.where(sd == 0, 0, m / sd)
+
def compute_energy(p, f, f1, f2):
+ if np.isnan(f1) or np.all(np.isnan(p)):
+ return np.nan
+ mask = (f > f1) & (f < f2)
+ df = f[1] - f[0]
+ return np.sum(p[mask]) * df
+
def find_theta_peak(p, f, f1, f2):
+ if np.all(np.isnan(p)):
+ return np.nan, np.nan
+ mask = (f > f1) & (f < f2)
+ p_m = p[mask]
+ f_m = f[mask]
+ peaks, _ = find_peaks(p_m)
+ idx = np.argmax(p_m[peaks])
+ return f_m[peaks[idx]], p_m[peaks[idx]]
+
def compute_half_width(p, f, m_p, m_f):
+ if np.isnan(m_p):
+ return np.nan, np.nan
+ m_p_half = m_p / 2
+ half_p = p - m_p_half
+ idx_f = np.where(f <= m_f)[0].max()
+ idxs_p1, = np.where(np.diff(half_p[:idx_f + 1] > 0) == 1)
+ if len(idxs_p1) == 0:
+ return np.nan, np.nan
+ m1 = idxs_p1.max()
+ idxs_p2, = np.where(np.diff(half_p[idx_f:] > 0) == 1)
+ m2 = idxs_p2.min() + idx_f
+ assert p[m1] < m_p_half < p[m1+1], (p[m1], m_p_half, p[m1+1])
+ assert p[m2] > m_p_half > p[m2+1], (p[m2], m_p_half, p[m2+1])
+
+ f1 = interp1d([half_p[m1], half_p[m1 + 1]], [f[m1], f[m1 + 1]])(0)
+ f2 = interp1d([half_p[m2], half_p[m2 + 1]], [f[m2], f[m2 + 1]])(0)
+ return f1, f2
+
# p = np.load('debug_p.npy')
+# f = np.load('debug_f.npy')
+# compute_half_width(p, f, 0.01038941, 30.30187709636872)
+
# plt.plot(f, p)
+# plt.xlim(29.9,30.6)
+
def compute_stim_peak(p, f, s_f):
+ if np.isnan(s_f):
+ return np.nan
+ return interp1d(f, p)(s_f)
+
def compute_spike_lfp_coherence(anas, sptr, NFFT):
+
+ sigs, freqs = el.sta.spike_field_coherence(anas, sptr, **{'nperseg': NFFT})
+ return sigs, freqs
+
def butter_bandpass(lowcut, highcut, fs, order=5):
+ nyq = 0.5 * fs
+ low = lowcut / nyq
+ high = highcut / nyq
+ b, a = scipy.signal.butter(order, [low, high], btype='band')
+ return b, a
+
+
+def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
+ b, a = butter_bandpass(lowcut, highcut, fs, order=order)
+ y = scipy.signal.filtfilt(b, a, data)
+ return y
+
+# def compute_spike_phase_func(lfp, times, return_degrees=False):
+# x_a = hilbert(lfp)
+# x_phase = np.angle(x_a)
+# if return_degrees:
+# x_phase = x_phase * 180 / np.pi
+# return interp1d(times, x_phase)
+
+
+def vonmises_kde(data, kappa=100, n_bins=100):
+ from scipy.special import i0
+ bins = np.linspace(-np.pi, np.pi, n_bins)
+ x = np.linspace(-np.pi, np.pi, n_bins)
+ # integrate vonmises kernels
+ kde = np.exp(kappa * np.cos(x[:, None] - data[None, :])).sum(1) / (2 * np.pi * i0(kappa))
+ kde /= np.trapz(kde, x=bins)
+ return bins, kde
+
+
+def spike_phase_score(phase_bins, density):
+ import math
+ import pycircstat as pc
+ ang = pc.mean(phase_bins, w=density)
+ vec_len = pc.resultant_vector_length(phase_bins, w=density)
+ # ci_lim = pc.mean_ci_limits(head_angle_bins, w=rate)
+ return ang, vec_len
+
def compute_clean_lfp(anas, width=500, threshold=2):
+ anas = np.array(anas)
+ idxs, = np.where(abs(anas) > threshold)
+ for idx in idxs:
+ anas[idx-width:idx+width] = 0 # TODO AR model prediction
+ return anas, idxs
+
+
+def compute_clean_spikes(spikes, idxs, times, width=500):
+
+ for idx in idxs:
+ t0 = times[idx-width]
+ stop = idx + width
+ if stop > len(times) - 1:
+ stop = len(times) - 1
+ t1 = times[stop]
+ mask = (spikes > t0) & (spikes < t1)
+ spikes = spikes[~mask]
+ spikes = spikes[spikes <= times[-1]]
+ return spikes
+
+
+def prepare_spike_lfp(anas, sptr, t_start, t_stop):
+
+ t_start = t_start * pq.s if t_start is not None else 0 * pq.s
+ sampling_rate = anas.sampling_rate
+ units = anas.units
+ times = anas.times
+ if t_start is not None and t_stop is not None:
+ t_stop = t_stop * pq.s
+ mask = (times > t_start) & (times < t_stop)
+ anas = np.array(anas)[mask,:]
+ times = times[mask]
+
+ # take best channel from other drive
+ best_channel = np.argmax(signaltonoise(anas))
+# best_channel = np.random.choice(anas.shape[1])
+
+ cleaned_anas, idxs = compute_clean_lfp(anas[:, best_channel])
+
+ cleaned_anas = neo.AnalogSignal(
+ signal=cleaned_anas * units, sampling_rate=sampling_rate, t_start=t_start
+ )
+
+ spike_units = sptr.units
+ spike_times = sptr.times
+ spike_times = compute_clean_spikes(spike_times, idxs, times)
+
+ sptr = neo.SpikeTrain(
+ spike_times[(spike_times > t_start) & (spike_times < times[-1])], units=spike_units,
+ t_start=t_start, t_stop=times[-1]
+ )
+
+ return cleaned_anas, sptr, best_channel
+
def compute_spike_phase_func(lfp, times, return_degrees=False):
+ from scipy.fftpack import next_fast_len
+ x_a = scipy.signal.hilbert(
+ lfp, next_fast_len(len(lfp)))[:len(lfp)]
+# x_a = hilbert(lfp)
+ x_phase = np.angle(x_a, deg=return_degrees)
+ return interp1d(times, x_phase)
+
+
+def compute_spike_phase(lfp, spikes, flim=[6,10]):
+
+ sample_rate = lfp.sampling_rate.magnitude
+
+ # sometimes the position is recorded after LFP recording is ended
+ times = np.arange(lfp.shape[0]) / sample_rate
+
+ spikes = np.array(spikes)
+ spikes = spikes[(spikes > times.min()) & (spikes < times.max())]
+
+ filtered_lfp = butter_bandpass_filter(
+ lfp.magnitude.ravel(), *flim, fs=sample_rate, order=3)
+
+ spike_phase_func = compute_spike_phase_func(filtered_lfp, times)
+
+ return spike_phase_func(spikes), filtered_lfp
+
plt.figure(figsize=(16,9))
+lfp = data_loader.lfp('1833-200619-2', 6)
+# lfp = data_loader.lfp('1834-220319-3', 6)
+# lfp = data_loader.lfp('1849-010319-4', 6)
+times = np.arange(lfp.shape[0]) / lfp.sampling_rate.magnitude
+clean_lfp, _ = compute_clean_lfp(lfp.magnitude[:, 0], threshold=2)
+plt.plot(times,lfp[:,0])
+plt.plot(times,clean_lfp)
+
+plt.figure(figsize=(16,9))
+plt.psd(lfp[:,0].ravel(), Fs=1000, NFFT=10000)
+plt.psd(clean_lfp, Fs=1000, NFFT=10000)
+plt.xlim(0,100)
+
# plt.figure(figsize=(16,9))
+
+# plt.plot(times,lfp[:,0])
+# # plt.plot(clean_lfp*100)
+# plt.plot(times[:-1], np.diff(lfp[:,0].magnitude.ravel()))
+# plt.xlim(64.5,65.5)
+# # plt.ylim(-250,250)
+
drive_0_channel_groups = [0, 1, 2, 3]
+drive_1_channel_groups = [4, 5, 6, 7]
+
# action_id_0, channel_0, unit_0 = '1833-200619-1', 6, 163
+# action_id_1, channel_1, unit_1 = '1833-200619-2', 6, 28
+action_id_0, channel_0, unit_0 = '1834-220319-3', 2, 46
+action_id_1, channel_1, unit_1 = '1834-220319-4', 2, 60
+
+# change data loader to get all LFPs and then selecte the best form the other
+lfp_0 = data_loader.lfp(action_id_0, channel_0)
+lfp_1 = data_loader.lfp(action_id_1, channel_1)
+
+# select best channel among other drive
+if channel_0 in drive_0_channel_groups:
+ lfps = []
+ for ch in drive_1_channel_groups:
+ lfps.append(data_loader.lfp(action_id_0, ch))
+elif channel_0 in drive_1_channel_groups:
+ lfps = []
+ for ch in drive_0_channel_groups:
+ lfps.append(data_loader.lfp(action_id_0, ch))
+
+# merge lfp of other drive into a single AnalogSignal
+lfp_arrays = np.hstack(lfps).as_array()
+lfp_0 = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,
+ t_start=lfps[0].t_start)
+
+# select best channel among other drive
+if channel_1 in drive_0_channel_groups:
+ lfps = []
+ for ch in drive_1_channel_groups:
+ lfps.append(data_loader.lfp(action_id_1, ch))
+elif channel_1 in drive_1_channel_groups:
+ lfps = []
+ for ch in drive_0_channel_groups:
+ lfps.append(data_loader.lfp(action_id_1, ch))
+# merge lfp of other tetrodes into a signle AnalogSignal
+lfp_arrays = np.hstack(lfps).as_array()
+lfp_1 = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,
+ t_start=lfps[0].t_start)
+
+
+sample_rate_0 = lfp_0.sampling_rate
+sample_rate_1 = lfp_1.sampling_rate
+
+lim_0 = get_lim(action_id_0)
+lim_1 = get_lim(action_id_1)
+
+sptrs_0 = data_loader.spike_trains(action_id_0, channel_0)
+
+sptrs_1 = data_loader.spike_trains(action_id_1, channel_1)
+
+cleaned_lfps_0, sptr_0, best_channel_0 = prepare_spike_lfp(lfp_0, sptrs_0[unit_0], *lim_0)
+
+cleaned_lfps_1, sptr_1, best_channel_1 = prepare_spike_lfp(lfp_1, sptrs_1[unit_1], *lim_1)
+
+coher_0, freq_0 = compute_spike_lfp_coherence(cleaned_lfps_0, sptr_0, 4096)
+
+coher_1, freq_1 = compute_spike_lfp_coherence(cleaned_lfps_1, sptr_1, 4096)
+
+spike_phase_0, filtered_lfp_0 = compute_spike_phase(cleaned_lfps_0, sptrs_0[unit_0], flim=[6,10])
+
+spike_phase_1, filtered_lfp_1 = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[6,10])
+
+# spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[10.5,11.5])
+spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[29.5,30.5])
+
+plt.figure()
+plt.plot(freq_0, coher_0.ravel())
+plt.plot(freq_1, coher_1.ravel())
+plt.xlim(0,20)
+
+plt.figure()
+bins_0, kde_0 = vonmises_kde(spike_phase_0, 100)
+ang_0, vec_len_0 = spike_phase_score(bins_0, kde_0)
+plt.polar(bins_0, kde_0, color='b')
+plt.polar([ang_0, ang_0], [0, vec_len_0], color='b')
+
+bins_1, kde_1 = vonmises_kde(spike_phase_1, 100)
+ang_1, vec_len_1 = spike_phase_score(bins_1, kde_1)
+plt.polar(bins_1, kde_1, color='r')
+plt.polar([ang_1, ang_1], [0, vec_len_1], color='r')
+
+bins_1_stim, kde_1_stim = vonmises_kde(spike_phase_1_stim, 100)
+ang_1_stim, vec_len_1_stim = spike_phase_score(bins_1_stim, kde_1_stim)
+plt.polar(bins_1_stim, kde_1_stim, color='k')
+plt.polar([ang_1_stim, ang_1_stim], [0, vec_len_1_stim], color='k')
+
+
# TODO fix artefact stuff from phase precession
+
NFFT = 8192
+theta_band_f1, theta_band_f2 = 6, 10
+drive_0_channel_groups = [0, 1, 2, 3]
+drive_1_channel_groups = [4, 5, 6, 7]
+
coherence_data, freqency_data = {}, {}
+theta_kde_data, theta_bins_data = {}, {}
+stim_kde_data, stim_bins_data = {}, {}
+
+def process(row):
+ action_id = row['action']
+ channel_group = row['channel_group']
+ unit_name = row['unit_name']
+
+ name = f'{action_id}_{channel_group}_{unit_name}'
+
+ # select best channel among other drive
+ if channel_group in drive_0_channel_groups:
+ lfps = []
+ for ch in drive_1_channel_groups:
+ lfps.append(data_loader.lfp(action_id, ch))
+ elif channel_group in drive_1_channel_groups:
+ lfps = []
+ for ch in drive_0_channel_groups:
+ lfps.append(data_loader.lfp(action_id, ch))
+
+ # merge lfp of other tetrodes into a signle AnalogSignal
+ lfp_arrays = np.hstack(lfps).as_array()
+ lfp = neo.AnalogSignal(signal=lfp_arrays * lfps[0].units, sampling_rate=lfps[0].sampling_rate,
+ t_start=lfps[0].t_start)
+ sptr = data_loader.spike_train(action_id, channel_group, unit_name)
+
+ lim = get_lim(action_id)
+
+ cleaned_lfp, sptr, best_channel = prepare_spike_lfp(lfp, sptr, *lim)
+
+ p_xys, freq = compute_spike_lfp_coherence(cleaned_lfp, sptr, NFFT=NFFT)
+
+ p_xy = p_xys.magnitude.ravel()
+ freq = freq.magnitude
+
+ theta_f, theta_p_max = find_theta_peak(p_xy, freq, theta_band_f1, theta_band_f2)
+
+ theta_energy = compute_energy(p_xy, freq, theta_band_f1, theta_band_f2) # theta band 6 - 10 Hz
+
+ theta_half_f1, theta_half_f2 = compute_half_width(p_xy, freq, theta_p_max, theta_f)
+
+ theta_half_width = theta_half_f2 - theta_half_f1
+
+ theta_half_energy = compute_energy(p_xy, freq, theta_half_f1, theta_half_f2) # theta band 6 - 10 Hz
+
+ theta_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[theta_band_f1, theta_band_f2])
+ theta_bins, theta_kde = vonmises_kde(theta_spike_phase)
+ theta_ang, theta_vec_len = spike_phase_score(theta_bins, theta_kde)
+ theta_kde_data.update({name: theta_kde})
+ theta_bins_data.update({name: theta_bins})
+
+ # stim
+
+ stim_freq = compute_stim_freq(action_id)
+
+ stim_p_max = compute_stim_peak(p_xy, freq, stim_freq)
+
+ stim_half_f1, stim_half_f2 = compute_half_width(p_xy, freq, stim_p_max, stim_freq)
+ stim_half_width = stim_half_f2 - stim_half_f1
+
+ stim_energy = compute_energy(p_xy, freq, stim_half_f1, stim_half_f2)
+
+ if np.isnan(stim_freq):
+ stim_spike_phase, stim_bins, stim_kde, stim_ang, stim_vec_len = [np.nan] * 5
+ else:
+ stim_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[stim_freq - .5, stim_freq + .5])
+ stim_bins, stim_kde = vonmises_kde(stim_spike_phase)
+ stim_ang, stim_vec_len = spike_phase_score(stim_bins, stim_kde)
+ stim_kde_data.update({name: stim_kde})
+ stim_bins_data.update({name: stim_bins})
+
+ coherence_data.update({name: p_xy})
+ freqency_data.update({name: freq})
+
+ result = pd.Series({
+ 'best_channel': best_channel,
+ 'theta_freq': theta_f,
+ 'theta_peak': theta_p_max,
+ 'theta_energy': theta_energy,
+ 'theta_half_f1': theta_half_f1,
+ 'theta_half_f2': theta_half_f2,
+ 'theta_half_width': theta_half_width,
+ 'theta_half_energy': theta_half_energy,
+ 'theta_ang': theta_ang,
+ 'theta_vec_len': theta_vec_len,
+ 'stim_freq': stim_freq,
+ 'stim_p_max': stim_p_max,
+ 'stim_half_f1': stim_half_f1,
+ 'stim_half_f2': stim_half_f2,
+ 'stim_half_width': stim_half_width,
+ 'stim_energy': stim_energy,
+ 'stim_ang': stim_ang,
+ 'stim_vec_len': stim_vec_len
+ })
+ return result
+
results = units.merge(
+ units.progress_apply(process, axis=1),
+ left_index=True, right_index=True)
+
# coher, freqs = {}, {}
+# for i, row in tqdm(units.iterrows(), total=len(units)):
+# action_id = row['action']
+# channel_group = row['channel_group']
+# unit_name = row['unit_name']
+
+# name = f'{action_id}_{channel_group}_{unit_name}'
+
+# lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response
+
+# sptr = data_loader.spike_train(action_id, channel_group, unit_name)
+
+# lim = get_lim(action_id)
+
+# p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)
+
+# snls = signaltonoise(clean_lfp)
+# best_channel = np.argmax(snls)
+# snl = snls[best_channel]
+# p_xy = p_xys[:,best_channel].magnitude
+# freq = freq.magnitude
+
+# coher.update({name: p_xy})
+# freqs.update({name: freq})
+
pd.DataFrame(coherence_data).to_feather(output / 'data' / 'coherence.feather')
+pd.DataFrame(freqency_data).to_feather(output / 'data' / 'freqs.feather')
+pd.DataFrame(theta_kde_data).to_feather(output / 'data' / 'theta_kde.feather')
+pd.DataFrame(theta_bins_data).to_feather(output / 'data' / 'theta_bins.feather')
+pd.DataFrame(stim_kde_data).to_feather(output / 'data' / 'stim_kde.feather')
+pd.DataFrame(stim_bins_data).to_feather(output / 'data' / 'stim_bins.feather')
+
action = project.require_action("stimulus-spike-lfp-response-other-drive")
+
action.modules['parameters'] = {
+ 'NFFT': NFFT,
+ 'theta_band_f1': theta_band_f1,
+ 'theta_band_f2': theta_band_f2
+}
+
action.data['results'] = 'results.csv'
+results.to_csv(action.data_path('results'), index=False)
+
copy_tree(output, str(action.data_path()))
+
septum_mec.analysis.registration.store_notebook(action, "10-calculate-stimulus-spike-lfp-response-other-drive.ipynb")
+
%load_ext autoreload
+%autoreload 2
+
import os
+import expipe
+import pathlib
+import numpy as np
+import spatial_maps.stats as stats
+import septum_mec
+import septum_mec.analysis.data_processing as dp
+import septum_mec.analysis.registration
+import head_direction.head as head
+import spatial_maps as sp
+import speed_cells.speed as spd
+import re
+import joblib
+import multiprocessing
+import shutil
+import psutil
+import pandas as pd
+import matplotlib.pyplot as plt
+import matplotlib
+import seaborn as sns
+from distutils.dir_util import copy_tree
+from neo import SpikeTrain
+import scipy
+from functools import reduce
+from tqdm.notebook import tqdm_notebook as tqdm
+tqdm.pandas()
+
+from spikewaveform.core import calculate_waveform_features_from_template, cluster_waveform_features
+
+from septum_mec.analysis.plotting import violinplot, despine
+
+from septum_mec.analysis.statistics import load_data_frames, make_paired_tables, make_statistics_table
+
#################################################
+
+# lfp_location = ''
+# lfp_location = '-other-tetrode'
+lfp_location = '-other-drive'
+
+##################################################
+
%matplotlib inline
+plt.rc('axes', titlesize=12)
+plt.rcParams.update({
+ 'font.size': 12,
+ 'figure.figsize': (6, 4),
+ 'figure.dpi': 150
+})
+
+output_path = pathlib.Path("output") / ("stimulus-spike-lfp-response" + lfp_location)
+(output_path / "statistics").mkdir(exist_ok=True, parents=True)
+(output_path / "figures").mkdir(exist_ok=True, parents=True)
+output_path.mkdir(exist_ok=True)
+
project_path = dp.project_path()
+project = expipe.get_project(project_path)
+actions = project.actions
+
data, labels, colors, queries = load_data_frames()
+
lfp_action = actions['stimulus-spike-lfp-response' + lfp_location]
+lfp_results = pd.read_csv(lfp_action.data_path('results'))
+
# lfp_results has old unit id's but correct on (action, unit_name, channel_group)
+lfp_results = lfp_results.drop('unit_id', axis=1)
+
data = data.merge(lfp_results, how='left')
+
data['stim_strength'] = data.stim_p_max / data.theta_peak
+
keys = [
+ 'theta_energy',
+ 'theta_peak',
+ 'theta_freq',
+ 'theta_half_width',
+ 'theta_vec_len',
+ 'theta_ang',
+ 'stim_energy',
+ 'stim_half_width',
+ 'stim_p_max',
+ 'stim_strength',
+ 'stim_vec_len',
+ 'stim_ang'
+]
+
results, labels = make_paired_tables(data, keys)
+
results['gridcell']['theta_peak']
+
xlabel = {
+ 'theta_energy': 'Theta coherence energy',
+ 'theta_peak': 'Theta peak coherence',
+ 'theta_freq': '(Hz)',
+ 'theta_half_width': '(Hz)',
+ 'theta_vec_len': 'a.u.',
+ 'theta_ang': 'rad'
+}
+for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
+ for key in xlabel:
+ fig = plt.figure(figsize=(3.7,2.2))
+ plt.suptitle(key + ' ' + cell_type)
+ legend_lines = []
+ for color, label in zip(colors, labels):
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, label=label))
+ sns.kdeplot(data=results[cell_type][key].loc[:,labels], cumulative=True, legend=False, palette=colors, common_norm=False)
+ plt.xlabel(xlabel[key])
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ despine()
+ figname = f'spike-lfp-coherence-histogram-{key}-{cell_type}'.replace(' ', '-')
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
xlabel = {
+ 'stim_energy': 'Coherence energy',
+ 'stim_half_width': '(Hz)',
+ 'stim_p_max': 'Peak coherence',
+ 'stim_strength': 'Ratio',
+ 'stim_vec_len': 'a.u.',
+ 'stim_ang': 'rad'
+}
+# key = 'theta_energy'
+# key = 'theta_peak'
+for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
+ for key in xlabel:
+ fig = plt.figure(figsize=(3.3,2.2))
+ plt.suptitle(key + ' ' + cell_type)
+ legend_lines = []
+ for color, label in zip(colors[1::2], labels[1::2]):
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, label=label))
+ sns.kdeplot(data=results[cell_type][key].loc[:,labels[1::2]], cumulative=True, legend=False, palette=colors[1::2], common_norm=False)
+ plt.xlabel(xlabel[key])
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+ plt.grid(False)
+ despine()
+ figname = f'spike-lfp-coherence-histogram-{key}-{cell_type}'.replace(' ', '-')
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
from septum_mec.analysis.statistics import VonMisesKDE
+
for paradigm in ['stim', 'theta']:
+ key = paradigm + '_vec_len'
+ for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
+ fig = plt.figure(figsize=(3.2,2.2))
+ plt.suptitle(key + ' ' + cell_type)
+ legend_lines = []
+ for color, query, label in zip(colors, queries, labels):
+ data_query = data.query(query + ' and ' + cell_type)
+ values = data_query[key].values
+ angles = data_query[paradigm + '_ang'].values
+ kde = VonMisesKDE(angles, weights=values, kappa=5)
+ bins = np.linspace(-np.pi, np.pi, 100)
+ plt.polar(bins, kde.evaluate(bins), color=color, lw=2)
+ plt.polar(angles, values, color=color, lw=1, ls='none', marker='.', markersize=2)
+# values.hist(
+# bins=bins[key], density=density, cumulative=cumulative, lw=lw,
+# histtype=histtype, color=color)
+ legend_lines.append(matplotlib.lines.Line2D([0], [0], color=color, lw=2, label=label))
+ plt.legend(
+ handles=legend_lines,
+ bbox_to_anchor=(1.04,1), borderaxespad=0, frameon=False)
+ plt.tight_layout()
+# plt.grid(False)
+ figname = f'spike-lfp-polar-plot-{paradigm}-{cell_type}'.replace(' ', '-')
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
stats = {}
+for cell_type, result in results.items():
+ stats[cell_type], _ = make_statistics_table(result, labels)
+
stats['gridcell']
+
for cell_type, stat in stats.items():
+ stat.to_latex(output_path / "statistics" / f"statistics_{cell_type}.tex")
+ stat.to_csv(output_path / "statistics" / f"statistics_{cell_type}.csv")
+
for cell_type, cell_results in results.items():
+ for key, result in cell_results.items():
+ result.to_latex(output_path / "statistics" / f"values_{cell_type}_{key}.tex")
+ result.to_csv(output_path / "statistics" / f"values_{cell_type}_{key}.csv")
+
from septum_mec.analysis.plotting import plot_bootstrap_timeseries
+
coher = pd.read_feather(output_path / 'data' / 'coherence.feather')
+freqs = pd.read_feather(output_path / 'data' / 'freqs.feather')
+
freq = freqs.T.iloc[0].values
+
+mask = (freq < 100)
+
for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
+ fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(5,2))
+ axs = axs.repeat(2)
+ for i, (ax, query) in enumerate(zip(axs.ravel(), queries)):
+ selection = [
+ f'{r.action}_{r.channel_group}_{r.unit_name}'
+ for i, r in data.query(query + ' and ' + cell_type).iterrows()]
+ values = coher.loc[mask, selection].dropna(axis=1).to_numpy()
+ values = 10 * np.log10(values)
+ plot_bootstrap_timeseries(freq[mask], values, ax=ax, lw=1, label=labels[i], color=colors[i])
+ # ax.set_title(titles[i])
+ ax.set_xlabel('Frequency Hz')
+ ax.legend(frameon=False)
+ ax.set_ylim(-30, 0)
+ axs[0].set_ylabel('Coherence')
+ despine()
+ figname = f'spike-lfp-coherence-{cell_type}'.replace(' ', '-')
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.png',
+ bbox_inches='tight', transparent=True)
+ fig.savefig(
+ output_path / 'figures' / f'{figname}.svg',
+ bbox_inches='tight', transparent=True)
+
nsi_vs_nsni = {}
+for key in keys:
+ df = pd.DataFrame()
+ dfs = [results[k][key].loc[:, ['entity', 'unit_idnum', 'Baseline I']].rename({'Baseline I': k}, axis=1) for k in ['ns_inhibited', 'ns_not_inhibited']]
+ df = pd.merge(*dfs, on=['entity', 'unit_idnum'], how='outer')
+ nsi_vs_nsni[key] = df
+
nsi_vs_nsni.keys()
+
nsi_vs_nsni['theta_energy']
+
from septum_mec.analysis.statistics import LMM
+
LMM(nsi_vs_nsni['theta_energy'], 'ns_inhibited', 'ns_not_inhibited', 'theta_energy')
+
stat, stat_vals = make_statistics_table(nsi_vs_nsni, ['ns_inhibited', 'ns_not_inhibited'], wilcoxon_test=False)
+
stat
+
stat.to_latex(output_path / "statistics" / f"statistics_nsi_vs_nsni.tex")
+stat.to_csv(output_path / "statistics" / f"statistics_nsi_vs_nsni.csv")
+
action = project.require_action("stimulus-spike-lfp-response" + lfp_location)
+
copy_tree(output_path, str(action.data_path()))
+
septum_mec.analysis.registration.store_notebook(action, "20_stimulus-spike-lfp-response.ipynb")
+
+
\n", + " | entity | \n", + "unit_idnum | \n", + "channel_group | \n", + "date | \n", + "Baseline I | \n", + "11 Hz | \n", + "Baseline II | \n", + "30 Hz | \n", + "
---|---|---|---|---|---|---|---|---|
51 | \n", + "1833 | \n", + "8 | \n", + "0 | \n", + "20719 | \n", + "0.214949 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
85 | \n", + "1833 | \n", + "13 | \n", + "0 | \n", + "20719 | \n", + "NaN | \n", + "0.052683 | \n", + "0.291401 | \n", + "0.055222 | \n", + "
86 | \n", + "1833 | \n", + "14 | \n", + "0 | \n", + "20719 | \n", + "NaN | \n", + "0.025182 | \n", + "0.451846 | \n", + "NaN | \n", + "
58 | \n", + "1833 | \n", + "23 | \n", + "0 | \n", + "200619 | \n", + "0.233113 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
127 | \n", + "1833 | \n", + "26 | \n", + "0 | \n", + "200619 | \n", + "NaN | \n", + "NaN | \n", + "0.202562 | \n", + "0.049574 | \n", + "
... | \n", + "... | \n", + "... | \n", + "... | \n", + "... | \n", + "... | \n", + "... | \n", + "... | \n", + "... | \n", + "
139 | \n", + "1849 | \n", + "835 | \n", + "4 | \n", + "150319 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "0.027863 | \n", + "
43 | \n", + "1849 | \n", + "851 | \n", + "5 | \n", + "60319 | \n", + "0.040082 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
65 | \n", + "1849 | \n", + "932 | \n", + "7 | \n", + "280219 | \n", + "0.025449 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
74 | \n", + "1849 | \n", + "937 | \n", + "7 | \n", + "280219 | \n", + "NaN | \n", + "0.135503 | \n", + "NaN | \n", + "NaN | \n", + "
105 | \n", + "1849 | \n", + "939 | \n", + "7 | \n", + "280219 | \n", + "NaN | \n", + "NaN | \n", + "0.435251 | \n", + "NaN | \n", + "
137 rows × 8 columns
\n", + "\n", + " | Theta energy | \n", + "Theta peak | \n", + "Theta freq | \n", + "Theta half width | \n", + "Theta vec len | \n", + "Theta ang | \n", + "Stim energy | \n", + "Stim half width | \n", + "Stim p max | \n", + "Stim strength | \n", + "Stim vec len | \n", + "Stim ang | \n", + "
---|---|---|---|---|---|---|---|---|---|---|---|---|
Baseline I | \n", + "2.0e-01 ± 2.6e-02 (63) | \n", + "1.7e-01 ± 1.9e-02 (63) | \n", + "7.7e+00 ± 7.5e-02 (63) | \n", + "6.4e-01 ± 3.8e-02 (63) | \n", + "2.0e-01 ± 1.5e-02 (63) | \n", + "3.7e+00 ± 1.0e-01 (63) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
11 Hz | \n", + "9.7e-02 ± 1.6e-02 (56) | \n", + "6.9e-02 ± 7.0e-03 (56) | \n", + "7.8e+00 ± 1.3e-01 (56) | \n", + "4.3e-01 ± 7.6e-02 (56) | \n", + "4.0e-02 ± 4.6e-03 (56) | \n", + "3.3e+00 ± 2.7e-01 (56) | \n", + "8.8e-02 ± 9.1e-03 (58) | \n", + "2.2e-01 ± 6.8e-03 (58) | \n", + "4.5e-01 ± 3.1e-02 (58) | \n", + "9.0e+00 ± 8.9e-01 (58) | \n", + "2.2e-01 ± 1.6e-02 (58) | \n", + "2.9e+00 ± 2.6e-01 (58) | \n", + "
Baseline II | \n", + "3.3e-01 ± 4.2e-02 (46) | \n", + "2.5e-01 ± 2.5e-02 (46) | \n", + "8.1e+00 ± 4.4e-02 (46) | \n", + "8.2e-01 ± 6.4e-02 (46) | \n", + "2.3e-01 ± 1.9e-02 (46) | \n", + "4.0e+00 ± 1.1e-01 (46) | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
30 Hz | \n", + "4.3e-02 ± 4.8e-03 (35) | \n", + "4.2e-02 ± 4.7e-03 (35) | \n", + "7.9e+00 ± 2.0e-01 (35) | \n", + "3.0e-01 ± 2.7e-02 (35) | \n", + "1.8e-02 ± 2.6e-03 (35) | \n", + "3.8e+00 ± 3.1e-01 (35) | \n", + "9.1e-02 ± 8.8e-03 (33) | \n", + "2.5e-01 ± 4.0e-03 (33) | \n", + "3.9e-01 ± 3.3e-02 (33) | \n", + "1.3e+01 ± 1.5e+00 (33) | \n", + "2.5e-01 ± 1.5e-02 (33) | \n", + "2.6e+00 ± 2.4e-01 (33) | \n", + "
LMM Baseline I - 11 Hz | \n", + "1.1e-01, 1.3e-01 | \n", + "5.3e-03, 1.2e-01 | \n", + "9.8e-01, -7.5e-03 | \n", + "7.8e-02, 2.7e-01 | \n", + "1.6e-03, 1.5e-01 | \n", + "2.0e-02, 4.3e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
LMM Baseline I - Baseline II | \n", + "1.9e-01, -5.3e-02 | \n", + "1.2e-01, -4.2e-02 | \n", + "7.6e-02, 2.2e-01 | \n", + "6.6e-01, -3.6e-02 | \n", + "2.6e-01, -3.2e-02 | \n", + "6.9e-01, 1.4e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
LMM Baseline I - 30 Hz | \n", + "5.5e-02, 1.6e-01 | \n", + "1.8e-02, 1.4e-01 | \n", + "3.7e-01, -2.3e-01 | \n", + "1.1e-11, 3.8e-01 | \n", + "9.7e-04, 1.8e-01 | \n", + "4.7e-01, 3.5e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
LMM 11 Hz - Baseline II | \n", + "1.2e-02, 2.3e-01 | \n", + "1.1e-02, 1.7e-01 | \n", + "7.1e-01, 1.4e-01 | \n", + "1.4e-12, 4.2e-01 | \n", + "6.3e-04, 1.6e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "
LMM 11 Hz - 30 Hz | \n", + "5.2e-02, -6.2e-02 | \n", + "4.5e-02, -3.2e-02 | \n", + "5.5e-01, 2.8e-01 | \n", + "3.3e-01, -1.3e-01 | \n", + "1.6e-02, -2.7e-02 | \n", + "4.0e-01, 2.8e-01 | \n", + "8.7e-01, 6.2e-03 | \n", + "2.7e-01, 2.8e-02 | \n", + "7.5e-01, -3.0e-02 | \n", + "1.1e-01, 4.5e+00 | \n", + "6.8e-01, 3.0e-02 | \n", + "4.9e-01, -4.2e-01 | \n", + "
LMM Baseline II - 30 Hz | \n", + "1.9e-02, 2.7e-01 | \n", + "2.5e-02, 1.9e-01 | \n", + "8.3e-01, 7.0e-02 | \n", + "6.4e-04, 4.8e-01 | \n", + "4.9e-04, 2.0e-01 | \n", + "6.1e-01, 3.0e-01 | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "NaN | \n", + "