{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "d47ae8a9", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import os\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline \n", "from ccf_utils import get_area_color\n", "import vbn_utils\n", "import decoding_utils as du\n", "import scipy.stats\n", "from vbn_utils import formatFigure" ] }, { "cell_type": "code", "execution_count": null, "id": "dea77128", "metadata": {}, "outputs": [], "source": [ "high_res = True\n", "if high_res:\n", " plt.rcParams['figure.dpi'] = 150\n", " plt.rcParams['savefig.dpi'] = 300\n", " plt.rcParams['font.size'] = 12\n", " plt.rcParams['pdf.fonttype'] = 42\n", "\n", " plt.rcParams['figure.facecolor'] = 'white'\n", " plt.rcParams['axes.facecolor'] = 'white'\n", " plt.rcParams['savefig.facecolor'] = 'white' # affects clipboard copy too\n", " plt.rcParams['savefig.transparent'] = False" ] }, { "cell_type": "code", "execution_count": null, "id": "f9979de0", "metadata": {}, "outputs": [], "source": [ "figure_save_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN Manuscript/VBN_revision_figures\"" ] }, { "cell_type": "code", "execution_count": null, "id": "8e178436", "metadata": {}, "outputs": [], "source": [ "#Paths to all of the useful supplemental tables and tensors\n", "active_tensor_file = \"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/supplemental_tables/vbnAllUnitSpikeTensor.hdf5\"\n", "stim_table_file = \"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/supplemental_tables/master_stim_table_no_filter.csv\"\n", "unit_table_file = \"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/supplemental_tables/master_units_with_responsiveness.csv\"\n", "\n", "sessions_table_file = \"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/vbn_s3_cache/visual-behavior-neuropixels-0.5.0/project_metadata/ecephys_sessions.csv\"" ] }, { "cell_type": "code", "execution_count": null, "id": "593b2a31", "metadata": {}, "outputs": [], "source": [ "structure_tree = pd.read_csv(\"/Volumes/programs/mindscope/workgroups/np-behavior/ccf_structure_tree_2017.csv\")" ] }, { "cell_type": "code", "execution_count": null, "id": "279855b5", "metadata": {}, "outputs": [], "source": [ "units = pd.read_csv(unit_table_file)\n", "units['cortical_layer'].replace('3-Feb', '2/3', inplace=True)\n", "units['cortical_layer'].replace('6a', '6', inplace=True)\n", "units['cortical_layer'].replace('6b', '6', inplace=True)" ] }, { "cell_type": "markdown", "id": "dbee6ea8", "metadata": {}, "source": [ "## Decoding dropouts and sufficiency tests" ] }, { "cell_type": "markdown", "id": "c2630adf", "metadata": {}, "source": [ "### Generated on HPC by 'run_decoding_dropouts.py'" ] }, { "cell_type": "code", "execution_count": null, "id": "b7fb5dfd", "metadata": {}, "outputs": [], "source": [ "change_decoding_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoding_dropouts/pooledChangeDecoding_basesub_active\"\n", "image_decoding_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoding_dropouts/pooledImageDecoding_basesub_active\"" ] }, { "cell_type": "code", "execution_count": null, "id": "3a8213ca", "metadata": {}, "outputs": [], "source": [ "decoding_results_base = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoding_dropouts\"\n", "label_to_dir = {'image': 'pooledImageDecoding_basesub_active',\n", " 'lick': 'pooledLickDecoding_basesub_active',\n", " 'change': 'pooledChangeDecoding_basesub_active',\n", " 'changeprechange': 'pooledChangePrechangeDecoding_basesub_active',\n", " 'reaction_time': 'pooledReactionTimeDecoding_basesub_active',\n", " 'visual_response': 'pooledVisualResponseDecoding_basesub_active',}" ] }, { "cell_type": "code", "execution_count": null, "id": "7d1619c1", "metadata": {}, "outputs": [], "source": [ "subset_files = []\n", "set_files = []\n", "tests = []\n", "labels = []\n", "subset_cell_types = []\n", "subset_layers = []\n", "subset_regions = []\n", "subset_clusters = []\n", "unitsamplesizes = []\n", "experience_levels = []\n", "\n", "for label in ['change', 'image']:\n", " decoding_dir = os.path.join(decoding_results_base, label_to_dir[label])\n", " decoding_files = os.listdir(decoding_dir)\n", " for file in decoding_files:\n", "\n", " if 'full' in file:\n", " continue\n", "\n", " set_info = file.split('_set_')[1].split('_subset_')[0]\n", " subset_info = file.split('subset_')[1]\n", " subset_info = subset_info.replace('VISp_VISl_VISal', 'VISp-VISl-VISal').replace('VISrl_VISpm_VISam', 'VISrl-VISpm-VISam')\n", "\n", " subset_regions.append(subset_info.split('_')[0])\n", " subset_layers.append(subset_info.split('_')[1])\n", " subset_cell_types.append(subset_info.split('_')[2])\n", " subset_clusters.append(subset_info.split('_')[3])\n", "\n", " if ('Familiar' in subset_info) or ('Novel' in subset_info):\n", " exp = 'Familiar' if 'Familiar' in subset_info else 'Novel'\n", " experience_levels.append(exp)\n", " tests.append(subset_info.split('_')[5])\n", " unitsamplesizes.append(subset_info.split('_')[6])\n", " else:\n", " experience_levels.append('all')\n", " tests.append(subset_info.split('_')[4])\n", " unitsamplesizes.append(subset_info.split('_')[5])\n", "\n", " set_file = [os.path.join(decoding_dir, f) for f in decoding_files if (f'_set_{set_info}' in f) and 'full' in f and f.split('_')[-4] == unitsamplesizes[-1]]\n", " if len(set_file) == 0:\n", " set_file = [None]\n", "\n", " subset_file = os.path.join(decoding_dir, file)\n", " if len(subset_file)>=260:\n", " subset_file = r\"\\\\?\\UNC\" + os.path.join(decoding_dir, subset_file)[1:]\n", "\n", " subset_files.append(subset_file)\n", " set_files.append(set_file[0])\n", " labels.append(label)\n", "\n", "data_table = pd.DataFrame({\n", " 'subset_file': subset_files,\n", " 'set_file': set_files,\n", " 'test': tests,\n", " 'label': labels,\n", " 'subset_region': subset_regions,\n", " 'subset_layer': subset_layers,\n", " 'subset_cell_type': subset_cell_types,\n", " 'subset_cluster': subset_clusters,\n", " 'unitsample_size': np.array(unitsamplesizes).astype(int),\n", " 'experience': experience_levels\n", "})\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6aeabd11", "metadata": {}, "outputs": [], "source": [ "from notebook_utils import DiD_test, bootstrapped_diff_ci" ] }, { "cell_type": "markdown", "id": "50f3ec0f", "metadata": {}, "source": [ "### Calculate dropout and sufficiency scores for each neuron subset, comparing performance to the full population with matched unitsamplesize" ] }, { "cell_type": "code", "execution_count": null, "id": "19f43f96", "metadata": {}, "outputs": [], "source": [ "binsize = 10\n", "time = np.arange(binsize, 750+binsize, binsize)\n", "\n", "subsets_to_plot = ['VISall_full_100', 'VISall_23_RS', 'VISall_4_RS', 'VISall_5_RS', 'VISall_6_RS', 'VISall_all_FS', 'VISall_all_SST', 'VISp-VISl-VISal_full_50', \n", " 'VISp-VISl-VISal_all_SST', 'VISrl-VISpm-VISam_all_SST', 'VISall_full_7', 'VISall_all_VIP', ]\n", "subset_labels = ['full pop' if 'full' in s else s for s in subsets_to_plot]\n", "colors = ['k', 'teal', 'teal', 'teal', 'teal', 'red', 'dodgerblue', 'k', 'dodgerblue', 'dodgerblue', 'k', 'orchid',]\n", "alphas = [1, 1, 0.75, 0.5, 0.25, 1, 1, 1, 1, 1, 1, 1]\n", "samplesizes = [40]*10 + [7]*2 \n", "axinds = [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2]\n", "\n", "plot_diff = False\n", "plot_chance = True\n", "\n", "index100 = list(time).index(100)\n", "for decoding_label in ['change', 'image']:\n", " for test in ['dropout', 'sufficiency']:\n", " if test == 'sufficiency':\n", " fig, ax = plt.subplots(1, 3, gridspec_kw={'width_ratios': [7, 4, 2]})\n", " else:\n", " fig, dax = plt.subplots()\n", " ax = [dax, dax, dax]\n", " for ind, (subset, color, alpha, samplesize) in enumerate(zip(subsets_to_plot, colors, alphas, samplesizes)):\n", " \n", " #always makes sense to use 100 sample size for dropouts\n", " if test == 'dropout':\n", " samplesize = 40\n", " if 'full' in subset and ind>0:\n", " continue\n", " \n", " subs_at_100 = []\n", " sets_at_100 = []\n", " for iexp, experience in enumerate(['Familiar', 'Novel']):\n", " if 'full' in subset:\n", " subset_info = data_table[\n", " (data_table['subset_region']==subset.split('_')[0]) & \n", " (data_table['unitsample_size']==samplesize) &\n", " (data_table['test']==test) &\n", " (data_table['label']==decoding_label) &\n", " (data_table['experience']==experience)\n", " ]\n", " subset_file = subset_info['set_file'].values\n", " subset_info = subset_info.iloc[[0]]\n", " else:\n", " subset_info = data_table[\n", " (data_table['subset_region']==subset.split('_')[0]) & \n", " (data_table['subset_layer']==subset.split('_')[1]) & \n", " (data_table['subset_cell_type']==subset.split('_')[2]) & \n", " (data_table['unitsample_size']==samplesize) &\n", " (data_table['test']==test) &\n", " (data_table['label']==decoding_label) &\n", " (data_table['experience']==experience)\n", " ]\n", " \n", " subset_file = subset_info['subset_file'].values\n", " if len(subset_file)>1:\n", " print(f'found more than one compatible file for subset: {subset}')\n", " break\n", "\n", " subset_data = np.load(subset_file[0])\n", " set_data = np.load(subset_info['set_file'].values[0])\n", "\n", " subset_at_100 = subset_data[:, index100]\n", " set_at_100 = set_data[:, index100]\n", "\n", " if plot_diff:\n", " mean_diff = subset_at_100.mean() - set_at_100.mean()\n", " ci95 = bootstrapped_diff_ci(subset_at_100, set_at_100)\n", " else:\n", " mean_diff = subset_at_100.mean()\n", " ci95 = np.percentile(subset_at_100, [2.5, 97.5])\n", "\n", " error = np.std(subset_at_100)\n", "\n", " subs_at_100.append(subset_at_100)\n", " sets_at_100.append(set_at_100)\n", "\n", " ax[axinds[ind]].errorbar([ind + iexp*0.25], [mean_diff], yerr=error, fmt='o', color=color, alpha=alpha)\n", " pval = scipy.stats.ranksums(subset_at_100, set_at_100, nan_policy='omit')\n", " print(f'{subset} sub vs set: {pval}')\n", "\n", " if plot_diff:\n", " DiD_pval = DiD_test(subs_at_100[0], sets_at_100[0], subs_at_100[1], sets_at_100[1])\n", " print(f'DiD p-value for {subset}: {DiD_pval}')\n", " if test == 'sufficiency':\n", " for ia, a in enumerate(ax):\n", " a_inds = np.where(np.array(axinds)==ia)[0]\n", " if len(a_inds)>0:\n", " a.set_xticks(np.arange(len(a_inds))+a_inds[0] + 0.125)\n", " a.set_xticklabels(np.array(subset_labels)[a_inds], rotation=90)\n", " a.set_xlim(a.get_xticks()[0]-0.25, a.get_xticks()[-1]+0.25)\n", "\n", " if plot_chance:\n", " chance = 0.125 if decoding_label=='image' else 0.5\n", " a.axhline(chance, ls='dotted', color='k')\n", " formatFigure(fig, a)\n", " else:\n", " dax.set_xticks([sind for sind, s in enumerate(subset_labels) if (sind==0 or 'full' not in s)])\n", " dax.set_xticklabels([s for sind, s in enumerate(subset_labels) if (sind==0 or 'full' not in s)], rotation=90)\n", " if plot_chance:\n", " chance = 0.125 if decoding_label=='image' else 0.5\n", " dax.axhline(chance, ls='dotted', color='k')\n", " formatFigure(fig, dax)\n", "\n", " \n", " delta_upper = '\\u0394' if plot_diff else ''\n", " ylabel = f'{delta_upper} accuracy: \\n{test} - full population' if plot_diff else f'{test}'\n", " ax[0].set_ylabel(ylabel)\n", " fig.suptitle(f'{test} - {decoding_label}')\n", " fig.tight_layout()\n", " " ] }, { "cell_type": "markdown", "id": "5c5fda77", "metadata": {}, "source": [ "### Look at time lag between familiar and novel image decoding for each subset" ] }, { "cell_type": "code", "execution_count": null, "id": "313ba5c9", "metadata": {}, "outputs": [], "source": [ "from notebook_utils import (\n", " get_sigmoidfit_midpoint\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "f656c0c0", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "warnings.filterwarnings(\n", " \"ignore\",\n", " message=\"overflow encountered in exp\",\n", " category=RuntimeWarning\n", ")\n", "\n", "plt.rcParams['font.size'] = 14\n", "binsize = 10\n", "time = np.arange(binsize, 750+binsize, binsize)\n", "\n", "subsets_to_plot = ['VISall_23_RS', 'VISall_4_RS', 'VISall_5_RS', 'VISall_6_RS', 'VISall_all_FS', 'VISall_all_SST', 'VISp-VISl-VISal_all_SST', 'VISrl-VISpm-VISam_all_SST', 'VISall_all_VIP',]\n", "colors = ['teal', 'teal', 'teal', 'teal', 'red', 'dodgerblue', 'dodgerblue', 'dodgerblue', 'orchid']\n", "alphas = [1, 0.75, 0.5, 0.25, 1, 1, 1, 1, 1]\n", "samplesizes = [40]*8 + [7] \n", "axind = [0,0,0,0,0,0,1,1,2]\n", "\n", "index100 = list(time).index(100)\n", "decoding_label = 'change'\n", "test = 'dropout'\n", "for decoding_label in ['image',]:\n", " for test in ['dropout', 'sufficiency']:\n", " all_fig, all_ax = plt.subplots(1,3, gridspec_kw={'width_ratios': [6,2,1]})\n", " all_fig.suptitle(f'{decoding_label} {test}')\n", " for ind, (subset, color, alpha, samplesize) in enumerate(zip(subsets_to_plot, colors, alphas, samplesizes)):\n", " \n", " if test == 'dropout':\n", " samplesize = 40\n", " \n", " subs_at_100 = []\n", " sets_at_100 = []\n", " exp_latencies = []\n", "\n", " fig, ax = plt.subplots()\n", " fig.suptitle(f'{subset} {test}')\n", " for iexp, experience in enumerate(['Familiar', 'Novel']):\n", " subset_info = data_table[\n", " (data_table['subset_region']==subset.split('_')[0]) & \n", " (data_table['subset_layer']==subset.split('_')[1]) & \n", " (data_table['subset_cell_type']==subset.split('_')[2]) & \n", " (data_table['unitsample_size']==samplesize) &\n", " (data_table['test']==test) &\n", " (data_table['label']==decoding_label) &\n", " (data_table['experience']==experience)\n", " ]\n", " \n", " subset_file = subset_info['subset_file'].values\n", " if len(subset_file)>1:\n", " print(f'found more than one compatible file for subset: {subset}')\n", " break\n", "\n", " subset_data = np.load(subset_file[0])\n", " set_data = np.load(subset_info['set_file'].values[0])\n", "\n", " subset_latencies = np.array([get_sigmoidfit_midpoint(time, sd) for sd in subset_data])\n", " exp_latencies.append(subset_latencies[:,0])\n", "\n", " mean_latency = np.nanmean(subset_latencies[:,0])\n", " \n", " error = np.nanstd(subset_latencies[:,0])\n", "\n", " ax.plot(time, subset_data.T, color=['b', 'r'][iexp], alpha=0.1)\n", " ax.plot(subset_latencies[:,0], subset_latencies[:,1], 'o', color=['b','r'][iexp], alpha=0.1)\n", " all_ax[axind[ind]].errorbar([ind + iexp*0.25], [mean_latency], yerr=error, fmt='o', color=color, alpha=alpha)\n", "\n", " diffs = exp_latencies[0] - exp_latencies[1]\n", " F0 = np.mean(diffs[~np.isnan(diffs)]<=0)\n", " lat_pval = 2 * min(F0, 1 - F0)\n", " \n", " print(f'latency pval for {subset}: {lat_pval}')\n", " ax.set_xlim(0, 100)\n", " ax.set_xlabel('Time from stimulus (ms)')\n", " ax.set_ylabel('Decoding accuracy')\n", " vbn_utils.formatFigure(fig, ax)\n", "\n", " for ia, a in enumerate(all_ax):\n", " a_inds = np.where(np.array(axind)==ia)[0]\n", " a.set_xticks(np.arange(len(a_inds))+a_inds[0] + 0.125)\n", " a.set_xticklabels(np.array(subsets_to_plot)[a_inds], rotation=90)\n", " a.set_xlim(a.get_xticks()[0]-0.25, a.get_xticks()[-1]+0.25)\n", " formatFigure(all_fig, a)\n", "\n", " all_ax[0].set_ylabel('Decoding latency (ms)')\n", " all_fig.suptitle(f'{test} - {decoding_label}')\n", " all_fig.tight_layout()\n" ] }, { "cell_type": "markdown", "id": "e80df0d9", "metadata": {}, "source": [ "## GLM predictions" ] }, { "cell_type": "code", "execution_count": null, "id": "0c9e7a14", "metadata": {}, "outputs": [], "source": [ "cluster_to_color_mapping = {\n", " \"On-trans\": \"#673092\",\n", " \"On-sust\": \"#0353a7\",\n", " \"Off-trans\": \"#0391d9\",\n", " \"Off-sust\": \"#57cfbd\",\n", " \"Stim-supp\": \"#05a05a\",\n", " \"Running\": \"#fea603\",\n", " \"Licking\": \"#fa6a50\",\n", " \"Lick ant\": \"#dc153c\",\n", " \"non-coding\": \"#888888\",\n", "}\n" ] }, { "cell_type": "markdown", "id": "5288dfb7", "metadata": {}, "source": [ "### Generated on HPC by 'run_GLM_prediction_psths.py'" ] }, { "cell_type": "code", "execution_count": null, "id": "5dfa3a32", "metadata": {}, "outputs": [], "source": [ "data_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_glm_prediction_psths\"\n", "unit_files = os.listdir(data_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "332782bd", "metadata": {}, "outputs": [], "source": [ "unitids = []\n", "num_hits = []\n", "num_misses = []\n", "num_nonchange_licks = []\n", "num_nonchange_no_licks = []\n", "predicted_psths = []\n", "predicted_sems = []\n", "psths = []\n", "psth_sems = []\n", "for unit_file in unit_files:\n", "\n", " unit_id = int(unit_file.split('_')[0])\n", "\n", " unit_data = np.load(os.path.join(data_dir, unit_file))\n", " unitids.append(unit_id)\n", " trial_counts = unit_data['trial_counts']\n", " num_hits.append(trial_counts[0])\n", " num_misses.append(trial_counts[1])\n", " num_nonchange_licks.append(trial_counts[2])\n", " num_nonchange_no_licks.append(trial_counts[3])\n", " predicted_psths.append(unit_data['prediction'])\n", " psths.append(unit_data['psth'])\n", " predicted_sems.append(unit_data['prediction_sem'])\n", " psth_sems.append(unit_data['psth_sem'])\n", "\n", "glm_results = pd.DataFrame({'unit_id': unitids, 'num_hits': num_hits, 'num_misses': num_misses, 'num_nonchange_licks': num_nonchange_licks, 'num_nonchange_no_licks': num_nonchange_no_licks, 'predicted_psth': predicted_psths, 'psth': psths, 'psth_sem':psth_sems, 'predicted_sems':predicted_sems})\n" ] }, { "cell_type": "code", "execution_count": null, "id": "65d7455a", "metadata": {}, "outputs": [], "source": [ "glm_results.to_csv(\"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/supplemental_tables/unit_glm_predicted_psths.csv\")" ] }, { "cell_type": "code", "execution_count": null, "id": "e88327fd", "metadata": {}, "outputs": [], "source": [ "from notebook_utils import variance_explained\n", "\n", "glm_results['variance_explained'] = glm_results.apply(variance_explained, axis=1)" ] }, { "cell_type": "code", "execution_count": null, "id": "f9c439d3", "metadata": {}, "outputs": [], "source": [ "glm_results = glm_results.merge(units[['unit_id', 'cluster_labels_new', 'SST', 'VIP', 'RS', 'FS', 'structure_acronym']], on='unit_id')" ] }, { "cell_type": "code", "execution_count": null, "id": "1c57c57d", "metadata": {}, "outputs": [], "source": [ "cluster_labels = ['On-trans', 'On-sust', 'Off-trans', 'Off-sust', 'Stim-supp', 'Lick ant', 'Licking', 'Running', 'action 4', 'action 5', 'action 6', 'non-coding']\n", "\n", "# example cells with representative variance explained\n", "exemplar_cell_ids = [\n", " 1126282659,\n", " 1177946509,\n", " 1179687443,\n", " 1094999579,\n", " 1100910805,\n", " 1094991649,\n", " 1095005619,\n", " 1100926769,\n", "]\n", "\n", "boxfig, boxax = plt.subplots()\n", "for cluster in np.arange(1,13):\n", "\n", " cluster_data = glm_results[(glm_results['cluster_labels_new']==cluster)]\n", " print(f'{cluster} {cluster_data[\"variance_explained\"].median()}')\n", "\n", " color = cluster_to_color_mapping.get(cluster_labels[cluster-1], '#888888')\n", " boxax.boxplot(cluster_data['variance_explained'], positions=[cluster], widths=0.5, showfliers=False, patch_artist=True, whis=[5,95],\n", " boxprops=dict(facecolor=color, color=color),\n", " whiskerprops=dict(color=color),\n", " capprops=dict(color=color),\n", " medianprops=dict(color='w', linewidth=2),\n", " )\n", "\n", "\n", " cluster_data = cluster_data[cluster_data['variance_explained'].between(0.7, 0.85)].sort_values(by='variance_explained', ascending=False) #get representative examples\n", " if cluster_data.empty:\n", " continue\n", " \n", " for ir, row in cluster_data.iloc[:10].iterrows():\n", " if not row['unit_id'] in exemplar_cell_ids:\n", " continue\n", " fig, ax = plt.subplots()\n", " fig.suptitle(f'{cluster_labels[cluster-1]}')\n", " fig.set_size_inches(8,5)\n", " ax.plot(row['predicted_psth']*40, color='k')\n", " ax.fill_between(np.arange(len(row['predicted_psth'])), row['predicted_psth']*40 + row['predicted_sems']*40, row['predicted_psth']*40 - row['predicted_sems']*40, alpha=0.3, color='k')\n", " ax.plot(row['psth']*40, color='orange')\n", " ax.fill_between(np.arange(len(row['psth'])), row['psth']*40 + row['psth_sem']*40, row['psth']*40 - row['psth_sem']*40, alpha=0.3, color='orange')\n", " ax.axvspan(40, 80, alpha=0.1, color='k', lw=0)\n", " ax.axvspan(120, 160, alpha=0.1, color='k', lw=0)\n", " ax.set_xticks(np.array([0, 10, 40, 50, 80, 90, 120, 130]) + 10)\n", " ax.set_xticklabels([0, 250]*4)\n", " formatFigure(fig,ax)\n", " [ax.text(40*ind + 20, ax.get_ylim()[1], ['hit', 'miss', 'false alarm', 'correct reject'][ind], ha='center', va='top') for ind in range(4)]\n", " ax.set_xlabel('Time from stimulus (ms)')\n", " ax.set_ylabel('Firing Rate (Hz)')\n", " fig.savefig(os.path.join(figure_save_dir, f'GLM_PSTH_variance_explained_{cluster_labels[cluster-1]}_{row[\"unit_id\"]}_test.pdf'))\n", "\n", "boxax.set_xticks(np.arange(1,13))\n", "boxax.set_xticklabels(cluster_labels, rotation=90)\n", "boxax.set_ylabel('Trial-averaged variance explained')\n", "boxax.set_xlim(0.5, 8.5)\n", "\n", "boxfig.savefig(os.path.join(figure_save_dir, f'GLM_PSTH_variance_explained_bycluster_boxplot_test.pdf'))\n" ] }, { "cell_type": "markdown", "id": "c6a2687c", "metadata": {}, "source": [ "## Examining feature weights on change decoders" ] }, { "cell_type": "code", "execution_count": null, "id": "ffcff8b4", "metadata": {}, "outputs": [], "source": [ "import h5py\n", "import math\n", "import warnings\n", "\n", "import sklearn\n", "from sklearn.svm import LinearSVC\n", "import matplotlib as mpl\n", "mpl.rcParams['pdf.fonttype'] = 42\n", "import decoding_utils as du" ] }, { "cell_type": "code", "execution_count": null, "id": "4399362e", "metadata": {}, "outputs": [], "source": [ "stim_table = pd.read_csv(stim_table_file)\n", "unit_table = pd.read_csv(unit_table_file)\n", "unitData = h5py.File(active_tensor_file)\n", "sessions = pd.read_csv(sessions_table_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "55cda0d4", "metadata": {}, "outputs": [], "source": [ "def decodeSingleClass(sp, y, nCrossVal, model, unitSampleSize, unitIDs):\n", " '''\n", " sp: units x trials\n", " y: labels \n", " nCrossVal: number of cross validation splits to run\n", " unitIDs: list of unit IDs for sp\n", " '''\n", " unitIDs = np.array(unitIDs)\n", " nUnits = sp.shape[0]\n", " warnings.filterwarnings('ignore')\n", " summary = {sampleSize: \n", " {metric: [] for metric in ('trainAccuracy',\n", " 'featureWeights',\n", " 'accuracy',\n", " 'prediction',\n", " 'confidence', \n", " 'balanced_accuracy',\n", " 'unit_ids')}\n", " for sampleSize in unitSampleSize}\n", " \n", " for sampleSize in unitSampleSize:\n", " if nUnits < sampleSize:\n", " continue\n", " if sampleSize>1:\n", " if sampleSize==nUnits:\n", " nSamples = 1\n", " unitSamples = [np.arange(nUnits)]\n", " else:\n", " # >99% chance each neuron is chosen at least once\n", " nSamples = int(math.ceil(math.log(0.01)/math.log(1-sampleSize/nUnits)))\n", " unitSamples = [np.random.choice(nUnits,sampleSize,replace=False) for _ in range(nSamples)]\n", " elif sampleSize==1:\n", " nSamples = nUnits\n", " unitSamples = [[i] for i in range(nUnits)]\n", " else: #if sampleSize<1, just run 1 iteration with all the units available\n", " nSamples = 1\n", " unitSamples = [np.arange(nUnits)]\n", "\n", " for metric in summary[sampleSize]:\n", " summary[sampleSize][metric].append([])\n", "\n", " for unitSamp in unitSamples:\n", " cv = du.trainDecoder(model,sp[unitSamp].T,y,nCrossVal)\n", " summary[sampleSize]['trainAccuracy'][-1].append(np.mean(cv['train_score']))\n", " summary[sampleSize]['featureWeights'][-1].append(np.mean(cv['coef'],axis=0).squeeze())\n", " summary[sampleSize]['accuracy'][-1].append(np.mean(cv['test_score']))\n", " summary[sampleSize]['prediction'][-1].append(cv['predict'])\n", " summary[sampleSize]['confidence'][-1].append(cv['decision_function'])\n", " summary[sampleSize]['balanced_accuracy'][-1].append(sklearn.metrics.balanced_accuracy_score(y.astype(bool), cv['predict'].astype(bool)))\n", " summary[sampleSize]['unit_ids'][-1].append(unitIDs[unitSamp])\n", "\n", " for metric in summary[sampleSize]:\n", " if metric == 'prediction':\n", " summary[sampleSize][metric][-1] = scipy.stats.mode(summary[sampleSize][metric][-1],axis=0)[0][0]\n", " elif metric not in ('featureWeights', 'unit_ids'):\n", " summary[sampleSize][metric][-1] = np.median(summary[sampleSize][metric][-1],axis=0)\n", " \n", " warnings.filterwarnings('default')\n", " \n", " return summary\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6c0f45b6", "metadata": {}, "outputs": [], "source": [ "regions = ('VISall', 'VISp', 'VISl', 'VISrl', 'VISal', 'VISpm', 'VISam', 'LGd', 'LP', 'SCMRN', 'Hipp')" ] }, { "cell_type": "markdown", "id": "a3c54a3c", "metadata": {}, "source": [ "### To generate flash metrics (skip to load pre-computed flash metrics)" ] }, { "cell_type": "code", "execution_count": null, "id": "4a1c3de2", "metadata": {}, "outputs": [], "source": [ "sessions = sessions[sessions['abnormal_histology'].isnull() & sessions['abnormal_activity'].isnull()]\n", "len(sessions)" ] }, { "cell_type": "code", "execution_count": null, "id": "28368165", "metadata": {}, "outputs": [], "source": [ "unitSampleSize = [40,]\n", "respWin = slice(20,100)\n", "\n", "change_decoder_results = {s: {r:[] for r in regions} for s in sessions['ecephys_session_id'].values}\n", "for count, (isess, session) in enumerate(sessions.iterrows()):\n", " if count%10==0:\n", " print(count)\n", "\n", " sessionId = session['ecephys_session_id']\n", "\n", " stim = stim_table[(stim_table['session_id']==sessionId) & stim_table['active']].reset_index()\n", "\n", " units = unit_table.set_index('unit_id').loc[unitData[str(sessionId)]['unitIds'][:]]\n", " spikes = unitData[str(sessionId)]['spikes']\n", " highQuality = du.apply_unit_quality_filter(units, no_abnorm=False)\n", "\n", " for region in regions:\n", " inRegion = du.getUnitsInRegion(units,region)\n", " final_unit_filter = highQuality & inRegion\n", "\n", " if np.sum(final_unit_filter)<10:\n", " continue\n", "\n", " sp = np.zeros((final_unit_filter.sum(),spikes.shape[1],spikes.shape[2]),dtype=bool)\n", " unitIDs = []\n", " for i,u in enumerate(np.where(final_unit_filter)[0]):\n", " sp[i]=spikes[u,:,:]\n", " unitIDs.append(unitData[str(sessionId)]['unitIds'][:][u])\n", "\n", " flash_response = sp[:, :, respWin].mean(axis=2)\n", "\n", " model = LinearSVC(C=1.0,max_iter=int(1e4), class_weight='balanced')\n", " change_results = decodeSingleClass(flash_response, stim['is_change'].values, 5, model, unitSampleSize, unitIDs)\n", " change_decoder_results[sessionId][region].append(change_results)" ] }, { "cell_type": "code", "execution_count": null, "id": "56852326", "metadata": {}, "outputs": [], "source": [ "dfs = []\n", "for sessionId, decoderdata in change_decoder_results.items():\n", " for region, regiondata in decoderdata.items():\n", " if len(regiondata)>0 and len(regiondata[0][40]['featureWeights'])>0:\n", " dfs.append(pd.DataFrame({'unit_id':np.array(regiondata[0][40]['unit_ids'][0]).flatten(), \n", " 'feature_weights': np.array(regiondata[0][40]['featureWeights'][0]).flatten(),\n", " 'region': region,\n", " 'sessionId': sessionId}))\n", "unit_decoder_df = pd.concat(dfs, ignore_index=True)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4f7e1fc5", "metadata": {}, "outputs": [], "source": [ "unit_decoder_df = unit_decoder_df.merge(unit_table[['unit_id', 'cluster_labels_new']], on='unit_id', )" ] }, { "cell_type": "code", "execution_count": null, "id": "2aecea03", "metadata": {}, "outputs": [], "source": [ "cluster_to_color_mapping = {\n", " \"On-trans\": \"#673092\",\n", " \"On-sust\": \"#0353a7\",\n", " \"Off-trans\": \"#0391d9\",\n", " \"Off-sust\": \"#57cfbd\",\n", " \"Stim-supp\": \"#05a05a\",\n", " \"Running\": \"#fea603\",\n", " \"Licking\": \"#fa6a50\",\n", " \"Lick ant\": \"#dc153c\",\n", " \"non-coding\": \"#888888\",\n", "}\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dfc434a6", "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.dpi'] = 300" ] }, { "cell_type": "code", "execution_count": null, "id": "480111d6", "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.facecolor'] = 'white'\n", "plt.rcParams['axes.facecolor'] = 'white'\n", "plt.rcParams['savefig.facecolor'] = 'white' # affects clipboard copy too\n", "plt.rcParams['savefig.transparent'] = False" ] }, { "cell_type": "code", "execution_count": null, "id": "68b62f35", "metadata": {}, "outputs": [], "source": [ "cluster_labels = ['On-trans', 'On-sust', 'Off-trans', 'Off-sust', 'Stim-supp', 'Lick ant', 'Licking', 'Running', 'action 4', 'action 5', 'action 6', 'non-coding']\n", "clusters_to_plot = list(np.arange(1,9)) + [12,]\n", "clusters_labels_to_plot = [cluster_labels[ind-1] for ind in clusters_to_plot]\n", "\n", "for region in regions + ['all',]:\n", " fig, ax = plt.subplots(1,2)\n", " fig.set_size_inches(12,6)\n", " fig.suptitle(region)\n", "\n", " if region == 'all':\n", " region_df = unit_decoder_df\n", " else:\n", " region_df = unit_decoder_df[unit_decoder_df['region'] == region]\n", "\n", " region_df_cluster_pivot = region_df.pivot_table(index=['sessionId', 'cluster_labels_new'], values='feature_weights')\n", " cluster_session_vals = []\n", " for c in clusters_to_plot:\n", " if (not c in region_df['cluster_labels_new'].unique()):\n", " cluster_session_vals.append([np.nan])\n", " else:\n", " cluster_session_vals.append(region_df_cluster_pivot.xs(c, level=1)['feature_weights'].values)\n", "\n", " cluster_session_vals = [c if len(c)>=3 else [np.nan] for c in cluster_session_vals]\n", "\n", " for ic, (c, clabel) in enumerate(zip(clusters_to_plot, clusters_labels_to_plot)):\n", " color = cluster_to_color_mapping.get(clabel, \"#888888\")\n", " ax[0].boxplot(cluster_session_vals[ic], positions=[ic], widths=0.5, showfliers=False, patch_artist=True,\n", " boxprops=dict(facecolor=color, color=color),\n", " whiskerprops=dict(color=color),\n", " capprops=dict(color=color),\n", " medianprops=dict(color='w', linewidth=2),\n", " )\n", " ax[0].set_xlabel('Cluster ID')\n", " ax[0].set_xticks(np.arange(len(clusters_to_plot)))\n", " ax[0].set_xticklabels(clusters_labels_to_plot, rotation=90)\n", " ax[0].set_ylabel('Feature Weight')\n", " ax[0].axhline(0, ls='dotted', c='k')\n", " vbn_utils.plot_comparison_matrix(*cluster_session_vals, ax=ax[1], labels=clusters_labels_to_plot, cmap='PiYG', colorbar=True, binarize=True)\n", "\n", " plt.tight_layout()" ] } ], "metadata": { "kernelspec": { "display_name": "vbn_manuscript", "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.20" } }, "nbformat": 4, "nbformat_minor": 5 }