-
Notifications
You must be signed in to change notification settings - Fork 0
/
pupil_saccade_removal.m
180 lines (127 loc) · 7.11 KB
/
pupil_saccade_removal.m
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
%% Check firing rates just after saccade removal
%
% compares firing rate during motion and stationary periods with and
% without the period just after saccade
%%
experiment_groups = 'visual_flow';
trial_types = {'RVT', 'RV', {'VT_RVT', 'VT_RV'}, {'V_RVT', 'V_RV'}};
duration_to_remove = 0.25; % amount of time in seconds to removal after a saccade
%%
ctl = RC2Analysis();
% probe recordings for this experiment group
probe_ids = ctl.get_probe_ids(experiment_groups);
probe_id_store = cell(1, length(trial_types));
prc_removed_motion = cell(1, length(trial_types));
prc_removed_stationary = cell(1, length(trial_types));
noncropped_motion_fr = cell(1, length(trial_types));
cropped_motion_fr = cell(1, length(trial_types));
noncropped_stationary_fr = cell(1, length(trial_types));
cropped_stationary_fr = cell(1, length(trial_types));
cluster_ids = cell(1, length(trial_types));
spike_class = cell(1, length(trial_types));
for probe_i = 1 : length(probe_ids)
% skip this mouse as the camera was started late
if strcmp(probe_ids{probe_i}, 'CAA-1114977_rec1_rec2_rec3')
continue
end
% get probe recording data
data = ctl.load_formatted_data(probe_ids{probe_i});
clusters = data.VISp_clusters();
% get table of saccades
sessions = data.motion_sessions;
saccade_tbl = ctl.load.camera0_saccades(sessions{1}.session_id);
% time of each saccade
saccade_times = sessions{1}.camera_t(saccade_tbl.saccade_frame);
for type_i = 1 : length(trial_types)
% get trials of correct type
trials = data.get_trials_with_trial_group_label(trial_types{type_i});
% setup storage arrays
prc_removed_motion{type_i}{probe_i} = nan(length(trials), 1);
prc_removed_stationary{type_i}{probe_i} = nan(length(trials), 1);
noncropped_motion_fr{type_i}{probe_i} = nan(length(trials), length(clusters));
cropped_motion_fr{type_i}{probe_i} = nan(length(trials), length(clusters));
noncropped_stationary_fr{type_i}{probe_i} = nan(length(trials), length(clusters));
cropped_stationary_fr{type_i}{probe_i} = nan(length(trials), length(clusters));
cluster_ids{type_i}{probe_i} = nan(length(clusters), 1);
spike_class{type_i}{probe_i} = cell(length(clusters), 1);
probe_id_store{type_i}{probe_i} = probe_ids{probe_i};
figure
for trial_i = 1 : length(trials)
fprintf('%i/%i, %i/%i\n', trial_i, length(trials), probe_i, length(probe_ids));
% get this trial and original trial
trial = trials{trial_i}.to_aligned;
% mask of motion in the trial
motion_mask = trial.motion_mask;
stationary_mask = trial.stationary_mask;
% find the saccades in this trial
timebase = trial.probe_t;
idx = saccade_times > timebase(1) & saccade_times < timebase(end);
trial_saccade_times = saccade_times(idx);
saccade_mask = false(size(timebase));
for sac_i = 1 : length(trial_saccade_times)
saccade_mask = saccade_mask | ...
(timebase > trial_saccade_times(sac_i) & timebase < ...
trial_saccade_times(sac_i) + duration_to_remove);
end
% plot the saccades on the trials
subplot(4, 5, trial_i)
plot(trial.velocity)
hold on;
plot(saccade_mask * 5)
box off
if iscell(trial_types{type_i})
str = trial_types{type_i}{1};
else
str = trial_types{type_i};
end
title(sprintf('%s, trial ID %i, %s', probe_ids{probe_i}, trial.trial_id, str), 'interpreter', 'none');
% percentage of data removed due to the motion energy
prc_removed_motion{type_i}{probe_i}(trial_i) = sum(motion_mask & saccade_mask) / sum(motion_mask);
prc_removed_stationary{type_i}{probe_i}(trial_i) = sum(stationary_mask & saccade_mask) / sum(stationary_mask);
% for each cluster
for clust_i = 1 : length(clusters)
% convolved firing rate during replay trial, aligned
fr = clusters(clust_i).fr.get_convolution(timebase);
% the original computed firing rate
noncropped_motion_fr{type_i}{probe_i}(trial_i, clust_i) = mean(fr(motion_mask));
noncropped_stationary_fr{type_i}{probe_i}(trial_i, clust_i) = mean(fr(stationary_mask));
% the average firing rate with 'camera motion' removed
cropped_motion_fr{type_i}{probe_i}(trial_i, clust_i) = mean(fr(motion_mask & ~saccade_mask));
cropped_stationary_fr{type_i}{probe_i}(trial_i, clust_i) = mean(fr(stationary_mask & ~saccade_mask));
% store info
cluster_ids{type_i}{probe_i}(clust_i) = clusters(clust_i).id;
spike_class{type_i}{probe_i}{clust_i} = clusters(clust_i).spiking_class;
end
end
end
end
%% print statistics
print_statistics(noncropped_motion_fr, cropped_motion_fr, trial_types, 'Motion')
print_statistics(noncropped_stationary_fr, cropped_stationary_fr, trial_types, 'Stationary')
function print_statistics(noncropped_fr, cropped_fr, trial_types, str)
delta_fr = cell(1, length(trial_types));
noncropped_median = cell(1, length(trial_types));
cropped_median = cell(1, length(trial_types));
for type_i = 1 : length(trial_types)
for probe_i = 1 : length(noncropped_fr)
trial_delta_fr = noncropped_fr{type_i}{probe_i} - cropped_fr{type_i}{probe_i};
delta_fr{type_i} = [delta_fr{type_i}, median(trial_delta_fr, 1)];
noncropped_median{type_i} = [noncropped_median{type_i}, median(noncropped_fr{type_i}{probe_i}, 1)];
cropped_median{type_i} = [cropped_median{type_i}, median(cropped_fr{type_i}{probe_i}, 1)];
end
end
for type_i = 1 : length(trial_types)
p_val_signrank = signrank(noncropped_median{type_i}(:), cropped_median{type_i});
p_val_sign = signtest(delta_fr{type_i});
% print info
if iscell(trial_types{type_i})
trial_str = strjoin(trial_types{type_i}, ',');
else
trial_str = trial_types{type_i};
end
fprintf('Period: %s\n', str);
fprintf('Trial type: %s\n', trial_str);
fprintf(' p-value Wilcoxon signrank: %.2f\n', p_val_signrank);
fprintf(' p-value sign-test: %.2f\n', p_val_sign);
end
end