forked from Masharipov/TMFC_toolbox
-
Notifications
You must be signed in to change notification settings - Fork 1
/
TMFC_command_window_example.m
370 lines (286 loc) · 14.2 KB
/
TMFC_command_window_example.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
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
clc
clear
close all
% BEFORE RUNNING THIS SCRIPT:
% 1) Set path to SPM12
% 2) Set path to TMFC_toolbox (Add with subfolders)
% 3) Change current working directory to: '...\TMFC_toolbox\examples'
cd(fileparts(matlab.desktop.editor.getActiveFilename)); % Set path to '...\TMFC_toolbox\examples'
%% Prepare example data and calculate basic first-level GLMs
data.SF = 1; % Scaling Factor (SF) for co-activations: SF = SD_oscill/SD_coact
data.SNR = 1; % Signal-to-noise ratio (SNR): SNR = SD_signal/SD_noise
data.STP_delay = 0.2; % Short-term synaptic plasticity (STP) delay, [s]
data.N = 20; % Sample size (Select 20 subjects out of 100 to reduce computations)
data.N_ROIs = 100; % Number of ROIs
data.dummy = 3; % Remove first M dummy scans
data.TR = 2; % Repetition time (TR), [s]
data.model = 'AR(1)'; % Autocorrelation modeling
% Set path for stat folder
spm_jobman('initcfg');
data.stat_path = spm_select(1,'dir','Select a folder for data extraction and statistical analysis');
% Set path for simulated BOLD time series *.mat file
data.sim_path = fullfile(pwd,'data','SIMULATED_BOLD_EVENT_RELATED_[2s_TR]_[1s_DUR]_[6s_ISI]_[40_TRIALS].mat');
% Set path for task design *.mat file (stimulus onset times, SOTs)
data.sots_path = fullfile(pwd,'data','TASK_DESIGN_EVENT_RELATED_[2s_TR]_[1s_DUR]_[6s_ISI]_[40_TRIALS].mat');
% Generate *.nii images and calculate GLMs
prepare_example_data(data)
% Change current directory to new TMFC project folder
cd(data.stat_path)
%% Setting up computation parameters
% Sequential or parallel computing (0 or 1)
tmfc.defaults.parallel = 1; % Parallel
% Store temporaty files during GLM estimation in RAM or on disk
tmfc.defaults.resmem = true; % RAM
% How much RAM can be used at the same time during GLM estimation
tmfc.defaults.maxmem = 2^32; % 4 GB
% Seed-to-voxel and ROI-to-ROI analyses
tmfc.defaults.analysis = 1;
%% Setting up paths
% The path where all results will be saved
tmfc.project_path = data.stat_path;
% Define paths to individual subject SPM.mat files
% tmfc.subjects(1).path = '...\Your_study\Subjects\sub_001\stat\Standard_GLM\SPM.mat';
% tmfc.subjects(2).path = '...\Your_study\Subjects\sub_002\stat\Standard_GLM\SPM.mat';
% tmfc.subjects(3).path = '...\Your_study\Subjects\sub_003\stat\Standard_GLM\SPM.mat';
% etc
% Alternativelly, use tmfc_select_subjects_GUI to select subjects
% Go to GLMs subfolder and select 20 subjects
SPM_check = 1; % Check SPM.mat files
[paths] = tmfc_select_subjects_GUI(SPM_check);
for i = 1:length(paths)
tmfc.subjects(i).path = paths{i};
end
clear SPM_check paths
%% Select ROIs
% Use tmfc_select_ROIs_GUI to select ROIs
%
% The tmfc_select_ROIs_GUI function creates group binary mask based on
% 1st-level masks (SPM.VM) and applies it to all selected ROIs. Empty ROIs
% will be removed. Masked ROIs will be limited to only voxels which have
% data for all subjects. The dimensions, orientation, and voxel sizes of
% the masked ROI images will be adjusted according to the group binary mask
%
% Go to ROI_masks subfolder and select 100 ROIs
[ROI_set] = tmfc_select_ROIs_GUI(tmfc);
tmfc.ROI_set(1) = ROI_set;
clear ROI_set
%% LSS regression
% Define conditions of interest
% tmfc.LSS.conditions(1).sess = 1;
% tmfc.LSS.conditions(1).number = 1;
% tmfc.LSS.conditions(2).sess = 1;
% tmfc.LSS.conditions(2).number = 2;
% Alternatively, use tmfc_LSS_GUI to select conditions of interest
[conditions] = tmfc_LSS_GUI(tmfc.subjects(1).path);
tmfc.LSS.conditions = conditions;
% Run LSS regression
start_sub = 1; % Start from the 1st subject
[sub_check] = tmfc_LSS(tmfc,start_sub);
clear conditions
%% BSC-LSS
% Extract and correlate mean beta series for conditions of interest
ROI_set_number = 1; % Select ROI set
[sub_check,contrasts] = tmfc_BSC(tmfc,ROI_set_number);
% Update contrasts info
% The tmfc_BSC function creates default contrasts for each
% condition of interest (i.e., Condition > Baseline)
tmfc.ROI_set(ROI_set_number).contrasts.BSC = contrasts;
% Define new contrasts:
tmfc.ROI_set(ROI_set_number).contrasts.BSC(3).title = 'TaskA_vs_TaskB';
tmfc.ROI_set(ROI_set_number).contrasts.BSC(4).title = 'TaskB_vs_TaskA';
tmfc.ROI_set(ROI_set_number).contrasts.BSC(3).weights = [1 -1];
tmfc.ROI_set(ROI_set_number).contrasts.BSC(4).weights = [-1 1];
% Calculate new contrasts
type = 3; % BSC-LSS
contrast_number = [3,4]; % Calculate contrasts #3 and #4
[sub_check] = tmfc_ROI_to_ROI_contrast(tmfc,type,contrast_number,ROI_set_number);
[sub_check] = tmfc_seed_to_voxel_contrast(tmfc,type,contrast_number,ROI_set_number);
% Load BSC-LSS matrices for the 'TaskA_vs_TaskB' contrast (contrast # 3)
for i = 1:data.N
M(i).paths = struct2array(load(fullfile(tmfc.project_path,'ROI_sets',tmfc.ROI_set(ROI_set_number).set_name,'BSC_LSS','ROI_to_ROI',...
['Subject_' num2str(i,'%04.f') '_Contrast_0003_[TaskA_vs_TaskB].mat'])));
end
matrices = cat(3,M(:).paths);
% Perform one-sample t-test (two-sided, FDR-correction)
contrast = 1; % A > B effect
alpha = 0.001/2; % alpha = 0.001 thredhold corrected for two-sided comparison
correction = 'FDR'; % False Discovery Rate (FDR) correction (Benjamini–Hochberg procedure)
[thresholded_1,pval,tval,conval_1] = tmfc_ttest(matrices,contrast,alpha,correction);
contrast = -1; % B > A effect
[thresholded_2] = tmfc_ttest(matrices,contrast,alpha,correction);
% Plot BSC-LSS results
f1 = figure(1); f1.Position = [382,422,1063,299];
try
sgtitle('BSC-LSS results');
catch
suptitle('BSC-LSS results');
end
subplot(1,3,1); imagesc(conval_1); title('Group mean'); axis square; colorbar; caxis(tmfc_axis(conval_1,1));
subplot(1,3,2); imagesc(thresholded_1); title('A>B (pFDR<0.001)'); axis square; colorbar;
subplot(1,3,3); imagesc(thresholded_2); title('B>A (pFDR<0.001)'); axis square; colorbar;
colormap(subplot(1,3,2),'parula')
colormap(subplot(1,3,3),'parula')
colormap(subplot(1,3,1),'redblue')
set(findall(gcf,'-property','FontSize'),'FontSize',16)
clear type contrasts contrast_number
%% FIR task regression (regress out co-activations and save residual time series)
% FIR window length in [s]
tmfc.FIR.window = 24;
% Number of FIR time bins
tmfc.FIR.bins = 24;
% Run FIR task regression
[sub_check] = tmfc_FIR(tmfc,start_sub);
%% LSS regression after FIR task regression (use residual time series)
% Define conditions of interest
tmfc.LSS_after_FIR.conditions = tmfc.LSS.conditions;
% Run LSS regression
[sub_check] = tmfc_LSS_after_FIR(tmfc,start_sub);
%% BSC-LSS after FIR task regression (use residual time series)
% Extract and correlate mean beta series for conditions of interest
ROI_set_number = 1; % Select ROI set
[sub_check,contrasts] = tmfc_BSC_after_FIR(tmfc,ROI_set_number);
% Update contrasts info
% The tmfc_BSC_after_FIR function creates default contrasts for each
% condition of interest (i.e., Condition > Baseline)
tmfc.ROI_set(ROI_set_number).contrasts.BSC_after_FIR = contrasts;
% Define new contrast:
tmfc.ROI_set(ROI_set_number).contrasts.BSC_after_FIR(3).title = 'TaskA_vs_TaskB';
tmfc.ROI_set(ROI_set_number).contrasts.BSC_after_FIR(4).title = 'TaskB_vs_TaskA';
tmfc.ROI_set(ROI_set_number).contrasts.BSC_after_FIR(3).weights = [1 -1];
tmfc.ROI_set(ROI_set_number).contrasts.BSC_after_FIR(4).weights = [-1 1];
% Calculate new contrast
type = 4; % BSC-LSS after FIR
contrast_number = [3,4]; % Calculate contrast #3 and #4
[sub_check] = tmfc_ROI_to_ROI_contrast(tmfc,type,contrast_number,ROI_set_number);
[sub_check] = tmfc_seed_to_voxel_contrast(tmfc,type,contrast_number,ROI_set_number);
% Load BSC-LSS (after FIR) matrices for the 'TaskA_vs_TaskB' contrast (contrast # 3)
clear M marices conval_1 thresholded_1
for i = 1:data.N
M(i).paths = struct2array(load(fullfile(tmfc.project_path,'ROI_sets',tmfc.ROI_set(ROI_set_number).set_name,'BSC_LSS_after_FIR','ROI_to_ROI',...
['Subject_' num2str(i,'%04.f') '_Contrast_0003_[TaskA_vs_TaskB].mat'])));
end
matrices = cat(3,M(:).paths);
% Perform one-sample t-test (two-sided, FDR-correction)
contrast = 1; % A > B effect
[thresholded_1,pval,tval,conval_1] = tmfc_ttest(matrices,contrast,alpha,correction);
contrast = -1; % B > A effect
[thresholded_2] = tmfc_ttest(matrices,contrast,alpha,correction);
% Plot BSC-LSS (after FIR) results
f2 = figure(2); f2.Position = [382,422,1063,299];
try
sgtitle('BSC-LSS (after FIR task regression) results');
catch
suptitle('BSC-LSS (after FIR task regression) results');
end
subplot(1,3,1); imagesc(conval_1); title('Group mean'); axis square; colorbar; caxis(tmfc_axis(conval_1,1));
subplot(1,3,2); imagesc(thresholded_1); title('A>B (pFDR<0.001)'); axis square; colorbar;
subplot(1,3,3); imagesc(thresholded_2); title('B>A (pFDR<0.001)'); axis square; colorbar;
colormap(subplot(1,3,2),'parula')
colormap(subplot(1,3,3),'parula')
colormap(subplot(1,3,1),'redblue')
set(findall(gcf,'-property','FontSize'),'FontSize',16)
clear type contrasts contrast_number
%% BGFC
% Calculate background functional connectivity (BGFC)
[sub_check] = tmfc_BGFC(tmfc,ROI_set_number,start_sub);
%% gPPI
% Define conditions of interest
[conditions] = tmfc_gPPI_GUI(tmfc.subjects(1).path);
tmfc.ROI_set(ROI_set_number).gPPI.conditions = conditions;
clear conditions
% VOI extraction
[sub_check] = tmfc_VOI(tmfc,ROI_set_number,start_sub);
% PPI calculation
[sub_check] = tmfc_PPI(tmfc,ROI_set_number,start_sub);
% gPPI calculation
[sub_check,contrasts] = tmfc_gPPI(tmfc,ROI_set_number,start_sub);
% Update contrasts info
% The tmfc_gPPI function creates default contrasts for each
% condition of interest (i.e., Condition > Baseline)
tmfc.ROI_set(ROI_set_number).contrasts.gPPI = contrasts;
% Define new contrasts:
tmfc.ROI_set(ROI_set_number).contrasts.gPPI(3).title = 'TaskA_vs_TaskB';
tmfc.ROI_set(ROI_set_number).contrasts.gPPI(4).title = 'TaskB_vs_TaskA';
tmfc.ROI_set(ROI_set_number).contrasts.gPPI(3).weights = [1 -1];
tmfc.ROI_set(ROI_set_number).contrasts.gPPI(4).weights = [-1 1];
% Calculate new contrasts
type = 1; % gPPI
contrast_number = [3,4]; % Calculate contrasts #3 and #4
[sub_check] = tmfc_ROI_to_ROI_contrast(tmfc,type,contrast_number,ROI_set_number);
[sub_check] = tmfc_seed_to_voxel_contrast(tmfc,type,contrast_number,ROI_set_number);
% Load gPPI matrices for the 'TaskA_vs_TaskB' contrast (contrast # 3)
clear M marices conval_1 thresholded_1
for i = 1:data.N
M(i).paths = struct2array(load(fullfile(tmfc.project_path,'ROI_sets',tmfc.ROI_set(ROI_set_number).set_name,'gPPI','ROI_to_ROI','symmetrical',...
['Subject_' num2str(i,'%04.f') '_Contrast_0003_[TaskA_vs_TaskB].mat'])));
end
matrices = cat(3,M(:).paths);
% Perform one-sample t-test (two-sided, FDR-correction)
contrast = 1; % A > B effect
[thresholded_1,pval,tval,conval_1] = tmfc_ttest(matrices,contrast,alpha,correction);
contrast = -1; % B > A effect
[thresholded_2] = tmfc_ttest(matrices,contrast,alpha,correction);
% Plot gPPI results
f3 = figure(3); f3.Position = [382,422,1063,299];
try
sgtitle('gPPI results');
catch
suptitle('gPPI results');
end
subplot(1,3,1); imagesc(conval_1); title('Group mean'); axis square; colorbar; caxis(tmfc_axis(conval_1,1));
subplot(1,3,2); imagesc(thresholded_1); title('A>B (pFDR<0.001)'); axis square; colorbar;
subplot(1,3,3); imagesc(thresholded_2); title('B>A (pFDR<0.001)'); axis square; colorbar;
colormap(subplot(1,3,2),'parula')
colormap(subplot(1,3,3),'parula')
colormap(subplot(1,3,1),'redblue')
set(findall(gcf,'-property','FontSize'),'FontSize',16)
clear type contrasts contrast_number
%% gPPI-FIR (gPPI model with psychological regressors defined by FIR functions)
% Define FIR parameters for gPPI-FIR
tmfc.ROI_set(ROI_set_number).gPPI_FIR.window = 24; % FIR window length in [s]
tmfc.ROI_set(ROI_set_number).gPPI_FIR.bins = 24; % Number of FIR time bins
% gPPI-FIR calculation
[sub_check,contrasts] = tmfc_gPPI_FIR(tmfc,ROI_set_number,start_sub);
% Update contrasts info
% The tmfc_gPPI_FIR function creates default contrasts for each
% condition of interest (i.e., Condition > Baseline)
tmfc.ROI_set(ROI_set_number).contrasts.gPPI_FIR = contrasts;
% Define new contrasts:
tmfc.ROI_set(ROI_set_number).contrasts.gPPI_FIR(3).title = 'TaskA_vs_TaskB';
tmfc.ROI_set(ROI_set_number).contrasts.gPPI_FIR(4).title = 'TaskB_vs_TaskA';
tmfc.ROI_set(ROI_set_number).contrasts.gPPI_FIR(3).weights = [1 -1];
tmfc.ROI_set(ROI_set_number).contrasts.gPPI_FIR(4).weights = [-1 1];
% Calculate new contrasts
type = 2; % gPPI-FIR
contrast_number = [3,4]; % Calculate contrasts #3 and #4
[sub_check] = tmfc_ROI_to_ROI_contrast(tmfc,type,contrast_number,ROI_set_number);
[sub_check] = tmfc_seed_to_voxel_contrast(tmfc,type,contrast_number,ROI_set_number);
% Load gPPI-FIR matrices for the 'TaskA_vs_TaskB' contrast (contrast # 3)
clear M marices conval_1 thresholded_1
for i = 1:data.N
M(i).paths = struct2array(load(fullfile(tmfc.project_path,'ROI_sets',tmfc.ROI_set(ROI_set_number).set_name,'gPPI_FIR','ROI_to_ROI','symmetrical',...
['Subject_' num2str(i,'%04.f') '_Contrast_0003_[TaskA_vs_TaskB].mat'])));
end
matrices = cat(3,M(:).paths);
% Perform one-sample t-test (two-sided, FDR-correction)
contrast = 1; % A > B effect
[thresholded_1,pval,tval,conval_1] = tmfc_ttest(matrices,contrast,alpha,correction);
contrast = -1; % B > A effect
[thresholded_2] = tmfc_ttest(matrices,contrast,alpha,correction);
% Plot gPPI-FIR results
f4 = figure(4); f4.Position = [382,422,1063,299];
try
sgtitle('gPPI-FIR results');
catch
suptitle('gPPI-FIR results');
end
subplot(1,3,1); imagesc(conval_1); title('Group mean'); axis square; colorbar; caxis(tmfc_axis(conval_1,1));
subplot(1,3,2); imagesc(thresholded_1); title('A>B (pFDR<0.001)'); axis square; colorbar;
subplot(1,3,3); imagesc(thresholded_2); title('B>A (pFDR<0.001)'); axis square; colorbar;
colormap(subplot(1,3,2),'parula')
colormap(subplot(1,3,3),'parula')
colormap(subplot(1,3,1),'redblue')
set(findall(gcf,'-property','FontSize'),'FontSize',16)
clear type contrasts contrast_number
%% Save TMFC project *.mat file
save(fullfile(data.stat_path,'TMFC_project.mat'),'tmfc');