-
Notifications
You must be signed in to change notification settings - Fork 1
/
onoffstepsanalyzer.py
278 lines (221 loc) · 9.73 KB
/
onoffstepsanalyzer.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 30 18:18:03 2017
@author: ycan
"""
import os
import warnings
import numpy as np
import matplotlib.pyplot as plt
import analysis_scripts as asc
import plotfuncs as plf
import iofuncs as iof
def stim_prefr_durations_frametimes(frametimes, nr_per=4):
"""
Returns the average duration for each period in the stimulus from
frametimes.
Parameters
--------
frametimes:
Array containing times where the frame was updated.
nr_per:
Number of periods in one trial. If preframes is not zero, this
will be four (grey-white-grey-black). If zero, it will be two
(white-black).
Returns
-----
trial_durs:
The duration for each period, shaped in a (nr_trials x nr_periods)
array. The last trial is discarded if it is not complete.
"""
# Calculate each period duration
durs = np.ediff1d(frametimes)
# Reshape into trial structure, discard incomplete trial at the end
durs = durs[:(durs.shape[0]//nr_per)*nr_per];
durs = durs.reshape((durs.shape[0]//nr_per, nr_per));
return durs
def onoffstepsanalyzer(exp_name, stim_nrs):
"""
Analyze onoffsteps data, plot and save it. Will make a directory
/data_analysis/<stimulus_name> and save svg [and pdf in subfolder.].
Parameters:
exp_name:
Experiment name.
stim_nr:
Order of the onoff steps stimulus.
"""
exp_dir = iof.exp_dir_fixer(exp_name)
exp_name = os.path.split(exp_dir)[-1]
if isinstance(stim_nrs, int):
stim_nrs = [stim_nrs]
for stim_nr in stim_nrs:
stim_nr = str(stim_nr)
stimname = iof.getstimname(exp_dir, stim_nr)
clusters, metadata = asc.read_spikesheet(exp_dir, cutoff=4)
clusterids = plf.clusters_to_ids(clusters)
parameters = asc.read_parameters(exp_dir, stim_nr)
refresh_rate = metadata['refresh_rate']
# Divide by the refresh rate to convert from number of
# frames to seconds
pars_stim_duration = parameters['Nframes']/refresh_rate
pars_preframe_duration = parameters.get('preframes',
0)/refresh_rate
if pars_preframe_duration == 0:
nopreframe = True
nr_periods = 2
else:
nopreframe = False
nr_periods = 4
# The first trial will be discarded by dropping the first four frames
# If we don't save the original and re-initialize for each cell,
# frametimings will get smaller over time.
frametimings_original = asc.readframetimes(exp_dir, stim_nr)
trial_durs = stim_prefr_durations_frametimes(
frametimings_original, nr_per=nr_periods)
avg_trial_durs = trial_durs.mean(axis=0)
if not nopreframe:
stim_duration = avg_trial_durs[1::2].mean()
preframe_duration = avg_trial_durs[::2].mean()
else:
stim_duration = avg_trial_durs.mean()
preframe_duration = 0
warnings.warn('On-off steps analysis with no preframes'
'is not tested, proceed with caution.')
contrast = parameters['contrast']
total_cycle = avg_trial_durs.sum()
# Set the bins to be 10 ms
tstep = 0.01
bins = int(total_cycle/tstep)+1
t = np.linspace(0, total_cycle, num=bins)
# Setup for onoff bias calculation
onbegin = preframe_duration
onend = onbegin+stim_duration
offbegin = onend+preframe_duration
offend = offbegin+stim_duration
# Determine the indices for each period
a = []
for i in [onbegin, onend, offbegin, offend]:
yo = np.asscalar(np.where(np.abs(t-i) < tstep/1.5)[0][-1])
a.append(yo)
# To exclude stimulus offset affecting the bias, use
# last 1 second of preframe period
prefs = []
for i in [onbegin-1, onbegin, offbegin-1, offbegin]:
yo = np.asscalar(np.where(np.abs(t-i) < tstep/1.5)[0][-1])
prefs.append(yo)
onper = slice(a[0], a[1])
offper = slice(a[2], a[3])
pref1 = slice(prefs[0], prefs[1])
pref2 = slice(prefs[2], prefs[3])
onoffbias = np.empty(clusters.shape[0])
baselines = np.empty(clusters.shape[0])
savedir = os.path.join(exp_dir, 'data_analysis', stimname)
os.makedirs(os.path.join(savedir, 'pdf'), exist_ok=True)
# Collect all firing rates in a list
all_frs = []
for i in range(len(clusters[:, 0])):
spikes = asc.read_raster(exp_dir, stim_nr,
clusters[i, 0], clusters[i, 1])
frametimings = frametimings_original
# Discard all the spikes that happen after the last frame
spikes = spikes[spikes < frametimings[-1]]
# Discard the first trial
spikes = spikes[spikes > frametimings[4]]
frametimings = frametimings[4:]
# Find which trial each spike belongs to, and subtract one
# to be able to use as indices
trial_indices = np.digitize(spikes, frametimings[::4])-1
rasterplot = []
# Iterate over all the trials, create an empty array for each
for j in range(int(np.ceil(frametimings.max()/total_cycle))):
rasterplot.append([])
# plt.eventplot requires a list containing spikes in each
# trial separately
for k in range(len(spikes)):
trial = trial_indices[k]
rasterplot[trial].append(spikes[k]-frametimings[::4][trial])
# Workaround for matplotlib issue #6412.
# https://github.com/matplotlib/matplotlib/issues/6412
# If a cell has no spikes for the first trial i.e. the first
# element of the list is empty, an error is raised due to
# a plt.eventplot bug.
if len(rasterplot[0]) == 0:
rasterplot[0] = [-1]
plt.figure(figsize=(9, 9))
ax1 = plt.subplot(211)
plt.eventplot(rasterplot, linewidth=.5, color='r')
# Set the axis so they align with the rectangles
plt.axis([0, total_cycle, -1, len(rasterplot)])
# Draw rectangles to represent different parts of the on off
# steps stimulus
plf.drawonoff(ax1, preframe_duration, stim_duration,
contrast=contrast)
plt.ylabel('Trial')
plt.gca().invert_yaxis()
ax1.set_xticks([])
plf.spineless(ax1)
# Collect all trials in one array to calculate firing rates
ras = np.array([])
for ii in range(len(rasterplot)):
ras = np.append(ras, rasterplot[ii])
# Sort into time bins and count how many spikes happened in each
fr = np.digitize(ras, t)
fr = np.bincount(fr)
# Normalize so that units are spikes/s
fr = fr * (bins/total_cycle) / (len(rasterplot)-1)
# Equalize the length of the two arrays for plotting.
# np.bincount(x) normally produces x.max()+1 bins
if fr.shape[0] == bins+1:
fr = fr[:-1]
# If there aren't any spikes at the last trial, the firing
# rates array is too short and plt.plot raises error.
while fr.shape[0] < bins:
fr = np.append(fr, 0)
prefr = np.append(fr[pref1], fr[pref2])
baseline = np.median(np.round(prefr))
fr_corr = fr - baseline
r_on = np.sum(fr_corr[onper])
r_off = np.sum(fr_corr[offper])
if r_on == 0 and r_off == 0:
bias = np.nan
else:
bias = (r_on-r_off)/(np.abs(r_on)+np.abs(r_off))
plt.suptitle(f'{exp_name}\n{stimname}'
f'\n{clusterids[i]} Rating: {clusters[i][2]}\n')
if fr.max() < 20:
bias = np.nan
onoffbias[i] = bias
baselines[i] = baseline
all_frs.append(fr)
ax2 = plt.subplot(212)
plt.plot(t, fr)
for eachslice in [onper, offper]:
ax2.fill_between(t[eachslice], fr[eachslice],
baseline, where=fr[eachslice] > baseline,
facecolor='lightgray')
plf.spineless(ax2)
plt.axis([0, total_cycle, fr.min(), fr.max()])
plt.title(f'Baseline: {baseline:2.0f} Hz Bias: {bias:0.2f}')
plt.xlabel('Time[s]')
plt.ylabel('Firing rate[spikes/s]')
# Save as svg for looking through data, pdf for
# inserting into presentations
plt.savefig(savedir+'/{:0>3}{:0>2}.svg'.format(clusters[i, 0],
clusters[i, 1]),
format='svg', bbox_inches='tight')
plt.savefig(os.path.join(savedir, 'pdf',
'{:0>3}'
'{:0>2}.pdf'.format(clusters[i, 0],
clusters[i, 1])),
format='pdf', bbox_inches='tight')
plt.close()
keystosave = ['clusters', 'total_cycle', 'bins', 'tstep',
'stimname', 'stim_duration', 'preframe_duration',
'contrast', 'all_frs', 't', 'exp_name', 'onoffbias',
'baselines']
data_in_dict = {}
for key in keystosave:
data_in_dict[key] = locals()[key]
np.savez(os.path.join(savedir, stim_nr + '_data'), **data_in_dict)
print(f'Analysis of {stimname} completed.')