DTI Analysis using rcamino for HCP data

All code for this document is located at here.

Resources and Goals

Much of this work has been adapted by the Tutorial for DTI analysis from ISMRM 2015: http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.ISMRM2015. Also, some of the model fitting (such as using gradient information) has been taken from http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.HCP.

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

Reading in the Data

First, we download the data from HCP. You must have your access keys set (see Getting Data from the Human Connectome Project (HCP)).

We will use the neurohcp package to download one subject data.

library(neurohcp)
hcp_id = "100307"
r = download_hcp_dir(
  paste0("HCP/", hcp_id, "/T1w/Diffusion"), 
  verbose = FALSE)
print(basename(r$output_files))
[1] "bvals"                   "bvecs"                  
[3] "data.nii.gz"             "grad_dev.nii.gz"        
[5] "nodif_brain_mask.nii.gz"

It contains 4 files:

  1. data.nii.gz - a 4D image of the DWI data.
  2. nodif_brain_mask.nii.gz - A brain mask of the DTI data
  3. bvals - a text file with the b-values
  4. bvecs - a text file with the b-vectors as the first 3 columns.

Creating

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)
camino_set_heap(heap_size = 10000)
outfiles = r$output_files
names(outfiles) = nii.stub(outfiles, bn = TRUE)
scheme_file = camino_fsl2scheme(
  bvecs = outfiles["bvecs"], bvals = outfiles["bvals"],
  bscale = 1)
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/fsl2scheme -bvecfile /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d23a16d742//bvecs -bvalfile /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d23a16d742//bvals -bscale 1   >  /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d27515acf9.scheme

The imaging scheme contains measurements at b=5, b=1000, b=2000, and b=3000 s / mm^2.

Subsetting data

By selecting a subset of the measurements, we can reduce processing time and memory requirements. Also, the high b-value shells aren’t optimal for estimating the diffusion tensor. So we’ll select data from the b=5 and b=1000 shells, which is still higher angular resolution than most DTI (90 directions).

If you have sufficient RAM, you can load the whole data set and extract a subset:

camino_ver = packageVersion("rcamino")
if (camino_ver < "0.5.2") {
  source("https://neuroconductor.org/neurocLite.R")
  neuroc_install("rcamino")  
}
sub_data_list = camino_subset_max_bval(
  infile = outfiles["data"],
  schemefile = scheme_file,
  max_bval = 1500,
  verbose = TRUE) 
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/split4dnii -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d23a16d742/data.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d2256a7b02_
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/image2voxel -imagelist '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d271b1d37e.txt' -outputfile '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d251b42f82.Bfloat' -outputdatatype float
sub_data = sub_data_list$image
sub_scheme = sub_data_list$scheme

This process may take a while, and using the RNifti package may be quicker.

nim = RNifti::readNifti(outfiles["data"])
sub_data = tempfile(fileext = ".nii.gz")
sub_scheme_res = camino_subset_max_bval_scheme(
  schemefile = scheme_file, max_bval = 1500,
  verbose = TRUE)
nim = nim[,,, sub_scheme$keep_files]
RNifti::writeNifti(image = nim, file = sub_data)
sub_scheme = sub_scheme_res$scheme
rm(list = "nim");
for (i in 1:10) gc();

Fit the diffusion tensor

# wdtfit caminoProc/hcp_b5_b1000.Bfloat caminoProc/hcp_b5_b1000.scheme \
# -brainmask 100307/T1w/Diffusion/nodif_brain_mask.nii.gz -outputfile caminoProc/wdt.Bdouble
# 
mod_file = camino_modelfit(
  infile = sub_data, scheme = sub_scheme, 
  mask = outfiles["nodif_brain_mask"], 
  gradadj = outfiles["grad_dev"],
  model = "ldt_wtd")
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/modelfit -inputfile '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d251b42f82.Bfloat' -outputfile '/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d24551fec4.Bdouble' -inputdatatype float -schemefile /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d22e0e2fe5.scheme -bgmask /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d23a16d742/nodif_brain_mask.nii.gz -maskdatatype float -model ldt_wtd -gradadj /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d23a16d742/grad_dev.nii.gz

Getting FA vlaues

# fa -inputfile caminoProc/wdt_dt.nii.gz -outputfile caminoProc/wdt_fa.nii.gz
fa_img = camino_fa_img(
  infile = mod_file,
  header = outfiles["nodif_brain_mask"],
  retimg = FALSE)
cat '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d24551fec4.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//RtmpP1xrb3/file100d216df835e.Bdouble'
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/voxel2image -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d216df835e.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d23a16d742/nodif_brain_mask.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d22c299c0d_ -components 1 -gzip

Visualizing FA images

We want to read the FA image into R:

fa_nii = readnii(fa_img)

In order to visualize the values, we are going to read in the mask so that we don’t visualize non-brain values:

mask = readnii(outfiles["nodif_brain_mask"])
hist(mask_vals(fa_nii, mask = mask), breaks = 1000)

Using ortho2, we can visualize these FA maps:

ortho2(fa_nii)

Getting MD vlaues

# md -inputfile caminoProc/wdt_dt.nii.gz -outputfile caminoProc/wdt_md.nii.gz
md_img = camino_md_img(
  infile = mod_file,
  header = outfiles["nodif_brain_mask"],
  retimg = FALSE)
cat '/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d24551fec4.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//RtmpP1xrb3/file100d2d5e73ba.Bdouble'
/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rcamino/camino/bin/voxel2image -inputfile /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d2d5e73ba.Bdouble -header /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/RtmpP1xrb3/file100d23a16d742/nodif_brain_mask.nii.gz -outputroot /var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpP1xrb3/file100d233ffbf31_ -components 1 -gzip

Visualizing MD images

We want to read the MD image into R:

md_nii = readnii(md_img)
hist(mask_vals(md_nii, mask = mask), breaks = 1000)

md2 = md_nii
md2[ md2 < 0] = 0
hist(mask_vals(md2, mask = mask), breaks = 1000)