An example of an fMRI analysis in FSL
John Muschelli
2021-02-16
All code for this document is located at here.
library(methods)
library(fslr)
library(neurobase)
library(extrantsr)
Attaching package: 'extrantsr'The following objects are masked from 'package:ANTsRCore':
origin, origin<-The following object is masked from 'package:neurobase':
zero_padThe following objects are masked from 'package:oro.nifti':
origin, origin<-In this tutorial we will discuss performing some preprocessing of a single subject functional MRI in FSL.
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) {
source("https://neuroconductor.org/neurocLite.R")
neuroc_install("kirby21.base")
}
if (!"kirby21.fmri" %in% packages) {
source("https://neuroconductor.org/neurocLite.R")
neuroc_install("kirby21.fmri")
}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).
library(kirby21.fmri)
library(kirby21.base)
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]
base_fname = nii.stub(fmri_fname, bn = TRUE)Parameter file
If you’d like to see the header information from the fMRI data, it is located by the following commands:
library(R.utils)
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)
info[11:23] [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).
Outline
The steps I will perform in this analysis:
- Calculation of Motion Parameters (
fslr::mcflirt) - Slice timing correction (
fslr::fsl_slicetimer), but we need to know how the scan was taken/slice order and repetition time (TR) - Motion Correction on Corrected Data (
fslr::mcflirt) - Coregistration of fMRI and a T1-weighted image (
fslr::flirt) - Registration to the Template space (
fslr::fnirt_with_affine) - De-meaning the data (fslr::fslmean(ts = TRUE))
- Skull stripping (fslr::fslbet)
- Registration to a template using the T1 and then transforming the fMRI with it
- Spatially smoothing the data (fslr:fslsmooth)
- Tissue-class segmentation (fslr::fast, ANTsR::atropos or extrantsr::otropos)?
- Bandpass/butterworth filtering (signal::butter, signal::buttord)
- 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)rm(list = "fmri") # just used for cleanup Stabilization of Signal
Volumes corresponding to the first 10 seconds of the rs-fMRI scan were dropped to allow for magnetization stabilization.
fmri = readnii(fmri_fname)
tr = 2 # 2 seconds
first_scan = floor(10.0 / tr) + 1 # 10 seconds "stabilization of signal"
sub_fmri = subset_4d(fmri, first_scan:ntim(fmri))Slice Timing Correction
library(fslr)
scorr_fname = paste0(base_fname, "_scorr.nii.gz")
if (!file.exists(scorr_fname)) {
scorr = fsl_slicetimer(file = sub_fmri,
outfile = scorr_fname,
tr = 2, direction = "z", acq_order = "contiguous",
indexing = "up")
} else {
scorr = readnii(scorr_fname)
}slicetimer -i "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpj90lVj/file10ebb50e05aa2.nii.gz" -v -r 2 -d 3 -o "113-01-fMRI_scorr";Slice Timing Correction
Similarly to the ANTsR processing, we can set opts = "-meanvol" so that the motion correction registers the images to the mean of the time series. The default otherwise is to register to the scan in the middle (closest to number of time points / 2). You can supercede this with either specifying -refvol for which time point to register to, or -reffile to specify a volume to register to. You can think of -meanvol as a wrapper for making the mean time series (using fslmaths(opt = "-Tmean")) and then passing that in as a reffile. We will do this explicitly in the code below. The -plots arguments outputs a .par file so that we can read the motion parameters after mcflirt is run.
moco_fname = paste0(base_fname,
"_motion_corr.nii.gz")
par_file = paste0(nii.stub(moco_fname), ".par")
avg_fname = paste0(base_fname,
"_avg.nii.gz")
if (!file.exists(avg_fname)) {
fsl_maths(file = sub_fmri,
outfile = avg_fname,
opts = "-Tmean")
}fslmaths "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpj90lVj/file10ebb65cd6e2b.nii.gz" -Tmean "113-01-fMRI_avg";[1] "113-01-fMRI_avg.nii.gz"if (!all(file.exists(c(moco_fname, par_file)))) {
moco_img = mcflirt(
file = sub_fmri,
outfile = moco_fname,
verbose = 2,
opts = paste0("-reffile ", avg_fname, " -plots")
)
} else {
moco_img = readnii(moco_fname)
}mcflirt -in "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpj90lVj/file10ebb11fbfd44.nii.gz" -reffile 113-01-fMRI_avg.nii.gz -plots -verbose 2 -o "113-01-fMRI_motion_corr";moco_params = readLines(par_file)
moco_params = strsplit(moco_params, split = " ")
moco_params = sapply(moco_params, function(x) {
as.numeric(x[ !(x %in% "")])
})
moco_params = t(moco_params)
colnames(moco_params) = paste0("MOCOparam", 1:ncol(moco_params))
head(moco_params) MOCOparam1 MOCOparam2 MOCOparam3 MOCOparam4 MOCOparam5 MOCOparam6
[1,] -0.000309931 0.001248200 -4.26752e-04 0.03832110 -0.558039 0.0768606
[2,] -0.000650414 0.000826515 -5.37836e-04 0.03833610 -0.607328 0.0589568
[3,] -0.000581280 0.001141430 -4.26752e-04 0.03831010 -0.566965 0.0704959
[4,] -0.000906044 0.000826516 -4.26752e-04 0.02266840 -0.577867 0.0465875
[5,] -0.000188335 0.001106950 -1.23971e-04 0.01542600 -0.556718 0.0412824
[6,] -0.000605871 0.000849419 -6.14711e-06 0.00108561 -0.577868 0.0439374Let’s Make a Matrix!
img_ts_to_matrix creates \(V\times T\) matrix, \(V\) voxels in mask, unlike ANTsR::timeseries2matrix. Therefore, we transpose the matrix so that it is consistent with the ANTsR tutorial for fMRI and so boldMatrix is \(T\times V\). We will get the average of the co-registered image using fslmaths. We wil use this average image to get a mask using the oMask function, which calls getMask, but for nifti objects. We will then zero out the average image using the mask image.
moco_avg_img = fslmaths(moco_fname, opts = "-Tmean")fslmaths "/Users/johnmuschelli/Dropbox/Projects/neuroc/fmri_analysis_fslr/113-01-fMRI_motion_corr.nii.gz" -Tmean "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpj90lVj/file10ebb7ed06e3f";maskImage = oMask(moco_avg_img,
mean(moco_avg_img),
Inf, cleanup = 2)
mask_fname = paste0(base_fname, "_mask.nii.gz")
writenii(maskImage, filename = mask_fname)
bet_mask = fslbet(moco_avg_img) > 0bet2 "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpj90lVj/file10ebb611db14.nii.gz" "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpj90lVj/file10ebb579b1353" bet_mask_fname = paste0(base_fname, "_bet_mask.nii.gz")
writenii(bet_mask, filename = bet_mask_fname)We can also look at the differences. We will assume the maskImage as the “gold standard” and bet_mask to be the “prediction”.
double_ortho(moco_avg_img, maskImage,
col.y = "white")double_ortho(moco_avg_img, bet_mask,
col.y = "white")ortho_diff( moco_avg_img, roi = maskImage, pred = bet_mask )Here we will create a matrix of time by voxels.
moco_avg_img[maskImage == 0] = 0
boldMatrix = img_ts_to_matrix(
moco_img)
boldMatrix = t(boldMatrix)
boldMatrix = boldMatrix[ , maskImage == 1]Calculation of DVARS
With this boldMatrix, we can calculate a series of information. For example, we can calculate DVARS based on the motion corrected data. We can also compare the DVARS to the DVARS calculated from the non-realigned data.
Here we show how you can still use ANTsR::computeDVARS to calculate DVARS. The first element of dvars is the mean of the my_dvars. By the definition of @power2012spurious, the first element of DVARS should be zero.
Here we will multiply the 3 first motion parameters (roll, pitch, yaw) by 50 to convert radians to millimeters by assuming a brain radius of 50 mm, as similar to @power2012spurious. The next 3 parameters are in terms of millimeters (x, y, z). We will plot each of the parameters on the same scale to look at the motion for each scan.
dvars = ANTsR::computeDVARS(boldMatrix)
dMatrix = apply(boldMatrix, 2, diff)
dMatrix = rbind(rep(0, ncol(dMatrix)), dMatrix)
my_dvars = sqrt(rowMeans(dMatrix^2))
head(cbind(dvars = dvars, my_dvars = my_dvars)) dvars my_dvars
[1,] 18861.75 0.00
[2,] 18335.35 18335.35
[3,] 18855.69 18855.69
[4,] 17292.62 17292.62
[5,] 18590.46 18590.46
[6,] 18375.67 18375.67print(mean(my_dvars))[1] 18861.75Similarly, we can calculate the marginal framewise displacement (FD). The rotation parameters are again in radians so we can translate these to millimeters based on a 50 mm radius of the head.
mp = moco_params
mp[, 1:3] = mp[, 1:3] * 50
mp = apply(mp, 2, diff)
mp = rbind(rep(0, 6), mp)
mp = abs(mp)
fd = rowSums(mp)mp = moco_params
mp[, 1:3] = mp[, 1:3] * 50
r = range(mp)
plot(mp[,1], type = "l", xlab = "Scan Number", main = "Motion Parameters",
ylab = "Displacement (mm)",
ylim = r * 1.25,
lwd = 2,
cex.main = 2,
cex.lab = 1.5,
cex.axis = 1.25)
for (i in 2:ncol(mp)) {
lines(mp[, i], col = i)
}rm(list = "mp")Heatmap of the values
We can look at the full trajectory of each voxel over each scan. We scaled the data (by column, which is voxel), which is somewhat equivalent to doing whole-brain z-score normalization of the fMRI.
We can find the index which has the highest mean value, which may indicate some motion artifact.
library(RColorBrewer)
library(matrixStats)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
r <- rf(32)
mat = scale(boldMatrix)
image(x = 1:nrow(mat),
y = 1:ncol(mat),
mat, useRaster = TRUE,
col = r,
xlab = "Scan Number", ylab = "Voxel",
main = paste0("Dimensions: ",
dim(mat)[1], "×", dim(mat)[2]),
cex.main = 2,
cex.lab = 1.5,
cex.axis = 1.25)
rmeans = rowMeans(mat)
bad_ind = which.max(rmeans)
print(bad_ind)[1] 174abline(v = bad_ind)sds = rowSds(mat)
print(which.max(sds))[1] 174rm(list = "mat")library(animation)
ani.options(autobrowse = FALSE)
gif_name = "bad_dimension.gif"
if (!file.exists(gif_name)) {
arr = as.array(moco_img)
pdim = pixdim(moco_img)
saveGIF({
for (i in seq(bad_ind - 1, bad_ind + 1)) {
ortho2(arr[,,,i], pdim = pdim, text = i)
}
}, movie.name = gif_name)
}
Much of the rest of the commands are again within R or already implemented in ANTSR. Analyses in FSL commonly use independent components analysis (ICA) with the melodic command. We will not cover that here and will cover it in a subsequent tutorial.
Session Info
devtools::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.0.2 (2020-06-22)
os macOS Catalina 10.15.7
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2021-02-16
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib source
abind 1.4-5 2016-07-21 [2] CRAN (R 4.0.0)
animation * 2.6 2018-12-11 [2] CRAN (R 4.0.0)
ANTsR * 0.5.6.1 2020-06-01 [2] Github (ANTsX/ANTsR@9c7c9b7)
ANTsRCore * 0.7.4.6 2020-07-07 [2] Github (muschellij2/ANTsRCore@61c37a1)
assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.0)
base64enc 0.1-3 2015-07-28 [2] CRAN (R 4.0.0)
bitops 1.0-6 2013-08-17 [2] CRAN (R 4.0.0)
cachem 1.0.4 2021-02-13 [1] CRAN (R 4.0.2)
callr 3.5.1 2020-10-13 [1] CRAN (R 4.0.2)
cli 2.3.0 2021-01-31 [1] CRAN (R 4.0.2)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.0.2)
colorout * 1.2-2 2020-06-01 [2] Github (jalvesaq/colorout@726d681)
colorspace 2.0-0 2020-11-11 [1] CRAN (R 4.0.2)
crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.2)
DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.2)
desc 1.2.0 2020-06-01 [2] Github (muschellij2/desc@b0c374f)
devtools 2.3.2 2020-09-18 [1] CRAN (R 4.0.2)
digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2)
dplyr * 1.0.4 2021-02-02 [1] CRAN (R 4.0.2)
ellipsis 0.3.1 2020-05-15 [2] CRAN (R 4.0.0)
evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.0)
extrantsr * 3.9.13.1 2020-09-03 [2] Github (muschellij2/extrantsr@00c75ad)
fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.0.2)
fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.2)
fslr * 2.25.0 2021-02-16 [1] local
generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.2)
ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 4.0.2)
git2r 0.28.0 2021-01-11 [1] Github (ropensci/git2r@4e342ca)
glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
gtable 0.3.0 2019-03-25 [2] CRAN (R 4.0.0)
highr 0.8 2019-03-20 [2] CRAN (R 4.0.0)
htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.2)
httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.2)
ITKR 0.5.3.3.0 2021-02-15 [1] Github (stnava/ITKR@ea0ac19)
kirby21.base * 1.7.4 2020-10-01 [1] local
kirby21.fmri * 1.7.0 2018-08-13 [2] CRAN (R 4.0.0)
kirby21.t1 * 1.7.3.2 2021-01-09 [1] local
knitr 1.31 2021-01-27 [1] CRAN (R 4.0.2)
lattice 0.20-41 2020-04-02 [2] CRAN (R 4.0.2)
lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2)
magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
matlabr 1.6.0 2020-07-01 [2] local
Matrix 1.3-2 2021-01-06 [1] CRAN (R 4.0.2)
matrixStats * 0.58.0 2021-01-29 [1] CRAN (R 4.0.2)
memoise 2.0.0 2021-01-26 [1] CRAN (R 4.0.2)
mgcv 1.8-33 2020-08-27 [1] CRAN (R 4.0.2)
munsell 0.5.0 2018-06-12 [2] CRAN (R 4.0.0)
neurobase * 1.31.0 2020-10-07 [1] local
neurohcp * 0.9.0 2020-10-19 [1] local
nlme 3.1-152 2021-02-04 [1] CRAN (R 4.0.2)
oro.nifti * 0.11.0 2020-09-04 [2] local
pillar 1.4.7 2020-11-20 [1] CRAN (R 4.0.2)
pkgbuild 1.2.0 2020-12-15 [1] CRAN (R 4.0.2)
pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.0)
pkgload 1.1.0 2020-05-29 [2] CRAN (R 4.0.0)
plyr 1.8.6 2020-03-03 [2] CRAN (R 4.0.0)
prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.0.0)
processx 3.4.5 2020-11-30 [1] CRAN (R 4.0.2)
ps 1.5.0 2020-12-05 [1] CRAN (R 4.0.2)
purrr 0.3.4 2020-04-17 [2] CRAN (R 4.0.0)
R.matlab 3.6.2 2018-09-27 [2] CRAN (R 4.0.0)
R.methodsS3 * 1.8.1 2020-08-26 [1] CRAN (R 4.0.2)
R.oo * 1.24.0 2020-08-26 [1] CRAN (R 4.0.2)
R.utils * 2.10.1 2020-08-26 [1] CRAN (R 4.0.2)
R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2)
RColorBrewer * 1.1-2 2014-12-07 [2] CRAN (R 4.0.0)
Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.2)
RcppEigen 0.3.3.9.1 2020-12-17 [1] CRAN (R 4.0.2)
remotes 2.2.0 2020-07-21 [2] CRAN (R 4.0.2)
reshape2 * 1.4.4 2020-04-09 [2] CRAN (R 4.0.0)
rlang 0.4.10 2020-12-30 [1] CRAN (R 4.0.2)
rmarkdown * 2.6 2020-12-14 [1] CRAN (R 4.0.2)
RNifti 1.3.0 2020-12-04 [1] CRAN (R 4.0.2)
rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.0.2)
rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.2)
scales 1.1.1 2020-05-11 [2] CRAN (R 4.0.0)
sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.0)
spm12r * 2.8.2 2021-01-11 [1] local
stapler 0.7.2 2020-07-09 [2] Github (muschellij2/stapler@79e23d2)
stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
stringr * 1.4.0 2019-02-10 [2] CRAN (R 4.0.0)
testthat 3.0.2 2021-02-14 [1] CRAN (R 4.0.2)
tibble 3.0.6 2021-01-29 [1] CRAN (R 4.0.2)
tidyselect 1.1.0 2020-05-11 [2] CRAN (R 4.0.0)
usethis 2.0.1 2021-02-10 [1] CRAN (R 4.0.2)
vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.2)
WhiteStripe 2.3.2 2019-10-01 [2] CRAN (R 4.0.0)
withr 2.4.1 2021-01-26 [1] CRAN (R 4.0.2)
xfun 0.21 2021-02-10 [1] CRAN (R 4.0.2)
yaml * 2.2.1 2020-02-01 [2] CRAN (R 4.0.0)
zoo * 1.8-8 2020-05-02 [2] CRAN (R 4.0.0)
[1] /Users/johnmuschelli/Library/R/4.0/library
[2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library