Resources and Goals

Much of this work has been adapted by the FSL guide for DTI: https://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 Location

The data located in this tutorial is located at http://camino.cs.ucl.ac.uk/uploads/Tutorials/example_dwi.zip. It contains 3 files:

1. 4Ddwi_b1000.nii.gz - a 4D image of the DWI data.
2. brain_mask.nii.gz - A brain mask of the DTI data
3. grad_dirs.txt - a 3 column text file with the b-vectors as the first 3 columns. The b-values are either 0 or 1000 depending on the b-vector.

Reading in the Data

First, we download the data into a temporary directory the unzip it:

tdir = tempdir()
tfile = file.path(tdir, "example_dwi.zip")
result = try({download.file("http://camino.cs.ucl.ac.uk/uploads/Tutorials/example_dwi.zip",destfile = tfile, mode = "wb")}, silent = TRUE)
if (inherits(result, "try-error")) {
  download.file("https://github.com/muschellij2/fslr/raw/master/inst/extdata/example_dwi.zip", destfile = tfile, mode = "wb")
}
out = unzip(zipfile = tfile, exdir = tdir, overwrite = TRUE)

library(fslr)

## Loading required package: oro.nifti

## oro.nifti 0.11.4

## Loading required package: neurobase

b_data_file = grep("[.]txt$", out, value = TRUE)
b_vecs = read.delim2(b_data_file, header = FALSE)
b_vecs = as.matrix(b_vecs)
class(b_vecs) = "numeric"

b_vals = rep(1000, nrow(b_vecs))
all_zero = apply(b_vecs  == 0, 1, all)
b_vals[all_zero] = 0

img = grep("4Ddwi_b1000", out, value = TRUE)
nim = readnii(img)
n_timepoints = dim(nim)[4]
stopifnot(nrow(b_vecs) == n_timepoints)

if (have.fsl()) {
  print(fsl_version())
}

## [1] "6.0.4"

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:

if (have.fsl()) {
  outfile = tempfile(fileext = ".nii.gz")
  ret = eddy_correct(infile = img, outfile = outfile, 
                     retimg = TRUE, reference_no = 0)
}

## eddy_correct "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpv3v5VZ/4Ddwi_b1000.nii.gz"  "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpv3v5VZ/file109232178d3e5" 0;

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.

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:

if (have.fsl()) {
  mask_fname = grep("mask", out, value = TRUE)
  res = dtifit(infile = ret, bvecs = b_vecs, 
               bvals = b_vals, mask = mask_fname)
}

## dtifit --data="/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpv3v5VZ/file1092358967c7d.nii.gz" --out="/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpv3v5VZ/file109237fc71ffa" --mask="/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmpv3v5VZ/brain_mask.nii.gz" --bvecs="/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpv3v5VZ/file109236a04979f.txt" --bvals="/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmpv3v5VZ/file10923537e7e19.txt" --verbose

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

if (have.fsl()) {
  res_imgs = lapply(res, readnii)
}

if (have.fsl()) {
  ortho2(res_imgs$FA)
}

if (have.fsl()) {
  ortho2(res_imgs$MD)
}

if (have.fsl()) {
  double_ortho(res_imgs$FA, res_imgs$MD)
}

if (have.fsl()) {
  mask = readnii(mask_fname)
  df = data.frame(FA = res_imgs$FA[ mask == 1], MD = res_imgs$MD[ mask == 1] )
  plot(df)
}

Fitting bedpostx

According to https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/UserGuide, the bedpostx function from FSL stands for “Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques. The X stands for modelling Crossing Fibres”. It also states “bedpostx takes about 15 hours to run”. We have implemented xfibres, the function bedpostx calls. Running it with the default options, the command would be:

xfibres(infile = outfile, 
        bvecs = b_vecs,
        bvals = b_vals,
        mask = mask_fname)

We are currently implementing probtracx2 and an overall wrapper for all these functions to work as a pipeline.