-
Notifications
You must be signed in to change notification settings - Fork 1
/
runNBS_old.m
executable file
·311 lines (247 loc) · 11.1 KB
/
runNBS_old.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
299
300
301
302
303
304
305
306
307
308
309
310
%% RUNNBS
% RUN NBS analyzes the output of hopf_particleswarm script. Specifically,
% runNBS applies the directed network-based statistic (NBS) to the
% effective connectivity matrices generated by hopf_particleswarm.
% The NBS utilizes an version of cluster-based statistical methods adapted
% to network topology. This allows the NBS to drastically reduce the
% multiple comparison problem when comparing links, thus dramatically
% improving its sensitivity to distributed effects. However, it is no more
% sensitive to focal effects than any other method.
%% SETUP
%% Set paths & directories
clear; clc; close all;
% Shuffle random seed. Necessary in array parallelization to avoid
% repeating same random seed across arrays.
rng('shuffle');
% Find general path (enclosing folder of current directory)
path{1} = strsplit(pwd, '/');
path{3,1} = strjoin(path{1}(1:end-1),'/');
path{4,1} = strjoin(path{1}, '/');
path{1,1} = strjoin(path{1}(1:end-2),'/');
path{2,1} = fullfile(path{1},'MATLAB');
% Set required subdirectories
path{5,1} = fullfile(path{3},'OCD','Results','EC');
path{6,1} = fullfile(path{3},'Atlases','AAL');
% Add relevant paths
addpath(genpath(fullfile(path{2}, 'BCT', 'NBS')));
addpath(fullfile(path{2},'spm12'));
addpath(genpath(fullfile(path{2}, 'mArrow3')));
addpath(genpath(fullfile(path{2}, 'permutationTest')));
addpath(genpath(fullfile(path{3}, 'Functions')));
addpath(genpath(fullfile(path{3}, 'LEICA', 'Functions')));
%% Settings
% Define files to load
fileName = 'LEICA90*.mat';
dirName = 'SubjectwithGroupPrior';
fileName = dir(fullfile(path{5}, dirName, fileName));
fileName = fileName.name;
% Analysis Settings
intercept = false; % determines whether to use intercept term in NBS
% Load network labels
label_AAL90 = load(fullfile(path{6},'AAL_labels.mat'));
label_AAL90 = string(label_AAL90.label90); % Convert labels to strings
label_AAL90 = strip(LR_version_symm(label_AAL90)); % Convert labels to symmetric format
% Generate MNI coordinates
coords_AAL90 = load(fullfile(path{6}, 'aal_cog.txt'), 'aal_cog');
coords_AAL90 = LR_version_symm(coords_AAL90);
origin = [65 45.5 35]; % center origin
MNIscale = 5.5/10; % scale MNI coordinates
sphereScale = 2.5; % scale sphere size
coords_AAL90 = MNIscale*coords_AAL90; % scale MNI coordinates
clear label90 MNIscale
% Set brain rendering parameters
cortex.file = fullfile(path{3},'OCD','Data','MNI152_T1_2mm_brain_mask.nii'); % file containing cortical atlas
cortex.color = [0.9 0.9 0.9]; % color for cortical rendering
cortex.transparency = 0.1; % set to 1 for opaque cortex
cortex.val = 0.3; % set isonormal line spacing
cortex.view = [-90 90]; % set camera angle
rdux = 0.7; % proportion of surface faces to keep (<1)
% set color index
cind = {[1 0 0], [0 0 1]; [1 0 1], [0 1 1]};
% Set number of brain images to generate
singleRender = 1;
%% Load and format EC data
% Load EC data
load(fullfile(path{5}, dirName, fileName), 'EC','I','N','condName','combs','mecDist','memberships','h','labels');
% Convert index format
ind = zeros(max(N.subjects), N.conditions);
n = 0;
for c = 1:N.conditions
ind(1:N.subjects(c),c) = I(n+(1:N.subjects(c)),c);
n = n + N.subjects(c);
end
% Set formatting indices
ind = logical(ind);
i = logical(I);
% Check if need intercept term
if intercept == true
I = horzcat(ones(size(I,1),1), I);
end
%% Run node strength analyses
% Declare strength array
stronk = nan(N.ROI, N.conditions, max(N.subjects), 2);
% Compute in-strength and out-strength
stronk(:,:,:,1) = squeeze(sum(EC, 2, 'omitnan')); % In-strength: sum over rows
stronk(:,:,:,2) = squeeze(sum(EC, 1, 'omitnan')); % Out-strength: sum over columns
% Run comparisions
strength.in = robustTests(stronk(:,1,:,1), stronk(:,2,:,1), N.ROI, 'p',0.05, 'testtype','permutation');
strength.out = robustTests(stronk(:,1,:,2), stronk(:,2,:,2), N.ROI, 'p',0.05, 'testtype','permutation');
% Tabulate FDR results
strength.summary = table(strength.in.FDR, strength.out.FDR, 'RowNames',label_AAL90, 'VariableNames',{'In','Out'});
% Format labels for strength results
instr = horzcat(squeeze(stronk(:,1,:,1))', squeeze(stronk(:,2,:,1))');
outstr = horzcat(squeeze(stronk(:,1,:,2))', squeeze(stronk(:,2,:,2))');
lbl = {repmat(label_AAL90, [2 1]), cat(1, repmat(string(labels.Properties.VariableNames{1}),[N.ROI 1]), repmat(string(labels.Properties.VariableNames{2}),[N.ROI 1]))};
cg = cat(1, repmat(["b";"r"],[N.ROI 1])); %, repmat("b",[N.ROI 1])); % cg = ["r", "b"];
[r, c] = find(strength.summary{:,:});
% Visualize strength results
S = figure('Position', [0 0 1280 1024]);
ax = subplot(3, nnz(strength.summary{:,:}), [1 nnz(strength.summary{:,:})]);
boxplot(ax, instr, lbl, 'PlotStyle','compact', 'Notch','on', 'ColorGroup',cg); hold on;
title('In-Strength');
scatter(r(c==1), (max(instr,[],'all','omitnan')+0.1).*ones(1,nnz(c==1)), 36, 'r', '*');
ax = subplot(3, nnz(strength.summary{:,:}), [1+nnz(strength.summary{:,:}) 2*nnz(strength.summary{:,:})]);
boxplot(ax, outstr, lbl, 'PlotStyle','compact', 'Notch','on', 'ColorGroup',cg); hold on;
title('Out-Strength');
scatter(r(c==2), (max(outstr,[],'all','omitnan')+0.1).*ones(1,nnz(c==2)), 36, 'r', '*');
for n = 1:nnz(strength.summary{:,:})
f = figure; hold on;
hg{1} = histogram(squeeze(stronk(r(n),1,:,c(n))), 'Normalization','probability');
hg{2} = histogram(squeeze(stronk(r(n),2,:,c(n))), 'Normalization','probability');
sz = min(hg{1}.BinWidth, hg{2}.BinWidth);
close(f);
figure(S);
ax = subplot(3, nnz(strength.summary{:,:}), 2*nnz(strength.summary{:,:})+n);
histogram(squeeze(stronk(r(n),1,:,c(n))), 'Normalization','Probability', 'BinWidth',sz, 'FaceAlpha',0.5); hold on;
histogram(squeeze(stronk(r(n),2,:,c(n))), 'Normalization','Probability', 'BinWidth',sz, 'FaceAlpha',0.5);
legend(labels.Properties.VariableNames);
title(['Modeled ', strength.summary.Properties.VariableNames{c(n)}, '-Strength of ', label_AAL90{r(n)}]);
end
clear n ax r c instr outstr lbl f cg
% Save as JPEG file
% saveas(S, fullfile(path{5}, dirName, strcat(fileName, '_strength')), 'jpeg');
%% Set parameter sweep for NBS
% Format EC data to NxNxM
C = EC; clear EC;
EC = zeros(N.ROI, N.ROI, sum(N.subjects));
for c = 1:N.conditions
EC(:,:,i(:,c)) = squeeze(C(:,:,c,ind(:,c)));
end
clear C c ind n i
% Set arrays for parameter sweeps
tstat = 3.5:0.5:6;
contrast = [-1,1; 1,-1];
strcont = {strjoin([labels.Properties.VariableNames(1), '<', labels.Properties.VariableNames(2)]), strjoin([labels.Properties.VariableNames(2), '<', labels.Properties.VariableNames(1)])};
% Set storage arrays
storarray = nan(length(tstat), size(contrast,1));
nbs = cell(length(tstat), size(contrast,1));
%% Analyze network using directed network-based statistic (NBS)
% Set directed NBS parameters
% C: connectivity matrices, arranged in a 3D (NxNxM) tensor
% GLM: structure defining the GLM for analyzing connectivity tensor
GLM.perms = 10000; % number of permutations (scalar)
GLM.X = I; % design matrix (observations x conditions)
GLM.test = 'ttest'; % type of statistical test to run ('ttest' or 'ftest')
STATS.size = 'Extent'; % method to measure size of cluster significance ('Intensity' or 'Extent')
STATS.alpha = 0.05; % significance level to test
% Sweep parameters
for c = 1:size(contrast,1)
GLM.contrast = contrast(c,:); % contrast vector (1 x conditions)
for t = 1:length(tstat)
STATS.thresh = tstat(t); % threshold test statistic
nbs{t, c} = NBSdirected(EC, GLM, STATS); % Run NBS
storarray(t, c) = ~isempty(nbs{t,c});
end
end
clear t c
%% Visualize NBS results
% locate significant components
j = h{:,'IC'};
if numel(j) > 1
i = j{1};
for k = 2:numel(j)
i = union(i, j{k});
end
i = unique(i);
else
i = j{:};
end
clear j
% Find average distance between groups
k(1) = combs(1,1); k(2) = combs(1,2);
scomb = nchoosek(1:sum(N.subjects(k)), 2);
scomb(scomb(:,2)<=N.subjects(k(1)),:) = [];
scomb(scomb(:,1)>N.subjects(k(1)),:) = [];
scomb(:,2) = scomb(:,2) - N.subjects(k(1));
d = nan(N.ROI, N.ROI, size(scomb,1));
for s = 1:size(scomb, 1)
d(:,:,s) = pdist2(squeeze(EC(:,:, scomb(s,1))), squeeze(EC(:,:, scomb(s,2))), mecDist);
end
[thresh, cont] = find(storarray);
if ~isempty(thresh) && singleRender
[thresh,j] = unique(thresh);
for t = 1:length(thresh)
map = nbs(thresh(t),:);
for m = 1:numel(map)
if iscell(map{m})
for n = 1:numel(map{m})
map{m}{n} = map{m}{n}.*mean(d, 3, 'omitnan')./10; % this would be an excellent place to apply recursive programming
end
else
map{m} = map{m}.*mean(d, 3, 'omitnan')./10;
end
end
F(t) = figure('Position', [0 0 1280 1024]); hold on;
plot_nodes_in_cortex(cortex, zscore(mean(memberships(:,i),2)), coords_AAL90, origin, sphereScale, [], map, cind(1,:), strcont, strength.summary, rdux);
end
elseif ~isempty(thresh)
imgs = cell(numel(unique(thresh)), numel(unique(cont)));
for t = 1:length(thresh)
ind = [thresh(t), cont(t)];
for f = 1:numel(nbs{ind(1), ind(2)})
F(f) = figure('Position', [0 0 1280 1024]);
% Plot significant connections
ax(1) = subplot(2, 4, [1 2]); colormap(ax(1),bone); hold on;
imagesc(full(nbs{ind(1), ind(2)}{f})); colorbar;
title(['NBS Network ', num2str(f)]);
yticks(1:N.ROI); yticklabels(label_AAL90);
xlim([1 N.ROI]); ylim([1 N.ROI]);
pbaspect([1 1 1]);
% Extract significant connections
[sconns(:,1), sconns(:,2)] = find(nbs{ind(1), ind(2)}{f});
% Highlight mean EC matrix
ax(2) = subplot(2, 4, [5 6]); colormap(ax(2),autumn); hold on;
imagesc(mean(d, 3, 'omitnan')); colorbar; hold on
scatter(sconns(:,2), sconns(:,1), 10, 'g', 's');
xlim([1 N.ROI]); ylim([1 N.ROI]);
title(['Mean Distance between ', condName{combs(1,1)}, ' and ', condName{combs(1,2)}, ' EC']);
yticks(1:N.ROI); yticklabels(label_AAL90); xticks([]);
pbaspect([1 1 1]);
% Render in SPM
map = nbs{ind(1),ind(2)}{f}.*mean(d, 3, 'omitnan')./10;
ax = subplot(2, 4, [3 4 7 8]); hold on
plot_nodes_in_cortex(cortex, zscore(mean(memberships(:,i),2)), coords_AAL90, origin, sphereScale, [], map, cind(1,:), strcont, strength.summary, rdux);
% plot_nodes_in_cortex(cortex, zscore(memberships(:,i)), coords_AAL90, origin, sphereScale, [], map, cind(1,:), strength.summary, rdux);
% plot_nodes_in_cortex(cortex, zscore(mean(memberships,2)), coords_AAL90, origin, sphereScale, [], map, rdux);
% Title figure
sgtitle(['Threshold ', num2str(tstat(ind(1))), ', ', strcont(ind(2),:)]);
% Save as JPEG file
% saveas(F(f), fullfile(path{5}, dirName, strcat(fileName(1:8), 'T', char(join(strsplit(num2str(tstat(ind(1))),'.'), '')), '_C', char(join(strsplit(num2str(contrast(ind(2),:))), '')), '_N', num2str(f))), 'jpeg');
clear sconns
end
% Place figure(s) in cell array
imgs{ind(1), ind(2)} = F;
% Save F
% savefig(F, fullfile(path{5}, dirName, strcat(fileName(1:8), 'T', char(join(strsplit(num2str(tstat(ind(1))),'.'), '')), '_C', char(join(strsplit(num2str(contrast(ind(2),:))), '')))));
clear F
end
end
clear K k n f t scomb ax ind nsig s thresh cont
%% Save results
% Save figure
if exist('F', 'var')
% savefig(F, fullfile(path{5}, dirName, strcat(fileName(1:end-4), '_NBS')), 'compact');
end
clear c C F
% Save results
% save(fullfile(path{5}, dirName, fileName(1:end-4)), 'strength','nbs','STATS','GLM','tstat','contrast','storarray', '-append');