+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [1]:
+
+
+
+
+
+
+%load_ext autoreload
+%autoreload 2
+
+
+
+
+
+In [2]:
+
+
+
+
+
+
+import os
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+import re
+import shutil
+import pandas as pd
+import scipy.stats
+
+import exdir
+import expipe
+from distutils.dir_util import copy_tree
+import septum_mec
+import spatial_maps as sp
+import head_direction.head as head
+import septum_mec.analysis.data_processing as dp
+import septum_mec.analysis.registration
+from septum_mec.analysis.plotting import violinplot
+
+from spike_statistics.core import permutation_resampling
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [3]:
+
+
+
+
+
+
+project_path = dp.project_path()
+project = expipe.get_project(project_path)
+actions = project.actions
+
+output_path = pathlib.Path("output") / "comparisons-gridcells"
+(output_path / "statistics").mkdir(exist_ok=True, parents=True)
+(output_path / "figures").mkdir(exist_ok=True, parents=True)
+
+
+
+
+
+Load cell statistics and shuffling quantiles¶
+
+
+
+
+
+In [4]:
+
+
+
+
+
+
+statistics_action = actions['calculate-statistics']
+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')
+statistics_results = pd.read_csv(statistics_action.data_path('results'))
+statistics = pd.merge(session_units, statistics_results, how='left')
+statistics.head()
+
+
+
+
+
+
+
+
+
+
+
+Out[4]:
+
+
+
+
+
+
+
+
+
+
+In [5]:
+
+
+
+
+
+
+shuffling = actions['shuffling']
+quantiles_95 = pd.read_csv(shuffling.data_path('quantiles_95'))
+quantiles_95.head()
+
+
+
+
+
+
+
+
+
+
+
+Out[5]:
+
+
+
+
+
+
+
+
+
+
+In [6]:
+
+
+
+
+
+
+action_columns = ['action', 'channel_group', 'unit_name']
+data = pd.merge(statistics, quantiles_95, on=action_columns, suffixes=("", "_threshold"))
+
+data['specificity'] = np.log10(data['in_field_mean_rate'] / data['out_field_mean_rate'])
+
+data.head()
+
+
+
+
+
+
+
+
+
+
+
+Out[6]:
+
+
+
+
+
+
+
+
+
+
+Statistics about all cell-sessions¶
+
+
+
+
+
+In [12]:
+
+
+
+
+
+
+data.groupby('stimulated').count()['action']
+
+
+
+
+
+
+
+
+
+
+
+Out[12]:
+
+
+
+
+
+
+
+
+
+
+
+Find all cells with gridness above threshold¶
+
+
+
+
+
+In [11]:
+
+
+
+
+
+
+query = 'gridness > gridness_threshold and information_rate > information_rate_threshold'
+sessions_above_threshold = data.query(query)
+print("Number of gridcells", len(sessions_above_threshold))
+print("Number of animals", len(sessions_above_threshold.groupby(['entity'])))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [13]:
+
+
+
+
+
+
+baseline = sessions_above_threshold.query('baseline')
+
+
+
+
+
+In [14]:
+
+
+
+
+
+
+gridcell_in_baseline = data[data.unit_id.isin(baseline.unit_id)]
+
+
+
+
+
+In [18]:
+
+
+
+
+
+
+baseline_i = gridcell_in_baseline.query('baseline and i')
+stimulated_11 = gridcell_in_baseline.query('frequency==11 and stim_location=="ms" and i')
+
+baseline_ii = gridcell_in_baseline.query('baseline and ii')
+stimulated_30 = gridcell_in_baseline.query('frequency==30 and stim_location=="ms" and ii')
+
+print("Number of gridcells in baseline i sessions", len(baseline_i))
+print("Number of gridcells in stimulated 11Hz ms sessions", len(stimulated_11))
+
+print("Number of gridcells in baseline ii sessions", len(baseline_ii))
+print("Number of gridcells in stimulated 30Hz ms sessions", len(stimulated_30))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+slice unique units¶
+
+
+
+
+
+In [19]:
+
+
+
+
+
+
+baseline_i = baseline_i.drop_duplicates('unit_id')
+stimulated_11 = stimulated_11.drop_duplicates('unit_id')
+baseline_ii = baseline_ii.drop_duplicates('unit_id')
+stimulated_30 = stimulated_30.drop_duplicates('unit_id')
+
+
+
+
+
+In [20]:
+
+
+
+
+
+
+print("Number of gridcells in baseline i sessions", len(baseline_i))
+print("Number of gridcells in stimulated 11Hz ms sessions", len(stimulated_11))
+
+print("Number of gridcells in baseline ii sessions", len(baseline_ii))
+print("Number of gridcells in stimulated 30Hz ms sessions", len(stimulated_30))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Calculate statistics¶
+
+
+
+
+
+In [21]:
+
+
+
+
+
+
+columns = [
+ 'average_rate', 'gridness', 'sparsity', 'selectivity', 'information_specificity',
+ 'max_rate', 'information_rate', 'interspike_interval_cv',
+ 'in_field_mean_rate', 'out_field_mean_rate',
+ 'burst_event_ratio', 'specificity', 'speed_score'
+]
+
+
+
+
+
+In [22]:
+
+
+
+
+
+
+gridcell_in_baseline.groupby('stimulated')[columns].mean()
+
+
+
+
+
+
+
+
+
+
+
+Out[22]:
+
+
+
+
+
+
+
+
+
+
+In [23]:
+
+
+
+
+
+
+gridcell_in_baseline.query('baseline')[columns].describe()
+
+
+
+
+
+
+
+
+
+
+
+Out[23]:
+
+
+
+
+
+
+
+
+
+
+In [24]:
+
+
+
+
+
+
+gridcell_in_baseline.query("stimulated")[columns].describe()
+
+
+
+
+
+
+
+
+
+
+
+Out[24]:
+
+
+
+
+
+
+
+
+
+
+Create nice table¶
+
+
+
+
+
+In [25]:
+
+
+
+
+
+
+def summarize(data):
+ return "{:.2f} ± {:.2f} ({})".format(data.mean(), data.sem(), sum(~np.isnan(data)))
+
+
+def MWU(column, stim, base):
+ '''
+ Mann Whitney U
+ '''
+ Uvalue, pvalue = scipy.stats.mannwhitneyu(
+ stim[column].dropna(),
+ base[column].dropna(),
+ alternative='two-sided')
+
+ return "{:.2f}, {:.3f}".format(Uvalue, pvalue)
+
+
+def PRS(column, stim, base):
+ '''
+ Permutation ReSampling
+ '''
+ pvalue, observed_diff, diffs = permutation_resampling(
+ stim[column].dropna(),
+ base[column].dropna(), statistic=np.median)
+
+ return "{:.2f}, {:.3f}".format(observed_diff, pvalue)
+
+
+def rename(name):
+ return name.replace("_field", "-field").replace("_", " ").capitalize()
+
+
+
+
+
+In [26]:
+
+
+
+
+
+
+_stim_data = gridcell_in_baseline.query('stimulated')
+_base_data = gridcell_in_baseline.query('baseline')
+
+result = pd.DataFrame()
+
+result['Stimulated'] = _stim_data[columns].agg(summarize)
+result['Baseline'] = _base_data[columns].agg(summarize)
+
+result.index = map(rename, result.index)
+
+result['MWU'] = list(map(lambda x: MWU(x, _stim_data, _base_data), columns))
+result['PRS'] = list(map(lambda x: PRS(x, _stim_data, _base_data), columns))
+
+result.to_latex(output_path / "statistics" / "statistics.tex")
+result.to_latex(output_path / "statistics" / "statistics.csv")
+result
+
+
+
+
+
+
+
+
+
+
+
+Out[26]:
+
+
+
+
+
+
+
+
+
+
+In [27]:
+
+
+
+
+
+
+_stim_data = stimulated_11
+_base_data = baseline_i
+
+result = pd.DataFrame()
+
+result['Stimulated'] = _stim_data[columns].agg(summarize)
+result['Baseline'] = _base_data[columns].agg(summarize)
+
+result.index = map(rename, result.index)
+
+result['MWU'] = list(map(lambda x: MWU(x, _stim_data, _base_data), columns))
+result['PRS'] = list(map(lambda x: PRS(x, _stim_data, _base_data), columns))
+
+
+result.to_latex(output_path / "statistics" / "statistics_11.tex")
+result.to_latex(output_path / "statistics" / "statistics_11.csv")
+result
+
+
+
+
+
+
+
+
+
+
+
+Out[27]:
+
+
+
+
+
+
+
+
+
+
+In [28]:
+
+
+
+
+
+
+_stim_data = stimulated_30
+_base_data = baseline_ii
+
+result = pd.DataFrame()
+
+result['Stimulated'] = _stim_data[columns].agg(summarize)
+result['Baseline'] = _base_data[columns].agg(summarize)
+
+result.index = map(rename, result.index)
+
+result['MWU'] = list(map(lambda x: MWU(x, _stim_data, _base_data), columns))
+result['PRS'] = list(map(lambda x: PRS(x, _stim_data, _base_data), columns))
+
+
+result.to_latex(output_path / "statistics" / "statistics_30.tex")
+result.to_latex(output_path / "statistics" / "statistics_30.csv")
+result
+
+
+
+
+
+
+
+
+
+
+
+Out[28]:
+
+
+
+
+
+
+
+
+
+
+In [29]:
+
+
+
+
+
+
+_stim_data = stimulated_30
+_base_data = stimulated_11
+
+result = pd.DataFrame()
+
+result['Stimulated 30Hz'] = _stim_data[columns].agg(summarize)
+result['Stimulated 11Hz'] = _base_data[columns].agg(summarize)
+
+result.index = map(rename, result.index)
+
+result['MWU'] = list(map(lambda x: MWU(x, _stim_data, _base_data), columns))
+result['PRS'] = list(map(lambda x: PRS(x, _stim_data, _base_data), columns))
+
+
+result.to_latex(output_path / "statistics" / "statistics_11_vs_30.tex")
+result.to_latex(output_path / "statistics" / "statistics_11_vs_30.csv")
+result
+
+
+
+
+
+
+
+
+
+
+
+Out[29]:
+
+
+
+
+
+
+
+
+
+
+In [30]:
+
+
+
+
+
+
+_stim_data = baseline_i
+_base_data = baseline_ii
+
+result = pd.DataFrame()
+
+result['Baseline i'] = _stim_data[columns].agg(summarize)
+result['Baseline ii'] = _base_data[columns].agg(summarize)
+
+result.index = map(rename, result.index)
+
+result['MWU'] = list(map(lambda x: MWU(x, _stim_data, _base_data), columns))
+result['PRS'] = list(map(lambda x: PRS(x, _stim_data, _base_data), columns))
+
+
+result.to_latex(output_path / "statistics" / "statistics_base_i_vs_base_ii.tex")
+result.to_latex(output_path / "statistics" / "statistics_base_i_vs_base_ii.csv")
+result
+
+
+
+
+
+
+
+
+
+
+
+Out[30]:
+
+
+
+
+
+
+
+
+
+
+Violinplot¶
+
+
+
+
+
+In [47]:
+
+
+
+
+
+
+%matplotlib inline
+plt.rc('axes', titlesize=12)
+plt.rcParams.update({
+ 'font.size': 12,
+ 'figure.figsize': (1.7, 3),
+ 'figure.dpi': 150
+})
+
+
+
+
+
+In [48]:
+
+
+
+
+
+
+stuff = {
+ '': {
+ 'base': gridcell_in_baseline.query('baseline'),
+ 'stim': gridcell_in_baseline.query('stimulated')
+ },
+ '_11': {
+ 'base': baseline_i,
+ 'stim': stimulated_11
+ },
+ '_30': {
+ 'base': baseline_ii,
+ 'stim': stimulated_30
+ }
+}
+
+
+
+
+
+Information rate¶
+
+
+
+
+
+In [49]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['information_rate'].to_numpy()
+ stimulated = data['stim']['information_rate'].to_numpy()
+ print(key)
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Spatial information")
+ plt.ylabel("bits/s")
+ plt.ylim(-0.2, 4)
+
+ plt.savefig(output_path / "figures" / f"spatial_information{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"spatial_information{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [50]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['specificity'].to_numpy()
+ stimulated = data['stim']['specificity'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Spatial specificity")
+ plt.ylabel("")
+ plt.ylim(-0.02, 1.25)
+ plt.savefig(output_path / "figures" / f"specificity{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"specificity{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [51]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['average_rate'].to_numpy()
+ stimulated = data['stim']['average_rate'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Average rate")
+ plt.ylabel("spikes/s")
+ plt.ylim(-0.2, 40)
+
+ plt.savefig(output_path / "figures" / f"average_rate{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"average_rate{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [52]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['max_rate'].to_numpy()
+ stimulated = data['stim']['max_rate'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Max rate")
+ plt.ylabel("spikes/s")
+ # plt.ylim(-0.2, 45)
+
+ plt.savefig(output_path / "figures" / f"max_rate{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"max_rate{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [53]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['interspike_interval_cv'].to_numpy()
+ stimulated = data['stim']['interspike_interval_cv'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("ISI CV")
+ plt.ylabel("Coefficient of variation")
+ # plt.ylim(0.9, 5)
+
+ plt.savefig(output_path / "figures" / f"isi_cv{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"isi_cv{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [54]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['in_field_mean_rate'].to_numpy()
+ stimulated = data['stim']['in_field_mean_rate'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("In-field rate")
+ plt.ylabel("spikes/s")
+ # plt.ylim(-0.1, 18)
+
+ plt.savefig(output_path / "figures" / f"in_field_mean_rate{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"in_field_mean_rate{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [55]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['out_field_mean_rate'].to_numpy()
+ stimulated = data['stim']['out_field_mean_rate'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Out-of-field rate")
+ plt.ylabel("spikes/s")
+ # plt.ylim(-0.2, 8)
+
+ plt.savefig(output_path / "figures" / f"out_field_mean_rate{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"out_field_mean_rate{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [56]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['burst_event_ratio'].to_numpy()
+ stimulated = data['stim']['burst_event_ratio'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Bursting ratio")
+ plt.ylabel("")
+ # plt.ylim(-0.02, 0.60)
+
+ plt.savefig(output_path / "figures" / f"burst_event_ratio{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"burst_event_ratio{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [57]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['max_field_mean_rate'].to_numpy()
+ stimulated = data['stim']['max_field_mean_rate'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Mean rate of max field")
+ plt.ylabel("(spikes/s)")
+ # plt.ylim(-0.5,25)
+
+ plt.savefig(output_path / "figures" / f"max_field_mean_rate{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"max_field_mean_rate{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [58]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['bursty_spike_ratio'].to_numpy()
+ stimulated = data['stim']['bursty_spike_ratio'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("ratio of spikes per burst")
+ plt.ylabel("")
+ # plt.ylim(-0.03,0.9)
+
+ plt.savefig(output_path / "figures" / f"bursty_spike_ratio{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"bursty_spike_ratio{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [59]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['gridness'].to_numpy()
+ stimulated = data['stim']['gridness'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Gridness")
+ plt.ylabel("Gridness")
+ plt.ylim(-0.005, 1.5)
+
+ plt.savefig(output_path / "figures" / f"gridness{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"gridness{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [60]:
+
+
+
+
+
+
+for key, data in stuff.items():
+ baseline = data['base']['speed_score'].to_numpy()
+ stimulated = data['stim']['speed_score'].to_numpy()
+ plt.figure()
+ violinplot(baseline, stimulated)
+ plt.title("Speed score")
+ plt.ylabel("Speed score")
+ # plt.ylim(-0.1, 0.5)
+
+ plt.savefig(output_path / "figures" / f"speed_score{key}.svg", bbox_inches="tight")
+ plt.savefig(output_path / "figures" / f"speed_score{key}.png", dpi=600, bbox_inches="tight")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Register in Expipe¶
+
+
+
+
+
+In [44]:
+
+
+
+
+
+
+action = project.require_action("comparisons-gridcells")
+
+
+
+
+
+In [45]:
+
+
+
+
+
+
+copy_tree(output_path, str(action.data_path()))
+
+
+
+
+
+
+
+
+
+
+
+Out[45]:
+
+
+
+
+
+
+
+
+
+
+
+In [46]:
+
+
+
+
+
+
+septum_mec.analysis.registration.store_notebook(action, "20_comparisons_gridcells.ipynb")
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+
+
+