An example of an fMRI analysis in ANTsR

All code for this document is located at here.

In this tutorial we will discuss performing some preprocessing of a single subject functional MRI in ANTsR.

1 Data Packages

For this analysis, I will use one subject from the Kirby 21 data set. The kirby21.base and kirby21.fmri packages are necessary for this analysis and have the data we will be working on. You need devtools to install these. Please refer to installing devtools for additional instructions or troubleshooting.

packages = installed.packages()
packages = packages[, "Package"]
if (!"kirby21.base" %in% packages) {
if (!"kirby21.fmri" %in% packages) {

2 Loading Data

We will use the get_image_filenames_df function to extract the filenames on our hard disk for the T1 image and the fMRI images (4D).

fnames = get_image_filenames_df(ids = 113, 
                    modalities = c("T1", "fMRI"), 
                    visits = c(1),
                    long = FALSE)
t1_fname = fnames$T1[1]
fmri_fname = fnames$fMRI[1]

2.1 Parameter file

If you’d like to see the header information from the fMRI data, it is located by the following commands:

par_file = system.file("visit_1/113/113-01-fMRI.par.gz", 
                       package = "kirby21.fmri")
# unzip it
con = gunzip(par_file, temporary = TRUE, 
             remove = FALSE, overwrite = TRUE)
info = readLines(con = con)
 [1] ".    Protocol name                      :   WIP Bold_Rest SENSE"
 [2] ".    Series Type                        :   Image   MRSERIES"   
 [3] ".    Acquisition nr                     :   11"                 
 [4] ".    Reconstruction nr                  :   1"                  
 [5] ".    Scan Duration [sec]                :   434"                
 [6] ".    Max. number of cardiac phases      :   1"                  
 [7] ".    Max. number of echoes              :   1"                  
 [8] ".    Max. number of slices/locations    :   37"                 
 [9] ".    Max. number of dynamics            :   210"                
[10] ".    Max. number of mixes               :   1"                  
[11] ".    Patient position                   :   Head First Supine"  
[12] ".    Preparation direction              :   Anterior-Posterior" 
[13] ".    Technique                          :   FEEPI"              

From the paper “Multi-parametric neuroimaging reproducibility: A 3-T resource study”, which this data is based on, it describes the fMRI sequence:

The sequence used for resting state functional connectivity MRI is typically identical to that used for BOLD functional MRI studies of task activation. Here, we used a 2D EPI sequence with SENSE partial-parallel imaging acceleration to obtain 3 × 3 mm (80 by 80 voxels) in-plane resolution in thirty-seven 3 mm transverse slices with 1 mm slice gap. An ascending slice order with TR/TE = 2000/30 ms, flip angle of 75°, and SENSE acceleration factor of 2 were used. SPIR was used for fat suppression. This study used an ascending slice acquisition order because a pilot studies revealed smaller motion induced artifacts with ascending slice order than with interleaved slice order. While using an ascending slice order, it was necessary to use a small slice gap to prevent cross talk between the slices. One 7-min run was recorded which provided 210 time points (discarding the first four volumes to achieve steady state).

3 Outline

The steps I will perform in this analysis:

  1. Calculation of Motion Parameters (fslr::mcflirt)
  2. Slice timing correction (fslr::fsl_slicetimer), but we need to know how the scan was taken/slice order and repetition time (TR)
  3. Motion Correction on Corrected Data (fslr::mcflirt)
  4. Coregistration of fMRI and a T1-weighted image (fslr::flirt)
  5. Registration to the Template space (fslr::fnirt_with_affine )
  6. Skull stripping (fslr::fslbet)
  7. Registration to a template using the T1 and then transforming the fMRI with it
  8. Spatially smoothing the data (fslr:fslsmooth)
  9. Tissue-class segmentation (fslr::fast, ANTsR::atropos or extrantsr::otropos)?
  10. Bandpass/butterworth filtering (signal::butter, signal::buttord)
  11. Get a connectivity matrix of certain regions, you need to specify an atlas.

Now we know that the head is first in (as usual) and the data was acquired in ascending order (i.e. bottom -> up) and the repetition time (TR) was 2 seconds The

fmri = readnii(fmri_fname)
ortho2(fmri, w = 1, add.orient = FALSE)