- Home
- Neuroconductor Tutorials
- DTI analysis rcamino hcp
DTI Analysis using rcamino for HCP data
John Muschelli
2021-02-16
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:
data.nii.gz
- a 4D image of the DWI data.nodif_brain_mask.nii.gz
- A brain mask of the DTI databvals
- a text file with the b-valuesbvecs
- 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)