From 6c23861ce3f3458fcae43953f33c27fe5950702b Mon Sep 17 00:00:00 2001 From: ycanerol Date: Mon, 12 Jun 2017 18:53:34 +0200 Subject: [PATCH] Add scripts to analyze mouse data from Dimos There are two differences between data from Fernando and Dimos; - Order of stimuli are different so stimulus type 2 and 3 refer to different stimuli in these experiments - Format of the .mat files are v7.3 (opens with h5py) and v7 (opens with scipy.io.loadmat() ) so they need different handling --- checkerflicker_analysis_dimos.py | 127 +++++++++++++++++++++++++++++ fff_analysis.py | 4 +- fff_analysis_dimos.py | 121 ++++++++++++++++++++++++++++ onoffplotter_dimos.py | 108 +++++++++++++++++++++++++ plotall.py | 10 ++- plotall_dimos.py | 133 +++++++++++++++++++++++++++++++ 6 files changed, 497 insertions(+), 6 deletions(-) create mode 100644 checkerflicker_analysis_dimos.py create mode 100644 fff_analysis_dimos.py create mode 100644 onoffplotter_dimos.py create mode 100644 plotall_dimos.py diff --git a/checkerflicker_analysis_dimos.py b/checkerflicker_analysis_dimos.py new file mode 100644 index 0000000..0d787ac --- /dev/null +++ b/checkerflicker_analysis_dimos.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jun 12 17:59:27 2017 + +@author: ycan +Analyze checkerflicker stimulus from Dimos +""" + +import sys +import numpy as np +from collections import Counter +import matplotlib.pyplot as plt +import scipy.io +# Custom packages +try: + import lnp_checkerflicker as lnpc + import lnp +except: + sys.path.append('/Users/ycan/Documents/official/gottingen/lab rotations/LR3 Gollisch/scripts/') + import lnp_checkerflicker as lnpc + import lnp + +main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ +LR3 Gollisch/data/Experiments/Mouse/2017_01_17/' + +cells_classified = main_dir+'sorting_info.mat' + +frames_path = 'frametimes/2_checkerflicker10x10bw2blinks_frametimings.mat' +stimulus_type = frames_path.split('/')[-1].split('_')[0] +stimulus_name = frames_path.split('/')[-1] + +sx = 60 +sy = 80 + +f = scipy.io.loadmat(cells_classified) +a = f['clusters'].T[0] +b = f['clusters'].T[1] +files = [] +for i in range(len(a)): + files.append('{}{:02.0f}'.format(int(a[i]), int(b[i]))) + +#files=['101'] + +first_run_flag = True + +for filename in files: + if first_run_flag: + f = scipy.io.loadmat(main_dir+frames_path) + ftimes = (np.array(f.get('ftimes'))/1000) + ftimes = ftimes.reshape((ftimes.size,)) + dt = np.average(np.ediff1d(ftimes)) + # Average difference between two frames in miliseconds + + total_frames = ftimes.shape[0] + ftimes = ftimes[:total_frames] + filter_length = 20 # Specified in nr of frames + + stimulus = np.load('/Users/ycan/Documents/official/gottingen/lab \ +rotations/LR3 Gollisch/data/checkerflickerstimulus.npy')[:, :, :total_frames] + + first_run_flag = False + + spike_path = main_dir+'rasters/'+str(stimulus_type)+'_SP_C'+filename+'.txt' + save_path = main_dir+'analyzed/'+str(stimulus_type)+'_SP_C'+filename + + spike_file = open(spike_path) + spike_times = np.array([float(line) for line in spike_file]) + spike_file.close() + + spike_counts = Counter(np.digitize(spike_times, ftimes)) + spikes = np.array([spike_counts[i] for i in range(total_frames)]) + + total_spikes = np.sum(spikes) + if total_spikes < 2: + continue + +# %% + sta_unscaled, max_i, temporal = lnpc.sta(spikes, + stimulus, + filter_length, + total_frames) + max_i = lnpc.check_max_i(sta_unscaled, max_i) + + stim_gaus = lnpc.stim_weighted(sta_unscaled, max_i, stimulus) + + sta_weighted, _ = lnp.sta(spikes, stim_gaus, filter_length, total_frames) + + w, v, _, _, _ = lnpc.stc(spikes, stim_gaus, + filter_length, total_frames, dt) + + bins = [] + spike_counts_in_bins = [] + for i in [sta_weighted, v[:, 0]]: + a, b = lnpc.nlt_recovery(spikes, stim_gaus, i, 60, dt) + bins.append(a) + spike_counts_in_bins.append(b) + + sta_weighted, bins[0], \ + spike_counts_in_bins[0], \ + _, _ = lnpc.onoffindex(sta_weighted, bins[0], + spike_counts_in_bins[0]) + + v[:, 0], bins[1],\ + spike_counts_in_bins[1],\ + peak, onoffindex = lnpc.onoffindex(v[:, 0], bins[1], + spike_counts_in_bins[1]) + +#%% + np.savez(save_path, + sta_unscaled=sta_unscaled, + sta_weighted=sta_weighted, + stimulus_type=stimulus_type, + total_frames=total_frames, + temporal=temporal, + v=v, + w=w, + max_i=max_i, + spike_path=spike_path, + filename=filename, + bins=bins, + spike_counts_in_bins=spike_counts_in_bins, + onoffindex=onoffindex, + total_spikes=total_spikes, + peak=peak, + stimulus_name=stimulus_name + ) diff --git a/fff_analysis.py b/fff_analysis.py index b72728f..6280e67 100644 --- a/fff_analysis.py +++ b/fff_analysis.py @@ -23,7 +23,7 @@ import lnp_checkerflicker as lnpc main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ -LR3 Gollisch/data/Experiments/Salamander/2014_02_25/' +LR3 Gollisch/data/Experiments/Salamander/2014_01_27/' stimulus_type = 2 # Change stimulus type: # full field flicker is 2 @@ -47,7 +47,7 @@ files.append('{}{:02.0f}'.format(a, int(b))) f.close() -#files = ['101'] +#files = ['502'] first_run_flag = True diff --git a/fff_analysis_dimos.py b/fff_analysis_dimos.py new file mode 100644 index 0000000..e2a2be6 --- /dev/null +++ b/fff_analysis_dimos.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jun 12 16:21:33 2017 + +@author: ycan + +Full field flicker analysis for data from Dimos + +""" + +import sys +import numpy as np +from collections import Counter +import matplotlib.pyplot as plt +import scipy.io +try: + import lnp + import lnp_checkerflicker as lnpc +except: + sys.path.append('/Users/ycan/Documents/official/gottingen/lab rotations/LR3 Gollisch/scripts/') + import lnp + import lnp_checkerflicker as lnpc + +main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ +LR3 Gollisch/data/Experiments/Mouse/2017_01_17/' + +cells_classified = main_dir+'sorting_info.mat' +#%% +stimulus_path = '/Users/ycan/Documents/official/gottingen/\ +lab rotations/LR3 Gollisch/data/fff2h.npy' +frames_path = 'frametimes/3_fff_gauss_2blinks_frametimings.mat' +stimulus_type = frames_path.split('/')[-1].split('_')[0] +stimulus_name = frames_path.split('/')[-1] + +f = scipy.io.loadmat(cells_classified) +a = f['clusters'].T[0] +b = f['clusters'].T[1] +files = [] +for i in range(len(a)): + files.append('{}{:02.0f}'.format(int(a[i]), int(b[i]))) + +#files = ['502'] +# %% +first_run_flag = True + +for filename in files: + if first_run_flag: + f = scipy.io.loadmat(main_dir+frames_path) + ftimes = (np.array(f.get('ftimes')/1000)) + ftimes = ftimes.reshape((ftimes.size,)) + dt = np.average(np.ediff1d(ftimes)) + # Average difference between two frames in miliseconds + + total_frames = ftimes.shape[0] + filter_length = 20 # Specified in nr of frames + + stimulus = np.load(stimulus_path)[:total_frames] + first_run_flag = False + + spike_path = main_dir+'rasters/'+str(stimulus_type)+'_SP_C'+filename+'.txt' + save_path = main_dir+'analyzed/'+str(stimulus_type)+'_SP_C'+filename + + spike_file = open(spike_path) + spike_times = np.array([float(line) for line in spike_file]) + spike_file.close() + + spike_counts = Counter(np.digitize(spike_times, ftimes)) + spikes = np.array([spike_counts[i] for i in range(total_frames)]) + + total_spikes = np.sum(spikes) + if total_spikes < 2: + continue + + # Bin spikes + + # Start the analysis + sta = lnp.sta(spikes, stimulus, filter_length, total_frames)[0] + # Use scaled STA + + generator = np.convolve(sta, stimulus, + mode='full')[:-filter_length+1] + + bins_sta, spikecount_sta = lnp.q_nlt_recovery(spikes, generator, + 60, dt) + + eigen_indices = [0] + bin_nr = 60 + + w, v, bins_stc, spikecount_stc, eigen_legends = lnp.stc(spikes, stimulus, + filter_length, + total_frames, dt, + eigen_indices, + bin_nr) + sta, bins_sta, \ + spikecount_sta, _, _ = lnpc.onoffindex(sta, bins_sta, + spikecount_sta) + + v[:, 0], bins_stc,\ + spikecount_stc, peak, onoffindex = lnpc.onoffindex(v[:, 0], bins_stc, + spikecount_stc) + +#%% + np.savez(save_path, + sta=sta, + stimulus_type=stimulus_type, + total_frames=total_frames, + spikecount_sta=spikecount_sta, + spikecount_stc=spikecount_stc, + bins_sta=bins_sta, + bins_stc=bins_stc, + filename=filename, + onoffindex=onoffindex, + v=v, + w=w, + spike_path=spike_path, + total_spikes=total_spikes, + peak=peak, + stimulus_name=stimulus_name + ) + diff --git a/onoffplotter_dimos.py b/onoffplotter_dimos.py new file mode 100644 index 0000000..9d06a72 --- /dev/null +++ b/onoffplotter_dimos.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Jun 6 15:27:42 2017 + +@author: ycan +""" +import os +import numpy as np +import matplotlib.pyplot as plt + +main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ +LR3 Gollisch/data/Experiments/Mouse/2017_01_17/analyzed/' + +exp_name = main_dir.split('/')[-4]+' '+main_dir.split('/')[-3] + +allfiles = os.listdir(main_dir) + +files_f = [] # Full field flicker +files_c = [] # Checkerflicker + +for i in allfiles: + if i[-4:] == '.npz': + if i[0] == str(3): + files_f.append(i) + elif i[0] == str(2): + files_c.append(i) +onoffindices_f = np.array([]) +onoffindices_c = np.array([]) + +spikenr_f = np.array([]) +spikenr_c = np.array([]) + +filenames_f = [] +filenames_c = [] + +exclude_spike_limit = 200 +excluded_f = 0 +excluded_c = 0 + +for i in files_f: + f = np.load(main_dir+i) + spikenr_f = np.append(spikenr_f, f['total_spikes']) + if spikenr_f[-1] < exclude_spike_limit: + spikenr_f = spikenr_f[:-1] + excluded_f += 1 + continue + onoffindices_f = np.append(onoffindices_f, f['onoffindex']) + filenames_f.append(str(f['filename'])) + +for i in files_c: + f = np.load(main_dir+i) + spikenr_c = np.append(spikenr_c, f['total_spikes']) + if spikenr_c[-1] < exclude_spike_limit: + spikenr_c = spikenr_c[:-1] + excluded_c += 1 + continue + onoffindices_c = np.append(onoffindices_c, f['onoffindex']) + filenames_c.append(str(f['filename'])) +# %% Only get cells that are in both sets +filenamesc_f = [] +spikenrc_f = np.array([]) +onoffindicesc_f = np.array([]) +filenamesc_c = [] +spikenrc_c = np.array([]) +onoffindicesc_c = np.array([]) + +for i in range(len(filenames_f)): + if filenames_f[i] in filenames_c: + filenamesc_f.append(filenames_f[i]) + spikenrc_f = np.append(spikenrc_f, spikenr_f[i]) + onoffindicesc_f = np.append(onoffindicesc_f, onoffindices_f[i]) +for i in range(len(filenames_c)): + if filenames_c[i] in filenames_f: + filenamesc_c.append(filenames_c[i]) + spikenrc_c = np.append(spikenrc_c, spikenr_c[i]) + onoffindicesc_c = np.append(onoffindicesc_c, onoffindices_c[i]) + +spikenr_f = spikenrc_f +onoffindices_f = onoffindicesc_f +spikenr_c = spikenrc_c +onoffindices_c = onoffindicesc_c + +# %% +plt.hist(spikenr_f, np.linspace(0, spikenr_f.max(), num=250)) +plt.title('Spike numbers for FFF') +plt.show() +plt.hist(spikenr_c, np.linspace(0, spikenr_c.max(), num=250)) +plt.title('Spike numbers for Checkerflicker') +plt.show() +# %% +plt.figure(figsize=(10, 6)) +plt.hist(onoffindices_f, bins=np.linspace(-1, 1, num=40), alpha=.6) +plt.hist(onoffindices_c, bins=np.linspace(-1, 1, num=40), alpha=.6) +plt.legend(['Full field', 'Checkerflicker']) +plt.title('Histogram of On Off indices \n{}'.format(exp_name)) +plt.xlabel('On-off index') +plt.ylabel('Frequency') +plt.show() + +plt.figure(figsize=(8, 8)) +plt.scatter(onoffindices_f, onoffindices_c) +plt.plot(np.linspace(-1, 1), np.linspace(-1, 1), '--') +plt.title('On-Off indices obtained from Full field vs Checkerflicker\n{}'.format(exp_name)) +plt.ylabel('Checkerflicker') +plt.xlabel('Full field flicker') +plt.axis('square') +plt.show() diff --git a/plotall.py b/plotall.py index 1dc4221..18ae046 100644 --- a/plotall.py +++ b/plotall.py @@ -10,7 +10,7 @@ import matplotlib.pyplot as plt main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ -LR3 Gollisch/data/Experiments/Salamander/2014_01_27/analyzed/' +LR3 Gollisch/data/Experiments/Salamander/2014_02_25/analyzed/' exp_name = main_dir.split('/')[-4]+' '+main_dir.split('/')[-3] @@ -56,7 +56,8 @@ ax = plt.subplot(3, 3, 2) plt.plot(f['bins_sta'], f['spikecount_sta'], '-') plt.plot(f['bins_stc'], f['spikecount_stc'], '-') - plt.text(.5, .99, 'On-Off Bias: {:2.2f}'.format(float(f['onoffindex'])), + plt.text(.5, .99, 'On-Off Bias: {:2.2f}\nTotal spikes: {}' + .format(float(f['onoffindex']), f['total_spikes']), horizontalalignment='center', verticalalignment='top', transform=ax.transAxes) @@ -86,7 +87,8 @@ ax = plt.subplot(3, 3, 5) for i in range(len(c['bins'])): plt.plot(c['bins'][i], c['spike_counts_in_bins'][i], '-') - plt.text(.5, .99, 'On-Off Bias: {:2.2f}'.format(float(c['onoffindex'])), + plt.text(.5, .99, 'On-Off Bias: {:2.2f}\nTotal spikes: {}' + .format(float(c['onoffindex']), c['total_spikes']), horizontalalignment='center', verticalalignment='top', transform=ax.transAxes) @@ -118,6 +120,6 @@ vmax=np.max(c['sta_unscaled'])) plt.title('Brightest pixel: {}'.format(c['max_i'])) plt.tight_layout(pad=5, h_pad=1, w_pad=1.8) - +# plt.show() plt.savefig(savepath, dpi=200, bbox_inches='tight') plt.close() diff --git a/plotall_dimos.py b/plotall_dimos.py new file mode 100644 index 0000000..981d036 --- /dev/null +++ b/plotall_dimos.py @@ -0,0 +1,133 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Jun 9 12:16:38 2017 + +@author: ycan + +Plotting for data from Dimos + +The stimulus names are different than Fernando's data. +2 is checkerflicker +3 is full field flicker + +""" +import os +import numpy as np +import matplotlib.pyplot as plt + +main_dir = '/Users/ycan/Documents/official/gottingen/lab rotations/\ +LR3 Gollisch/data/Experiments/Mouse/2017_01_17/analyzed/' + +exp_name = main_dir.split('/')[-4]+' '+main_dir.split('/')[-3] + +allfiles = os.listdir(main_dir) + +files_f = [] # Full field flicker +files_c = [] # Checkerflicker + +for i in allfiles: + if i[-4:] == '.npz': + if i[0] == str(2): + files_f.append(i.split('C')[-1].split('.')[0]) + elif i[0] == str(3): + files_c.append(i.split('C')[-1].split('.')[0]) + +files = [i for i in files_c if i in files_f] + +for i in files: + # Changed this part because of stimulus order difference + fname_c = main_dir+'2_SP_C'+i+'.npz' + fname_f = main_dir+'3_SP_C'+i+'.npz' + + f = np.load(fname_f) + c = np.load(fname_c) + + savepath = '/'.join(main_dir.split('/')[:-1])+'/SP_C'+i + + # %% plot all + plt.figure(figsize=(12, 12), dpi=200) + plt.suptitle([' '.join(str(c['spike_path']) + .split('rasters')[0].split('Experiments')[1] + .split('/'))+str(i)]) + plt.subplot(3, 3, 1) + plt.plot(f['sta']) + plt.plot(f['v'][:, 0]) + plt.title('Filters') + plt.axvline(f['peak'], linewidth=1, color='r', linestyle='dashed') + plt.legend(['STA', 'Eigenvalue 0', 'Peak'], fontsize='small') + plt.xticks(np.linspace(0, 20, 20/2+1)) + plt.ylabel('Full field flicker\n$\\regular_{Linear\,\,\,output}$', + fontsize=16) + plt.xlabel('Time') + + ax = plt.subplot(3, 3, 2) + plt.plot(f['bins_sta'], f['spikecount_sta'], '-') + plt.plot(f['bins_stc'], f['spikecount_stc'], '-') + plt.text(.5, .99, 'On-Off Bias: {:2.2f}\nTotal spikes: {}' + .format(float(f['onoffindex']), f['total_spikes']), + horizontalalignment='center', + verticalalignment='top', + transform=ax.transAxes) + plt.title('Non-linearities') + plt.ylabel('Firing rate') + plt.xlabel('Linear output') + + plt.subplot(3, 3, 3) + plt.plot(f['w'], 'o') + plt.title('Eigenvalues of covariance matrix') + plt.xticks(np.linspace(0, 20, 20/2+1)) + plt.xlabel('Eigenvalue index') + plt.ylabel('Variance') + + plt.subplot(3, 3, 4) + plt.plot(c['sta_weighted']) + plt.plot(c['v'][:, 0]) + plt.plot(c['temporal']) + plt.axvline(c['peak'], linewidth=1, color='r', linestyle='dashed') + plt.title('Filters') + plt.ylabel('Checkerflicker\n$\\regular_{Linear\,\,\,output}$', fontsize=16) + plt.xlabel('Time') + plt.xticks(np.linspace(0, 20, 20/2+1)) + plt.legend(['Weighted stimulus', 'Eigenvalue 0', 'Brightest pixel', + 'Peak'], fontsize='small') + + ax = plt.subplot(3, 3, 5) + for i in range(len(c['bins'])): + plt.plot(c['bins'][i], c['spike_counts_in_bins'][i], '-') + plt.text(.5, .99, 'On-Off Bias: {:2.2f}\nTotal spikes: {}' + .format(float(c['onoffindex']), c['total_spikes']), + horizontalalignment='center', + verticalalignment='top', + transform=ax.transAxes) + plt.title('Non-linearities') + plt.xlabel('Linear output') + plt.ylabel('Firing rate') + + plt.subplot(3, 3, 6) + plt.plot(c['w'], 'o') + plt.title('Eigenvalues of covariance matrix') + plt.xticks(np.linspace(0, 20, 20/2+1)) + plt.xlabel('Eigenvalue index') + plt.ylabel('Variance') + + plt.subplot(3, 3, 7) + plt.imshow(c['sta_unscaled'][:, :, c['max_i'][2]].reshape((60, 80,)), + cmap='Greys', + vmin=np.min(c['sta_unscaled']), + vmax=np.max(c['sta_unscaled'])) + plt.title('Receptive field') + + plt.subplot(3, 3, 8) + f_size = 5 + plt.imshow(c['sta_unscaled'][c['max_i'][0]-f_size:c['max_i'][0]+f_size+1, + c['max_i'][1]-f_size:c['max_i'][1]+f_size+1, + int(c['max_i'][2])], + cmap='Greys', + vmin=np.min(c['sta_unscaled']), + vmax=np.max(c['sta_unscaled'])) + plt.title('Brightest pixel: {}'.format(c['max_i'])) + plt.tight_layout(pad=5, h_pad=1, w_pad=1.8) +# plt.show() + plt.savefig(savepath, dpi=200, bbox_inches='tight') + plt.close()