# code for Figure 4 panels # import libs import matplotlib from matplotlib import pyplot as plt from matplotlib.ticker import PercentFormatter import numpy as np import pandas import seaborn as sns from scipy import stats from importlib import reload import pickle import spatint_utils import os # reload modules reload(spatint_utils) # define variables trn_red, lgn_green, _ = spatint_utils.get_colors() spatint_utils.plot_params() class Fig4: """Class to for plotting panels for Fig. 4""" def __init__(self): """Init class""" parentdir = os. path. dirname(os. getcwd()) filename = parentdir + '/data/params_mouse.yaml' # read trn size tuning dataframe self.trn_sztun_df = pandas.read_pickle( filepath_or_buffer= parentdir + '/data/trn_sztun_df.pkl') # read dict for trn_sztun_ex with open(parentdir + '/data/trn_sztun_ex_dict.pkl', 'rb') as f: self.trn_sztun_ex_dict = pickle.load(f) # read trn retinotopy dataframe self.trn_retino_df = pandas.read_pickle( filepath_or_buffer=parentdir + '/data/trn_retino_df.pkl') # get data for trn retinotopy example series self.trn_retino_ex = self.trn_retino_df[(self.trn_retino_df.m == 'BL6_2018_0003') & (self.trn_retino_df.s == 7)] # read trn/lgn rf area dataframe self.rf_area_df = pandas.read_pickle( filepath_or_buffer=parentdir + '/data/rf_area_df.pkl') # read lgn size tuning dataframe self.lgn_sztun_df = pandas.read_pickle( filepath_or_buffer=parentdir + '/data/lgn_sztun_df.pkl') def exrfs(self): """Plot example receptive fields (Fig. 4d,h) Returns ------- axs: list list with axes for example rfs """ # define keys for example rfs exrfs = [ # trn example rfs {'m': 'BL6_2018_0003', 's': 7, 'e': 1, 'u': 25}, {'m': 'BL6_2018_0003', 's': 4, 'e': 1, 'u': 9}, {'m': 'PVCre_2018_0009', 's': 4, 'e': 10, 'u': 10}, {'m': 'BL6_2018_0003', 's': 2, 'e': 1, 'u': 17}, {'m': 'BL6_2018_0003', 's': 3, 'e': 1, 'u': 12}, # dlGN example rfs {'m': 'Ntsr1Cre_2019_0003', 's': 4, 'e': 1, 'u': 7}, {'m': 'Ntsr1Cre_2018_0003', 's': 2, 'e': 1, 'u': 23}, {'m': 'Ntsr1Cre_2018_0003', 's': 2, 'e': 1, 'u': 29}] # init figure f, axs = plt.subplots(1, 8, figsize=(7, 1.5)) for i, (ax, exrf) in enumerate(zip(axs, exrfs)): # get target row in df ex_row = self.rf_area_df[(self.rf_area_df.m == exrf['m']) & (self.rf_area_df.s == exrf['s']) & (self.rf_area_df.e == exrf['e']) & (self.rf_area_df.u == exrf['u'])] if i == len(exrfs) - 1: scalebar = {'dist': 5, 'width': 1, 'length': 20, 'col': 'w'} else: scalebar = {'width': None, 'dist': None, 'length': None, 'col': None} # plot rfs spatint_utils.plot_rf(means=ex_row.mean_fr.values[0], stim_param_values=ex_row.ti_axes.values[0], grat_width=ex_row.grat_width, grat_height=ex_row.grat_height, ax=ax, scalebar=scalebar) # add labels if i == 0: ax.set_ylabel('visTRN') elif i == 5: ax.set_ylabel('dLGN') f = plt.gcf() f.tight_layout() # describe two example rfs self.desc_exrfs(exrfs) return axs def desc_exrfs(self, exrfs): """Describe two example RFs for Fig 4d Parameters ----- exrfs: list example keys """ labels = ['small', 'large'] for i, label in zip(np.array([3, 4]), labels): exrf = exrfs[i] ex_area = self.rf_area_df[(self.rf_area_df.m == exrf['m']) & (self.rf_area_df.s == exrf['s']) & (self.rf_area_df.e == exrf['e']) & (self.rf_area_df.u == exrf['u'])] print('%s area= %0.3f \n' '%s rsq = %0.3f \n' % (label, ex_area.area, label, ex_area.add_rsq)) def retino_ex(self, figsize=(2.5, 2.5), axs=None): """Plot RFs map of example trn recording (Fig. 4e) Parameters ------- figsize: tuple Figure size (width, height) axs: list two axes used for plot and colorbar Returns ------- axs: list two axes with plot and colorbar """ if axs is None: # create figure f = plt.figure(figsize=figsize) axs = [] # axis for plot l = 0.15 b = 0.15 w = 0.6 h = 0.6 axs.append(f.add_axes([l, b, w, h])) # axis for colorbar l = 0.8 w_cbar = 0.05 axs.append(f.add_axes([l, b, w_cbar, h])) # print number of units in series exseries = self.trn_retino_ex print('number of units in example series:', len(exseries)) # define colors vmin = 0.3 vmax = 0.95 colors = np.linspace(vmin, vmax, len(exseries)) mymap = plt.cm.get_cmap("Reds") my_colors = mymap(colors) # iterate over units for urowi, (_, urow) in enumerate(exseries.iterrows()): # compute parameters params = urow.params angle = params[2] * 180 / np.pi fitpars_deg = spatint_utils.degdiff(180, angle, 180) # calculate ellipse x, y = spatint_utils.calculate_ellipse(params[0], params[1], params[4], params[5], fitpars_deg) # plot axs[0].plot(x, y, c=my_colors[urowi]) # layout axs[0].set_xlim(-20, 70) axs[0].set_ylim(-20, 70) axs[0].set_xticks((0, 30, 60)) axs[0].set_yticks((0, 30, 60)) axs[0].spines['bottom'].set_bounds(-10, 70) axs[0].spines['left'].set_bounds(-10, 70) axs[0].set_ylabel('Elevation ($\degree$)') axs[0].set_xlabel('Azimuth ($\degree$)') # draw colorbar sm = matplotlib.colors.LinearSegmentedColormap.from_list('Reds', my_colors) vmin_plot = 0.15 vmax_plot = 1 norm = matplotlib.colors.Normalize(vmin=vmin_plot, vmax=vmax_plot) cbar = matplotlib.colorbar.ColorbarBase(axs[1], cmap=sm, norm=norm) cbar.set_label('Depth (μm)', rotation=270) cbar.ax.invert_yaxis() barunits = (vmax - vmin) / (exseries.depth.iloc[-1] - exseries.depth.iloc[0]) ticklabels = np.array([3050, 3100, 3150]) ticks = ((ticklabels - exseries.depth.iloc[0]) * barunits) + vmin cbar.set_ticks(ticks) cbar.set_ticklabels(ticklabels) f = plt.gcf() f.tight_layout() return axs def trn_retino(self, figsize=(6, 2.5), axs=None): """Plot RFs map of example trn recording (Fig. 4f-g) Parameters ------- figsize: tuple Figure size (width, height) axs: list two axes, one per dimension (azim, elev) Returns ------- axs: list two axes """ if axs is None: # create figure f, axs = plt.subplots(1, 2, figsize=figsize) # get data trn_retino_df = self.trn_retino_df # plot azimuth against depth axs[0].scatter(trn_retino_df.azim, trn_retino_df.depth, facecolors='none', edgecolors='k', linewidth=0.5) # print n n = len(trn_retino_df.azim) print('n azimuth = %d' % n) # plot regression line from ancova model xeval = np.arange(np.min(trn_retino_df.azim), np.max(trn_retino_df.azim), 1) # parameters for the model copied from R intercept = 3044.511585 vis_angle = -1.191659 azim = 60.728882 # for category 1: azim 0: elev interaction = -1.838795 ymodel_azim = intercept + xeval * vis_angle + 1 * azim + 1 * xeval * interaction axs[0].plot(xeval, ymodel_azim, 'r') # layout axs[0].set_xticks((0, 30, 60)) axs[0].set_yticks((2500, 2900, 3300)) axs[0].invert_yaxis() axs[0].set_ylabel('Depth (μm)') axs[0].set_xlabel('Azimuth ($\degree$)') axs[0].spines['bottom'].set_bounds(-20, 70) axs[0].spines['left'].set_bounds(2500, 3300) # plot elevation against depth axs[1].scatter(trn_retino_df.elev, trn_retino_df.depth, facecolors='none', edgecolors='k', linewidth=0.5) # print n n = len(trn_retino_df.elev) print('n elev = %d' % n) # plot regression line from model xeval = np.arange(np.min(trn_retino_df.elev), np.max(trn_retino_df.elev), 1) ymodel_elev = intercept + xeval * vis_angle + 0 * azim + 0 * xeval * interaction axs[1].plot(xeval, ymodel_elev, 'r') # layout axs[1].set_xticks((0, 30, 60)) axs[1].set_yticks((2500, 2900, 3300)) axs[1].invert_yaxis() axs[1].set_ylabel('') axs[1].set_yticklabels([]) axs[1].set_xlabel('Elevation ($\degree$)') axs[1].spines['bottom'].set_bounds(-20, 70) axs[1].spines['left'].set_bounds(2500, 3300) f = plt.gcf() f.tight_layout() return axs def rf_area(self, figsize=(2, 2), ax=None): """Create violin plot for comparison of TRN and LGN RF sizes (Fig. 4i) Parameters ------- figsize: tuple figure size (width, height) ax: instance of matplotlib.axes class axis to use for plotting Returns ------- ax: mpl axis Axis with plot """ if ax is None: # make figure f, ax = plt.subplots(figsize=figsize) # get data rf_area_df = self.rf_area_df # split data trn_area = rf_area_df['area'][rf_area_df['region'] == 'PGN'].array lgn_area = rf_area_df['area'][rf_area_df['region'] == 'LGN'].array # plot data sns.violinplot(data=[np.log(trn_area), np.log(lgn_area)], palette=[trn_red, lgn_green], ax=ax, linewidth=1, inner=None) # plot mean ax.plot([0, 1], [np.log(trn_area.mean()), np.log(lgn_area.mean())], linestyle='', c='k', marker='.') # format plot ylabels = np.array([10, 100, 1000]) ax.set_yticks(np.log(ylabels)) ax.set_yticklabels(ylabels) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ylims = ax.get_ylim() ax.spines['bottom'].set_bounds(0, 1) ax.spines['left'].set_bounds(ylims) plt.gca().get_xticklabels()[0].set_color(trn_red) plt.gca().get_xticklabels()[1].set_color(lgn_green) ax.set_xticklabels(['visTRN', 'dLGN']) ax.set_ylabel('RF area (deg$^2$)') ax.grid(False) # mannwhitneyUtest u_stat, parea = stats.mannwhitneyu(trn_area, lgn_area) # test for differences in variance # (with center = median, levene's test = brown-forsythe test) f_stat, pvar = stats.levene(trn_area, lgn_area, center='median') # ratio ratio = trn_area.mean() / lgn_area.mean() # print stats print('dispersion stats: Brown–Forsythe test\n' 'Fstat = %0.3f \n' 'pval = 10**%0.3f\n\n' 'central tendency stats\n' 'Ustat = %0.3f\n' 'pval area = 10**%0.3f \n' 'N area visTRN = %d \n' 'N area dLGN = %d\n' 'visTRN mean area +- sem = %0.3f (+- %0.3f)\n' 'dLGN mean area +- sem = %0.3f (+- %0.3f)\n' 'visTRN rfs are on average %0.3f x larger than dLGN rfs' % (f_stat, np.log10(pvar), u_stat, np.log10(parea), len(trn_area), len(lgn_area), trn_area.mean(), stats.sem(trn_area), lgn_area.mean(), stats.sem(lgn_area), ratio)) f = plt.gcf() f.tight_layout() return ax def norm_szcurves(self, figsize=(3, 3), ax=None, eval_range=76, thres=128, mark_ex=True, lw=0.5, xticks=(0, 25, 50, 75), colormap='Greys'): """Plot normalized fitted size tuning curves for visTRN population (Fig. 4l) Parameters ------- figsize: tuple Figure size (width, height) ax: instance of matplotlib.axes class Axis to use for plotting eval_range: int Range over which to evaluate model thres: int Lower threshold for darkness of line mark_ex: bool If true plots example neuron in different color lw: float Linewidth xticks: tuple xticks colormap: string Colormap Returns ------- ax: mpl axis Axis with plot """ if ax is None: # create figure f, ax = plt.subplots(figsize=figsize) # define range over which to evaluate model x_eval = range(eval_range) # define colormap cmap = plt.cm.get_cmap(colormap) # plot all tuning curves for row in self.trn_sztun_df.itertuples(): # get y data params = row.tun_pars y = spatint_utils.rog_offset(x_eval, *params) # subtract offset y_sub = y - y[0] # normalize y_norm = y_sub / np.nanmax(y_sub) # define color si = row.si_76 col = int(np.round((1 - si) * 255)) col = np.max((col, thres)) # plot ax.plot(x_eval, y_norm, c=cmap(col), lw=lw) if mark_ex: # plot example session in red params = self.trn_sztun_ex_dict['tun_pars'] y = spatint_utils.rog_offset(x_eval, *params) y_sub = y - y[0] y_norm = y_sub / np.nanmax(y_sub) ax.plot(x_eval, y_norm, c=trn_red, lw=lw) # layout ax.set_xlabel('Diameter ($\degree$)') ax.set_ylabel('Normalized firing rate') ax.set_xticks(xticks) ax.set_yticks((0, 0.5, 1)) ax.spines['bottom'].set_bounds(0, 75) ax.spines['left'].set_bounds(0, 1) return ax def ex_sztun_curve(self, figsize=(4, 2), axs=None): """Plot example visTRN size-tuning curve and raster (Fig. 4jk) Parameters ------- figsize: tuple figure size (width, height) axs: instance of matplotlib.axes class axis to use for plotting Returns ------- ax: mpl axis Axis with plot """ if axs is None: # create figure f, axs = plt.subplots(1, 2, figsize=figsize) # get data ex_data = self.trn_sztun_ex_dict # plot raster spatint_utils.plot_raster(raster=ex_data['rasters'][ex_data['u']], tranges=ex_data['tranges'], opto=ex_data['opto'], ax=axs[0]) # plot curve spatint_utils.plot_tun(means=ex_data['tun_mean'], sems=ex_data['tun_sem'], spons=ex_data['tun_spon_mean'], xs=ex_data['ti_axes'], params=ex_data['tun_pars'], ax=axs[1]) # format layout axs[1].set_xticks((0, 25, 50, 75)) axs[1].set_yticks((0, 10, 20)) axs[1].spines['bottom'].set_bounds(0, 75) f = plt.gcf() f.tight_layout() # print info for example cell sz_ex_si = self.trn_sztun_ex_dict['si_76'] sz_ex_rfcs = self.trn_sztun_ex_dict['rfcs_76'] print('size tuning example cell: \n' 'Preferred size = %0.3f \n' 'SI = %0.3f' % (sz_ex_rfcs, sz_ex_si)) def si(self, figsize=(2, 2), ax=None): """Plot histogram of suppression indices for visTRN and dLGN population (Fig. 4m) Parameters ------- figsize: tuple figure size (width, height) ax: instance of matplotlib.axes class axis to use for plotting Returns ------- ax: mpl axis Axis with plot """ if ax is None: # create figure f, ax = plt.subplots(figsize=figsize) # get data trn_sis = self.trn_sztun_df.si_76.to_numpy() lgn_sis = self.lgn_sztun_df.si_76.str[0].to_numpy() # assert that no data is beyond limits assert (~np.any(trn_sis < 0)) & (~np.any(trn_sis > 1)), 'Datapoints beyond bounds' assert (~np.any(lgn_sis < 0)) & (~np.any(lgn_sis > 1)), 'Datapoints beyond bounds' # plot si_bins = np.arange(0, 1.05, 0.05) # define bins trn_red_trsp = (*trn_red[0:3], 0.5) # define color for trn lgn_green_trsp = (*lgn_green[0:3], 0.5) # define color for dlgn ax.hist(trn_sis, bins=si_bins, weights=np.ones(len(trn_sis)) / len(trn_sis), lw=0, fc=trn_red_trsp) ax.hist(lgn_sis, bins=si_bins, weights=np.ones(len(lgn_sis)) / len(lgn_sis), lw=0, fc=lgn_green_trsp) # layout ax.yaxis.set_major_formatter(PercentFormatter(1)) ax.set_xlabel('Suppression index') ax.set_ylabel('neurons') ax.set_xticks((0, 0.5, 1)) ax.set_yticks((0, 0.25, 0.5)) ax.spines['bottom'].set_bounds(0, 1) ax.set_xlim((0, 1)) f = plt.gcf() f.tight_layout() # compute and print stats for trn n_trn = len(trn_sis) si_mean_trn = np.mean(trn_sis) si_sem_trn = stats.sem(trn_sis) si_med_trn = np.median(trn_sis) lower05_trn = (len(trn_sis[trn_sis < 0.05]) / len(trn_sis) * 100) # percentage of cells with si smaller 0.05 print('trn size tuning population: \n' 'n = %d\n' 'mean si +/- sem = %0.3f (+- %0.3f) \n' 'median si = %0.3f \n' '%0.3f percent of pgn cells have si < 0.05\n' % (n_trn, si_mean_trn, si_sem_trn, si_med_trn, lower05_trn)) # compute and print stats for dlgn n_lgn = len(lgn_sis) si_mean_lgn = np.mean(lgn_sis) si_sem_lgn = stats.sem(lgn_sis) si_med_lgn = np.median(lgn_sis) print('lgn size tuning population:\n' 'n = %d\n' 'mean si +/- sem = %0.3f (+- %0.3f)\n' 'median si = %0.3f\n' % (n_lgn, si_mean_lgn, si_sem_lgn, si_med_lgn)) # compare the two ustat_sis, p_sis = stats.mannwhitneyu(trn_sis, lgn_sis) print('mannwhitneyu test to compare si in dlgn and trn:\n' 'Ustat: %0.3f\n' 'pvalue: 10**%0.3f\n' % (ustat_sis, np.log10(p_sis))) return ax