123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395 |
- import numpy as np
- import matplotlib.pyplot as plt
- import matplotlib
- import pandas as pd
- idx = pd.IndexSlice
- import pickle as pkl
- import os
- import spatint_utils as su
- datapath = os.getcwd()+'/../data/'
- with open(datapath+'lgn_spat_profile.pkl','rb') as r:
- lgn_spat_profile = pkl.load(r)
-
- with open(datapath+'lgn_ori_examples.pkl','rb') as r:
- lgn_ori_examples = pkl.load(r)
-
- with open(datapath + 'v1_raw_muaRFs.pkl','rb') as r:
- v1_raw_muaRFs = pkl.load(r)
-
- with open(datapath + 'v1_muaRFs.pkl','rb') as r:
- v1_muaRFs = pkl.load(r)
-
- with open(datapath + 'lgn_raw_muaRFs.pkl','rb') as r:
- lgn_raw_muaRFs = pkl.load(r)
-
- with open(datapath + 'lgn_muaRFs.pkl','rb') as r:
- lgn_muaRFs = pkl.load(r)
- def agg_overlap_bins(data, start, end, nbins, binsize, agg_op, boot_value, filtr = False):
- """computes aggregate values for bins along a given axis, aggregate operation is defined by agg_op"""
-
- bins_data = []
- print('Averaging and bootstrapping in overlapping bins')
- for bin_start in np.linspace(start, end, nbins):
- bin_data = data.loc[(data['d']>bin_start) & (data['d']<(bin_start+binsize))]
-
- if filtr != False and len(bin_data)<filtr:
- continue
- booties = []
- for i in range(1000):
- bt_bin = bin_data[boot_value][np.random.randint(0,len(bin_data),len(bin_data))]
- booties.append(np.mean(bt_bin))
- booties = np.sort(booties)
- lower = booties[24]
- upper = booties[974]
- boot_std = np.std(booties)
- bin_agg = bin_data.agg(agg_op)
- bin_agg['lower'] = lower
- bin_agg['upper'] = upper
- bin_agg['boot'] = boot_std
- bin_agg['d'] = bin_start
- bins_data.append(bin_agg)
- return pd.concat(bins_data, axis =1)
- #bin_mean = agg_overlap_bins(lgn_spat_profile,0,100,31,15,agg_op='mean',boot_value='gain',filtr=5)
- def plot_fold_change(ax=None, data=lgn_spat_profile):#, bin_mean=bin_mean):
- """Plot mean fold change values of all units as scatter plot"""
-
- if ax is None:
- fig, ax = plt.subplots()
-
- # Clip low and high values for better visibility
- data.loc[data['gain']>2,'gain'] = 2.35
- data.loc[data['gain']<-2,'gain'] = -2.35
-
- ax.plot(data['d'],data['gain'],c='k', alpha=0.7, linestyle='', fillstyle='none', clip_on = False, marker='o',mec='k',markersize=1)
- examples = data.reindex([('Ntsr1Cre_2015_0080',3,4027),
- ('Ntsr1Cre_2018_0003',2,22,3),
- ('Ntsr1Cre_2015_0080',4,3049)])
- ax.plot([-100,100],[0,0],c='k', lw=0.35)
- ax.plot(bin_mean.loc['d'][:15]+7.5,bin_mean.loc['gain'][:15], lw=3,alpha=0.5, c='#1090cfff')
- ax.plot(bin_mean.loc["d"][8:15]+7.5,bin_mean.loc['gain'][8:15], lw=3,alpha=1, c ='#1090cfff')
- ax.set(xlim=[0,60], ylim=[-2.5,2.5])
- ax.set_xlabel('Distance of dLGN RFs to V1 RFs at injection site $(^{\circ})$',labelpad=5)
- ax.set_xticks([0,20,40,60])
- ax.set_ylabel('Fold change',labelpad=0)
- ax.minorticks_off()
- ax.set_yticks(list(np.linspace(-2,2,5))+[-2.35,2.35])
- ax.set_yticklabels([str("{0:.2f}".format(round(ytick,2)))for ytick in np.geomspace(0.25,4,5)]+['< 0.25','> 4.00'])
- ax.spines['left'].set_bounds(-2,2)
- ax.spines['right'].set_visible(False)
- ax.spines['top'].set_visible(False)
-
- def plot_modulation_histogram(ax = None, data = lgn_spat_profile):
-
- if ax is None:
- fig = plt.figure()
- fig.set_figheight(2.576)
- fig.set_figwidth(4.342)
- ax = fig.add_axes((0.15,0.25,0.75,0.75))
-
-
- bins = np.linspace(0,60,13)
- data = data.sort_values('d')
- data['bin'] = data.apply(lambda x: np.searchsorted(bins, x['d'], side='right')-1, axis=1)
- ratios = data.groupby('bin').apply(lambda x: pd.DataFrame({'sup' :[len(x[x['modlab']==-1])/len(x)],
- 'fac' :[len(x[x['modlab']== 1])/len(x)],
- 'snull':[len(x[(x['modlab']==0) & (x['gain']<0)])/len(x)],
- 'fnull' :[len(x[(x['modlab']==0) & (x['gain']>0)])/len(x)]}))
- ratios.index = ratios.index.set_names('dummy', level=1)
- ratios = ratios.reset_index("dummy").drop(columns='dummy')
- ratios['binpos'] = bins[ratios.index]+2.5
- ax.bar(ratios['binpos'],-ratios['sup'], bottom = -ratios['snull'], color = 'blue', width=4)
- ax.bar(ratios['binpos'], ratios['fac'], bottom = ratios['fnull'], color = 'orange',width = 4)
- ax.bar(ratios['binpos'], -ratios['snull'], bottom = 0, color='lightskyblue', width=4)
- ax.bar(ratios['binpos'], ratios['fnull'], bottom = 0, color = 'moccasin', width=4)
- ax.set_xlabel('Distance of dLGN RF to \n mean V1 RF at injection site ($\degree$)')
- ax.spines['top'].set_visible(False)
- ax.spines['right'].set_visible(False)
- ax.spines['left'].set_bounds(high=0.75,low=-0.75)
- ax.plot([0,70],[0,0],c='k',lw=0.35)
- ax.set(xlim=[2.5,60],ylim=[-0.9,0.9])
- ax.set_yticks([-0.75,-0.25,0,0.25,0.75])
- ax.set_yticklabels(['75','25','','25','75'])
- ax.set_xticks([0,20,40,60])
- ax.set_ylabel('Proportion (%)')
- ax.yaxis.set_label_coords(-0.1,0.55)
- ax.annotate('enhanced', xy=(1,0.15), xytext=(1,0.15),xycoords='data',rotation='vertical',color='darkorange')
- ax.annotate('suppressed', xy=(1,-0.8), xytext=(1,-0.8),xycoords='data',rotation='vertical',color='blue')
-
- def plot_ori_examples(ax = None, data = lgn_ori_examples):
-
- def compute_ori_curve(df,opto,xs):
-
- return su.sum_of_gaussians(xs,*df['tun_pars'][opto])
-
- xs = np.linspace(0,360,361)
- ctrl_curves = data.apply(compute_ori_curve, axis=1, opto=0,xs=xs)
- opto_curves = data.apply(compute_ori_curve, axis=1, opto=1,xs=xs)
-
- if ax is None:
- gridspec_kw = {'left':0.175,'right':0.95,'bottom':0.225,'top':0.95,'hspace':0.1}
- fig, axes = plt.subplots(1,3,figsize=(6.327,1.971), gridspec_kw = gridspec_kw)
-
- for i,ax in enumerate(axes):
- ax.plot(xs, ctrl_curves.iloc[i],c='k',ms=2)
- ax.plot(xs, opto_curves.iloc[i],c='#1090cfff',ms=2)
- ax.plot(np.linspace(0,330,12),data.iloc[i]['tun_mean'][:,0],mfc='k',ms=8,ls='',marker='.',mew=0)
- ax.plot(np.linspace(0,330,12),data.iloc[i]['tun_mean'][:,1],mfc='#1090cfff',ms=8,ls='',marker='.',mew=0)
- ax.errorbar(np.linspace(0,330,12), data.iloc[i]['tun_mean'][:,0], yerr=data.iloc[i]['tun_sem'][:,0],fmt='none', ecolor='k', elinewidth=1.5)
- ax.errorbar(np.linspace(0,330,12), data.iloc[i]['tun_mean'][:,1], yerr=data.iloc[i]['tun_sem'][:,1],fmt='none', ecolor='#1090cfff', elinewidth=1.5)
- ax.plot(np.linspace(0,330,12), 12*[data.iloc[i]['tun_spon_mean'][0]],c='k',lw=0.5)
- ax.plot(np.linspace(0,330,12), 12*[data.iloc[i]['tun_spon_mean'][1]],c='#1090cfff',lw=0.5)
- ax.spines['top'].set_visible(False)
- ax.spines['right'].set_visible(False)
- ax.spines['bottom'].set_visible(True)
- ax.spines['bottom'].set_bounds(0,360)
- xticks = [0,180,360]
- ax.set_xticks(xticks)
- ax.set_yticks([])
- ax.set_xticklabels([])
-
- axes[0].spines['left'].set_visible(True)
- axes[0].spines['bottom'].set_visible(True)
- axes[0].set(ylim=[0,100])
- axes[0].spines['left'].set_bounds(0,80)
- axes[0].set_yticks([0,40,80])
-
- axes[1].set(ylim=[0,62.5])
- axes[1].spines['left'].set_bounds(0,50)
- axes[1].set_yticks([0,25,50])
- axes[2].set(ylim=[0,62.5])
- axes[2].spines['left'].set_bounds(0,50)
- axes[2].set_yticks([0,25,50])
- axes[0].set_ylabel('Firing rate (sp/s)')
- axes[0].yaxis.set_label_coords(-0.33,0.4)
- axes[0].set_xlabel('Direction ($\degree$)',labelpad=0)
- axes[0].set_xticks([0,180,360])
- axes[0].set_xticklabels(['0','180','360'])
- axes[0].xaxis.set_tick_params(length=1)
-
-
- def plot_v1_raw_muaRFs(ax = None, data = v1_raw_muaRFs):
-
- if ax is None:
- gridspec_kw = {'left':0.1,'right':0.975,'bottom':0.15,'top':0.9}
- fig,ax = plt.subplots(3,1,gridspec_kw=gridspec_kw)
- fig.set_figheight(1.984)
- fig.set_figwidth(1.332)
- cmap = matplotlib.cm.get_cmap('Greens')
- cmap = matplotlib.colors.LinearSegmentedColormap.from_list('custom',
- [(0,'white'),
- (1,'green')])
- for axi, (index, row) in zip(ax,data.iterrows()):
- left = row['ti_axes'][:,0].min()-2.4
- right = row['ti_axes'][:,0].max()+2.5
- bottom = row['ti_axes'][:,1].min()-2.5
- top = row['ti_axes'][:,1].max()+2.5
- extent = (left,right,bottom,top)
- axi.imshow(np.rot90(row['rfs'].mean(axis=(0,3))),cmap='gray', extent = extent)
- plot_ellipse(row,axi,cmap, linewidth=1.5)
- axi.axis('scaled')
- axi.set(xlim=[left,right],ylim=[bottom,top])
- axi.set_xticks([])
- axi.set_yticks([])
- else:
- xticks = [left,right]
- yticks = [bottom,top]
- axi.set_xticks(xticks)
- axi.set_yticks([])
- ax[0].set_yticks(yticks)
- axi.set_xticklabels([str("{0:.0f}".format(round(xtick,0))) +'$\degree$' for xtick in xticks])
- ax[0].set_yticklabels([str("{0:.0f}".format(round(ytick,0))) +'$\degree$' for ytick in yticks],
- rotation=90)
- plt.setp(ax[0].yaxis.get_majorticklabels(), va='center')
- plt.setp(axi.xaxis.get_majorticklabels(), ha='center')
- axi.tick_params(length=2)
- ax[0].tick_params(axis='y', pad=0, length=2)
- axi.tick_params(axis='x', pad=1)
-
- def plot_v1_rf_ellipses(ax = None, data = v1_muaRFs):
-
-
- if ax is None:
- fig = plt.figure(figsize=(7,5))
- ax = fig.add_axes((0.15,0.15,0.7,0.7))
-
- cmap = matplotlib.cm.get_cmap('Greens')
- data = data.reset_index('ch')
-
- one_series = data.loc['Ntsr1Cre_2015_0080',2]
- one_series = one_series[one_series['sigma_x_mix']<15]
- maxchan,minchan = one_series['ch'].max(), one_series['ch'].min()
- V1_x = data.groupby('m').agg('mean').loc['Ntsr1Cre_2015_0080']['x_mix']
- V1_y = data.groupby('m').agg('mean').loc['Ntsr1Cre_2015_0080']['y_mix']
- one_series.apply(lambda x: plot_ellipse(x,ax=ax,cmap=cmap,linewidth=1,maxchan=maxchan),axis=1)
- ax.scatter(V1_x,V1_y,marker='+',s=100,lw=2,c='k',zorder=3)
- ax.set(xlim=[-20,80],ylim=[-35,50])
- ax.set_xlabel('Azimuth ($\degree$)'
- ,labelpad=0)
- ax.set_ylabel('Elevation ($\degree$)'
- ,labelpad=0)
- xticks = [tick.get_text() for tick in ax.get_xticklabels()]
- ax.spines['top'].set_visible(False)
- ax.spines['right'].set_visible(False)
-
- def plot_lgn_raw_muaRFs(ax = None, data = lgn_raw_muaRFs):
-
- if ax is None:
- gridspec_kw = {'left':0.1,'right':0.975,'bottom':0.1,'top':0.93}
- fig,ax = plt.subplots(5,2,gridspec_kw=gridspec_kw)
- fig.set_figheight(3.552)
- fig.set_figwidth(2.190)
- ax_column1 = ax[:,0]
- ax_column2 = ax[:,1]
-
- cmap = matplotlib.colors.LinearSegmentedColormap.from_list('custom',
- [(0,'xkcd:off white'),
- (1,'xkcd:neon purple')])
- cmap2 = matplotlib.colors.LinearSegmentedColormap.from_list('custom',
- [(0,'xkcd:off white'),
- (1,'xkcd:saffron')]) #deep tea
- cols = {5:'#00555a',9:'blue'}
- ser1 = data.loc['Ntsr1Cre_2015_0080',9].iloc[14:19]
- ser2 = data.loc['Ntsr1Cre_2015_0080',5].iloc[7:12]
- for axi, (index, row) in zip(ax_column1,ser1.iterrows()):
- left = row['ti_axes'][:,0].min()-2.4
- right = row['ti_axes'][:,0].max()+2.5
- bottom = row['ti_axes'][:,1].min()-2.5
- top = row['ti_axes'][:,1].max()+2.5
- extent = (left,right,bottom,top)
- axi.imshow(np.rot90(row['rfs'].mean(axis=(0,3))),cmap='gray', extent = extent)
- plot_ellipse(row,axi,cmap, linewidth =1.5)
- axi.axis('scaled')
- axi.set(xlim=[left,right],ylim=[bottom,top])
- axi.set_xticks([])
- axi.set_yticks([])
- else:
- xticks = [left,right]
- yticks = [bottom,top]
- axi.set_xticks(xticks)
- axi.set_yticks([])
- ax_column1[0].set_yticks(yticks)
- axi.set_xticklabels([str("{0:.0f}".format(round(xtick,0))) +'$\degree$' for xtick in xticks])
- ax_column1[0].set_yticklabels([str("{0:.0f}".format(round(ytick,0))) +'$\degree$' for ytick in yticks],
- rotation=90)
- plt.setp(ax_column1[0].yaxis.get_majorticklabels(), va='center')
- plt.setp(axi.xaxis.get_majorticklabels(), ha='center')
- ax_column1[0].tick_params(axis='y', pad=0)
- for axi, (index, row) in zip(ax_column2,ser2.iterrows()):
- left = row['ti_axes'][:,0].min()-2.4
- right = row['ti_axes'][:,0].max()+2.5
- bottom = row['ti_axes'][:,1].min()-2.5
- top = row['ti_axes'][:,1].max()+2.5
- extent = (left,right,bottom,top)
- axi.imshow(np.rot90(row['rfs'].mean(axis=(0,3))),cmap='gray', extent = extent)
- plot_ellipse(row,axi,cmap2, linewidth=1.5)
- axi.axis('scaled')
- axi.set(xlim=[left,right],ylim=[bottom,top])
- axi.set_xticks([])
- axi.set_yticks([])
- else:
- xticks = [left,right]
- yticks = [bottom,top]
- axi.set_xticks(xticks)
- axi.set_yticks([])
- ax_column2[0].set_yticks(yticks)
- axi.set_xticklabels([str("{0:.0f}".format(round(xtick,0))) +'$\degree$' for xtick in xticks]
- )
- ax_column2[0].set_yticklabels([str("{0:.0f}".format(round(ytick,0))) +'$\degree$' for ytick in yticks],
- rotation=90)
- plt.setp(ax_column2[0].yaxis.get_majorticklabels(), va='center')
- plt.setp(axi.xaxis.get_majorticklabels(), ha='center')
- ax_column2[0].tick_params(axis='y', pad=0)
-
- def plot_lgn_rf_ellipses(ax = None, data = lgn_muaRFs, v1data=v1_muaRFs):
-
- if ax is None:
- fig = plt.figure(figsize=(2.995,3.531))
- ax = fig.add_axes((0.275,0.2,0.65,0.725))
- cmap1 = matplotlib.colors.LinearSegmentedColormap.from_list('custom',
- [(0,'xkcd:neon purple'),
- (1,'xkcd:neon purple')])
- cmap2 = matplotlib.colors.LinearSegmentedColormap.from_list('custom',
- [(0,'xkcd:saffron'),
- (1,'xkcd:saffron')])
- cmaps = {9:cmap1,5:cmap2}
- two_series = data.loc[idx['Ntsr1Cre_2015_0080',(5,9),:]]
- V1_x = v1data.groupby('m').agg('mean').loc['Ntsr1Cre_2015_0080']['x_mix']
- V1_y = v1data.groupby('m').agg('mean').loc['Ntsr1Cre_2015_0080']['y_mix']
- two_series.groupby('s').apply(lambda x: plot_RF_ellipses(x, ax, cmaps[x.index.unique(level=1)[0]],linewidth=0.5))
- ax.scatter(V1_x,V1_y,marker='+',s=100,lw=2,c='k',zorder=3)
- ax.set(xlim=[-5,70],ylim=[-25,60])
- ax.set_xlabel('Azimuth ($\degree$)'
- ,labelpad=0)
- ax.set_ylabel('Elevation ($\degree$)'
- ,labelpad=0)
- xticks = [tick.get_text() for tick in ax.get_xticklabels()]
- xticks = [0,40,80]
- yticks=[-20,0,30,60]
- xticklabels = [str(xtick) for xtick in xticks]
- yticklabels = [str(ytick) for ytick in yticks]
- ax.yaxis.set_tick_params(pad=0)
- ax.set_xticks(xticks)
- ax.set_xticklabels(xticklabels)
- ax.set_yticks(yticks)
- ax.set_yticklabels(yticklabels)
- ax.spines['top'].set_visible(False)
- ax.spines['right'].set_visible(False)
- ax.spines['left'].set_bounds(-20,60)
- ax.spines['bottom'].set_bounds(0,80)
- norm= matplotlib.colors.Normalize(vmin=0.5,vmax=1)
-
- def plot_RF_ellipses(data, ax, cmap,linewidth=1):
- """Plotting RF ellipses from fitted parameters."""
- data.apply(lambda x: plot_ellipse(x,ax,cmap,linewidth=linewidth),axis=1)
-
- def plot_ellipse(data, ax, cmap, linewidth=1, maxchan=60):
- """Plot from single df row"""
- ddiff = su.degdiff(180,(data['theta_mix']*180/np.pi),180)
- params = data[['x_mix',
- 'y_mix',
- 'sigma_x_mix',
- 'sigma_y_mix']]
- x_ellipse, y_ellipse = su.calculate_ellipse(*params,ddiff)
- col = cmap(data['rsq_mix'])
- ax.plot(x_ellipse,y_ellipse, lw=linewidth,c=col)
- ax.axis('scaled')
|