diannepat/bip repository overview

⁠BIP PROJECT

Dianne Patterson, University of Arizona, SLHS Dept.
Created Feb 23, 2017
Updated April 23, 2019

BIP (bidirectional iterative parcellation) uses FSL DWI processing to not only characterize each tract, but also the connected grey matter at each end of a tract:

Patterson, D. K., Van Petten, C., Beeson, P., Rapcsak, S. Z., & Plante, E. (2014). Bidirectional iterative parcellation of diffusion weighted imaging data: separating cortical regions connected by the arcuate fasciculus and extreme capsule. NeuroImage, 102 Pt 2, 704–716. http://doi.org/10.1016/j.neuroimage.2014.08.032⁠

Special thanks to Patricia Klobusiakova without whose eagle eye and sharp mind, there would be many more errors in this set of scripts.

⁠REVISION CONTROL

The bip project is under version control using git.
git clone https://[email protected]/dpat/bipbids.git⁠ to create a copy of the bipbids directory.

⁠DIRECTORY DESCRIPTION

The bip project contains bash scripts, lists, and standard space image masks used to run bip.

The scripts primarily run FSL 5.0.9 commands (these should be backward compatible to ~5.0.7). Everything was written to run in an Ubuntu 14 docker container using BIDS compliant file structure and naming: http://bids.neuroimaging.io/⁠

⁠REF subdirectory:

This subdirectory contains non-scripts organized into two subdirectories.
LISTS contains lists of defined tracts and rois to facilitate iteration.
IMAGES contains binary image masks for particular tracts (endpoint rois and termination masks).

⁠DOCKER

The bipbids directory contains a Dockerfile and run.py at the top level. If you have docker running on Mac or Linux, the easiest way to get the current container is to pull it from dockerhub like this:

Alternatively, you can build the container from the downloaded bipbids directory containing the Dockerfile:

⁠docker_wrap.sh

To facilitate running bip on your local machine, use docker_wrap.sh.

If docker_wrap.sh is in your path, the above call will show you help.

By examining the contents of docker_wrap.sh, you can see examples of the different calls that can be made to the docker container.

You must be in your main bids directory for docker_wrap.sh to work correctly.

Always choose a stage to run: setup, prep, bip or stats.
If you choose bip, you must specify a tract to run.
Additional arguments are assumed to be subject numbers. All of the steps run at the participant level, except stats, which runs at the group level.

⁠SINGULARITY

To facilitate running bip on the HPC, you must build a singularity container. The singularity recipe is included. The singularity recipe builds off the docker container and then adds CUDA 8 binaries for running eddy, bedpostX and probtrackx2. These are MUCH faster (when they work). They do work on the University of Arizona Ocelote computer with the P100 Nvidia GPUs.

To use a singularity container, you must first build it, which means you need access to a linux machine where you have root privileges (i.e., NOT the HPC).

Upload the file Singularity and build the container:

Add the container to the machine where you wish to run and run it with a singularity command like this:

This command says to run bip.simg (the singularity container) on subjects under the present working dir in Data and write to Data/derivatives. We are running for one subject, sub-327, we are running the bip stage of processing on the arc_l tract. We do want to use the gpu. Look at the example run commands in docker_wrap.sh for other examples. Run commands are a bit different for singularity and docker, but a lot of the bids structures are the same.

N.B. Testing indicates that the results of probtrackx2 and probtrackx2_gpu are similar, but somewhat different...so DO NOT mix the two approaches! Pick one and stick to it.

⁠Sample Dataset

https://arizona.box.com/v/bip⁠

⁠PARTICIPANT LEVEL STEPS

examples
Run setup for all subjects:

Run setup for two subjects: sub-001 and sub-329:

Run prep for all subjects:

Run prep for one subject: sub-001:

Run bip on the left arcuate for all subjects:

Run bip on the right arcuate for sub-001:

Available tracts for bip are: arc_l arc_r aslant_l aslant_r b3tob3_ih cb2th_lr cb2th_rl cst_l cst_r echo extcap_l extcap_r fr2th_l fr2th_r ilf_l ilf_r iof_l iof_r mdlf_l mdlf_r par2th_l par2th_r echo slf2_l slf2_r tmp2th_l tmp2th_r unc_l unc_r vof_l vof_r vwfa2vwfa_ih w5tow5_ih

