-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathspin_receptor_map.m
235 lines (163 loc) · 7.4 KB
/
spin_receptor_map.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
%% spin test (spin receptor map - maps obtained from Vasa code) - Figure 4a/b iv
% ref: Vasa et al Cerebral Cortex 2018
% (https://doi.org/10.1093/cercor/bhx249)
clear all; close all;
basedir = '/Users/sps253/Documents/energy_landscape';
cd(basedir);
savedir = fullfile(basedir,'results','example');mkdir(savedir); % set save directory
%%
split='main'
load(fullfile(['data/',split,'.mat']))
numClusters=4; T=0.001;
load(['subjenergy_split',num2str(split),'_k',num2str(numClusters),'_T',num2str(T),'.mat'], 'centroids', 'E_weighted_PLavg');
load(['Partition_bp',num2str(split),'_k',num2str(numClusters),'.mat'], 'clusterNames');
%% load SC - for spin test to work, need to add missing region back in
if nparc == 461
load data/Lausanne463_HCP_DTI_count.mat vol_normalized_sc% connectivity
% load ls463_DTI_fiber_consensus_HCP.mat connectivity %consensus matrix (less sparse)
A = centroids(:,1:13,:);
B = centroids(:,244,:);
C = centroids(:,14:end,:);
% D = zeros(30,1,4); sans brainstem
centroids = horzcat(A,B,C);
nparc=462;
connectivity = vol_normalized_sc;
connectivity(463,:)=[]; %sans brainstem
connectivity(:,463)=[];
sc = connectivity;
end
%%
c = 0;
Anorm = NORMALIZE(sc,c);
% define x0 and xf, initial and final states as cluster centroids for each
% state transition
Xf_ind = repmat(1:numClusters,[1 numClusters]); % final state order
Xo_ind = repelem(1:numClusters,numClusters); % paired with different initial states, use reshape(x,[numClusters numClusters])' to get matrix
onDiag = (1:numClusters) + (numClusters*(0:(numClusters-1)));
offDiag = 1:(numClusters^2); offDiag(onDiag) = []; % isolate off diagonal from linearized transition probabilities
if nparc == 454
load data/5HTvecs_sch454.mat mean5HT2A_sch454
HT = mean5HT2A_sch454;
elseif nparc == 232
load data/5HTvecs_sch232.mat mean5HT2A_sch232
HT = mean5HT2A_sch232;
elseif nparc == 462
load data/5HTvecs_ls463.mat mean5HT2A_ls463
HT = mean5HT2A_ls463;
end %weight towards 5HT2a
norm_HT = (HT/max(HT))';
%% spin test
if nparc ==232
load(['data/SpinTests/rotated_maps/rotated_Schaefer_200.mat'])
load('data/sch232_to_yeo.csv');
network7labels=sch232_to_yeo';
nosub_HT=norm_HT;
nosub_HT(network7labels>7)=[];
spun_2a = [];
for i=1:size(perm_id,2)
for r=1:size(perm_id,1)
spun_2a(r,i) = nosub_HT(perm_id(r,i));
end
end
a = spun_2a; % cortical assignments
b = reshape(repelem(norm_HT(201:232),10000),10000,32)'; % subcortical assignments
perm_map_2a = vertcat(a,b);
end
if nparc ==462
load(['data/SpinTests/rotated_maps/rotated_Lausanne_463.mat'])
load('data/Lausanne_463_subnetworks.mat');
network7labels=subnetworks_reorder;
nosub_HT=norm_HT;
nosub_HT(network7labels>7)=[];
perm_id = perm_id+1;
spun_2a = [];
for i=1:size(perm_id,2)
for r=1:size(perm_id,1)
spun_2a(r,i) = nosub_HT(perm_id(r,i));
end
end
a = spun_2a(1:223,:); %right hemi cortical assignments
b = reshape(repelem(norm_HT(224:230),10000),10000,7)'; %right hemi subcortical assignments
c = spun_2a(224:end,:); %left hemi cortical assignments
d = reshape(repelem(norm_HT(456:463),10000),10000,8)'; %left hemi subcortical assignments
perm_map_2a = vertcat(a,b,c,d);
perm_map_2a(463,:) = []; %sans brainstem
end
%%
nperms=size(perm_map_2a,2);
x0 = squeeze(mean(centroids(nsubjs+1:nsubjs*2,:,Xo_ind)));
xf = squeeze(mean(centroids(nsubjs+1:nsubjs*2,:,Xf_ind)));
NullTransitionEnergy = zeros(nperms,numClusters^2);
for P = 1:nperms
disp(['Perm ',num2str(P)])
InputVector = perm_map_2a(:,P)'; %randomly reshuffle input vector
B = InputVector .*eye(nparc) + eye(nparc); % construct input matrix allowing input only into selected regions in InputVector
E_weighted_RAND = zeros(1,numClusters^2);
for transition = 1:numClusters^2
[x, u] = MIN_ENG_CONT(Anorm, T, B, x0(:,transition), xf(:,transition), 0);
E_weighted_RAND(transition) = sum(sum(u.^2))*T/1001; % integrate over inputs
end
NullTransitionEnergy(P,:) = E_weighted_RAND;
end
%% compute pvalues
E_matrix = reshape(E_weighted_PLavg, [numClusters numClusters])';
[nperms,num_transitions] = size(NullTransitionEnergy);
NullTransitionEnergyMatrix = zeros(numClusters,numClusters,nperms);
for P = 1:nperms
NullTransitionEnergyMatrix(:,:,P) = reshape(NullTransitionEnergy(P,:),[numClusters numClusters])';
end
% sig_thresh = 0.05 / num_transitions; % bonferroni correct over all transitions
pvals_onetail = mean(E_matrix > NullTransitionEnergyMatrix,3); % get p-value as % of times real energy is higher than energy in null matrix for a transition
fdrpv1t = mafdr(reshape(pvals_onetail,1,numClusters^2),'BHFDR',1);
fdrpv1t = reshape(fdrpv1t,numClusters,numClusters);
%% visualize
randdist = mean(NullTransitionEnergyMatrix,3);
truedist = reshape(E_weighted_PLavg, [numClusters numClusters])';
maxVal = max(max([randdist,truedist])); % sync color scales
minVal = min(min([randdist,truedist]));
clusterNames={'MS-1a';'MS-1b';'MS-2a';'MS-2b'};
figure;
subplot(1,3,1);
imagesc(randdist);
xticks(1:numClusters); yticks(1:numClusters);
xticklabels(clusterNames); xtickangle(90); yticklabels(clusterNames); axis square;
COLOR_TICK_LABELS(true,true,numClusters);
ylabel('Current State'); xlabel('Next State');
title('Random');
set(gca,'FontSize',18);
set(gca,'TickLength',[0 0]);
set(gca,'Fontname','arial');
caxis([minVal maxVal]); colorbar
subplot(1,3,2);
imagesc(truedist);
xticks(1:numClusters); yticks(1:numClusters);
xticklabels(clusterNames); xtickangle(90); yticklabels(clusterNames); axis square;
COLOR_TICK_LABELS(true,true,numClusters);
ylabel('Current State'); xlabel('Next State');
title('True');
set(gca,'FontSize',18);
set(gca,'TickLength',[0 0]);
set(gca,'Fontname','arial');
caxis([minVal maxVal]); colorbar
subplot(1,3,3);
LSDMinusPLTP = fdrpv1t;%(grpDiff); %((grpAvgPL - grpAvgLSD)); %switching order for manuscript figures
imagesc(LSDMinusPLTP); colormap('viridis');
xticks(1:numClusters); xticklabels(clusterNames); xtickangle(90);
yticks(1:numClusters); yticklabels(clusterNames); axis square
ylabel('Initial State'); xlabel('Final State');
sig_thresh = 0.05; %/numClusters^2; % bonferroni correction, for two-tailed p-values so only
[y,x] = find(fdrpv1t < sig_thresh); %this was used with two-tailed perm test-> p___.*~eye(numClusters)
text(x-.15,y+.18,'**','Color','w','Fontsize', 36);
[y2,x2] = find(pvals_onetail < sig_thresh); %this was used with two-tailed perm test-> p___.*~eye(numClusters)
text(x2-.15,y2+.18,'*','Color','w','Fontsize', 36);
u_caxis_bound = .5;%max(max(LSDMinusPLTP));
l_caxis_bound = 0;%min(min(LSDMinusPLTP));
h = colorbar; ylabel(h,'corrected p'); caxis([l_caxis_bound u_caxis_bound]); h.Ticks = [l_caxis_bound (u_caxis_bound+l_caxis_bound)/2 u_caxis_bound];
h.TickLabels = [round(l_caxis_bound,2,'significant') round((l_caxis_bound+u_caxis_bound)/2,1,'significant') round(u_caxis_bound,1,'significant')];
COLOR_TICK_LABELS(true,true,numClusters);
title('Shuffled > True');
set(gca,'FontSize',18);
set(gca,'TickLength',[0 0]);
set(gca,'Fontname','arial');
%% save
save(fullfile(savedir,['receptorSPIN__split',num2str(split),'_k',num2str(numClusters),'_T',num2str(T),'.mat']),'E_weighted_PLavg','NullTransitionEnergy','pvals_onetail','fdrpv1t','nperms');