-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexvivoRGC_RF_Analysis_plotting.m
298 lines (239 loc) · 8.08 KB
/
exvivoRGC_RF_Analysis_plotting.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
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
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
%Tested on Matlab 2019b
%requires:
% panel (from Matlab file exchange)
% othercolor (from Matlab file exchange)
load("./exvivo_RFs.mat");
addpath("./utils");
%%
fprintf("Flipping ON rfs to OFF\n");
spatial_unsigned_all = flip_on_rfs(spatial_all);
clear spatial_all; %save RAM
%% Fig 3a. Binned Spatial RFs from one retina (retina 5)
f = figure;
idx_ret = find(ret_number == 5);
idx_map = bin_2d(centroids(idx_ret,:), 'binsize', 300, 'mincount',5);
counts = plot_binned_rfs(idx_map, spatial_unsigned_all(:,:,idx_ret));
sgtitle(sprintf("Fig 3a: Cell count = %.1f ± %.1f", mean(counts), std(counts)));
f.Position = [0 0 750 750]; %ensure figure is square for proper alignment
%% Extended Data Fig 5. Binned Spatial RFs from all sopsin retinas
f = figure;
idx_ret = find(sopsin_gt);
idx_map = bin_2d(centroids(idx_ret,:), 'binsize', 300, 'mincount',5);
counts = plot_binned_rfs(idx_map, spatial_unsigned_all(:,:,idx_ret));
sgtitle(sprintf("Extended Fig 5: Cell count = %.1f ± %.1f", mean(counts), std(counts)));
f.Position = [0 0 750 750]; %ensure figure is square for proper alignment
%% Fig 3b. Mean RFs from below, at and above horizon
figure;
cond = sopsin_gt; %only use sopsin retinas
as = [-1500, -900, 500]; bs = [-1000, -700, 600]; %bin limits in um
for i = 1:3
subplot(2,3,3+i);
a = as(i); b = bs(i);
idx = centroids(:,2) > a & centroids(:,2) < b & cond;
spatial = squeeze(mean(spatial_unsigned_all(:,:,idx), 3));
plot_spatial(spatial, sprintf("N = %d", nnz(idx)));
title(sprintf("[%d, %d] µm", a, b));
end
sgtitle("Fig 3b");
%% Fig 3c. Radial RF profiles (same bins as above)
fprintf("computing radial profiles of all rfs, ~2mins\n");
radial_all = zeros(size(spatial_unsigned_all,3), 50);
for i = 1:size(spatial_unsigned_all,3)
radial_all(i,:) = radial_profile(spatial_unsigned_all(:,:,i));
end
%%
for i = 1:3
a = as(i); b = bs(i);
idx = centroids(:,2) > a & centroids(:,2) < b & cond;
radials = radial_all(idx,:);
avg_rad(i,:) = mean(radials,1);
caption{i} = sprintf("[%d, %d] µm, n = %d", a, b, nnz(idx));
end
figure;
subplot(2,1,2);
plot(linspace(0,500,size(radials,2)),avg_rad', 'LineWidth', 3);
yline(0, 'LineStyle', '--');
legend(caption,'Location', 'southeastoutside');
xlabel("Radius (um)");
ylabel("RF Profile Weight");
ylim([-.1 .1]);
yticks([-.1 0 .1]);
xticks([0 500]);
xlim([0 500]);
colorv = [144 25 28]/255;
colorh = [184 72 38]/255;
colord = [231 138 36]/255;
colororder([colord;colorh;colorv]);
sgtitle("Fig 3c");
%% Fig 3d-i RF Properties binned in 1D or 2D
% and extended Data Fig. 7
clear cond;
cond{1} = find(sopsin_gt);
cond{2} = find(~sopsin_gt);
for c = 1:2
idx = cond{c};
idx_map = bin_2d(centroids(idx,:), 'binsize', 100);
colors = {flipud(othercolor('PRGn7')), flipud(othercolor('PRGn7')), othercolor('BrBG11')};
titles = ["Relative Surround Strength", "Center Size", "Vertical Surround Asymmetry"];
params{1} = scRatio(idx);
params{2} = centSize(idx);
params{3} = rfAsymm(idx);
figure;
if c == 1; sgtitle("Fig 3d-i"); else sgtitle("Extended Fig. 7"); end
for i = 1:3
% 2D Bins
subplot(2,3,i);
output = propertyMap(idx_map, params{i});
nz = output(output~=0);
mini = quantile(nz(:), .05); %remove extremes from dynamic range
maxi = quantile(nz(:), .95);
title(titles(i));
imAlpha = ones(size(output));
imAlpha(isnan(output)) = 0;
hold on;
if i~=3
imagesc(output, 'AlphaData', imAlpha, [mini, maxi]);
else
imagesc(output, 'AlphaData', imAlpha, [-1,1]);
end
xticks([]); yticks([]);
colorbar;
ax = imgca();
axis image; axis xy;
colormap(ax, colors{i});
mid = size(output,1)/2;
plot(mid,mid, '+');
% 1D Bins (per retina)
subplot(2,3,3+i);
pts = plot_trend_by_group(centroids(idx,:), params{i}, ret_number(idx));
ylabel(titles(i));
if i~=3; legend('hide'); end
dors = pts(:,1) < 0;
[~,pval] = kstest2(pts(dors,2), pts(~dors,2));
title(sprintf("Pvalue: %e", pval));
end
end
%% Fig. 3j (Example surround orientations from retina 9)
figure;
colors = [colorv;colorh;colord];
for i = 1:3
subplot(3,1,i);
orts = [gprops_all(idxs{i}).orientation_centers_of_mass_data_masks];
alpha = wrapTo180(sopsin_dorsal(9)-90);
orts = wrapToPi(-deg2rad(alpha) - orts);
histogram(orts, 20,'Normalization', 'probability', 'FaceColor', colors(i,:), 'EdgeAlpha', 0);
ylim([0 .5]);
hold on;
xlim([-pi pi]);
xticks(linspace(-pi,pi, 5));
xticklabels(["T", "D", "N", "V", "T"]);
end
sgtitle("Fig. 3j");
%% Fig. 3k (Surround asymmetry in sinusoidal projection of visual space)
alpha0 = deg2rad(22); %final coordinates of optical axis = resting position of mouse eye
theta0 = deg2rad(64);
cents = polar_to_visual(polar_coords, alpha0, theta0);
minlim = min(cents, [], 'all'); maxlim = max(cents, [], 'all');
lims = [minlim maxlim];
binsize = (maxlim - minlim)/40; %divide into 40 bins
idx_map = bin_2d(cents, 'binsize', binsize, 'lims', lims, 'mincount', 0);
map = propertyMap(idx_map, rfAsymm);
figure;
colormap(othercolor('BrBG11'));
imAlpha = ones(size(map));
imAlpha(isnan(map)) = 0;
image('CData', map, 'XData', lims, 'YData', lims,'CDataMapping','scaled', 'AlphaData', imAlpha);
caxis([-1 1]);
colorbar;
hold on;
plotGrid();
axis equal
% Also plot the rim of a retina
phiRim = zeros(1,100)'; %latitude = 0 defines the equator in retinal coords
lambdaRim = linspace(-pi,pi,100)'; % make a line passing through all longitudes
cent_rim = polar_to_visual([lambdaRim phiRim], alpha0, theta0);
x = cent_rim(:,1); y = cent_rim(:,2);
kink = find(abs(diff(x)) > 2);
x = circshift(x,-kink);
y = circshift(y, -kink);
plot(x,y, 'LineWidth', 2, 'Color', 'k');
title("Surround Asymmetry in Visual Coordinates");
xticks([-pi 0 pi]);
xticklabels({'-180', 'Front', '180'});
yticks([-pi/2 0 pi/2]);
yticklabels({'ground', 'horizon', 'sky'});
xlabel("Azimuth");
ylabel("Elevation");
%% Fig 4a
k = numel(unique(group_idx));
temps = temporal_all./max(abs(temporal_all), [], 2); %normalize
nCellsPerClust = zeros(k,1);
for i = 1:k
nCellsPerClust(i) = nnz(group_idx == i);
end
inds = zeros(length(group_idx), 1);
idxPic = zeros(length(group_idx), 1);
ccmap = jet(k);
cnt = 0;
for cluster_id = 1:k
cnt1 = cnt + sum(group_idx == cluster_id);
inds(cnt+1:cnt1) = find(group_idx == cluster_id);
idxPic(cnt+1:cnt1) = cluster_id;
cnt = cnt1;
end
figure;
ax_cluster = subplot(1, 12, 1);
imagesc(idxPic);
ylabel('Cluster ID');
xlabel('');
set(gca, 'YTickLabels', {''}, 'YTick', []);
colormap(ax_cluster, jet)
ax_data = subplot(1, 12, 2:12);
imagesc(temps(inds,:), [-1.5,1.5]);
colormap(flipud(othercolor('RdBu11')));
xticks([1 size(temps,2)]);
xticklabels({0 num2str(25*size(temps,2))});
title('Temporal RFs');
xlabel('Time (ms)');
sgtitle("Fig 4a");
figure;
for id = 1:k
subplot(5,2,id);
idx = group_idx == id;
traces = temps(idx,:);
avg = mean(traces, 1)';
stds = std(traces, 0, 1)';
hold on;
plot(avg-stds, 'Color', 'k', 'LineWidth', 1);
plot(avg+stds, 'Color', 'k', 'LineWidth', 1);
plot(avg, 'LineWidth', 2, 'Color', "k");
title(sprintf("cluster %d, n = %d", id,nnz(group_idx == id)));
ylim([-1.2 1.2]);
xlim([1, size(temps, 2)]);
xticks([1 size(temps,2)]);
xticklabels({0 num2str(25*size(temps,2))});
end
sgtitle("Fig 4a");
%% Fig 4c-e 1D Bins of properties, within each cluster
idx = find(sopsin_gt);
titles = ["Relative Surround Strength", "Center Size", "Vertical Surround Asymmetry"];
params{1} = scRatio(idx);
params{2} = centSize(idx);
params{3} = rfAsymm(idx);
figure;
for i = 1:3
subplot(1,3,i);
pts = plot_trend_by_group(centroids(idx,:), params{i}, group_idx(idx));
ylabel(titles(i));
if i~=3; legend('hide'); end
dors = pts(:,1) < 0;
[~,pval] = kstest2(pts(dors,2), pts(~dors,2));
title(sprintf("Pvalue: %e", pval));
end
%% Helper functions
function unsigned = flip_on_rfs(signed)
unsigned = signed;
mid = ceil(size(signed,1)/2);
onrfs = signed(mid,mid,:) > 0;
unsigned(:,:,onrfs) = -signed(:,:,onrfs);
end