{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import os\n", "from matplotlib import pyplot as plt\n", "from ccf_utils import get_area_color\n", "import vbn_utils\n", "import decoding_utils as du\n", "import scipy.stats\n", "%matplotlib inline \n" ] }, { "cell_type": "code", "execution_count": null, "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", "passive_tensor_file = \"/Volumes/programs/mindscope/workgroups/np-behavior/vbn_data_release/supplemental_tables/vbnAllUnitSpikeTensor_passive.hdf5\"\n", "\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, "metadata": {}, "outputs": [], "source": [ "structure_tree = pd.read_csv(\"/Volumes/programs/mindscope/workgroups/np-behavior/ccf_structure_tree_2017.csv\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generated on HPC by 'run_pooled_decoding.py'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "decoding_results_base = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_decoding_from_sensory_action_clusters\"\n", "label_to_dir = {'image': 'pooledImageDecoding',\n", " 'lick': 'pooledLickDecoding',\n", " 'change': 'pooledChangeDecoding',\n", " 'changeprechange': 'pooledChangePrechangeDecoding',\n", " 'reaction_time': 'pooledReactionTimeDecoding',\n", " 'visual_response': 'pooledVisualResponseDecoding',}" ] }, { "cell_type": "code", "execution_count": null, "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'\n", " plt.rcParams['savefig.transparent'] = False" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from notebook_utils import (\n", " plot_decoding_results,\n", " plot_facemap_decoding,\n", " get_latency2,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Session decoding for LP and VISp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "session_change_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_decoding_from_sensory_action_clusters/sessionChangeDecoding_basesub\"\n", "unit_sample_size = 10\n", "\n", "qualifying_files = [f for f in os.listdir(session_change_dir) if f.endswith(f'[{unit_sample_size}].npy')]\n", "qualifying_sessions = [f.split('_')[1] for f in qualifying_files if 'LP' in f or 'VISp' in f]\n", "sessions_with_LP_and_VISp = []\n", "for session in qualifying_sessions:\n", " session_files = [f for f in qualifying_files if f.split('_')[1] == session and ('LP' in f or 'VISp' in f)]\n", " if len([s for s in session_files if 'LP' in s])>0 and len([s for s in session_files if 'VISp' in s])>0:\n", " sessions_with_LP_and_VISp.append(session_files)\n", "\n", "sessions_with_LP_and_VISp = np.unique(sessions_with_LP_and_VISp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "len(sessions_with_LP_and_VISp)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "summary_dict = {'LP':[], 'VISp':[]}\n", "for session_file in sessions_with_LP_and_VISp:\n", " full_path = os.path.join(session_change_dir, session_file)\n", " data = np.load(full_path)\n", " \n", " area = 'LP' if 'LP' in session_file else 'VISp'\n", " print(f'{area} : {data.shape}')\n", " summary_dict[area].append(np.mean(data, axis=0))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "binsize = 10\n", "time = np.arange(binsize, 750+binsize, binsize)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lp_lats = get_latency2(np.stack(summary_dict['LP']), 0.5, time, slice(0,10), upsample_factor=10, norm=True, threshold=0.5, sub='chance')[1]\n", "visp_lats = get_latency2(np.stack(summary_dict['VISp']), 0.5, time, slice(0,10), upsample_factor=10, norm=True, threshold=0.5, sub='chance')[1]\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(lp_lats, visp_lats, 'ko')\n", "ax.plot([45,65], [45,65], 'k--')\n", "ax.set_aspect('equal')\n", "ax.set_xlabel('LP latency')\n", "ax.set_ylabel('VISp latency')\n", "vbn_utils.formatFigure(fig, ax)\n", "\n", "nonans = ~np.isnan(lp_lats) & ~np.isnan(visp_lats)\n", "p = scipy.stats.wilcoxon(lp_lats[nonans], visp_lats[nonans])[1]\n", "ax.text(50, 45, 'p = {:.2e}'.format(p))\n", "\n", "print(np.sum(nonans))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['font.size'] = 18\n", "binsize = 10\n", "time = np.arange(binsize, 750+binsize, binsize)\n", "lp_array = np.stack(summary_dict['LP'])\n", "visp_array = np.stack(summary_dict['VISp'])\n", "above_thresh_lp = np.max(lp_array[:,:10], axis=1) > 0.6\n", "above_thresh_visp = np.max(visp_array[:,:10], axis=1) > 0.6\n", "\n", "above_thresh = above_thresh_lp & above_thresh_visp\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(time, (lp_array[above_thresh]).T, color='r', alpha=0.1)\n", "ax.plot(time, (visp_array[above_thresh]).T, color='b', alpha=0.1)\n", "ax.plot(time, (lp_array[above_thresh]).mean(axis=0), color='r', label='LP')\n", "ax.plot(time, (visp_array[above_thresh]).mean(axis=0), color='b', label='VISp')\n", "ax.legend()\n", "ax.set_xlim(0,100)\n", "ax.set_ylim(0.4, 0.8)\n", "\n", "\n", "fig, ax = plt.subplots()\n", "vbn_utils.mean_sem_plot(lp_array[above_thresh], ax, time, color=get_area_color('LP', structure_tree), label='LP')\n", "vbn_utils.mean_sem_plot(visp_array[above_thresh], ax, time, color=get_area_color('VISp', structure_tree), label='VISp')\n", "ax.set_xlim(0,120)\n", "ax.set_ylim(0.4, 0.75)\n", "ax.axhline(0.5, color='k', ls='dotted')\n", "ax.set_xlabel('Time from stimulus (ms)')\n", "ax.set_ylabel('Decoding Accuracy')\n", "ax.legend(frameon=False)\n", "vbn_utils.formatFigure(fig, ax)\n", "print(np.sum(above_thresh), np.sum(above_thresh))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timing of image, change and lick decoding across areas" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "decoding_results_base = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_decoding_from_sensory_action_clusters/with_unitsamp_replacement_and_same_splits_for_sessionunits\"\n", "\n", "basesub = '_basesub'\n", "label_to_dir = {'image': f'pooledImageDecoding{basesub}',\n", " 'lick': f'pooledLickDecoding{basesub}',\n", " 'change': f'pooledChangeDecoding{basesub}',\n", " 'changeprechange': f'pooledChangePrechangeDecoding{basesub}',\n", " 'reaction_time': f'pooledReactionTimeDecoding{basesub}',\n", " 'visual_response': f'pooledVisualResponseDecoding{basesub}',}\n", "\n", "\n", "labels = ['change','image', 'lick', 'reaction_time', 'visual_response']\n", "labels = ['change', 'image', 'lick']\n", "clusters = ['sensory', 'action',]\n", "\n", "areas = ['LGd', 'LP', 'VISall', 'SCMRN', 'Hipp', 'VISp', 'VISl', 'VISal', 'VISrl', 'VISam', 'VISpm',]\n", "\n", "unitSampleSizes = [50,100,]\n", "nPsuedoFlashes = [100,]\n", "nUnitSamples = [1000,]\n", "binsize = 10\n", "\n", "def binsize_filter(filename, binssize):\n", " if binsize==5:\n", " return '5binsize' in filename\n", " else:\n", " return '5binsize' not in filename\n", "\n", "decoding_dict_basesub = {cond: {l: {a: {c: {uss: {npf: {nus: [] for nus in nUnitSamples} for npf in nPsuedoFlashes} \n", " for uss in unitSampleSizes} for c in clusters} \n", " for a in areas} for l in labels} \n", " for cond in ['active', 'passive']}\n", "for cond in ['active', 'passive']:\n", " for label in labels:\n", " dir = label_to_dir[label] + '_' + cond\n", " decoding_results_dir = os.path.join(decoding_results_base, dir)\n", " decoding_results = os.listdir(decoding_results_dir)\n", "\n", " for area in areas:\n", "\n", " area_data_files = [f for f in decoding_results if (f.split('_')[1]==area) and binsize_filter(f, binsize)]\n", " for cluster in clusters:\n", " \n", " area_cluster_data = [f for f in area_data_files if f.split('_')[2]==cluster]\n", " \n", " for uss in unitSampleSizes:\n", " area_cluster_uss = [f for f in area_cluster_data if f.split('_')[3]==str(uss)]\n", "\n", " for npf in nPsuedoFlashes:\n", " area_cluster_uss_npf = [f for f in area_cluster_uss if f.split('_')[4]==str(npf)]\n", "\n", " for nus in nUnitSamples:\n", " area_cluster_uss_npf_nus = [f for f in area_cluster_uss_npf if f.replace('.npy', '').split('_')[5]==str(nus) in f]\n", " if len(area_cluster_uss_npf_nus) == 1:\n", " decoding_dict_basesub[cond][label][area][cluster][uss][npf][nus] = np.load(os.path.join(decoding_results_dir, area_cluster_uss_npf_nus[0]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load facemap data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import h5py\n", "\n", "f = h5py.File(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_video_analysis/facemapData/1044385384_524761_20200819.behavior_proc.hdf5\", 'r')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_video_analysis/facemapDecoding_basesub\"\n", "session_files = os.listdir(data_dir)\n", "\n", "all_session_data = []\n", "for session_file in session_files:\n", " session_data = np.load(os.path.join(data_dir, session_file), allow_pickle=True)\n", " session_data = session_data.item()\n", " all_session_data.append(session_data['balancedAccuracy']['non-change lick'])\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "all_session_data = [a for a in all_session_data if len(a)>0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculate decoding latencies" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# suppress warnings\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from vbn_utils import formatFigure\n", "plt.rcParams['font.size'] = 14\n", "neuron_num = 100\n", "bootstrap_num = 1000\n", "areas = ['LGd', 'LP', 'VISall', 'VISp', 'VISl', 'VISal', 'VISrl', 'VISam', 'VISpm', 'SCMRN', 'Hipp', ]\n", "\n", "latency_summary = {area: {cluster:{} for cluster in ['sensory', 'action']} for area in areas + ['facemap']}\n", "time = np.arange(binsize, 750+binsize, binsize)\n", "latency_threshold = 0.5\n", "for label in ['image', 'change', 'lick',]:\n", "\n", " fig, ax = plt.subplots(1,2)\n", " fig.suptitle(label)\n", " fig.set_size_inches(8,3.8)\n", "\n", " cluster = 'sensory' if label != 'lick' else 'action'\n", " lims = [0,int(100/binsize)] if label != 'lick' else [0, int(750/binsize)]\n", " axlims = np.array(lims)*binsize\n", " chance = 1/8 if label=='image' else 0.5\n", " for i, area in enumerate(areas):\n", " if len(decoding_dict_basesub['active'][label][area][cluster][neuron_num][100][bootstrap_num]) == 0:\n", " continue\n", "\n", " tup, latencies, mean_latency, low_error, hi_error = get_latency2(decoding_dict_basesub['active'][label][area][cluster][neuron_num][100][bootstrap_num], \n", " chance, time, slice(*lims), threshold=latency_threshold, upsample_factor=10, sub='min', norm=True)\n", "\n", " latencies = latencies.astype(float)\n", " \n", " if not area in ['VISp', 'VISl', 'VISal', 'VISrl', 'VISam', 'VISpm']:\n", " plot_decoding_results(decoding_dict_basesub, label, area, cluster, neuron_num, 100, bootstrap_num, time,\n", " bootstrap_iterations=100, condition='active', ax=ax[0], norm=False, norm_slice=slice(*lims),\n", " plotlabel=area, lw=2)\n", " \n", " \n", " print(f'{area} {cluster} {label} latency: {np.nanmean(latencies)}')\n", " \n", " latency_summary[area][cluster][label] = {'mean': mean_latency, 'low': low_error, 'high': hi_error, 'latencies': latencies}\n", " ax[1].errorbar(i, mean_latency, yerr=[[low_error], [hi_error]], color=get_area_color(area, structure_tree), fmt='o')\n", "\n", " if label == 'lick':\n", " frameInterval = 1/60\n", " facemap_time = np.arange(0,0.75+frameInterval/2,frameInterval)*1000\n", " tup, latencies, mean_latency, low_error, hi_error = get_latency2(np.array(all_session_data), 0.5, facemap_time, slice(0,46), threshold=latency_threshold, upsample_factor=17)\n", " latency_summary['facemap']['action']['lick'] = {'mean': mean_latency, 'low': low_error, 'high': hi_error}\n", " ax[1].errorbar(len(areas), mean_latency, yerr=[[low_error], [hi_error]], color='k', fmt='o')\n", " ax1_x_upper_lim = len(areas)+0.5\n", " ax1_y_lims = (50, 250)\n", " print(f'facemap {label} latency: {mean_latency}')\n", " else:\n", " ax1_x_upper_lim = len(areas)-0.5\n", " ax1_y_lims = (28, 80)\n", "\n", "\n", " if label == 'image':\n", " ax0_y_lims = (0.08, 1)\n", " else:\n", " ax0_y_lims = (0.45, 1) \n", " \n", " \n", "\n", " ax[1].set_xlim(-0.5, ax1_x_upper_lim)\n", " ax[1].set_ylim(ax1_y_lims)\n", " ax[1].set_xticks(np.arange(len(areas)))\n", " ax[1].set_xticklabels(areas, rotation=45)\n", " ax[0].set_xlim(*axlims)\n", " ax[0].set_ylim(*ax0_y_lims)\n", " ax[0].axhline(chance, color='k', ls='--', alpha=0.5)\n", "\n", " ax[1].set_ylabel('Time to half max (ms)')\n", " ax[0].set_ylabel('Decoding Accuracy')\n", "\n", " xlabel = label if label != 'lick' else 'change'\n", " formatFigure(fig, ax[0], xLabel=f'Time from {xlabel} onset (ms)')\n", " formatFigure(fig, ax[1])\n", " plt.tight_layout()\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run stats on decoding latencies for VISp/LP change and VIS/SCMRN lick" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generated on HPC by 'run_decoder_area_comparison_monte_carlo.py'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Load monte carlo null distribution for LP-VISp\n", "null_data_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoder_area_comparison_nulls/pooledChangeDecoding_basesub_active\"\n", "null_files = os.listdir(null_data_dir)\n", "\n", "num_permutations = 1000\n", "null_values = []\n", "for perm in range(num_permutations):\n", " v1_null_file = [f for f in null_files if f'pop1_{perm}_' in f]\n", " v1_null_data = np.load(os.path.join(null_data_dir, v1_null_file[0]))\n", " v1_null_latencies = get_latency2(v1_null_data, 0.5, time, slice(*[0,int(100/binsize)]), threshold=0.5, upsample_factor=10, sub='min', norm=True)[1]\n", "\n", " lp_null_file = [f for f in null_files if f'pop2_{perm}_' in f]\n", " lp_null_data = np.load(os.path.join(null_data_dir, lp_null_file[0]))\n", " lp_null_latencies = get_latency2(lp_null_data, 0.5, time, slice(*[0,int(100/binsize)]), threshold=0.5, upsample_factor=10, sub='min', norm=True)[1]\n", "\n", " null_values.append(np.nanmean(v1_null_latencies - lp_null_latencies))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "v1_lats = latency_summary['VISp']['sensory']['change']['latencies']\n", "lp_lats = latency_summary['LP']['sensory']['change']['latencies']\n", "\n", "obs = np.nanmean(v1_lats - lp_lats)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pvalue = (1 + np.sum(np.abs(null_values) >= np.abs(obs))) / (num_permutations+1)\n", "print(pvalue)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(null_values, align='mid')\n", "plt.axvline(obs, color='k')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Load monte carlo null distribution for SCMRN/VISall action neurons (lick decoding)\n", "null_data_dir = \"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoder_area_comparison_nulls/pooledLickDecoding_basesub_active\"\n", "null_files = os.listdir(null_data_dir)\n", "\n", "num_permutations = 1000\n", "null_values = []\n", "for perm in range(num_permutations):\n", " visall_null_file = [f for f in null_files if f'pop1_{perm}_' in f]\n", " visall_null_data = np.load(os.path.join(null_data_dir, visall_null_file[0]))\n", " visall_null_latencies = get_latency2(visall_null_data, 0.5, time, slice(*[0, int(750/binsize)]), threshold=0.5, upsample_factor=10, sub='min', norm=True)[1]\n", "\n", " scmrn_null_file = [f for f in null_files if f'pop2_{perm}_' in f]\n", " scmrn_null_data = np.load(os.path.join(null_data_dir, scmrn_null_file[0]))\n", " scmrn_null_latencies = get_latency2(scmrn_null_data, 0.5, time, slice(*[0, int(750/binsize)]), threshold=0.5, upsample_factor=10, sub='min', norm=True)[1]\n", "\n", " null_values.append(np.nanmean(visall_null_latencies - scmrn_null_latencies))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "visall_lats = latency_summary['VISall']['action']['lick']['latencies']\n", "scmrn_lats = latency_summary['SCMRN']['action']['lick']['latencies']\n", "\n", "obs = np.nanmean(visall_lats - scmrn_lats)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pvalue = (1 + np.sum(np.abs(null_values) >= np.abs(obs))) / (num_permutations+1)\n", "print(pvalue)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['font.size'] = 12\n", "plt.hist(null_values, align='mid', color='gray', label='null distribution')\n", "plt.axvline(obs, color='r', label='observed difference')\n", "plt.xlabel('$\\\\Delta$ lick decoding latency (VIS - SC/MRN) (ms)')\n", "plt.ylabel('Count')\n", "plt.xlim(-65, 65)\n", "plt.legend(frameon=False, loc='upper left')\n", "plt.tight_layout()\n", "vbn_utils.formatFigure(plt.gcf(), plt.gca())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load opto data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from analysis_utils import fitCurve, calc_gompertz\n", "\n", "opto_results = (np.load(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_opto/mouse_opto_fa_rates.npy\"), \n", " np.load(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_opto/mouse_opto_hit_rates.npy\"),\n", " np.load(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_opto/opto_time.npy\"))\n", "opto_hit_rates = opto_results[1]\n", "opto_time = opto_results[2]\n", "opto_time[-1] = 150\n", "fig, ax = plt.subplots()\n", "for rates, color in zip([opto_results[0], opto_results[1]], ['r', 'b']):\n", " sigmoid_fits = [fitCurve(calc_gompertz, opto_time[:-1], rate[:-1], \n", " [rate.max(), 0.1, 75, rate.min()],\n", " bounds=([rate.max()*0.9, -np.inf, -np.inf, -np.inf], [rate.max()*1.1, np.inf, np.inf, np.inf]))\n", " for rate in list(rates) + [rates.mean(axis=0)]]\n", " sigmoids = [calc_gompertz(np.arange(-10, 140), *sigmoid_fit) for sigmoid_fit in sigmoid_fits]\n", " ax.plot(opto_time, rates.mean(axis=0), color + 'o', ms=8)\n", " ax.plot(np.arange(-10, 140), sigmoids[-1], color)\n", "\n", " \n", " mean = np.mean(rates, axis=0)\n", " sem = np.std(rates, axis=0)/(len(rates)**0.5)\n", " ax.errorbar(opto_time, mean, yerr = sem, color=color, linestyle='None')\n", " ax.plot(opto_time[-1], mean[-1], 'o', color=color, mfc='w', ms=8)\n", "\n", "ax.set_xlabel('Laser onset relative to change (ms)')\n", "ax.set_ylabel('Hit rate')\n", "ax.set_xticks(np.arange(0, 200, 50))\n", "ax.set_xticklabels(list(np.arange(0, 200, 50)[:-1]) + ['no opto'])\n", "\n", "opto_latencies = []\n", "for rate, sigmoid in zip(opto_hit_rates, sigmoids):\n", " opto_midpoint = sigmoid.min() + (sigmoid.max() - sigmoid.min())*0.2\n", " opto_latencies.append(np.argmin(np.abs(sigmoid - opto_midpoint)))\n", "\n", "ax.axvline(np.median(opto_latencies), color='b', linestyle='--')\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax_dummy = plt.subplots()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stim_table = pd.read_csv(stim_table_file)\n", "units = pd.read_csv(unit_table_file)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "session_list = [str(s) for s in stim_table['ecephys_session_id'].unique()]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "unit_filter = du.getUnitsInRegion(units, 'VISall', rs=True) & du.apply_unit_quality_filter(units) & du.get_units_in_cluster(units, *np.arange(6))\n", "unit_ids = units.loc[unit_filter]['unit_id'].values\n", "session_list = [str(s) for s in units.loc[unit_filter]['ecephys_session_id'].unique()]\n", "psth_data = vbn_utils.unit_averaged_psth(active_tensor_file, stim_table, session_list, unit_ids, *['engaged', 'is_change', 'lickbout_for_flash_during_response_window'])\n", "passive_psth_data = vbn_utils.unit_averaged_psth(passive_tensor_file, stim_table, session_list, unit_ids, *['engaged', 'is_change', 'lickbout_for_flash_during_response_window'])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "all_psth_mean = np.concatenate(psth_data[0]).mean(axis=0)\n", "passive_all_psth_mean = np.concatenate(passive_psth_data[0]).mean(axis=0)\n", "\n", "all_psth_mean = all_psth_mean - all_psth_mean[:100].min()\n", "passive_all_psth_mean = passive_all_psth_mean - passive_all_psth_mean[:100].min()\n", "\n", "max_val = np.max([all_psth_mean.max(), passive_all_psth_mean.max()])\n", "\n", "all_psth_norm = all_psth_mean/max_val\n", "passive_all_psth_norm = passive_all_psth_mean/max_val\n", "\n", "fig, ax = plt.subplots()\n", "vbn_utils.mean_sem_plot(np.concatenate(psth_data[0]), color='b', ax=ax)\n", "vbn_utils.mean_sem_plot(np.concatenate(passive_psth_data[0]), color='gray', ax=ax)\n", "\n", "plt.xlim(50, 250)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## aggregate timing plots\n", "\n", "fig, ax = plt.subplots()\n", "opto_sigmoid = calc_gompertz(np.arange(-10, 250), *sigmoid_fits[-1])\n", "image_mean = plot_decoding_results(decoding_dict_basesub, 'image', 'LGd', 'sensory', 100, 100, 100, time,\n", " bootstrap_iterations=100, ax=ax_dummy, norm=False, norm_slice=slice(*[0,int(150/binsize)]), return_mean=True)\n", "change_mean = plot_decoding_results(decoding_dict_basesub, 'change', 'VISall', 'sensory', 100, 100, 100, time,\n", " bootstrap_iterations=100, ax=ax_dummy, norm=False, norm_slice=slice(*[0,int(150/binsize)]), return_mean=True)\n", "scmrn_lick_mean = plot_decoding_results(decoding_dict_basesub, 'lick', 'SCMRN', 'action', 100, 100, 100, time,\n", " bootstrap_iterations=100, ax=ax_dummy, norm=False, norm_slice=slice(*[0,int(750/binsize)]), return_mean=True)\n", "\n", "facemap_mean = plot_facemap_decoding(all_session_data, facemap_time, ax_dummy, bootstrap_iterations=100, color='k', return_mean=True)\n", "\n", "\n", "normed = [(a[:lim] - a[:lim].min())/(max(a[:lim]-a[:lim].min())) for a, lim in zip([opto_sigmoid, change_mean, \n", " scmrn_lick_mean, image_mean, facemap_mean], [None, 50, 50, 50, 30])]\n", "ax.plot(np.arange(-10, 250), normed[0], 'b', label='opto hit rates')\n", "ax.plot(time[:25], normed[1][:25], 'r', label='VISall change decoding')\n", "ax.plot(time[:25], normed[2][:25], color=get_area_color('MB', structure_tree), label='SCMRN lick decoding')\n", "ax.plot(time[:25], normed[3][:25], color=get_area_color('LGd', structure_tree), label='LGd image decoding')\n", "ax.plot(facemap_time[:16], normed[4][:16], 'k', label='facemap decoding')\n", "\n", "ax.set_xlim(0, 250)\n", "formatFigure(fig, ax, xLabel='Time from change onset (ms)', yLabel='Normalized')\n", "ax.plot([20, 100], [1.05, 1.05], 'gray', lw=5)\n", "\n", "fig.savefig(os.path.join(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN Manuscript/Fig3 timing_figure\", 'timing_aggregate_test.pdf'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from notebook_utils import set_spine_linewidth" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "areas = ['LGd', 'LP', 'VISall', 'VISp', 'VISl', 'VISal', 'VISrl', 'VISam', 'VISpm', 'SCMRN', ]\n", "area_labels = ['LGd', 'LP', 'VISall', 'VISp', 'VISl', 'VISal', 'VISrl', 'VISam', 'VISpm', 'SCm/MRN', ]\n", "\n", "fig, ax = plt.subplots()\n", "fig.set_size_inches(12,3)\n", "fighist, axhist = plt.subplots(2,1)\n", "for ia, area in enumerate(areas):\n", " for il, label in enumerate(['image', 'change']):\n", " if not label in latency_summary[area]['sensory']:\n", " continue\n", " mean = latency_summary[area]['sensory'][label]['mean']\n", " hi_error = latency_summary[area]['sensory'][label]['high']\n", " low_error = latency_summary[area]['sensory'][label]['low']\n", " if label=='change':\n", " fill_color = 'w'\n", " edge_color = get_area_color(area, structure_tree)\n", " size = 12\n", " else:\n", " fill_color = get_area_color(area, structure_tree)\n", " edge_color = 'none'\n", " size = 12\n", " ax.errorbar(mean, ia, xerr=[[low_error], [hi_error]], color=get_area_color(area, structure_tree), fmt='o', mfc=fill_color, ms=size, mec=edge_color, alpha=0.75, markeredgewidth=2)\n", " axhist[il].hist(latency_summary[area]['sensory'][label]['latencies'], histtype='step')\n", "\n", "ax.set_xlabel('Time from image onset to half max (ms)')\n", "ax.set_yticks(np.arange(len(areas)))\n", "ax.set_yticklabels(area_labels)\n", "\n", "ax.set_xlim(30,80)\n", "ax.set_xlabel('Time from image onset to half max (ms)')\n", "ax.set_yticks(np.arange(len(areas)))\n", "ax.set_yticklabels(area_labels)\n", "ax.set_ylim(-0.6, len(areas))\n", "vbn_utils.formatFigure(fig, ax)\n", "set_spine_linewidth(ax, linewidth=2)\n", "for ia, _ in enumerate(areas):\n", " if not ia%2==0:\n", " ax.axhspan(ia-0.5, ia+0.5, facecolor='lightgray', alpha=0.3, lw=0)\n", "\n", "fig, ax = plt.subplots()\n", "fig.set_size_inches(4,3)\n", "fighist, axhist = plt.subplots()\n", "areas = ['VISall', 'SCMRN']\n", "for ys, area in zip([2,9], areas):\n", " if not 'lick' in latency_summary[area]['action']:\n", " continue\n", " mean = latency_summary[area]['action']['lick']['mean']\n", " hi_error = latency_summary[area]['action']['lick']['high']\n", " low_error = latency_summary[area]['action']['lick']['low']\n", " color = 'k' if area == 'facemap' else get_area_color(area, structure_tree)\n", " ax.errorbar(mean, ys, xerr=[[low_error], [hi_error]], color=color, fmt='>', ms=12, mec='none', markeredgewidth=2)\n", " axhist.hist(latency_summary[area]['action']['lick']['latencies'], histtype='step')\n", "\n", "ax.set_xlim(100,250)\n", "ax.set_ylim(0, 10)\n", "ax.set_ylim(-0.6, 10)\n", "ax.set_yticks([2,9])\n", "ax.set_yticklabels(['VISall', 'SCm/MRN'])\n", "vbn_utils.formatFigure(fig, ax)\n", "ax.yaxis.set_visible(False)\n", "ax.spines['left'].set_visible(False) # Hides the left spine\n", "\n", "set_spine_linewidth(ax, linewidth=2)\n", "\n", "for ia in range(10):\n", " if not ia%2==0:\n", " ax.axhspan(ia-0.5, ia+0.5, facecolor='lightgray', alpha=0.3, lw=0)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_rt2(row):\n", " \n", " trial_id = row['behavior_trial_id']\n", " session_id = row['ecephys_session_id']\n", " first_lick_in_trial = stim_table[(stim_table['behavior_trial_id']==trial_id)&(stim_table['ecephys_session_id']==session_id)]['lick_time'].values\n", " if len(first_lick_in_trial)==0:\n", " rt = np.nan \n", " else:\n", " rt = np.nanmin(first_lick_in_trial) - row['start_time']\n", "\n", " return rt\n", "\n", "lick_time_from_flash2 = stim_table[stim_table['is_change']].apply(get_rt2, axis=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reaction_times = lick_time_from_flash2\n", "fig, ax = plt.subplots()\n", "ax.hist(reaction_times*1000, bins=np.arange(0, 1.01, 0.01)*1000, color='k', linewidth=2, alpha=0.3, ec='none', density=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## comparing decoding across action/sensory clusters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from vbn_utils import formatFigure\n", "plt.rcParams['font.size'] = 16\n", "neuron_num = 100\n", "num_bootstraps = 1000\n", "areas = ['VISall', 'SCMRN',]\n", "\n", "time = np.arange(binsize, 750+binsize, binsize)\n", "latency_threshold = 0.5\n", "for label in ['image', 'change', 'lick',]:\n", "\n", " fig, ax = plt.subplots(1,2)\n", " fig.suptitle(label)\n", " fig.set_size_inches(12,5)\n", " for cluster, ls in zip(['sensory', 'action'], ['-', '--']):\n", "\n", "\n", " lims = [0,int(300/binsize)] if label != 'lick' else [0, int(750/binsize)]\n", " axlims = np.array(lims)*binsize\n", " chance = 1/8 if label=='image' else 0.5\n", " for i, area in enumerate(areas):\n", " if len(decoding_dict_basesub['active'][label][area][cluster][neuron_num][100][num_bootstraps]) == 0:\n", " continue\n", "\n", " tup, latencies, mean_latency, low_error, hi_error = get_latency2(decoding_dict_basesub['active'][label][area][cluster][neuron_num][100][num_bootstraps], \n", " chance, time, slice(*lims), threshold=latency_threshold, upsample_factor=10)\n", " \n", " \n", " plot_decoding_results(decoding_dict_basesub, label, area, cluster, neuron_num, 100, num_bootstraps, time,\n", " bootstrap_iterations=100, condition='active', ax=ax[0], norm=False, norm_slice=slice(*lims),\n", " plotlabel=area, ls=ls)\n", " \n", " \n", " print(f'{area} {cluster} {label} latency: {latencies.mean()}')\n", " \n", " ax[1].errorbar(i, mean_latency, yerr=[[low_error], [hi_error]], color=get_area_color(area, structure_tree), fmt='o')\n", "\n", " if label == 'lick':\n", " frameInterval = 1/60\n", " facemap_time = np.arange(0,0.75+frameInterval/2,frameInterval)*1000\n", " ax[1].errorbar(len(areas), mean_latency, yerr=[[low_error], [hi_error]], color='k', fmt='o')\n", " ax1_x_upper_lim = len(areas)+0.5\n", " ax1_y_lims = (50, 250)\n", " print(f'facemap {label} latency: {mean_latency}')\n", " else:\n", " ax1_x_upper_lim = len(areas)-0.5\n", " ax1_y_lims = (28, 80)\n", "\n", "\n", " if label == 'image':\n", " ax0_y_lims = (0.08, 1)\n", " else:\n", " ax0_y_lims = (0.45, 1) \n", " \n", " \n", " ax[1].set_xlim(-0.5, ax1_x_upper_lim)\n", " ax[1].set_ylim(ax1_y_lims)\n", " ax[1].set_xticks(np.arange(len(areas)))\n", " ax[1].set_xticklabels(areas, rotation=45)\n", " ax[0].set_xlim(*axlims)\n", " ax[0].set_ylim(*ax0_y_lims)\n", "\n", " ax[0].set_ylabel('Decoding Accuracy')\n", "\n", " xlabel = label if label != 'lick' else 'change'\n", " formatFigure(fig, ax[0], xLabel=f'Time from {xlabel} onset (ms)')\n", " formatFigure(fig, ax[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Decoding dropouts and sufficiency tests" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generated on HPC by 'run_decoding_dropouts.py'" ] }, { "cell_type": "code", "execution_count": null, "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, "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, "metadata": {}, "outputs": [], "source": [ "decoding_dir = os.path.join(decoding_results_base, label_to_dir['change'])\n", "subset_files = [f for f in decoding_files if 'subset' in f and 'dropout' in f]\n", "subsets = [f.split('subset_')[1].split(f'_{test}')[0] for f in subset_files]\n", "subsets" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "time = np.arange(binsize, 750+binsize, binsize)\n", "\n", "decoding_dir = os.path.join(decoding_results_base, label_to_dir['change'])\n", "decoding_files = os.listdir(decoding_dir)\n", "subset_files = [f for f in decoding_files if 'subset' in f and 'dropout' in f]\n", "subsets = [f.split('subset_')[1].split(f'_dropout')[0] for f in subset_files]\n", "\n", "decoding_labels = ['change', 'image']\n", "\n", "for unitsubset in subsets:\n", " for decoding_label in decoding_labels:\n", " fig, axes = plt.subplots(1,2)\n", " fig.set_size_inches(12,4)\n", "\n", " for ax, test in zip(axes, ['dropout', 'sufficiency']):\n", "\n", " this_decoding_dir = os.path.join(decoding_results_base, label_to_dir[decoding_label])\n", " these_decoding_files = os.listdir(this_decoding_dir)\n", " subset_file = [f for f in these_decoding_files if f'subset_{unitsubset}' in f and test in f]\n", "\n", " if len(subset_file)==0:\n", " continue\n", "\n", " unitset = subset_file[0].split('set_')[1].split('_sub')[0]\n", " subset_data = np.load(os.path.join(this_decoding_dir, subset_file[0]))\n", "\n", " fig.suptitle(f\"{decoding_label} \\n Subset: {unitsubset}, Full: {unitset}\")\n", "\n", " full_file = [f for f in these_decoding_files if unitset in f and '_full_' in f]\n", "\n", " if len(full_file)==0:\n", " continue\n", "\n", " full_data = np.load(os.path.join(this_decoding_dir, full_file[0]))\n", "\n", " \n", " ax.plot(time, subset_data.mean(axis=0), label= test.capitalize())\n", " ax.plot(time, full_data.mean(axis=0), label='Full')\n", " ax.legend()\n", " ax.set_title(f'{test.capitalize()}')\n", " ax.set_xlim(0, 300)\n", " \n", " plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fam = np.load(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoding_dropouts/pooledImageDecoding_basesub_active/pooledImageDecoding_set_VISall_all_all_sensory_Familiar_full_100_100_100_10binsize.npy\")\n", "nov = np.load(\"/Volumes/programs/mindscope/workgroups/np-behavior/VBN_revision_decoding_dropouts/pooledImageDecoding_basesub_active/pooledImageDecoding_set_VISall_all_all_sensory_Novel_full_100_100_100_10binsize.npy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(fam.mean(axis=0))\n", "plt.plot(nov.mean(axis=0))\n", "plt.xlim(0,15)" ] } ], "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": 2 }