1 | function [maps_out,fdrth,H] = fmri_analysis(DATA,varargin) |
---|
2 | % |
---|
3 | % This function does a quick fMRI analysis in SPM style. No motion |
---|
4 | % correction or slice timing at this time... |
---|
5 | % |
---|
6 | |
---|
7 | % Defaults |
---|
8 | TR = []; |
---|
9 | contr = [1]; |
---|
10 | onset = []; |
---|
11 | durat = []; |
---|
12 | smooth_kernel = []; |
---|
13 | d_cutoff = []; % Don't temporally filter data by default |
---|
14 | qFDR = 0.05; % Default threshold for FDR correction |
---|
15 | globalNorm = false; % Don't do global normalization by default |
---|
16 | omitVols = []; % Don't omit any volumes from the analysis by default |
---|
17 | show_wbar = true; % Show waitbar by default |
---|
18 | mask = []; % Don't mask by default |
---|
19 | regr = []; % Don't use additional regressors by default |
---|
20 | |
---|
21 | if rem(length(varargin),2)~=0 |
---|
22 | error('parameters/values missing!') |
---|
23 | end |
---|
24 | |
---|
25 | % Parse varargin |
---|
26 | for ii=1:2:length(varargin) |
---|
27 | param = lower(varargin{ii}); |
---|
28 | value = varargin{ii+1}; |
---|
29 | switch param |
---|
30 | case 'tr' |
---|
31 | TR = value; |
---|
32 | case {'onset','onsets'} |
---|
33 | onset = value; |
---|
34 | onset = onset(:); |
---|
35 | case {'durat','durations','duration'} |
---|
36 | durat = value; |
---|
37 | durat = durat(:); |
---|
38 | case 'smooth' |
---|
39 | smooth_kernel = value; |
---|
40 | case 'contrast' |
---|
41 | contr = value; |
---|
42 | case 'fdrth' |
---|
43 | qFDR = value; |
---|
44 | case {'hipass','hicutoff'} |
---|
45 | d_cutoff = value; |
---|
46 | case 'omitvolumes' |
---|
47 | omitVols = value; |
---|
48 | case 'wbar' |
---|
49 | if value==1 |
---|
50 | show_wbar = true; |
---|
51 | else |
---|
52 | show_wbar = false; |
---|
53 | end |
---|
54 | case 'mask' |
---|
55 | if isnumeric(value) || islogical(value) |
---|
56 | mask = value; |
---|
57 | else |
---|
58 | error('Invalid mask!') |
---|
59 | end |
---|
60 | case 'regressor' |
---|
61 | regr = value; |
---|
62 | end |
---|
63 | end |
---|
64 | |
---|
65 | % Check that we have something to fit... |
---|
66 | if isempty(onset) && isempty(regr) |
---|
67 | warning('No onsets or regressors defined. Only mean will be fitted!') |
---|
68 | end |
---|
69 | |
---|
70 | % Check that TR has been given |
---|
71 | if isempty(TR) |
---|
72 | error('Repetition time (TR) has not been defined!') |
---|
73 | end |
---|
74 | |
---|
75 | % Check onsets and durations |
---|
76 | if ~isempty(onset) |
---|
77 | if isempty(durat) |
---|
78 | error('Durations for onsets have not been defined!') |
---|
79 | end |
---|
80 | if length(durat)==1 |
---|
81 | durat = repmat(durat,1,length(onset)); |
---|
82 | durat = durat(:); |
---|
83 | elseif length(durat)~=length(onset) |
---|
84 | error('Mismatch in the number of elements in onset and duration!') |
---|
85 | end |
---|
86 | end |
---|
87 | |
---|
88 | % Check data |
---|
89 | if isstruct(DATA) && isfield(DATA,'FTDATA') |
---|
90 | data = DATA.FTDATA; |
---|
91 | elseif isnumeric(DATA) && ndims(DATA)>2 |
---|
92 | data=DATA; |
---|
93 | else |
---|
94 | error('Input data has to be a 3D/4D numeric matrix or Aedes data structure!') |
---|
95 | end |
---|
96 | |
---|
97 | % Permute 3D matrix to 4D |
---|
98 | if ndims(data)==3 |
---|
99 | data = permute(data,[1 2 4 3]); |
---|
100 | end |
---|
101 | |
---|
102 | % Initialize mask |
---|
103 | if isempty(mask) |
---|
104 | mask = true(size(data,1),size(data,2),size(data,3)); |
---|
105 | end |
---|
106 | |
---|
107 | % Check regressors |
---|
108 | if ~isempty(regr) |
---|
109 | if size(regr,1)~=size(data,4) |
---|
110 | error('The lengths of the regressors do not match data length!') |
---|
111 | end |
---|
112 | regr = detrend(regr); |
---|
113 | end |
---|
114 | |
---|
115 | % Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!! |
---|
116 | bf = [ |
---|
117 | 0 |
---|
118 | 0.000008876945146 |
---|
119 | 0.000387621003097 |
---|
120 | 0.003012441669394 |
---|
121 | 0.011548214895067 |
---|
122 | 0.030056522080168 |
---|
123 | 0.061233628208021 |
---|
124 | 0.105349995559045 |
---|
125 | 0.160158500633323 |
---|
126 | 0.221528149631163 |
---|
127 | 0.284405279593501 |
---|
128 | 0.343761753554314 |
---|
129 | 0.395323573307322 |
---|
130 | 0.436002661868442 |
---|
131 | 0.464045998242252 |
---|
132 | 0.478965423658889 |
---|
133 | 0.481327079129854 |
---|
134 | 0.472473786978796 |
---|
135 | 0.454237528221769 |
---|
136 | 0.428680153860941 |
---|
137 | 0.397883140401152 |
---|
138 | 0.363793656845973 |
---|
139 | 0.328124931280142 |
---|
140 | 0.292303457117353 |
---|
141 | 0.257453130446147 |
---|
142 | 0.224406054242129 |
---|
143 | 0.193730677047637 |
---|
144 | 0.165769521647024 |
---|
145 | 0.140680556701329 |
---|
146 | 0.118477987483032 |
---|
147 | 0.099069734765398 |
---|
148 | 0.082290068881472 |
---|
149 | 0.067926763712502 |
---|
150 | 0.055742762013807 |
---|
151 | 0.045492744488986 |
---|
152 | 0.036935220196777 |
---|
153 | 0.029840852217657 |
---|
154 | 0.023997740489123 |
---|
155 | 0.019214335904674 |
---|
156 | 0.015320580887648 |
---|
157 | 0.012167779438633 |
---|
158 | 0.009627606069476 |
---|
159 | 0.007590575511228 |
---|
160 | 0.005964217696074 |
---|
161 | 0.004671136982604 |
---|
162 | 0.003647081075751 |
---|
163 | 0.002839102788446 |
---|
164 | 0.002203865389816 |
---|
165 | 0.001706118264261 |
---|
166 | 0.001317352436400 |
---|
167 | 0.001014633776781 |
---|
168 | 0.000779604140824 |
---|
169 | 0.000597636251764 |
---|
170 | 0.000457125954290 |
---|
171 | 0.000348904854607 |
---|
172 | 0.000265756797153 |
---|
173 | 0.000202022712087 |
---|
174 | 0.000153279812313 |
---|
175 | 0.000116082720651 |
---|
176 | 0.000087755728843 |
---|
177 | 0.000066226941806 |
---|
178 | 0.000049896490408 |
---|
179 | 0.000037532277385 |
---|
180 | ]; |
---|
181 | |
---|
182 | fprintf(1,'\n******************************************\n'); |
---|
183 | fprintf(1,'Starting fMRI analysis.\n'); |
---|
184 | if isstruct(DATA) && isfield(DATA,'HDR') && isfield(DATA.HDR,'fpath') |
---|
185 | filename = [DATA.HDR.fpath,DATA.HDR.fname]; |
---|
186 | fprintf(1,'File: %s\n\n',filename) |
---|
187 | else |
---|
188 | fprintf(1,'\n'); |
---|
189 | end |
---|
190 | |
---|
191 | % Omit volumes from data and stimulus function if requested |
---|
192 | if ~isempty(omitVols) |
---|
193 | |
---|
194 | fprintf(1,'Skipping requested volumes...\n'); |
---|
195 | if ~isempty(onset) |
---|
196 | % Calculate new onsets and durations |
---|
197 | ton = onset; |
---|
198 | tof = onset+durat+1; |
---|
199 | tmp=zeros(1,size(data,4)); |
---|
200 | tmp(ton)=1;tmp(tof)=-1; |
---|
201 | sf=cumsum(tmp); |
---|
202 | sf(omitVols)=[]; |
---|
203 | tmp=diff([0 sf]); |
---|
204 | onset = find(tmp==1); |
---|
205 | durat = find(tmp==-1); |
---|
206 | if isempty(durat) |
---|
207 | % Block stays up... |
---|
208 | durat=length(tmp)-onset(end)-1; |
---|
209 | elseif length(durat) < length(onset) |
---|
210 | durat(1:end-1) = durat(1:end-1)-onset(1:end-1)-1; |
---|
211 | durat(end) = length(tmp)-onset(end)-1; |
---|
212 | else |
---|
213 | durat = durat-onset-1; |
---|
214 | end |
---|
215 | onset=onset(:); |
---|
216 | durat=durat(:); |
---|
217 | |
---|
218 | fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']); |
---|
219 | fprintf(1,['New duration: ',num2str(durat(:)'),'\n']); |
---|
220 | end |
---|
221 | |
---|
222 | if ~isempty(regr) |
---|
223 | regr(omitVols,:)=[]; |
---|
224 | end |
---|
225 | data(:,:,:,omitVols) = []; |
---|
226 | |
---|
227 | end |
---|
228 | |
---|
229 | % Create stimulus function (32 bin offset) |
---|
230 | k = size(data,4); % Number of scans |
---|
231 | T = 16; |
---|
232 | dt = TR/T; |
---|
233 | if ~isempty(onset) |
---|
234 | u = onset.^0; |
---|
235 | if ~any(durat) |
---|
236 | u = u/dt; |
---|
237 | end |
---|
238 | ton = round(onset*TR/dt) + 32; % onsets |
---|
239 | tof = round(durat*TR/dt) + ton + 1; % offset |
---|
240 | sf = zeros((k*T + 128),size(u,2)); |
---|
241 | ton = max(ton,1); |
---|
242 | tof = max(tof,1); |
---|
243 | for j = 1:length(ton) |
---|
244 | if numel(sf)>ton(j), |
---|
245 | sf(ton(j),:) = sf(ton(j),:) + u(j,:); |
---|
246 | end; |
---|
247 | if numel(sf)>tof(j), |
---|
248 | sf(tof(j),:) = sf(tof(j),:) - u(j,:); |
---|
249 | end; |
---|
250 | end |
---|
251 | sf = cumsum(sf); % integrate |
---|
252 | |
---|
253 | % Convolve stimulus with the HRF |
---|
254 | conv_sf = conv(sf,bf); |
---|
255 | |
---|
256 | % Resample the convolved stimulus |
---|
257 | RR = conv_sf([0:(k-1)]*16+1+32); |
---|
258 | RR=RR(:); |
---|
259 | else |
---|
260 | RR=[]; |
---|
261 | end |
---|
262 | |
---|
263 | |
---|
264 | % Smooth data if requested |
---|
265 | if ~isempty(smooth_kernel) && ~ismember(smooth_kernel,[1 1 1],'rows') |
---|
266 | if show_wbar |
---|
267 | wbh = aedes_calc_wait('Smoothing data...'); |
---|
268 | drawnow |
---|
269 | end |
---|
270 | if all(smooth_kernel) |
---|
271 | tic |
---|
272 | smooth_data = fmri_smooth(data,smooth_kernel); |
---|
273 | toc |
---|
274 | else |
---|
275 | fprintf(1,'fMRI analysis warning: Could not smooth data!\n'); |
---|
276 | smooth_data = data; |
---|
277 | end |
---|
278 | if show_wbar |
---|
279 | close(wbh); |
---|
280 | end |
---|
281 | else |
---|
282 | smooth_data = data; |
---|
283 | end |
---|
284 | |
---|
285 | |
---|
286 | if ~isempty(d_cutoff) |
---|
287 | % Create filtering matrix |
---|
288 | nCosine = ceil((2*k*TR)/(d_cutoff + 1)); |
---|
289 | S = sqrt(2/k)*cos([1:nCosine]'*pi*([1:k]/k)).'; |
---|
290 | KKT = eye(size(S,1))-2*S*S'+S*S'*S*S'; |
---|
291 | else |
---|
292 | KKT = eye(k); |
---|
293 | S = 0; |
---|
294 | end |
---|
295 | H = [RR regr ones(k,1)]; |
---|
296 | H = H-S*(S'*H); % Filter design matrix |
---|
297 | nParam = size(H,2); |
---|
298 | maps_out = struct('pmap',[],'tmap',[]); |
---|
299 | maps_out.pmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3),nParam); |
---|
300 | maps_out.tmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3)); |
---|
301 | |
---|
302 | % Calculate parametric map(s) |
---|
303 | nPlanes = size(smooth_data,3); |
---|
304 | nCols = size(smooth_data,2); |
---|
305 | nRows = size(smooth_data,1); |
---|
306 | if length(contr)<nParam |
---|
307 | contr(nParam)=0; % Pad with zeros |
---|
308 | end |
---|
309 | c = repmat(contr,nRows,1); |
---|
310 | HTH = pinv(H'*H); |
---|
311 | R = eye(size(H,1))-H*HTH*H'; |
---|
312 | |
---|
313 | if show_wbar |
---|
314 | wbh = aedes_wbar(0,sprintf('Estimating parameters. Processing plane 0/%d',nPlanes)); |
---|
315 | drawnow |
---|
316 | end |
---|
317 | |
---|
318 | % Process data in columns |
---|
319 | for ii=1:nPlanes |
---|
320 | |
---|
321 | %fprintf(1,'Processing plane %d/%d\n',ii,nPlanes); |
---|
322 | for kk=1:nCols |
---|
323 | if show_wbar |
---|
324 | aedes_wbar(ii/nPlanes,wbh,sprintf('Estimating parameters. Processing plane %d/%d, column %d/%d',ii,nPlanes,kk,nCols)); |
---|
325 | end |
---|
326 | col_data = squeeze(smooth_data(:,kk,ii,:)).'; |
---|
327 | col_data = col_data-S*(S'*col_data); |
---|
328 | th=H\col_data; |
---|
329 | r = col_data-H*th; |
---|
330 | rr=diag(r'*r); |
---|
331 | rr=rr(:); |
---|
332 | sig2 = rr./trace(R*KKT); |
---|
333 | sig2 = repmat(sig2,1,nParam); |
---|
334 | T = diag(c*th)./sqrt(c.*sig2*HTH*H'*KKT*H*HTH*contr'); |
---|
335 | T(find(mask(:,kk,ii)==0))=0; |
---|
336 | maps_out.pmap(:,kk,ii,:) = th.'; |
---|
337 | maps_out.tmap(:,kk,ii)=T.'; |
---|
338 | end |
---|
339 | end |
---|
340 | if show_wbar |
---|
341 | close(wbh); |
---|
342 | end |
---|
343 | |
---|
344 | if show_wbar |
---|
345 | wbh = aedes_calc_wait('Calculating threshold...'); |
---|
346 | drawnow |
---|
347 | end |
---|
348 | |
---|
349 | % Calculate effective degrees of freedom |
---|
350 | dof = (trace(R*KKT).^2)/trace(R*KKT*R*KKT); |
---|
351 | |
---|
352 | % p-values |
---|
353 | pval_map = 1-tdist(maps_out.tmap,dof); |
---|
354 | |
---|
355 | % Perform FDR (False Discovery Rate) correction |
---|
356 | cV = 1; |
---|
357 | |
---|
358 | pValues = pval_map(:); |
---|
359 | tValues = maps_out.tmap(:); |
---|
360 | |
---|
361 | [pValuesSorted,sortInd] = sort(pValues); |
---|
362 | tValuesSorted = tValues(sortInd); |
---|
363 | nP = length(pValues); |
---|
364 | |
---|
365 | pFDR = [1:nP]'/nP*qFDR/cV; % FDR-correction |
---|
366 | thresFDRind = find(pValuesSorted<=pFDR,1,'last'); |
---|
367 | if ~isempty(thresFDRind) |
---|
368 | thresFDR = tValuesSorted(thresFDRind); |
---|
369 | else |
---|
370 | thresFDR = []; |
---|
371 | end |
---|
372 | |
---|
373 | if show_wbar |
---|
374 | close(wbh); |
---|
375 | end |
---|
376 | |
---|
377 | if nargout>=2 |
---|
378 | fdrth = thresFDR; |
---|
379 | end |
---|
380 | |
---|
381 | if ~isempty(thresFDR) |
---|
382 | fprintf(1,['FDR threshold at p<',num2str(qFDR),': %.3f\n'],thresFDR) |
---|
383 | else |
---|
384 | fprintf(1,['No significant voxels over FDR threshold at p<',num2str(qFDR),'!\n']) |
---|
385 | end |
---|
386 | |
---|
387 | fprintf(1,'******************************************\n'); |
---|
388 | |
---|
389 | |
---|