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

@@ -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
+             )

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
 

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
+             )
+

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()