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