diannepat/bip2 repository overview

⁠BIP2 Docker Image

Dianne Patterson, University of Arizona, SLHS Dept.
Created March 24, 2020 Updated April 5, 2020

⁠What is BIP

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.

⁠What's New

• BIP2 is an enhanced version of BIP project (diannepat/bip):
• BIP2 now handles lesion masks for segmentation and for probtrackX.
• BIP2 no longer does defacing nor does it expect defacing. You may pass it data that has been defaced, or data that has not been defaced.
• BIP2 creates separate derivatives directories for:
  • dwi (fsl_dwi_proc),
  • anatomical processing (fsl_anat_proc),
  • the actual bip iteration results (bip)
  • If you use diannepat/dwiqc to get quality control metrics for your dwi data, you'll find it uses the same directory structure and they interoperate to avoid redundant processing or files.
• BIP2 now saves original copies of anatomical files and lesion masks in sourcedata not sourcedir, as per the BIDS specification.
• BIP2 checks for appropriate files at each step (setup, prep and bip) and exits with an error log if there is a problem.
• PE (phase encoding direction) is calculated for you if pe.txt does not exist.
• index.txt is created if no file like index*.txt is found.
• setup will generate a file, mask_check.txt, which can be run on a machine with your processed data and fsleyes so you can check your masks before continuing.
• BIP2 will set the permissions to allow writes to the bids anat dir, and, if necessary, the dwi and reverse PE images (if they have an odd number of slices they must be cropped to have an even number of slices and rewritten to the same directory.)
• Fixed bug that kept iterating when the volume of the result was 0.
• Added cst_l and cst_r (corticospinal tracts) back into options.

⁠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, python scripts, lists, and standard space image masks used to run bip. There is a Singularity recipe for building an HPC compliant version of BIP with GPU processing for BedpostX and Probtrack2X.

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 (on Bitbucket, see above) containing the Dockerfile:

⁠bip_wrap.sh

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

If bip_wrap.sh is in your path, the above call will show you help. By examining the contents of bip_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 bip_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 bip2.sif (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 run commands in bip_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://osf.io/fbh7q/⁠

⁠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: bip_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
|   |   `--sub-001_label-lesion_roi.nii.gz  (optional if you have a lesion mask)
|   |-- 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
|   `--sub-001_label-lesion_roi.nii.gz  (optional if you have a lesion mask)
`-- 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

The following will be computed if text files are not found: 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. In addition, the lesion mask will be added to the stop mask used by probtrackX (so, tracking will not be allowed through the lesion).

⁠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 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⁠

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, dwi mask, and if you have an fmap directory, magnitude mask. setup.sh uses a lesion mask if one is present and properly named (e.g., sub-001_label-lesion_roi.nii.gz). setup.sh generated mask_check.txt which you can run at the command line to evaluate all of your masks in fsleyes.

⁠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 (A DWI image containing more directions take longer to process). The shortest time I've seen is about 5 hours, unless you use the GPU processing implemented in Singularity. 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 derivatives/bip directory. bip runs 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 it takes some work. I'd suggest downloading the bitbucket directory and then working to generate and test masks.

⁠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 (optional, will be generated for you if not present) pe.txt (optional, will be generated for you if not present)

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.