In [1]:
%load_ext autoreload
%autoreload 2
In [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 tqdm.notebook import tqdm_notebook as tqdm
tqdm.pandas()
from spike_statistics.core import permutation_resampling
from spikewaveform.core import calculate_waveform_features_from_template, cluster_waveform_features
from septum_mec.analysis.plotting import violinplot, despine
In [3]:
%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"
(output_path / "statistics").mkdir(exist_ok=True, parents=True)
(output_path / "figures").mkdir(exist_ok=True, parents=True)
output_path.mkdir(exist_ok=True)
In [4]:
data_loader = dp.Data()
actions = data_loader.actions
project = data_loader.project
In [5]:
identification_action = actions['identify-neurons']
sessions = pd.read_csv(identification_action.data_path('sessions'))
units = pd.read_csv(identification_action.data_path('units'))
session_units = pd.merge(sessions, units, on='action')
In [6]:
lfp_action = actions['stimulus-spike-lfp-response']
lfp_results = pd.read_csv(lfp_action.data_path('results'))
In [7]:
# lfp_results has old unit id's but correct on (action, unit_name, channel_group)
lfp_results = lfp_results.drop('unit_id', axis=1)
In [8]:
stim_action = actions['stimulus-response']
stim_results = pd.read_csv(stim_action.data_path('results'))
In [9]:
# lfp_results has old unit id's but correct on (action, unit_name, channel_group)
stim_results = stim_results.drop('unit_id', axis=1)
In [10]:
statistics_action = actions['calculate-statistics']
shuffling = actions['shuffling']
statistics_results = pd.read_csv(statistics_action.data_path('results'))
statistics_results = session_units.merge(statistics_results, how='left')
quantiles_95 = pd.read_csv(shuffling.data_path('quantiles_95'))
action_columns = ['action', 'channel_group', 'unit_name']
data = pd.merge(statistics_results, quantiles_95, on=action_columns, suffixes=("", "_threshold"))
In [11]:
data['unit_day'] = data.apply(lambda x: str(x.unit_idnum) + '_' + x.action.split('-')[1], axis=1)
In [12]:
data = data.merge(lfp_results, how='left')
In [13]:
data = data.merge(stim_results, how='left')
In [14]:
waveform_action = actions['waveform-analysis']
waveform_results = pd.read_csv(waveform_action.data_path('results')).drop('template', axis=1)
In [15]:
data = data.merge(waveform_results, how='left')
In [16]:
colors = ['#1b9e77','#d95f02','#7570b3','#e7298a']
labels = ['Baseline I', '11 Hz', 'Baseline II', '30 Hz']
queries = ['baseline and Hz11', 'frequency==11', 'baseline and Hz30', 'frequency==30']
In [17]:
data.bs = data.bs.astype(bool)
In [18]:
data.loc[data.eval('t_i_peak == t_i_peak and not bs'), 'ns_inhibited'] = True
data.ns_inhibited.fillna(False, inplace=True)
data.loc[data.eval('t_i_peak != t_i_peak and not bs'), 'ns_not_inhibited'] = True
data.ns_not_inhibited.fillna(False, inplace=True)
In [19]:
# make baseline for inhibited vs not inhibited
data.loc[data.unit_id.isin(data.query('ns_inhibited').unit_id.values), 'ns_inhibited'] = True
data.loc[data.unit_id.isin(data.query('ns_not_inhibited').unit_id.values), 'ns_not_inhibited'] = True
In [20]:
query = (
'gridness > gridness_threshold and '
'information_rate > information_rate_threshold and '
'gridness > .2 and '
'average_rate < 25'
)
sessions_above_threshold = data.query(query)
print("Number of sessions above threshold", len(sessions_above_threshold))
print("Number of animals", len(sessions_above_threshold.groupby(['entity'])))
In [21]:
gridcell_sessions = data[data.unit_day.isin(sessions_above_threshold.unit_day.values)]
print("Number of gridcells", gridcell_sessions.unit_idnum.nunique())
print("Number of gridcell recordings", len(gridcell_sessions))
print("Number of animals", len(gridcell_sessions.groupby(['entity'])))
In [22]:
gridcell_sessions
Out[22]:
In [23]:
data.loc[:,'gridcell'] = np.nan
data['gridcell'] = data.isin(gridcell_sessions)
data.loc[data.eval('not gridcell and bs'), 'bs_not_gridcell'] = True
data.bs_not_gridcell.fillna(False, inplace=True)
In [24]:
data.query('baseline and Hz11 and gridcell').head()
Out[24]:
In [25]:
density = True
cumulative = True
histtype = 'step'
lw = 2
bins = {
'theta_energy': np.arange(0, 2.1, .03),
'theta_peak': np.arange(0, 1, .03),
'theta_freq': np.arange(4, 12, .1),
'theta_half_width': np.arange(0, 5, .1)
}
xlabel = {
'theta_energy': 'Theta coherence energy',
'theta_peak': 'Theta peak coherence',
'theta_freq': '(Hz)',
'theta_half_width': '(Hz)',
}
results = {}
for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
results[cell_type] = {}
for key in bins:
results[cell_type][key] = pd.DataFrame()
fig = plt.figure(figsize=(3.5,2.2))
plt.suptitle(key + ' ' + cell_type)
legend_lines = []
for color, query, label in zip(colors, queries, labels):
values = data.query(query + ' and ' + cell_type)[key]
results[cell_type][key] = pd.concat([
results[cell_type][key],
values.rename(label).reset_index(drop=True)], axis=1)
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=lw, label=label))
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)
plt.xlim(-0.05, bins[key].max() - bins[key].max()*0.02)
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)
In [26]:
data['stim_strength'] = data.stim_p_max / data.theta_peak
In [27]:
density = True
cumulative = True
histtype = 'step'
lw = 2
bins = {
'stim_energy': np.arange(0, .4, .01),
'stim_half_width': np.arange(0, 1.5, .01),
'stim_p_max': np.arange(0, 1, .01),
'stim_strength': np.arange(0, 50, 1)
}
xlabel = {
'stim_energy': 'Coherence energy',
'stim_half_width': '(Hz)',
'stim_p_max': 'Peak coherence',
'stim_strength': 'Ratio',
}
# key = 'theta_energy'
# key = 'theta_peak'
for cell_type in ['gridcell', 'ns_inhibited', 'ns_not_inhibited']:
for key in bins:
results[cell_type][key] = pd.DataFrame()
fig = plt.figure(figsize=(3.2,2.2))
plt.suptitle(key + ' ' + cell_type)
legend_lines = []
for color, query, label in zip(colors[1::2], queries[1::2], labels[1::2]):
values = data.query(query + ' and ' + cell_type)[key]
results[cell_type][key] = pd.concat([
results[cell_type][key],
values.rename(label).reset_index(drop=True)], axis=1)
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=lw, label=label))
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)
plt.xlim(-0.05, bins[key].max() - bins[key].max()*0.02)
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)
stats¶
In [28]:
def summarize(data):
return "{:.2f} ± {:.2f} ({})".format(data.mean(), data.sem(), sum(~np.isnan(data)))
def MWU(df, keys):
'''
Mann Whitney U
'''
Uvalue, pvalue = scipy.stats.mannwhitneyu(
df[keys[0]].dropna(),
df[keys[1]].dropna(),
alternative='two-sided')
return "{:.2f}, {:.3f}".format(Uvalue, pvalue)
def PRS(df, keys):
'''
Permutation ReSampling
'''
pvalue, observed_diff, diffs = permutation_resampling(
df[keys[0]].dropna(),
df[keys[1]].dropna(), statistic=np.median)
return "{:.2f}, {:.3f}".format(observed_diff, pvalue)
def rename(name):
return name.replace("_field", "-field").replace("_", " ").capitalize()
In [ ]:
stats = {}
for cell_type in results:
stat = pd.DataFrame()
for key, df in results[cell_type].items():
Key = rename(key)
stat[Key] = df.agg(summarize)
stat[Key] = df.agg(summarize)
for i, c1 in enumerate(df.columns):
for c2 in df.columns[i+1:]:
stat.loc[f'MWU {c1} {c2}', Key] = MWU(df, [c1, c2])
stat.loc[f'PRS {c1} {c2}', Key] = PRS(df, [c1, c2])
stats[cell_type] = stat
In [ ]:
stats['gridcell']
In [ ]:
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")
psd plots¶
In [ ]:
from septum_mec.analysis.plotting import plot_bootstrap_timeseries
In [ ]:
coher = pd.read_feather(output_path / 'data' / 'coherence.feather')
freqs = pd.read_feather(output_path / 'data' / 'freqs.feather')
In [ ]:
freq = freqs.T.iloc[0].values
mask = (freq < 100)
In [ ]:
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)
Store results in Expipe action¶
In [ ]:
action = project.require_action("stimulus-spike-lfp-response")
In [ ]:
copy_tree(output_path, str(action.data_path()))
In [ ]:
septum_mec.analysis.registration.store_notebook(action, "20_stimulus-spike-lfp-response.ipynb")
In [ ]: