Scripts for processing brain MRI data

This library is a collection of Bash scripts, Slurm scripts, and R functions written for the processing of diffusion weighted imaging (DWI) and resting-state fMRI (rs-fMRI) data. The scripts start with just the raw DICOM images and perform steps up to network creation, based on the results from probtrackx2 for DWI and other methods for rs-fMRI. The Slurm scripts are appropriate for use on a high performance computing (HPC) system/cluster; I use the Texas Advanced Computing Center (TACC) systems.

The code has been written to work with projects following the Brain Imaging Data Structure (BIDS). That is, if you only supply a tar.gz file with DICOM images, the appropriate BIDS-compliant directories will be created. Note that this only applies to inputs or minimal processing (e.g., conversion to NIfTI); while there are some BIDS Derivatives proposals, nothing has found consistent use (to my knowledge). So I place some output directories which will be described later.

Table of Contents

• Requirements
  • Files/Formats and Directories
    • DICOM files
    • Basic QC
    • Parcellations
  • Software
• Installation
  • dcmtk
  • dcm2niix
  • jo
  • jq
  • Xvfb
  • FSL
  • Freesurfer
• Processing Steps
  • DWI
• Variables
• Known Issues
  • Slice acquisition times
    • GE
  • fslroi
  • FS to DTI registration script

Requirements

Files/Formats and Directories

The project directory is where all scripts will be run from; its variable is ${projdir} and is automatically set as the working directory. It is the top-level directory for your project; i.e., all relevant data should be accessible from here.

DICOM files

When running the initial scripts, you will need to provide one of the following:

1. A sourcedata directory, which contains BIDS-compliant subject directories and the DICOM files in ${target}_dicom.tar.gz within that directory tree. In this case, you would not use the --tgz option to dti_dicom2nifti_bet.sh.

   Here, the ${target} variable should follow the BIDS spec; for example (with optional information in square brackets):

   sub-<studyID>[_ses-<sessionID>][_acq-<acquisition>]_dwi_dicom.tar.gz
   sub-<studyID>[_ses-<sessionID>]_task-rest[_acq-<acquisition>][_run-<runID>]_dicom.tar.gz

2. A .tar.gz that you provide as input to the initial script, dti_dicom2nifti_bet.sh. This .tar.gz file MUST be either directly in ${projdir} or you must provide the full path. This will be renamed to ${target}_dicom.tar.gz (see above) and placed under the sourcedata directory.

Basic QC

To perform the most basic QC check — checking that image dimensions match — there must be a simple text file for each acquisition.

For example, if your DWI acquisition matrix is 96 x 96, there are 65 slices, and 32 diffusion-weighted volumes + 1 non-diffusion weighted volume:

cat ${projdir}/data/sizes/dwi_size.txt
96 96 65 33

If the files do not exist, the initial script will exit with an error. If the files do exist, and the dimensions do not match, the data will be moved to the ${projdir}/unusable directory.

Parcellations

Subject-specific parcellations will be used as the sources/targets of the network (at least for DTI tractography). The results from Freesurfer's recon-all should be in ${projdir}/freesurfer.

Software

In addition to good-quality T1-weighted and DWI data, some software requirements are:

• dcmtk for reading from the DICOM headers
• A recent version of dcm2niix
  • The version I used at the time of writing is v1.0.20181125 GCC4.8.5 (64-bit Linux)
• jo for writing out JSON files containing the parameters used for each tool. For example, it will record the -f value used with bet.
• jq also processes JSON files.
• FSL version >= 6.0.0
• The ImageMagick suite
  • Available in the repositories for Red Hat-based systems (RHEL, CentOS, Scientific Linux)
• Xvfb, the X Virtual Frame Buffer
• eddyqc (bundled with FSL since v6.0.0)
• Freesurfer version >= 5.3.0
  • Required for parcellation, the results of which will be used in the tractography step
• R is required for eddy-related movement QC. R is available on all major OS's. Necessary packages include:
  • optparse
  • data.table
  • ggplot2
  • gridExtra
