-
Notifications
You must be signed in to change notification settings - Fork 1
/
mprf_main.m
226 lines (185 loc) · 9.14 KB
/
mprf_main.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
function results = mprf_main(subjID, opt)
% Wrapper script containing MEG and MRI preprocessing and analyses subfunctions
% involved in the MEG Retinotopy project.
%
% %%% WORKFLOW %%%
% 0. Load paths and define parameters
%
% 1. MEG data preprocessing:
% 1.0: Define preprocessing options (if not defined as input variable)
% 1.1: Preprocess MEG data from raw (if requested)
% 1.2: Load MEG stim
% 1.3: Load MEG Gain matrix
%
% 2. MRI data preprocessing: (if requested)
% 2.1: Get pRF parameters in mrVista voxel space and smooth parameters
% 2.2: Project pRF parameters to freesurfer anatomy, get visual ROIs
% 2.3: Project pRF parameters from Freesurfer to Brainstorm surface
% (2.4): if perturbing pRF parameters on surface, we do so here.
%
% 3. Forward model:
% 3.1: Predict response for MEG stimulus at cortical surface level
% (could be Brainstorm or FreeSurfer surface)
% 3.2: Predict response for MEG stimulus at MEG sensor level
% (multiply with gain matrix)
% 3.3: Computing phase referenced amplitude from preprocessed MEG data
% and predicted MEG responses
% 3.4: Comparing predicted MEG time series and phase-referenced MEG
% steady-state responses
%
% %%% DEPENDENCIES %%%
% A. MRI preprocessing steps:
% A1: FreeSurfer's auto-segmentation of structural MRI data (v5.3)
% A2: Preprocessing of raw functional MRI data into preprocessed nifti's (i.e. MRI
% distortion correction w/ fsl top-up)
% A3: Running pRF model in mrVista (voxel space) to get Gray retinotopy modelfits
% A4: Apply probabilistic visual ROI atals by Wang et al. (2015) to
% subject's anatomy in FreeSurfer directory
%
% B. External MATLAB Toolboxes:
% - Brainstorm (https://github.com/brainstorm-tools/brainstorm3)
% last version commit 49a4f9b6 (June 24, 2021)
% - FieldTrip (github.com/fieldtrip/fieldtrip)
% last version commit 5cf4cf4ce (March 9, 2020)
% - VistaSoft (github.com/vistalab/vistasoft)
% last run version commit d6d0207 (March 17, 2021)
% - meg_utils (github.com/WinawerLab/meg_utils)
% last run version commit f0c35f (December 28, 2020)
%
% Add all to path with the ToolboxToolbox:
% tbUse('retmeg')
%
% Written by Eline Kupers (NYU) and Akhil Edadan (UU) - 2019-2021
%
%% 0. Load paths
% Set options if not defined (see getOpts function for more options)
if ~exist('opt', 'var') || isempty(opt)
opt = getOpts('saveFig', true,'verbose', true, 'fullSizeMesh', true, ...
'perturbOrigPRFs', false, 'addOffsetParam', false, ...
'refitGainParam', false);
end
% Load paths with data files for this subject
dirPth = loadPaths(subjID);
% Go back to root
cd(mprf_rootPath)
fprintf('(%s): Starting analysis of subject %s, using %s\n', mfilename, subjID, regexprep(opt.subfolder,'/',' '));
%% 1. MEG data preprocessing
% Start MEG data processing (if requested)
if opt.doMEGPreproc
if opt.verbose; fprintf('(%s): Preprocess MEG data..\n', mfilename); end
% 1.1 Get preprocessed data from raw MEG data (.sqd) to preprocessed MEG data
% (matfile, MEG sensors x epochs x time points x run nr)
[data, conditions, opt] = preprocessMEGRetinotopyData(subjID, dirPth, opt);
data = data.data;
% 1.2 Get MEG stimulus (binarized and reduced to epochs x 10201 pixels)
stim = loadStim(subjID, dirPth, opt);
stim = stim.MEG;
% 1.3 Get Gain matrix (produced via brainstorm)
gainMtx = loadGainMtx(subjID, dirPth, opt);
else % If you want to skip preprocessing, load data structs
load(fullfile(dirPth.meg.processedDataPth, 'epoched_data_hp_preproc_denoised.mat'), 'data');
load(fullfile(dirPth.meg.processedDataPth, 'megStimConditions.mat'), 'triggers');
load(fullfile(dirPth.meg.processedDataPth, 'meg_stimulus.mat'), 'meg_stim');
gainMtx = loadGainMtx(subjID, dirPth, opt);
% remove similar-named fields
data = data.data;
stim = meg_stim;
conditions = triggers;
end
% Get MEG data struct and fill with loaded data
meg = struct();
meg.data = data;
meg.stim = stim;
meg.stim.conditions = conditions;
meg.gain = gainMtx;
% Save some memory
clear gainMtx data stim conditions
%% 2. MRI data preprocessing
% This section will:
% 2.1 Smoothing pRF params in voxel space (and recompute beta's)
% 2.2 Export pRF params to FreeSurfer mesh + get Wang et al visual rois
% 2.3 Export pRF to Brainstorm mesh (if opt.fullSizeMesh = false)
% Get MRI data struct
mri = struct();
% What pRF retinotopy maps are used on what size mesh?
if opt.mri.useBensonMaps % Use benson maps or not
[mri.data, mri.prfSurfPath] = loadBensonRetinotopyMaps(subjID, dirPth, opt);
elseif opt.mri.useNYU3TAveMaps
mri.prfSurfPath = loadNYU3TRetinotopyMaps(subjID, dirPth, opt);
elseif opt.fullSizeMesh % if using full gain matrix, we need pRF params on FreeSurfer surface
mri.prfSurfPath = dirPth.fmri.saveDataPth_prfFS;
else % if using downsampled gain matrix, we need pRF params on Brainstorm surface
mri.prfSurfPath = dirPth.fmri.saveDataPth_prfBS;
end
% Start MRI data processing (if requested)
if opt.doMRIPreproc && ~loadBensonRetinotopyMaps && ~opt.mri.useNYU3TAveMaps
if opt.verbose; fprintf('(%s): Preprocess MRI data..\n', mfilename); end
% 2.1 Smoothing pRF params in voxel space (and then recompute beta
% values with smoothed parameters). We use mrVISTA for this.
% mrVista gray nodes (i.e. voxels) --> mrVista gray nodes (i.e. voxels)
mprf_pRF_sm(dirPth, opt);
% 2.2 Get smoothed pRF params on FreeSurfer mid gray surface.
% --> mrVista gray nodes (i.e. voxels) to FreeSurfer vertices
mprf_pRF_sm_FS(dirPth,opt);
% 2.3 Get smoothed pRF params and ROIs on Brainstorm pial surface.
% --> Freesurfer vertices to downsampled Brainstorm vertices
mprf_pRF_sm_FS_BS(subjID, dirPth,opt);
if opt.verbose
% Get summary figures for 2.1 pRF parameters before/after smoothing
mprf_pRF_sm_fig(dirPth, opt);
close all;
% Get summary figures for 2.2 pRF parameters on FreeSurfer surface
mprf_pRF_sm_FS_fig(dirPth,opt);
close all;
% Get summary figures for 2.3 pRF parameters on Brainstorm surface
mprf_pRF_sm_FS_BS_fig(dirPth,opt);
close all;
end
end
% 2.4 If perturbing original pRF parameters on the cortical surface:
if opt.vary.perturbOrigPRFs
if opt.verbose; fprintf('(%s): Perturb local pRFs on cortex..\n', mfilename); end
mprf_perturbOrigPRFs(mri.prfSurfPath, opt)
end
%% 3. Forward model
% 3.1 Predict response to MEG stimulus on mesh vertices (could be BS or FS)
% INPUTS (1) path to pRF parameters on surface (string)
% (2) MEG stimulus (struct with x, y, images, etc)
% OUTPUTS (1) predicted responses on surface (epochs x vertices)
predSurfResponse = mprf_MEGPredictionFromSurfaceWrapper(mri.prfSurfPath, meg.stim, dirPth, opt);
results.predSurfResponse = predSurfResponse;
%% 3.2 Predicted response for MEG stimulus at MEG sensor level (weighting
% predicted surface responses with gain matrix)
% INPUTS (1) predicted surface responses (on BS or FS surface)
% (epochs x vertices)
% (2) gain matrix (sensors x vertices)
%
% OUTPUTS (1) predicted MEG responses (epochs x sensors)
predMEGResponse = mprf_MEGPredictionSensorsWrapper(predSurfResponse, meg.gain, dirPth, opt);
results.predMEGResponse = predMEGResponse;
%% 3.3 Computing phase referenced amplitude from preprocessed MEG data
% and predicted MEG responses from cortical surface
% INPUTS (1) preprocessed MEG data (time x epochs x run x sensors)
% (2) predicted MEG responses (epochs x sensors)
% OUTPUTS (1) phase-referenced MEG time series
% (sensors x run group x epochs)
% (2) bestBetas (1 x run group x sensors)
% (3) bestRefPhase (1 x run group x sensors)
% (4) bestOffsets (1 x run group x sensors)
[phaseRefMEGResponse, bestBetas, bestRefPhase, bestOffsets] = mprf_MEGPhaseReferenceDataWrapper(meg.data, predMEGResponse, dirPth, opt);
results.phaseRefMEGResponse = phaseRefMEGResponse;
results.bestBetas = bestBetas;
results.bestRefPhase = bestRefPhase;
results.bestOffsets = bestOffsets;
%% 3.4 Comparing predicted MEG time series and phase-referenced MEG steady-state responses
% INPUTS (1) Phase referenced MEG time series (sensors x run groups x epochs)
% (2) predicted MEG sensor responses from MRI prfs
% (epochs x sensors)
% OUTPUTS (1) modelfit to mean phase-referenced MEG data, scaled by
% gain (and if requested, with offset) (epochs x sensors)
% (2) average variance explained per MEG sensor (1 x sensors)
[predMEGResponseScaled,meanVarExpl] = mprf_CompareMEGDataToPRFPredictionWrapper(phaseRefMEGResponse, predMEGResponse, bestBetas, bestOffsets, dirPth, opt);
results.predMEGResponseScaled = predMEGResponseScaled;
results.meanVarExpl = meanVarExpl;
fprintf('(%s) finished without error!\n', mfilename)
end