DTI Analysis in Camino using rcamino

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 reconstruction: http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.DTI. We will show you a few steps that have been implemented in rcamino: camino_pointset2scheme, camino_modelfit, camino_fa, camino_md, and camino_dteig.

Data Location

The data located in this tutorial is located at http://cmic.cs.ucl.ac.uk/camino//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

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")
download.file("http://cmic.cs.ucl.ac.uk/camino//uploads/Tutorials/example_dwi.zip",
              destfile = tfile)
files = unzip(zipfile = tfile, exdir = tdir, overwrite = TRUE)

Making b-vectors and b-values

As dtifit requires the b-values and b-vectors to be separated, and this data has b-values of \(1000\) when the b-vectors is not zero. This is very important and you must know where your b-values and b-vectors are when doing your analyses and what units they are in.

library(rcamino)
b_data_file = grep("[.]txt$", files, value = TRUE)
scheme_file = camino_pointset2scheme(infile = b_data_file,
                                bvalue = 1e9)
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/pointset2scheme -inputfile '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/grad_dirs.txt' -bvalue 1000000000 -outputfile /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e4f033bb7.scheme

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.

img_fname = grep("4Ddwi_b1000", files, value = TRUE)
img = neurobase::readnii(img_fname)
ntim(img)
[1] 33
grads = readLines(b_data_file)
length(grads)
[1] 33
# cleanup
rm(list= "img"); gc()
            used   (Mb) gc trigger   (Mb) limit (Mb)   max used   (Mb)
Ncells   2777394  148.4    4907837  262.2         NA    4907837  262.2
Vcells 311498041 2376.6  839477736 6404.8      65536 1048353822 7998.4

Running Image Conversion

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:

float_fname = camino_image2voxel(infile = img_fname, 
                                outputdatatype = "float")
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/image2voxel -inputfile '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/4Ddwi_b1000.nii.gz' -outputfile '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e3161603a.Bfloat' -outputdatatype float

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.

Fit the diffusion tensor

mask_fname = grep("mask", files, value = TRUE)
model_fname = camino_modelfit(
  infile = float_fname,
  scheme = scheme_file,
  mask = mask_fname,
  outputdatatype = "double"
  )
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/modelfit -inputfile '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e3161603a.Bfloat' -outputfile '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e23738c2b.Bdouble' -inputdatatype float -schemefile /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e4f033bb7.scheme -bgmask /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/brain_mask.nii.gz -maskdatatype float -model dt

Getting FA vlaues

fa_fname = camino_fa(infile = model_fname)
cat '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e23738c2b.Bdouble' |  /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/fa -inputmodel dt -outputdatatype double > '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e6e247981.Bdouble'

Converting FA values back into an image

library(neurobase)
fa_img_name = camino_voxel2image(infile = fa_fname, 
                            header = img_fname, 
                            gzip = TRUE, 
                            components = 1)
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/voxel2image -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e6e247981.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/4Ddwi_b1000.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e24308354_ -components 1 -gzip 
fa_img = readnii(fa_img_name)

Converting with piping

We can chain Camino commands using the magrittr pipe operation (%>%):

library(magrittr)

Attaching package: 'magrittr'
The following object is masked from 'package:R.utils':

    extract
The following object is masked from 'package:R.oo':

    equals
fa_img2 = model_fname %>% 
  camino_fa() %>% 
  camino_voxel2image(header = img_fname, gzip = TRUE, components = 1) %>% 
  readnii
cat '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e23738c2b.Bdouble' |  /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/fa -inputmodel dt -outputdatatype double > '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e1c674c44.Bdouble'
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/voxel2image -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e1c674c44.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/4Ddwi_b1000.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e41c012ed_ -components 1 -gzip 
all.equal(fa_img2, fa_img2)
[1] TRUE

Visualizing FA images

Using ortho2, we can visualize these FA maps:

ortho2(fa_img)

Getting MD vlaues

Similar to getting FA maps, we can get mean diffusivity (MD) maps, read them into R, and visualize them using ortho2:

md_img = model_fname %>% 
  camino_md() %>% 
  camino_voxel2image(header = img_fname, gzip = TRUE, components = 1) %>% 
  readnii
cat '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e23738c2b.Bdouble' |  /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/md -inputmodel dt -outputdatatype double > '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e61a0ae7b.Bdouble'
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/voxel2image -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e61a0ae7b.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/4Ddwi_b1000.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e7c1f3f4f_ -components 1 -gzip 
ortho2(md_img)

Export DTs to NIfTI

Using camino_dt2nii, we can export the diffusion tensors into NIfTI files. We see the result is the filenames of the NIfTI files, and that they all exist (otherwise there’d be an errors.)

nifti_dt = camino_dt2nii(
  infile = model_fname, 
  inputmodel = "dt",
  header = img_fname, 
  gzip = TRUE
)
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/dt2nii -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/file1422e23738c2b.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpnacR6C/4Ddwi_b1000.nii.gz -inputmodel dt -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e5a2024ae_ -gzip 
stopifnot(all(file.exists(nifti_dt)))
print(nifti_dt)
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e5a2024ae_exitcode.nii.gz"
[2] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e5a2024ae_lns0.nii.gz"    
[3] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpnacR6C/file1422e5a2024ae_dt.nii.gz"

We can read these DT images into R again using readnii, but we must set drop_dim = FALSE for diffusion tensor images because the pixel dimensions are zero and readnii assumes you want to drop “empty” dimensions

dt_imgs = lapply(nifti_dt, readnii, drop_dim = FALSE)