ih=interhemispheric

⁠GROUP LEVEL STEP

Concatenate statistics from individual subjects: docker_wrap.sh stats

⁠OVERVIEW OF STEPS

⁠GETTING READY

Your file naming and directory structure are presumed to be bids compliant:

Example With fieldmaps:

|-- sub-001
|   |-- anat
|   |   |-- sub-001_T1w.json
|   |   `-- sub-001_T1w.nii.gz
|   |-- dwi
|   |   |-- sub-001_acq-AP_dwi.bval
|   |   |-- sub-001_acq-AP_dwi.bvec
|   |   |-- sub-001_acq-AP_dwi.json
|   |   `-- sub-001_acq-AP_dwi.nii.gz
|   `-- fmap
|       |-- sub-001_dir-PA_epi.json
|       |-- sub-001_dir-PA_epi.nii.gz
|       |-- sub-001_magnitude1.json
|       |-- sub-001_magnitude1.nii.gz
|       |-- sub-001_magnitude2.json
|       |-- sub-001_magnitude2.nii.gz
|       |-- sub-001_phasediff.json
|       `-- sub-001_phasediff.nii.gz

Example Without Fieldmaps

sub-001
|-- anat
|   |-- sub-001_T1w.json
|   `-- sub-001_T1w.nii.gz
`-- dwi
    |-- sub-001_dwi.bval
    |-- sub-001_dwi.bvec
    |-- sub-001_dwi.json
    `-- sub-001_dwi.nii.gz

In other words the names of the dwi files must at least contain subject and dwi, but additional information may be added (e.g., acqp-AP)

You will need at least a structural image and a dwi image with bvals and bvecs. They should be named and organized as above.

Optionally, if you have the fieldmaps (magnitude images, phase image and reverse phase encoded B0 image), place these in the fmap directory. These fieldmaps can be used with topup and eddy to improve distortion correction. Currently this handles the standard Siemens fieldmaps.

To take advantage of fieldmap processing, you will need the following text files in your bids_dir. The content of these files is explained more below. See sub-001 dataset for examples. acqp_??.txt (e.g., acqp_AP.txt)
acqparams.txt
index_??.txt (e.g., index_AP.txt)
pe.txt

⁠Lesion Mask

Also include a lesion mask if you have one (1's in the lesion, 0's elsewhere). If setup.sh finds a lesion image named as follows and in the anat directory, e.g. anat/sub-001_label-lesion_roi.nii.gz, the lesion image will be used to improve the registration of the structural image into MNI space.

⁠RUNNING THE SCRIPTS

Your file naming and directory structure are presumed to be bids compliant (see above). The T1w image will be defaced to facilitate data sharing. The docker container expects you to choose either setup, prep or bip. If you choose bip, then you must also select a tract from the available options.

⁠setup

The setup step runs setup.sh to do defacing, cropping and skull stripping. After this step you can manually check and edit masks. I particularly recommend that you check the masks for the B0 image.

optiBET.sh is used to improve the brain masks: Lutkenhoff, E. S., Rosenberg, M., Chiang, J., Zhang, K., Pickard, J. D., Owen, A. M., & Monti, M. M. (2014). Optimized Brain Extraction for Pathological Brains (optiBET). PLoS ONE, 9(12), e115551–13. http://doi.org/10.1371/journal.pone.0115551⁠

The function Deface works on the T1w image: reorient, bias-correct, crop, brain mask and deface. Original identifiable images are stored in sourcedir.

PrepDWI creates the B0 image (nodif and the default mask). If you have the fmap directory, then this includes the function Top which does topup and applytopup. This takes ~25 minutes to run on a 2013 mac pro. topup calculates the distortion and movement in the DWI images. Top looks for ${bids_dir}/acqparams.txt

B0 images are identified by examining the bvals file. If you do not have fieldmaps, then one B0 image is identified and used. If you have fieldmpas, then two B0 images are used .for blip 1: the first and the last found in the file.

applytopup applies the distortion and movement corrections to the B0 images. applytopup works like this: we have 2 images: blip1 consists of 2 B0 blip up volumes. Call these volumes 1 and 2. blip2 consists of 2 B0 blip down volumes. Call these volumes 3 and 4. applytopup combines the first blip1 volume (1) with the first blip2 volume (3) to create a volume (1+3). Then the 2nd blip1 image (2) is combined with the second blip2 image (4) to create a volume (2+4). These two new volumes are stored in ${output_dir}/${subj}/dwi/b0_hifi.nii.gz.

⁠Manual Checks

setup.sh creates all the masks that you might wish to check: T1w brain mask, T1w defacing mask, dwi mask, and if you have an fmap directory, magnitude mask. If you need lesion_masks, they will be added to the directory at this endpoint (details to be worked out).

⁠prep

prep.sh will apply your corrected masks, generate registrations, run dtifit and BedpostX.

prep.sh corrects the dwi image (with eddy_correct or eddy depending on your optional files), then runs epi_reg to register your DWI B0 image into structural space (~ 5 minutes). Finally, prep.sh will run dtifit and bedpostX (BedpostX can take days to run. It benefits from a good graphics card and more CPUs. More directions take longer to process). The shortest time I've seen is about 5 hours. If you have the fmap directory, then prep.sh expects parameter files in the bids_dir: acqp_AP.txt acqparams.txt index_AP.txt pe.txt These are explained in more detail below.

⁠bip

bip.sh expects a tract. It will create the tract subdirectory in the output dwi dir, run probtrackx2 iteratively until the seed stabilizes, and calculate some statistics on the final standard space output. Specifically, it'll generate a tsv file in the stats dir for the tract (volume, scalar dti measures), and a second tsv file for the endpoints (volume and center of gravity). tracts are trimmed before performing statistics so they do not overlap with GM or CSF.

Available tracts are: arc_l arc_r aslant_l aslant_r b3tob3_ih cb2th_lr cb2th_rl cst_l cst_r extcap_l extcap_r fr2th_l fr2th_r ilf_l ilf_r iof_l iof_r mdlf_l mdlf_r par2th_l par2th_r slf2_l slf2_r tmp2th_l tmp2th_r unc_l unc_r vof_l vof_r vwfa2vwfa_ih w5tow5_ih

ih=interhemispheric, arc=arcuate, aslant=aslant, b3tob3 is brocas to brocas (interhamispheric), cb2th_lr=left cerebellum to right thalamus, cb2th_rl=right cerebellum to left thalamus, cst=corticospinal tract, extcap=extreme capsule, fr2th=frontal lobe to thalamus, ilf=inferior longitudinal fasciculus, iof=inferior occipital fasciculus, mdlf=middle longitudinal fasciculus, par2th=parietal lobe to thalamus, slf2=superior longitudinal fasciculus part 2, tmp2th=temporal lobe to thalamus, unc=uncinate, vof=vertical occipital fasciculus, vwfa2vwfa_ih= visual word form area to visual word form area (interhemispheric), w5tow5=posterior language area to posterior language area (interhemispheric).

People disagree about tract names and trajectories. The particulars of each available tract have been defined based on the literature and extensive experimentation. You may or may not agree with these choices. It is possible to create your own masks...but the documentation is not here yet to make that easy.

⁠stats

When you have finished running all the tracts of interest, you can concatenate the statistics for each subject into two files: one for endpoints and one for the tracts

⁠FIELDMAP AND TOPUP PARAMETER FILES

If using field maps and a reverse-phase encode image, certain parameters need to be set:

⁠FIELD MAPS

https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FUGUE/Guide#SIEMENS_data⁠

Acquiring and Using Fieldmaps: https://lcni.uoregon.edu/kb-articles/kb-0003⁠

⁠DIFFERENCE IN ECHO TIMES OF FIELDMAPS

This value is needed for fsl_prepare_fieldmap. This can be calculated automatically from the json files associated with the fieldmap images and is called ECHO_DIFF in the prep.sh script. By default this is 2.46 for Siemens, and refers to characteristics of the fieldmaps.

To calculate the difference in echo times, you need the TE for the first and 2nd echo: TE=4.92 (first magnitude image; e.g., sub-001_magnitude1.nii.gz); TE=7.38 (phasediff and 2nd magnitude image; sub-001_phasediff.nii.gz and sub-001_magnitude2.nii.gz); 1 echo (First mag map); 2 echoes (phase map and 2nd mag); 7.38-4.92=2.46. prep.sh looks up these values in the json files.

⁠PHASE ENCODE DIRECTION OF DWI

epi_reg requires a phase encode direction (--pedir) for the primary dwi image if you have fmap. We read ${bids_dir}/pe.txt to get this value (e.g., -y). Unfortunately, the JSON files use i,j,k and FSL uses x,y,z to record this value. "PhaseEncodingDirection": "j-" in the json file corresponds to -y in FSL's terminology (-y=AP).

⁠EFFECTIVE ECHO SPACING OF DWI

To run epi_reg with field maps, we need the effective echo spacing of the dwi image being processed. This is called ECHO_SPACE in the prep.sh script. The echo spacing should be listed in the parameters of the scan (e.g., 0.94 ms for my 32 direction B=1000 dti scans). Effective echo spacing=echo spacing divided by the grappa (acceleration, "ParallelReductionFactorInPlane") factor (2 in this case) and then divide by 1000 to get units in seconds instead of ms.

e.g., (0.94/2)/1000=0.00047

Effective Echo Spacing is sometimes reported directly in the BIDS json file, e.g., in sub-001_acq-AP_dwi.json "EffectiveEchoSpacing": 0.00047001,

⁠TOTAL READOUT TIME OF DWI

Total Readout time is used by topup, and appears in the acqp and acqparams files. There are several ways to compute total readout time.

The Topup Users guide says:
"If your readout time is identical for all acquisitions you don't necessarily have to specify a valid value in this column (you can e.g. just set it to 1), but if you do specify correct values the estimated field will be correctly scaled in Hz, which may be a useful sanity check."
https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide/⁠

FSL calculates total readout time as the echo spacing (dwell time) * the number of phase encoding steps. If you use grappa (may be called "ParallelReductionFactorInPlane" or "AccelFactPE" in the JSON file), divide the number of phase encoding steps (128 for sub-001) by the grappa factor (2 for sub-001). "PhaseEncodingLines": 128 "ParallelReductionFactorInPlane": 2 (n.b. sometimes called "AccelFactPE": 2) "DerivedVendorReportedEchoSpacing": 0.000940019,

e.g., (128/2)*.000940019 sec =0.0602 sec

From the JSON file, we see that dcm2niix has performed a similar (but apparently slightly different) calculation for sub-001: "TotalReadoutTime": 0.0596912

I use the calculated value (but either is probably fine)

0 -1 0 0.0602 (acqp_A-P.txt...this line is repeated 32 times, once for each DWI volume)

acqparams.txt: Used for topup. Specifies the values for the 4 blip_down and blip_up B0 volumes

index_A-P.txt: Used by eddy to index all the values in acqp_A-P.txt (all 1's)

⁠REVERSE PHASE ENCODE DWI DATA

DWI data are acquired in one direction (e.g., A->P). Because the echo time for DWI images is so long, this means there is tremendous distortion. If you look at A->P images in axial view, you'll see stretchy eyeballs. You can acquire a second image (1-2 B0s) in the opposite phase-encode direction (P->A), this will take ~ 45 seconds. This second image experiences distortion in the opposite direction of the original A->P images. If you look at the P->A image in axial view, you'll see it has squishy eyeballs. BIDS recommends saving this reverse phase encode image with the fieldmaps: e.g., sub-001_dir-PA.nii.gz

Default Siemens A-P: stretchy eyeballs, (negative phase encode, blip down -1) Reverse phase encode B0 P-A: squishy eyeballs, (positive phase encode, blip_up 1)

In FSL we can use topup to calculate the distortion parameters based on both images. We can then use applytopup to apply these parameters to another image. https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/ApplyTopupUsersGuide⁠

N.B. My experiments suggest it is better to correct the A-P data (less noise, fewer spikes and holes in weird places, better eyeball reconstruction)

Finally, we feed the topup results into eddy for correcting eddy currents and movement. https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy⁠

Place the following parameter files describing your primary dwi image in your bids_dir for topup processing: acqp.txt acqparams.txt index.txt pe.txt (actually this is for eddy...but don't forget it!)

The acqp file is used by eddy, as described here: http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY/UsersGuide⁠ This is a text-file describing the acquisition parameters for the different images in --imain. The format of this file is identical to that used by topup (though the parameter is called --datain there). The value is repeated 32 times for the 32 volumes. The -1 indicated the phase encode direction for y in A-P. The last value ** Total Readout time of DWIs ** is described above.