• (Optional) A CUDA-capable GPU (for eddy_cuda and bedpostx_gpu)
  • If you don't have this, you can use eddy_openmp and the regular bedpostx instead; you'll have to change the code accordingly.
  • To run bedpostx on a SLURM system without a GPU, you will need the launcher utility

Installation

You can simply clone the repository and add it to your search path. Example:

git clone https://github.com/cwatson/mri_library.git
echo "export PATH=PATH:${PWD}/mri_library/bin" >> ~/.bash_profile

If you are on a cluster or other system that uses, e.g., .bashrc instead, you may have to add the path manually.

dcmtk

For CentOS 7, at least, this is available in the nux-dextop repository. If you don't already have this repo, run the following as root:

yum -y install http://li.nux.ro/download/nux/dextop/el7/x86_64/nux-dextop-release-0-5.el7.nux.noarch.rpm
yum install dcmtk\*

The version I use at the time of writing is v3.6.0.

If you are on Debian: sudo apt install dcmtk.

For other systems, you will probably have to install from source. See the DCMTK site for more information.

dcm2niix

I prefer to clone the Github repository.

cd /usr/local
git clone git://github.com/rordenlab/dcm2niix.git
cd dcm2niix
mkdir build && cd build
cmake ..
make install

jo

This utility can generate JSON from the command line. To install it, follow the instructions on the jo repository page. You will also need automake and autoconf.

cd /usr/local
git clone git://github.com/jpmens/jo.git
cd jo
autoreconf -i
./configure
make check
make install

If you are running CentOS 6, you will have to install the autoconf268 package, and then call autoreconf268 instead.

jq

This should be in repositories for all major Linux OS's. For both CentOS 6 and CentOS 7, it is in the epel repository, with versions v1.3.2 and v1.5.1, respectively.

Xvfb

There should be a package readily available on most systems. On CentOS, it is called xorg-x11-server-Xvfb.

FSL

To install the latest version of FSL (which is v6.0.0 as of October 2018), you simply run their installer. This requires that you already have an older version of FSL on your system.

cd ${FSLDIR}
python fslinstaller.py

Freesurfer

You can find download and install instructions at the Freesurfer wiki.

Processing Steps

DWI

The scripts will perform the following steps. Freesurfer's recon-all should be run before this (or before step 5, at least).

1. Run dti_dicom2nift_bet.sh to extract DICOM files and convert to NIfTI using dcm2niix, skullstrip, and create images for QC purposes.

   1. Renames and moves the tgz file to the correct, BIDS-compatible filename and location (if necessary)
   2. Extracts the DICOM files from the tgz file and convert to NIfTI (using dcm2niix)
   3. Moves the nii.gz, bvecs, bvals, and json files to the appropriate subject directory under rawdata.
   4. Checks the image dimensions against the study's dimensions (specified by the user in a text file) via qc_basic.sh. If this fails, the subject's data are moved to a directory called unusable.
   5. If there are multiple b0 volumes, they will be averaged when creating nodif.nii.gz
   6. Skullstrips the data using bet.
   7. Runs dti_qc_bet.sh to generate screenshots of the skullstripped data. See example images.

   dti_dicom2nifti_bet.sh -s s001 --acq multishell --tgz s001_dicom.tar.gz
   # Creates the files:
   ${projdir}/rawdir/sub-s001/dwi/sub-s001_acq-multishell_dwi.{nii.gz,bval,bvec,json}
   # Creates the directory:
   ${projdir}/tractography/sub-s001/dwi/qc/bet/

