-
Notifications
You must be signed in to change notification settings - Fork 0
/
visualize_divergence_mouse_coronal_neurons.m
125 lines (110 loc) · 3.44 KB
/
visualize_divergence_mouse_coronal_neurons.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
% Compute overlap of samples per organ by computing the Jensen-Shannon
% divergence
clear variables
close all
%% Data
resultsdatadir = [LUMCDATADIR 'tSNE_ABA/normalized/results_mouse_genes_selection_cell_types/'];
datadir = [LUMCDATADIR 'tSNE_ABA/AllenMouseBrain_coronal/'];
% Load mapped data
load([resultsdatadir 'MappedMouseGenesSelection2Dneuronsid10'],'mappedX');
% Load voxel labels
load([datadir 'voxelLabels_coronal']);
mainStructures = {'Isocortex', 'OLF', 'HPF' 'CTXsp', 'STR', 'PAL', 'CB', 'TH', 'HY', 'MB', 'P', 'MY'};
load([datadir 'voxels']);
nlabels = max(voxelLabels);
%% Determine histogram size
mappedX = mappedX(:,1:2);
minbin = min(mappedX);
maxbin = max(mappedX);
nbins = 40;
%% Compute divergence
D = zeros(nlabels,nlabels);
for lnr1 = 1:nlabels
for lnr2 = 1:nlabels
if lnr2 >= lnr1-1
P = mappedX(voxelLabels==lnr1,:);
Q = mappedX(voxelLabels==lnr2,:);
d = jsdiv_pts(P,Q,nbins,minbin,maxbin);
D(lnr1,lnr2) = d;
else
D(lnr1,lnr2) = NaN;
% D(lnr1,lnr2) = 0;
end
end
end
binwidth = (maxbin-minbin)./nbins;
% save([resultsdatadir 'jsdiv_mouse_coronal_SVD'],'D','minbin','maxbin','nbins','mainStructures');
%% Show divergencies
Ddisp = [D zeros(size(D,1),1); nan(1,size(D,2)-1) 0 0];
fs = 15;
figure(1)
clf
hsurf = surface(Ddisp);
set(gcf,'Renderer','Painters')
set(gca,'FontSize',fs)
% set(gca,'XTick',1:nlabels,'XTickLabel',mainStructures);
% set(gca,'YTick',1:nlabels,'YTickLabel',mainStructures);
axis equal off
% set(gca,'YLim',[0.5 nlabels+0.5])
% rotateXLabels(gca,45)
cm = hot;
colormap(cm);
set(gcf,'Color','White')
rotate(hsurf,[0 0 1],-45)
textxloc = ((1:size(D,1))-size(D,1)/2)*sqrt(2)+size(D,1)/2-sqrt(2)/2+1;
textyloc = zeros(1,size(D,1))+size(D,1)/2+1.2+sqrt(2)/2;
for lnr = 1:size(D,1)
text(textxloc(lnr),textyloc(lnr),mainStructures{lnr},'FontSize',fs,'Rotation',90)
end
h = colorbar('horiz');
set(h,'FontSize',fs)
set(gca,'XLim',[textxloc(1)-sqrt(2) textxloc(end)+sqrt(2)])
set(gca,'YLim',[textxloc(1)-sqrt(2) textyloc(end)+4])
% orient landscape
set(gcf, 'InvertHardcopy', 'off','PaperPositionMode','auto');
saveas(1,'figure3d.eps','epsc2');
% tightInset = get(gca, 'TightInset');
% position(1) = tightInset(1);
% position(2) = tightInset(2);
% position(3) = 1 - tightInset(1) - tightInset(3);
% position(4) = 1 - tightInset(2) - tightInset(4);
% set(gca, 'Position', position);
% %% Plot clusters
% lnr1 = 1;
% lnr2 = 5;
%
% figure(2)
% clf
% hold on
% plot(mappedX(voxelLabels==lnr1,1),mappedX(voxelLabels==lnr1,2),'.');
% plot(mappedX(voxelLabels==lnr2,1),mappedX(voxelLabels==lnr2,2),'r.');
% lx = minbin(1):binwidth(1):maxbin(1);
% ly = minbin(2):binwidth(2):maxbin(2);
% plot([lx;lx],[ly(1)*ones(size(lx));ly(end)*ones(size(lx))],'k')
% plot([lx(1)*ones(size(ly));lx(end)*ones(size(ly))],[ly;ly],'k')
%
% %%
% figure(3)
% clf
% hold on
% sl1 = find(voxelLabels==lnr1);
% sl2 = find(voxelLabels==lnr2);
% for snr = sl1
% if ~isempty(voxel.color_HEX{snr})
% c = rgbconv(voxel.color_HEX{snr});
% else
% c = [0 0 0];
% end
% plot(mappedX(snr,1),mappedX(snr,2),'.','Color',c);
% end
% for snr = sl2
% if ~isempty(voxel.color_HEX{snr})
% c = rgbconv(voxel.color_HEX{snr});
% else
% c = [0 0 0];
% end
% plot(mappedX(snr,1),mappedX(snr,2),'.','Color',c);
% end
%
% plot([lx;lx],[ly(1)*ones(size(lx));ly(end)*ones(size(lx))],'k')
% plot([lx(1)*ones(size(ly));lx(end)*ones(size(ly))],[ly;ly],'k')