DTI Analysis in fslr

All code for this document is located at here.

Resources and Goals

Much of this work has been adapted by the FSL guide for DTI: http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/UserGuide. We will show you a few steps that have been implemented in fslr: eddy_correct and dtifit. Although xfibres has been adapted for fslr, which is the backend for bedpostx from FSL, it takes a long time to run and the results are not seen in this vignette, though code is given to illustrate how it would be run.

Data Packages

For this analysis, I will use one subject from the Kirby 21 data set. The kirby21.base and kirby21.dti 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.dti" %in% packages) {

Loading Data

We will use the get_image_filenames_df function to extract the filenames on our hard disk for the T1 image.

[1] TRUE
fnames = get_image_filenames_df(ids = 113, 
                    modalities = c("T1", "DTI"), 
                    visits = c(1),
                    long = FALSE)
t1_fname = fnames$T1[1]
dti_fname = fnames$DTI[1]
base_fname = nii.stub(dti_fname, bn = TRUE)
dti_data = get_dti_info_filenames(
  ids = 113, 
  visits = c(1))
b_fname = dti_data$fname[ dti_data$type %in% "b"]
grad_fname = dti_data$fname[ dti_data$type %in% "grad"]

Making b-vectors and b-values

As dtifit requires the b-values and b-vectors to be separated, we will take a look at this data.

b_vals = readLines(b_fname)
b_vals = as.numeric(b_vals)
b0 = which(b_vals == 0)[1]

b_vecs = read.delim(grad_fname, header = FALSE)
b_vecs$V4 = NULL
colnames(b_vecs) = c("x", "y", "z")

Printing out FSL Version

[1] "5.0.10"

Checking our data

Here we ensure that the number of b-values/b-vectors is the same as the number of time points in the 4D image.

n_timepoints = fslval(dti_fname, "dim4")
fslval "/Users/johnmuschelli/Library/R/4.0/library/kirby21.dti/visit_1/113/113-01-DTI.nii.gz" dim4
stopifnot(nrow(b_vecs) == n_timepoints)
stopifnot(length(b_vals) == n_timepoints)

Running eddy_correct

Here, we will run an eddy current correction using FSL’s eddy_correct through fslr. We will save the result in a temporary file (outfile), but also return the result as a nifti object ret, as retimg = TRUE. We will use the first volume as the reference as is the default in FSL. Note FSL is zero-indexed so the first volume is the zero-ith index:

eddy_fname = paste0(base_fname, "_eddy.nii.gz")
if (!file.exists(eddy_fname)) {
  eddy = eddy_correct(
    infile = dti_fname, 
    outfile = eddy_fname, 
    retimg = TRUE, 
    reference_no = 0)
} else {
  eddy = readnii(eddy_fname)
eddy_correct "/Users/johnmuschelli/Library/R/4.0/library/kirby21.dti/visit_1/113/113-01-DTI.nii.gz"  "113-01-DTI_eddy" 0;

Let’s look at the eddy current-corrected (left) and the non-corrected data (right). Here we will look at the image where the b-value is equal to zero.

dti = readnii(dti_fname)
eddy0 = extrantsr::subset_4d(eddy, b0)
dti0 = extrantsr::subset_4d(dti, b0)
double_ortho(robust_window(eddy0), robust_window(dti0))

Note, from here on forward we will use either the filename for the output of the eddy current correction or the eddy-current-corrected nifti object.

Getting a brain mask

Let’s get a brain mask from the eddy-corrected data:

mask_fname = paste0(base_fname, "_mask.nii.gz")
if (!file.exists(mask_fname)) {
  fsl_bet(infile = dti0, outfile = mask_fname)
bet2 "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmp253ES1/file367e7e76fe3c.nii.gz" "./113-01-DTI_mask" 
[1] "113-01-DTI_mask.nii.gz"
mask = readnii(mask_fname)
ortho2(robust_window(dti0), mask, col.y = alpha("red", 0.5))

Running DTI Fitting as a cursor

Now that we have eddy current corrected our data, we can pass that result into dtifit to get FA maps and such:

outprefix = base_fname
suffixes = c(paste0("V", 1:3),
             paste0("L", 1:3),
             "MD", "FA", "MO", "S0")
outfiles = paste0(outprefix, "_", 
                  suffixes, get.imgext())
names(outfiles) = suffixes
if (!all(file.exists(outfiles))) {
  res = dtifit(infile = eddy_fname, 
               bvecs = as.matrix(b_vecs),
               bvals = b_vals, 
               mask = mask_fname,
               outprefix = outprefix)
dtifit --data="/Users/johnmuschelli/Dropbox/Projects/neuroc/DTI_analysis_fslr/113-01-DTI_eddy.nii.gz" --out="113-01-DTI" --mask="/Users/johnmuschelli/Dropbox/Projects/neuroc/DTI_analysis_fslr/113-01-DTI_mask.nii.gz" --bvecs="/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmp253ES1/file367e4220886f.txt" --bvals="/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmp253ES1/file367e65dd5153.txt" --verbose 

By default, the result of dtifit is the filenames of the resultant images. Here we will read in those images:

res_imgs = lapply(outfiles[c("FA", "MD")], readnii)

Plotting an FA map

Using the ortho2 function, you can plot the fractional anisotropy (FA) map