2. (manual) Check the quality of bet by viewing the images in ${resdir}/qc/bet.

   1. Re-run Step 1 if the skullstrip wasn't acceptable. Do this by passing the --rerun and -t|--threshold options to dti_dicom2nifti_bet.sh.

   For example,

   eog tractography/sub-s001/dwi/qc/bet/*.png
   dti_dicom2nifti_bet.sh -s s001 --acq multishell --rerun -t 0.4

3. Run eddy via dti_eddy.sh.

   1. Also calculates eddy-specific QC metrics via eddy_quad from eddyqc.

   dti_eddy.sh -s s001 -acq multishell

4. If you do not have a GPU but do have a SLURM scheduler, run BEDPOSTX via dti_bedpostx_run.sh.

   1. You will then also need to run dti_bedpostx_postproc.sh.
   2. If you do have a GPU, you can run bedpostx_gpu on ${projdir}/${resdir}.
   3. If you have a system with an SGE scheduler, you can run bedpostx normally.

5. Run the setup script dti_reg_FS_to_diff.sh to register the Freesurfer parcellation to diffusion space.

6. Check the quality of the registration/parcellation by viewing the images produced from Step 5. See an example image.

7. Run the tractography via dti_probtrackx2_run.sh. You can choose to use the GPU version, or you can run it in parallel (which requires GNU parallel.

8. Create connectivity matrices in which the entries are the estimated streamline counts between all region pairs, using the R script fsl_fdt_matrix.R. This creates a file called fdt_network_matrix (assuming you haven't already run probtrackx2 in network mode).

9. Create connectivity matrices in which the entries are the mean of some microstructural measure (e.g., FA) along the top N% of streamlines. The default is to use the top 10%.

Variables

The following variables are exported by setup_vars.sh and are used to create the appropriate directory and filename targets that conform to the BIDS standard See the code block following the list for examples.

projdir

The project's top-level directory. All preprocessing scripts must be called from this directory.

target

The character string for the subject (plus the session and acquisition, if applicable).
Directories and filenames will both be generated using this variable.

srcdir

The directory that holds the "source" data. Here, "source" indicates data that came directly from the scanner; i.e., DICOM files, PAR/REC files (for Philips data), Siemens mosaic, etc.

rawdir

The directory that will hold the "raw" data. Here, "raw" simply indicates that no preprocessing has been applied to them.

fs_sub_dir

The name of the directories containing Freesurfer results. Since Freesurfer's SUBJECTS_DIR is "flat" (i.e., all data must be in a single directory), if the study is longitudinal then the _ses-${session} will be in the directory names.

resdir

The directory where results will be stored. For DWI data, this will be ${projdir}/tractography (containing the results from both bedpostx and probtrackx2).

For example, in a theoretical test-retest study, the outputs might look like the following. The directories live directly under ${projdir}.

${target}       sub-s001_ses-retest_acq-highres_T1w
                sub-s001_ses-retest_acq-multishell_dwi
                sub-s001_ses-retest_task-rest_bold

${srcdir}       sourcedata/sub-s001/ses-retest/{anat,dwi,func}
                    `- ${target}_dicom.tar.gz
${rawdir}       rawdata/sub-s001/ses-retest/dwi
                    `- ${target}.{bvals,bvecs,json,nii.gz}
${resdir}       tractography/sub-s001/ses-retest/dwi
${fs_sub_dir}   freesurfer/sub-s001_ses-retest

Known Issues

Slice acquisition times

It seems that GE and Philips do not record the SliceTiming information which is necessary for eddy's slice-to-volume motion correction. There isn't a foolproof way of getting this information, aside from at the scanner console itself.

GE

From this thread, helpful DICOM tags are:

• 0020,1002 ImagesInAcquisition -- should be equal for all in a single acquisition
• 0020,9057 InStackPositionNumber -- the slice number for each volume
• 0020,0013 InstanceNumber -- slice number in the whole acquisition
• 0020,1041 SliceLocation

That thread links to another thread which references some more DICOM tags, and links to a Github repo that may be able to find out this information.

fslroi

In FSL v6.0.0, fslroi has undesired behavior; for some files, it remapped voxel values so that the nodif images were completely wrong. A temporary workaround is to use fslmaths to change the image type of dwi_orig.nii.gz to float. This bug should be fixed in v6.0.1.

FS to DTI registration script

The QC script for Freesurfer to diffusion space relies on xvfb-run. The argument -n allows you to specify a specific server number, but sometimes this still results in errors. The SLURM script (dti_reg.launcher.slurm) will automatically increment this from 0 for each subject, but if you have more than 100 subjects processed at once, there might be an error.