% Looking at the fMRI data using pca_image input_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; mask_file=input_file; mask_thresh=fmri_mask_thresh(mask_file); pca_image(input_file, [], 4, mask_file, mask_thresh); colormap(spectral); exclude=[1 2]; pca_image(input_file, exclude, 4, mask_file, mask_thresh); X_remove=[ones(120,1) (1:120)']; X_interest=repmat(eye(12),10,1); pca_image(input_file, exclude, 4, mask_file, mask_thresh, X_remove, X_interest); % Plotting the hemodynamic response function (hrf) using fmridesign hrf_parameters=[5.4 5.2 10.8 7.35 0.35] time=(0:240)/10; hrf0=fmridesign(time,0,[1 0],[],hrf_parameters); plot(time,squeeze(hrf0.X(:,1,1)),'LineWidth',2) xlabel('time (seconds)') ylabel('hrf') title('Glover hrf model') % Making the design matrices using fmridesign frametimes=(0:119)*3; slicetimes=[0.14 0.98 0.26 1.10 0.38 1.22 0.50 1.34 0.62 1.46 0.74 1.58 0.86]; eventid=kron(ones(10,1),[1; 2]); eventimes=(0:19)'*18+9; duration=ones(20,1)*9; height=ones(20,1); events=[eventid eventimes duration height] S=kron(ones(10,1),kron([0 0; 1 0; 0 0; 0 1],ones(3,1))); X_cache=fmridesign(frametimes,slicetimes,events,[],hrf_parameters); X_cache=fmridesign(frametimes,slicetimes, [] , S,hrf_parameters); plot(squeeze(X_cache.X(:,:,1,4)),'LineWidth',2) legend('Hot','Warm') xlabel('frame number') ylabel('response') title('Hot and Warm responses') % Analysing one run with fmrilm contrast=[1 0; 0 1; 1 -1]; exclude=[1 2]; which_stats=[1 1 1 1 1 0 0 0 1]; input_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['g:/keith/results/ha_100326_hot'; 'g:/keith/results/ha_100326_wrm'; 'g:/keith/results/ha_100326_hmw'] [df1 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); % Visualizing the results using glass_brain t_file='g:/keith/results/ha_100326_hmw_mag_t.mnc'; m=fmris_read_image(t_file,4,1); imagesc(m.data',[-6 6]); colorbar; axis xy; colormap(spectral); mask_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; mask_thresh=fmri_mask_thresh(mask_file); a=fmris_read_image(mask_file,4,1); imagesc((m.data.*(a.data>mask_thresh))',[-6 6]); colorbar; axis xy; glass_brain(t_file,3,mask_file,mask_thresh); % F-tests contrast=[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1] which_stats=[0 0 0 1]; output_file_base='g:/keith/results/ha_100326_drift'; fwhm_rho='g:/keith/results/ha_100326_hot_rho.mnc'; fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats,fwhm_rho); clf; m=fmris_read_image('g:/keith/results/ha_100326_drift_mag_F.mnc',4,1); imagesc(m.data'); colorbar; axis xy; colormap(spectral); % A linear effect of temperature temperature=[45.5 35.0 49.0 38.5 42.0 49.0 35.0 42.0 38.5 45.5 ... 38.5 49.0 35.0 45.5 42.0 45.5 38.5 42.0 35.0 49.0]'; Temperature=35:3.5:49; subplot(2,2,1) plot(Temperature,1:5,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(a) Linear temperature effect'); subplot(2,2,2) plot(Temperature,10-((1:5)-4).^2,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(b) Quadratic temperature effect'); subplot(2,2,3) plot(Temperature,10-((1:5)-4).^2+((1:5)-3).^3,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(c) Cubic temperature effect'); subplot(2,2,4) plot(Temperature,[1 4 3 5 3.5],'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title(['(d) Quartic = arbitrary temperature effect']); events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature] contrast=[0 1; 1 49; 1 35; 0 14]; events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature; zeros(20,1)+3 eventimes duration temperature.^2] contrast=[0 0 1]; contrast=[1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1]; contrast=[eye(5) ones(5,4)]; which_stats=[0 0 0 1]; contrast=[.8 -.2 -.2 -.2 -.2; -.2 .8 -.2 -.2 -.2; -.2 -.2 .8 -.2 -.2; -.2 -.2 -.2 .8 -.2; -.2 -.2 -.2 -.2 .8]; contrast=[eye(5)-ones(5)/5 ones(5,4)]; which_stats=[0 0 0 1]; contrast=[-7.0 -3.5 0 3.5 7.0]; which_stats=[1 1 1]; % Combining runs/sessions/subjects with multistat contrast=[1 -1]; which_stats=[1 1 1]; input_file='g:/keith/data/brian_ha_19971124_1_093923_mri_MC.mnc'; output_file_base='g:/keith/results/ha_093923_hmw'; [df2 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='g:/keith/data/brian_ha_19971124_1_101410_mri_MC.mnc'; output_file_base='g:/keith/results/ha_101410_hmw'; [df3 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='g:/keith/data/brian_ha_19971124_1_102703_mri_MC.mnc'; output_file_base='g:/keith/results/ha_102703_hmw'; [df4 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); X=[1 1 1 1]'; contrast=[1]; which_stats=[1 1 1]; input_files_effect=['g:/keith/results/ha_093923_hmw_mag_ef.mnc'; 'g:/keith/results/ha_100326_hmw_mag_ef.mnc'; 'g:/keith/results/ha_101410_hmw_mag_ef.mnc'; 'g:/keith/results/ha_102703_hmw_mag_ef.mnc']; input_files_sdeffect=['g:/keith/results/ha_093923_hmw_mag_sd.mnc'; 'g:/keith/results/ha_100326_hmw_mag_sd.mnc'; 'g:/keith/results/ha_101410_hmw_mag_sd.mnc'; 'g:/keith/results/ha_102703_hmw_mag_sd.mnc']; output_file_base='g:/keith/results/ha_multi_hmw_fixed' input_files_df=[df1 df2 df3 df4] clf; m=fmris_read_image('g:/keith/results/ha_100326_hot_fwhm.mnc',4,1); imagesc(m.data',[0 20]); colorbar; axis xy; input_files_fwhm=8; df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,Inf) m=fmris_read_image('g:/keith/results/ha_multi_hmw_fixed_t.mnc',4,1); imagesc(m.data',[-6 6]); colorbar; axis xy; colormap(spectral); % Fixed and random effects which_stats=[1 1 1 1 0 0 0 0 1]; output_file_base='g:/keith/results/ha_multi_hmw' [df, df_resid]=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,20) m=fmris_read_image('g:/keith/results/ha_multi_hmw_sdratio.mnc',4,1); imagesc(m.data'); colorbar; axis xy; m=fmris_read_image('g:/keith/results/ha_multi_hmw_fwhm.mnc',4,1); blured_fwhm=gauss_blur('g:/keith/results/ha_multi_hmw_fwhm.mnc',10); mb=fmris_read_image(blured_fwhm,4,1); imagesc([m.data' mb.data'],[0 20]); colorbar; axis xy; colormap(spectral); % Thresholding the tstat image with stat_threshold and fdr_threshold stat_threshold(1000000,26000,8,112,0.05) stat_threshold(1000000,26000,8,112,0.05,3) stat_threshold(1000000,26000,8,[112; 112],0.05,3) stat_threshold(1000000,26000,8,[112 0; 3 112],0.05,3) mask_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; input_file='g:/keith/results/ha_multi_hmw_t.mnc'; fdr_threshold(input_file,[],mask_file,mask_thresh,112,0.05) % Locating peaks and clusters with locmax stat_file='g:/keith/results/ha_multi_hmw_t.mnc'; mask_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; lm=locmax(stat_file, 4.86, mask_file,mask_thresh) pval=stat_threshold(1000000,26000,8,112,lm(:,1)) % Producing an SPM style summary with stat_summary fwhm_file='g:/keith/results/ha_multi_hmw_fwhm.mnc'; summmary=stat_summary(stat_file, fwhm_file, [112 0; 3 112], ... mask_file, mask_thresh, 0.001); % Confidence regions for the spatial location of local maxima using conf_region conf_region(stat_file,4.86,mask_file,mask_thresh) glass_brain('g:/keith/results/ha_multi_hmw_t_95cr.mnc', ... 4.86,mask_file,mask_thresh); % Conjunctions which_stats=[0 0 0 0 1]; multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,20) stat_summary('g:/keith/results/ha_multi_hmw_conj.mnc', fwhm_file, ... [112 0; 3 112], mask_file, mask_thresh, 0.001, 1, 4); stat_threshold(1000000,26000,8,112,0.05,0.001,0.05,4) % Extracting values from a minc file using extract voxel=[62 66 4] effect=extract(voxel,'g:/keith/results/ha_multi_hmw_ef.mnc') sdeffect=extract(voxel,'g:/keith/results/ha_multi_hmw_sd.mnc') ref_data=squeeze(extract(voxel,'g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc')); ef_hot=extract(voxel,'g:/keith/results/ha_100326_hot_mag_ef.mnc') ef_wrm=extract(voxel,'g:/keith/results/ha_100326_wrm_mag_ef.mnc') fitted=mean(ref_data)+ef_hot*X_cache(:,1,1,5)+ef_wrm*X_cache(:,2,1,5); clf; plot(frametimes,[ref_data fitted],'LineWidth',2); ylim([930 960]); legend('Reference data','Fitted values'); xlabel('time (seconds) from start of epoch'); ylabel('fMRI response'); % Estimating the time course of the response eventid=kron(ones(10,1),(1:12)'); eventimes=frametimes'; duration=ones(120,1)*3; height=ones(120,1); events=[eventid eventimes duration height] X_bases=fmridesign(frametimes,slicetimes,events,[],zeros(1,5)); plot(squeeze(X_bases(:,1:4,1,5)),'LineWidth',2) legend('Time01','Time02','Time03','Time04') xlabel('frame number') ylabel('response') title('Slice 5, first 4 times only') contrast=[eye(12)-ones(12)/12]; exclude=[1 2]; which_stats=[0 1 1 1]; input_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['g:/keith/results/ha_100326_time01'; 'g:/keith/results/ha_100326_time02'; 'g:/keith/results/ha_100326_time03'; 'g:/keith/results/ha_100326_time04'; 'g:/keith/results/ha_100326_time05'; 'g:/keith/results/ha_100326_time06'; 'g:/keith/results/ha_100326_time07'; 'g:/keith/results/ha_100326_time08'; 'g:/keith/results/ha_100326_time09'; 'g:/keith/results/ha_100326_time10'; 'g:/keith/results/ha_100326_time11'; 'g:/keith/results/ha_100326_time12']; fmrilm(input_file,output_file_base,X_bases,contrast,exclude,which_stats,fwhm_rho); stat_threshold(1000000,26000,8,[11 103]) lm=locmax([output_file_base(1,:) '_mag_F.mnc'],5.18) voxel=[63 66 4] values=extract(voxel,output_file_base,'_mag_ef.mnc') sd=extract(voxel,output_file_base,'_mag_sd.mnc') b_hot=extract(voxel,'g:/keith/results/ha_100326_hot_mag_ef.mnc') b_wrm=extract(voxel,'g:/keith/results/ha_100326_wrm_mag_ef.mnc') time=(1:360)/10; X_hrf=fmridesign(time,0,[1 9 9 1],[],hrf_parameters); hrf=squeeze(X_hrf(:,1,1)); plot((0:12)*3+slicetimes(voxel(3)),values([1:12 1]),'k', ... [0:11; 0:11]*3+slicetimes(voxel(3)), [values+sd; values-sd],'g', ... time,[zeros(1,90) ones(1,90) zeros(1,180)],'r', ... time,[zeros(1,270) ones(1,90)],'b', ... time,hrf*b_hot+hrf([181:360 1:180])*b_wrm,'g','LineWidth',2); legend('Estimated response','Modeled response'); text(10,0.5,'Hot') text(28,0.5,'Warm') xlabel('time (seconds) from start of epoch'); ylabel('fMRI response'); % Estimating the delay contrast=[1 0; 0 1; 1 -1] which_stats=[1 1 1] input_file='g:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['g:/keith/results/ha_100326_hot'; 'g:/keith/results/ha_100326_wrm'; 'g:/keith/results/ha_100326_hmw'] n_poly=3; fwhm_rho='g:/keith/results/ha_100326_hot_rho.mnc'; confounds=[]; contrast_is_delay=[1 1 1]; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, fwhm_rho, n_poly, confounds, contrast_is_delay); t=fmris_read_image('g:/keith/results/ha_100326_hot_mag_t.mnc',4,1); d=fmris_read_image('g:/keith/results/ha_100326_hot_del_ef.mnc',4,1); s=fmris_read_image('g:/keith/results/ha_100326_hot_del_sd.mnc',4,1); a=fmris_read_image(input_file,4,1); imagesc([(d.data.*(t.data>4.86)+(a.data4.86)+(a.datamask_thresh))'); colorbar; axis xy; % Higher order autoregressive models which_stats=[1 1 1 0 1 0 0 1]; contrast=[1 -1]; output_file_base='g:/keith/results/ha_100326_hmw_ar2'; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, 15, 3, [], 0, [], 'spectral', 2); m=fmris_read_image('g:/keith/results/ha_100326_hmw_ar2_A.mnc',4,1:2); mm=reshape(permute(m.data,[2 1 3 4]),m.dim(1),prod(size(m.data))/m.dim(1)); imagesc(mm); colorbar; axis xy;