Simple Correlation Analysis


The following steps are only suggestive. A better alternative is to follow Example 9 in "afni_proc.py -help".

This type of correlation analysis typically looks for the correlation between the time series at the seed region and the rest of the brain. The seed region time series may be the signal under one specific task/condition or resting state.

The following example of running correlation analysis demonstrates how such analysis is done with AFNI. The steps are built based on discussion on AFNI message board and through personal communications (Thanks especially go to Giuseppe Pagnoni for his contribution and insight).

(1) Remove the baseline, drifting effect, head motion and any effects of no interest from the original time series. The purpose of this step is to make the seed region's time series immune to any undesirable effects other than the resting state (or task effect). This can be done with program 3dSynthesize (and 3dcalc), but it has to take some output from 3dDeconvolve as input. More specifically it needs two inputs, one from option -cbucket in 3dDeconvolve while the other from option -x1D. If you haven't done so in the previous analysis, you may have to re-run 3dDeconvolve with these options. The following command isolates the effects of no interest:

3dSynthesize -cbucket CbucketFileFrom3dDeconvolve+orig -matrix x1DFrom3dDeconvolve -select ...  prefix EffectsOfNoInterest

Then you may want to  remove the effects of no interest from the original signal by using 3dcalc:

3dcalc -a InputFor3Deconvolve -b EffectsOfNoInterest -expr 'a-b' -prefix CleanData

(2) Convert the subject's EPI data in time series to common space

adwarp -apar anat_subj01+tlrc -dpar 'CleanData_subj01+orig' -dxyz 3 -prefix EPI_subj01

where anat_subj01+tlrc is an anatomical data of the subject.

(3) Use 3dmaskave (or 3dmaskdump/3dROIstat) to extract average time series of the ROI or time series at a voxel with a peak t value (one voxel - (46, -29, 21) - in this case)

3dmaskdump -noijk -dbox 46 -29 21 EPI_subj01+tlrc > EPI_subj01.1D

Notes: 

A. -noijk avoids writing out 3 coordinate index numbers at the beginning to the output file.
B. Regarding -dbox, be careful with different coordinate systems. Here is some detail about the issue. Also check out '3dmaskdump -help'.
C. Use option -mask in 3dmaskave/3dmaskdump/3dROIstat if you want to use an ROI as the seed.

The output EPI_subj01.1D is a one-row text file. Convert the one-row time series to one column:

1dtranspose EPI_subj01.1D Seed_subj01.1D

(4) Correlation analysis (with quadratic fitting for baseline and trend), 

3dfim+ -input EPI_subj01+tlrc \
 -polort 2 \                             # Quadratic fitting of the baseline  + drifting 
 -ort_file covariates.1D \         # optional, e.g., physiological regressors
 -ideal_file Seed_subj01.1D \  # Seed 
 -out Correlation \                  
 -bucket Corr_subj01

3dDeconvolve is favorable especially if some time points need to be censored out or if regressors of no interest (e.g., head motion) are involved in the analysis,

3dDeconvolve -input EPI_subj01+tlrc \
 -polort 2 \
 -censor ... \
 -num_stimts 2 \
 -stim_file 1 Seed_subj01.1D -stim_label 1 SeedCorr \
 -stim_file 2 HeadMotion.1D'[0] -stim_label 2 Roll \
 -stim_file 3 HeadMotion.1D'[1] -stim_label 3 Pitch \

 ...

 -tout -rout -fitts fit_subj01 \
 -bucket Corr_subj01

However 3dDeconvolve can only spill out R2, not correlation coefficient R itself. So if you want R, you would have to take square root and find out its sign based on the sign of its corresponding beta value. 

3dcalc -a Corr_subj01+tlrc'[SubbrickForR2]' -b Corr_subj01+tlrc'[SubbrickForBeta]'-expr 'ispositive(b)*sqrt(a)-isnegative(b)*sqrt(a)'-prefix Corr_subj01R

(5) Group analysis

a. Convert correlation coefficients to Gaussian

Correlation coefficients range from -1 to 1. To be able to run group analysis, Fisher's Z transformation formula can be used to reduce skewness and make the sampling distribution more normal when sample size is big enough: z =  (1/2) * ln((1+r)/(1-r)), where z is approximately normally distributed with mean r, and standard error 1/(n-3)0.5 (n: sample size).

3dcalc -a Corr_subj01R+tlrc -expr 'log((1+a)/(1-a))/2' -prefix Corr_subj01_Z

or more concisely,

3dcalc -a Corr_subj01R+tlrc -expr 'atanh(a)' -prefix Corr_subj01_Z

 

b. Group analysis with  Z scores or regressor coefficients of all subjects

 

3dttest -prefix Grp_Result -base1 0 \

-set2 \
Corr_subj01_Z+tlrc \
Corr_subj02_Z+tlrc \
Corr_subj03_Z+tlrc \
Corr_subj04_Z+tlrc \
Corr_subj05_Z+tlrc \
Corr_subj06_Z+tlrc \
Corr_subj07_Z+tlrc \
Corr_subj08_Z+tlrc \
Corr_subj09_Z+tlrc \
Corr_subj10_Z+tlrc \

Corr_subj11_Z+tlrc \

Corr_subj12_Z+tlrc 

 

Keep in mind a high t value only indicates r is with some significance level different from 0, does not necessarily mean a high correlation. And if you want the correlation coefficient at group level, you can convert the z-score (sub-brick #0 in the output from the group analysis) to r (the reverse conversion of step (5) above):

3dcalc -a Grp_Result+tlrc'[0]' -expr '(exp(2*a)-1)/(exp(2*a)+1)' -prefix GrpR

or,

3dcalc -a Grp_Result+tlrc'[0]' -expr 'atan(a)' -prefix GrpR