PET CBF example

We use FMRISTAT to analyze some non-kinetic PET CBF data. In this study there were 13 subjects who were given a right-hand finger tapping task at four different frequencies: 1Hz, 2Hz, 3Hz and 4Hz (labelled h1, h2, h3 and h4 in the file names). These are to be compared to a baseline task of rest (labelled hb).

Example_pet contains all the matlab commands for the example used here. The first step is to set up a matrix containing all the file names for the analysis. This could be done using an editor. For the example here, there are 13 subjects and 5 scans per subject:

input_files=[ 'C:/keith/fMRI/manou/cbf_non_kin/subj01_-h1_tal_200008161003.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj01_-h2_tal_200008161025.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj01_-h3_tal_200008161043.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj01_-h4_tal_200008161101.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj01_-hb_tal_200008160942.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj02_-h1_tal_200012051211.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj02_-h2_tal_200012051246.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj02_-h3_tal_200012051305.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj02_-h4_tal_200012051229.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj02_-hb_tal_200012051154.mnc'; . . . (10 other subjects) . . 'C:/keith/fMRI/manou/cbf_non_kin/subj13_-h1_tal_200008281408.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj13_-h2_tal_200008281350.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj13_-h3_tal_200008281333.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj13_-h4_tal_200008281319.mnc'; 'C:/keith/fMRI/manou/cbf_non_kin/subj13_-hb_tal_200008281427.mnc'] The first thing might be to look at one slice of the data (here slice 21, zero based), laid out in a task x subjects matrix. To do this, use view_slices with a last parameter layout is a matrix which lays out the slices (in the order in input_files) where you want them in the array (0 leaves a gap in the array):

n=65; layout=reshape(1:n,5,13) layout = 1 6 11 16 21 26 31 36 41 46 51 56 61 2 7 12 17 22 27 32 37 42 47 52 57 62 3 8 13 18 23 28 33 38 43 48 53 58 63 4 9 14 19 24 29 34 39 44 49 54 59 64 5 10 15 20 25 30 35 40 45 50 55 60 65 clf; view_slices(input_files,[],[],21,1:n,[0 40000],layout) This will put the first 5 slices into the first column, etc. to give

It is immediately obvious that the data need normalizing so that all images have roughly the same average intensity (here 1). To do this use normalize on each volume, putting the results into a new directory whose names are kept in input_files_norm, which can be created by hand and read in as above, or done using matlab commands as follows:

input_files_norm=[input_files(:,1:32) repmat('normalized/',n,1) input_files(:,33:size(input_files,2))] for i=1:n normalize(input_files(i,:),input_files_norm(i,:)); end clf; view_slices(input_files_norm,[],[],21,1:n,[0 1.5],layout) Now the data look much better:

Next we create the design matrix for the linear model. This is the tricky bit, which usually causes the most headaches. We need 65 rows, one for each scan, and columns for an effect for each subject, followed by columns for the effect of each task (h1,h2,h3,h4), leaving out the baseline (hb). Again you can type this in by hand:

X=[ 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . . . (other 10 subjects) . . 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0] or using some cunning advanced matlab commands: X=[kron(eye(13),ones(5,1)) kron(ones(13,1),eye(5))]; X=X(:,1:17); We next specify the contrast we are interested in. We want an overall mean (to act as a mask later on), then we want to look at each task relative to baseline. To do this we make up a matrix with the same number of columns as X, and a row for each contrast in those columns:

contrast=[mean(X,1); zeros(4,13) eye(4)] contrast = Columns 1 through 10 0.0769 0.0769 0.0769 0.0769 0.0769 0.0769 0.0769 0.0769 0.0769 0.0769 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 11 through 17 0.0769 0.0769 0.0769 0.2000 0.2000 0.2000 0.2000 0 0 0 1.0000 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 1.0000 Now we must select a matrix of names for the output bases, one row for each row of contrast. Again you can type this by hand or use some cunning matlab:

output=[repmat('c:/keith/fMRI/manou/cbf_non_kin_smoother_',5,1) num2str((0:4)')] output = c:/keith/fMRI/manou/cbf_non_kin_smoother_0 c:/keith/fMRI/manou/cbf_non_kin_smoother_1 c:/keith/fMRI/manou/cbf_non_kin_smoother_2 c:/keith/fMRI/manou/cbf_non_kin_smoother_3 c:/keith/fMRI/manou/cbf_non_kin_smoother_4 Then we specify which statistics we want as output. Let's ask for effects (_ef), sds (_sd), T statistics (_t) and fwhm images (_fwhm):

which_stats='_t _ef _sd _fwhm' Now we are ready for the final run. Sd input files, and their df's are left empty, so multistat will smooth the sample sd to achieve 100 (default) df for the t statistics:

multistat(input_files_norm,[],[],[],X,contrast,output,which_stats); Multistat makes a first pass through the data to work out the sd, effective smoothness of the data (fwhm_data = 8.3371mm) and smoothing of the sd needed to achieve 100 df (fwhm_varatio = 4.6505mm). It then passes through a second time to calculate the which_stats.

To look at the results, we can choose the overall average effect as a mask, which we can threshold at say 0.65 (half its maximum). Then we can again use view_slices to take a look at every 5th slice through the volume of say task 3. For the effect of task 3 - baseline (in units of normalized cbf):

mask='c:/keith/fMRI/manou/cbf_non_kin_smoother_0_ef.mnc' clf; view_slices('c:/keith/fMRI/manou/cbf_non_kin_smoother_3_ef.mnc',mask,0.65,[0:5:79])

Standard error, in units of normalized cbf:

clf; view_slices('c:/keith/fMRI/manou/cbf_non_kin_smoother_3_sd.mnc',mask,0.65,[0:5:79])

T statistic = effect / standard error:

clf; view_slices('c:/keith/fMRI/manou/cbf_non_kin_smoother_3_t.mnc',mask,0.65,[0:5:79])

Effective FWHM, in mm:

clf; view_slices('c:/keith/fMRI/manou/cbf_non_kin_smoother_0_fwhm.mnc',mask,0.65,[0:5:79])

Or we might want to compare the t statistics on all tasks at a range of slices (64 65,66) that cover the region of greatest activation (rows are tasks, columns are slices):

m=['c:/keith/fMRI/manou/cbf_non_kin_smoother_1_t.mnc'; 'c:/keith/fMRI/manou/cbf_non_kin_smoother_2_t.mnc'; 'c:/keith/fMRI/manou/cbf_non_kin_smoother_3_t.mnc'; 'c:/keith/fMRI/manou/cbf_non_kin_smoother_4_t.mnc'] clf; view_slices(m,mask,0.65,64:66,1:4)

Finally we can use stat_summary to get locations of clusters, local maxima, and their P-values:

stat_summary('c:/keith/fMRI/manou/cbf_non_kin_smoother_3_t.mnc', ... 'c:/keith/fMRI/manou/cbf_non_kin_smoother_0_fwhm.mnc',[],mask,0.65); clus vol resel Pval (one) 1 3056 4.25 0 ( 0) 2 1217 2.9 0.002 ( 0) 3 256 0.61 0.473 (0.025) 4 228 0.53 0.583 (0.034) 5 73 0.19 0.989 (0.176) 6 62 0.13 0.998 (0.249) 7 52 0.08 1 (0.363) 8 31 0.05 1 (0.484) 9 17 0.04 1 (0.508) 10 17 0.03 1 (0.617) 11 14 0.02 1 (0.638) clus peak Pval (one) Qval (i j k) ( x y z ) 1 6.42 0 ( 0) 0 (37 62 64) (-36.2 -19.4 58.5) 1 6.26 0 ( 0) 0 (38 61 66) (-34.8 -21.2 61.5) 2 6.12 0.001 ( 0) 0 (77 43 12) ( 17.4 -52.1 -19.5) 3 4.61 0.483 (0.019) 0.008 (61 70 62) ( -4 -5.7 55.5) 4 4.01 3.376 (0.122) 0.044 (70 37 19) ( 8 -62.4 -9) 4 3.98 3.605 ( 0.13) 0.047 (69 37 18) ( 6.7 -62.4 -10.5) 4 3.97 3.787 (0.137) 0.049 (70 38 20) ( 8 -60.7 -7.5) 6 3.72 7.498 (0.269) 0.094 (86 42 6) ( 29.5 -53.8 -28.5) 5 3.67 8.653 ( 0.31) 0.107 (29 58 59) (-46.9 -26.3 51) 11 3.43 15.522 (0.553) 0.185 (38 51 62) (-34.8 -38.4 55.5) 7 3.43 15.619 (0.557) 0.186 (57 74 64) ( -9.4 1.2 58.5) 8 3.37 18.151 (0.646) 0.215 (67 74 62) ( 4 1.2 55.5) 9 3.33 20.117 (0.715) 0.239 (64 35 14) ( 0 -65.9 -16.5) 9 3.32 20.295 (0.721) 0.241 (63 34 14) ( -1.3 -67.6 -16.5) 10 3.31 20.734 (0.736) 0.246 (30 58 66) (-45.6 -26.3 61.5) 7 3.26 23.278 (0.825) 0.274 (55 72 65) (-12.1 -2.2 60) As usual it will also produce a glass brain and a 3D picture of the significant clusters that you can rotate:

The above analysis assumed that scans were independent and had equal sd on a given subject. If we are not sure about this, we can do a "safer" analysis of say task 3 - baseline by keeping only the task3 and baseline scans, and eliminating the rest. This reduces the residual degrees of freedom from (#subjects -1)x(#tasks - 1)=48 to (#subjects - 1)=12, so more smoothing of the sd is required. The matlab commands are identical:

keep=[3 5 8 10 13 15 18 20 23 25 28 30 33 35 38 40 43 45 48 50 53 55 58 60 63 65] input_files_norm3b=input_files_norm(keep,:) X1=[kron(eye(13),ones(2,1)) kron(ones(13,1),eye(2))]; X1=X1(:,1:14) contrast1=[zeros(1,13) 1] output3b='c:/keith/fMRI/manou/cbf_non_kin_smooth_3b' multistat(input_files_norm3b,[],[],[],X1,contrast1,output3b,which_stats); This time the sd is smoothed by fwhm_varatio = 11.2977mm to achieve 100 df. We can compare the two analyses as follows:

m=['c:/keith/fMRI/manou/cbf_non_kin_smoother_3_t.mnc'; 'c:/keith/fMRI/manou/cbf_non_kin_smooth_3b_t.mnc '] clf; view_slices(m,mask,0.65,64:66,1:2,[-3 6],[1 2 3 4; 5 6 7 8])

The top row is the first analysis (assuming independent scans with equal sd), the second row is the "safer" analysis which doesn't make this assumption. The difference is imperceptible in this case.

Finally here is a list of clusters, peaks and their P-values:

stat_summary('c:/keith/fMRI/manou/cbf_non_kin_smooth_3b_t.mnc', ... 'c:/keith/fMRI/manou/cbf_non_kin_smooth_3b_fwhm.mnc',[],mask,0.65); clus vol resel Pval (one) 1 3315 4.51 0 ( 0) 2 844 1.4 0.064 (0.003) 3 311 0.64 0.408 (0.025) 4 138 0.22 0.949 (0.139) 5 86 0.22 0.951 (0.141) 6 10 0.04 1 (0.522) 7 3 0.01 1 (0.804) 8 3 0.01 1 (0.809) 9 3 0.01 1 (0.834) clus peak Pval (one) Qval (i j k) ( x y z ) 1 6.38 0 ( 0) 0 (38 61 64) (-34.8 -21.2 58.5) 1 6.28 0 ( 0) 0 (37 62 62) (-36.2 -19.4 55.5) 2 5.33 0.022 (0.001) 0.001 (77 43 12) ( 17.4 -52.1 -19.5) 3 4.26 1.299 (0.057) 0.024 (61 70 61) ( -4 -5.7 54) 5 3.69 6.718 ( 0.29) 0.11 (69 37 18) ( 6.7 -62.4 -10.5) 4 3.68 6.981 (0.301) 0.114 (55 72 65) (-12.1 -2.2 60) 4 3.58 8.959 (0.386) 0.144 (55 71 66) (-12.1 -4 61.5) 3 3.47 11.708 (0.502) 0.187 (67 74 62) ( 4 1.2 55.5) 3 3.35 15.917 ( 0.68) 0.249 (65 73 63) ( 1.3 -0.5 57) 4 3.34 15.975 (0.683) 0.249 (57 74 64) ( -9.4 1.2 58.5) 6 3.33 16.64 (0.711) 0.259 (60 53 9) ( -5.4 -34.9 -24) 8 3.26 19.324 (0.824) 0.296 (64 35 14) ( 0 -65.9 -16.5) 9 3.21 21.872 (0.931) 0.333 (85 42 5) ( 28.1 -53.8 -30) 7 3.18 23.242 (0.988) 0.35 (63 34 14) ( -1.3 -67.6 -16.5)