FDT & TBSS Practical

Practical Overview

• Diffusion data
  Fitting diffusion tensors
  
• Bedpostx output
  Mapping Uncertainty
  
• Tractography
  Generating connectivity distributions from single voxels
  
• Connectivity-based segmentation
  Using connectivity to classify voxels, and running probabilistic tractography from seed masks
  
• Crossing-Fibre Tractography
  Generating connectivity distributions using multi-fibre directions
  
• Anatomical priors in tractography
  Using exclusion/inclusion/termination masks to inform tractography
  
• TBSS - Tract-Based Spatial Statistics
  Doing whole-brain voxelwise stats on DTI data

Diffusion Data and DTI

• cd ~/fsl_course_data/fdt/subj1

  List this directory and you should see a standard FDT directory.

• Bring up fslview and open data - you will need to reset the maximum display range to around 500. This is the diffusion data before correction for eddy currents. It is a 4D image. Turn on movie mode. Notice the movement and distortions.

• Turn off movie mode and have a look through the different volumes (different "timepoints" in the 4D data). See if you can work out which volumes had no diffusion gradient applied (e.g. "0" and "1", can you find more?), and, for some others, in which direction the diffusion gradients were applied. (You can check if you were right - e.g for diffusion encoding direction for volume 4 (the fifth volume), type
  cat bvecs | awk '{print $5}'
  Note - fslview starts counting from 0, and awk starts from 1 ! The awk command is extracting the 5th vector from the bvec file, which contains the normalized vectors giving the gradient directions in x, y and z.

• Extracting nodif and generating a brain mask.

  From the raw data you need to identify a volume without diffusion weighting (there may be more than one - you just need to find the first, for example). This will be used for several purposes, including being used for generating a brain mask and in registration to other images. If you look at bvals you can find the first zero entry - that signifies a volume with no diffusion weighting ("b=0"). For example,
  cat bvals
  tells you that the first timepoint is such a volume. You would also need to know this number ("0") to tell the head-motion/eddy-current correction which image to use as the reference volume.

  You can now extract this as a standalone 3D image, using
  fslroi data nodif 0 1
  generating new image nodif

  Now use the BET GUI (Bet or Bet_gui) to brain-extract this image, saving the brain as nodif_brain and turning on the binary mask output option. Have a look in FSLView - is the extraction good? It probably is a bit too small in places - this is common when using the default BET options - try rerunning with a smaller "fractional intensity threshold" until you're happy.

• Correcting for eddy current and head motion.

  You also need to run eddy_correct to align all volumes to each other. This attempts to correct for subject head motion and the effects of eddy currents. It uses the original images, not the brain extracted ones. We will skip this step here as it takes 5-15 minutes, and assume that our data is already reasonably aligned.

• DTIfit - fitting diffusion tensors to the data.

  Open the FDT GUI by typing Fdt or (Fdt_gui on a Mac), and select dtifit from the main menu.
  
  

• You can run DTIfit on a standard FDT directory by simply selecting the input directory (subj1) and pressing "Go".
  
  
• Alternatively, you can select the input files manually by clicking "Specify input files manually"
  • Then select the following files for the relevant tabs.
    • data: data
    • mask: nodif_brain_mask
    • gradient directions: bvecs
    • bvals: bvals
    • output : dti
  and press Go
  
  
• Note, we could also carry out dti fitting from the command line
  dtifit -k data -m nodif_brain_mask -r bvecs -b bvals -o dti
  
  
• Run dtifit using any of these three ways
  
  
• Load the anisotropy map into fslview
  fslview dti_FA &
  
  
• Note the areas with high diffusion anisotropy (white matter). Are there any areas where the FA is greater than 1? Why might this happen?
  answer
  
  
• Add the principal eigenvector map to your display
  file -> Add -> dti_V1
  • At first this looks slightly messy!
  • Select (i) -> DTI Display -> Lines
  • Zoom in and out of your data. You should see clear white matter pathways through the vector field.
  • Now select (i) -> DTI Display -> Display as: RGB. Diffusion direction is now coded by colour. To clarify the display, we can modulate the colour intensity with the FA map so anisotropic voxels appear bright. Select (i) -> Modulation -> dti_FA
  • You can work out which colours mean which directions by selecting (i) -> DTI Display -> Display as: Lines (RGB). The vectors are now coded by colour. You can turn back on the visibility of dti_FA so that the vectors are overlayed on the FA map.
  
  
• Now load the Mean Diffusivity (MD) map into fslview.
  fslview dti_MD &
  Observe how uniform the MD appears in both gray and white matter areas. Higher intensities in the CSF-filled areas indicate larger diffusivity -on average-, due to the absence of barriers to water molecule motions. Now load the principal DTI eigenvalue dti_L1. Larger values exist in white matter regions compared to gray matter. Next load the secondary eigenvalue dti_L2. We can now observe that L2 is lower in the white matter. Try also dti_L3. The contrast between white and gray matter is even larger, with much lower values being observed in the former. On the contrary, L1, L2 and L3 do not vary that much in the gray matter (or in the CSF). These differences between L1, L2 and L3 give rise to diffusion anisotropy and lead to high FA values in the white matter and low FA in the gray matter. Keep in mind the diffusion tensor ellipsoid for a pictorial representation of the DTI eigenvalues.
  
  
• What about dti_V2 and dti_V3? These vectors along with dti_V1 define a local coordinate system for each voxel, so that V1 points along the PDD (principal diffusion direction). L1, L2 and L3 are the diffusivities along directions V1, V2 and V3, respectively. Only V1, though, can be assosiated with the underlying fibre orientations, as shown by the colour maps before. V2 and V3 are just perpendicular to this orientation.
  

Bedpostx output - single fibre orientation

cd ../subj1.bedpostX

This is what a bedpostx directory looks like if you use only a single fibre orientation per voxel (i.e. bedpostx -n 1 ). You will see a two-fibre bedpostx directory later.

• As it takes several hours to run, bedpostx has already been run for you on the data in subj1. This contains all the files essential for tractography.

• For now, ignore all the files beginning with MASK. They will be used later.

• Open mean_f1samples in fslview (you may need to adjust the maximum display intensity to 1). This is the mean of the distribution on the anisotropy parameter in the partial volume model - it should look very similar to the FA map you saw earlier.

• Add dyads1 and select (i) -> DTI Display -> Lines . This is the mean of the distribution of fibre orientations - it should look very similar to the dti_V1 you saw earlier. As before with V1, you can color code dyads1.

• Now add dyads1_dispersion. You need to adjust the Min and Max display intensity. Set Min to 0.005 and Max to 0.03. Also change the colour of the image, selecting (i) -> Lookup Table Options -> Red-Yellow. This map encodes for each voxel the uncertainty (i.e. spread of the distribution) around the mean fibre orientation dyads1. As expected, highly anisotropic voxels in white matter have very low dispersion index, while gray matter and CSF voxels have orders of magnitude higher dispersion. Notice that some voxels in white matter have also high dispersion. These are voxels where a model with a single fibre compartment is not good, indicating a complex fibre architecture (i.e. crossing fibres). An indicative example is the centrum semiovale. Set Y=48 and look at the coronal view. Can you locate it?

• Apart from the mean of the fibre orientation distribution, bedpostx returns samples from this distribution that are utilized in probabilistic tractography. These are contained in merged_th1samples and merged_ph1samples. Fibre orientations are described in spherical coordinates* and these files contain the angle values for each orientation sample. Viewing the merged_ files with fslview is not particularly informative, normally dyads are used for visualization.

  ^* Spherical polar coordinates

FDT - Tractography Practical - single fibre per voxel

• Simple tracking - connectivity from a single seed point.
  • In fslview, open the FA map that you created using DTIFIT. Find the co-ordinates of a voxel in the internal capsule.
    
  • Bring up Fdt - it should default to "ProbTrackx" with "single voxel" as the selected seed mode.
  • Select subj1.bedpostX as the Bedpostx directory and enter your selected voxel as the seed.
  • Choose an output name, and press "Go".
  • Overlay the connectivity distribution on top of the FA map - note the pathway's name ends with the coordinates of the seed voxel (e.g. basename_X_Y_Z.nii.gz). Adjust the min and max display thresholds to look at the relative probabilities of paths. Change the color look-up-table to Red-Yellow, by (i) -> Lookup Table Options -> Red-Yellow. Note the narrowness of the distribution - there is very low uncertainty in the Internal Capsule.
  
  
  
• Tracking from a voxel in a non-diffusion space.
  • Open subj1.bedpostX/struct_brain into fslview and find an interesting seed location.
  • in Fdt, remaining in "single voxel" seed mode, select "seed space is not diffusion" - This should bring up two new dialogues, specifying the Seed space.
  • The "Seed space reference image" is the image you selected the seed point from (struct_brain).
  • The seed to diff space registration transform can be found in the xfms directory (again, the registration transforms have already been computed for you): In this case str2diff.mat
  • Run this tractography and load the output into FSLView on top of the structural image.
  
  
  
• Try tracking from voxels in MNI152 standard space - you can also specify the seed points in standard space mm (e.g. FMRI activation blobs output from FEAT). Either linear or nonlinear transformations to standard space can be utilized. In the former case (linear), use standard2diff.mat as a transform and standard as a reference image. In the latter (nonlinear), use standard2diff_warp and diff2standard_warp as transforms and standard as a reference image. Try the nonlinear option.
  
  
• Don't try too many tracts now though, because you will see a big difference when it comes to using two fibres per voxel.

FDT - Segmentation Practical

• A connectivity-based seed classification analysis has been run for you and is saved in a sub-directory called THAL2CTX_right within your subj1.bedpostX directory.
  
  
• This analysis was carried out in standard space. cd back into subj1.bedpostX. Open the standard brain in fslview and overlay the seed mask, MASK_average_thal_right - you should see that the seed is thalamus.
  
  
• In the subj1.bedpostX directory, there is a file called THAL2CTX_right/targets.txt. If you cat this file, you will see a list of target masks used in this analysis. Overlay some of these masks in your open fslview window - they are different cortical zones (you may need to set the intensity display range to say 0:1 to view these).
  
  
• The objective of this type of analysis is to ask the following question: "For each location in the seed mask, what is the relative probability of connection to each of my target masks". The output is therefore a single image for each target mask in which only voxels within the seed mask contain data. At each voxel within the seed mask the voxel value is the number of samples which reached this target mask from this seed voxel. These outputs can be found in THAL2CTX_right in the subj1.bedpostX directory.
  
  
• Now overlay in FSLView, in turn, each of the seeds_to... files in THAL2CTX_right. Compare the spatial distributions and probabilities of connection from different thalamic voxels to different cortical target zones.
  
  
• Run find_the_biggest on the outputs of seeds to targets to generate a hard segmentation of the thalamus and overlay this segmentation onto the standard brain. Something like:
  find_the_biggest seeds_to_* biggest_segmentation
  The different numbers in the output image (biggest_segmentation) each relate to a different target mask. You may want to change the colour look-up-table to something like Random-Rainbow.

Crossing fibre tractography

cd ../subj1_2fibres.bedpostX

This is a two-fibre bedpostX analysis output directory.

• Just like before, mean_f1samples tells you the estimated proportion of the diffusion signal that can be accounted for by the first fibre orientation. Now though, mean_f2samples tells you the same for the second orientation. dyads1 and dyads2 are the mean estimates for the orientation vectors for these. Load these four images into fslview.
  
  
• First turn off dyads1 and dyads2 and change mean_f2samples to a different colour ( (i) -> look-up table). Threshold it at 0.05.
  
  
• You should see that in much of the white matter, there are two fibre populations, but that this does not depend strongly on the threshold (e.g. try 0.01) - this is because ARD (automatic relevance determination, part of the Bayesian analysis of the crossing fibres) has removed voxels where there is little evidence for a second population.
  
  
• Now turn back on dyads1 and dyads2 and change them both to lines mode ( (i)->DTI display options, display as ->lines). The first fibre will be in red, and the second in blue. You will note that there are blue lines everywhere (even where there is no evidence for a second fibre in mean_f2samples), but you should notice the blue lines are well ordered in the regions where 2 fibres are supported (where mean_f2samples survives the threshold), but random elsewhere.

• Back in your shell, now type: maskdyads dyads2 mean_f2samples
  This makes a new image of dyads2, but only in regions where 2 fibres are supported. Load this image up into your fslview. It will be called dyads2_thr0.05 If you display only mean_f1samples, dyads1 and dyads2_thr0.05 and turn everything else off, you should be able to see some beautiful crossing fibre architecture!

• Now try to visualise these crossing fibres in standard space. For that you will have to transform the dyads, applying the registration transformation in xfms/diff2standard_warp, using the program vecreg. Type:
  vecreg -i dyads1 -o dyads1_to_standard -r standard -w xfms/diff2standard_warp
  Do the equivalent for dyads2_thr0.05. Open fslview with a standard brain and the two registered dyads.
  
  
• Try some of the same tractography you have already done in the previous section but, this time, specify subj1_2fibres.bedpostX as your subject directory. You should find that you can see many pathways that were previously invisible. For example, see if you can find the medial portion of the SLF in the 2-fibre data (see example in lecture).
  
  
• Try tracking from cortical grey matter - this is likely to produce variable results due to high uncertainty on fibre direction in the cortex (try moving into the white matter adjacent to the grey matter area of interest).

Anatomical priors in tractography

Here we will make use of target masks to inform tractography. Provided a seed mask in the lateral geniculate nucleus of the thalamus (LGN), we will construct the pathway connecting LGN to the ipsilateral primary visual cortex (V1) that passes via Meyer's loop.

• Run tractography from the LGN mask (file called SEED_LGN_LEFT.nii.gz). You will need to turn the seed mode from "single voxel" to "single mask". Note that this mask is in standard space (i.e. use the transformation as in the example above). You may want to reduce the number of samples to 500 to speed up tractography (under Options tab).
  
  
• Use a waypoint mask in V1 (TARGET_V1_LEFT.nii.gz) to isolate only those tracts that reach V1.
  
  
• Use V1 also as a termination mask to avoid tracts that reach V1 and then continue to other parts of the brain.
  
  
• There are still some artefactual tracts. Draw an exclusion mask in fslview and include it in the tractography. (or use the file EXCLUSION.nii.gz).

TBSS - Tract-Based Spatial Statistics

We will now run a TBSS analysis of a small dataset - 3 young people (mean age 25) and 3 older people (mean age 65). We will attempt to find where on the tract skeleton the two groups of subjects have significantly different FA.

0 - creating FA data from a diffusion study

We have already created the FA images from the 6 subjects, using dtifit.

cd ~/fsl_course_data/tbss

Remember that all of the main scripts (until you reach the randomise stage near the end) need to be run from within this directory.

LOOK AT YOUR DATA:

• Open one image with fslview to get a feel for the raw data quality and resolution. You may need to adjust the intensity display range. Look at the image histogram.
• Run slicesdir *.nii.gz and open the resulting web page report; this is a quick-and-dirty way of checking through all the original images.

1 - preparing your FA data for TBSS

tbss_1_preproc *nii.gz

2 - registering all the FA data

cp precomputed_registration/* FA

This puts the necessary registration output files into FA, as if you had run the registration yourself.

3 - post-registration processing

Before reading this section, start the next script running so that you can read what's going on while it runs (it will take about 5 minutes):

tbss_3_postreg -S

The previous script only got as far as registering all subjects to the chosen template. The tbss_3_postreg script applies these registrations to take all subjects into 1x1x1mm standard space.

Next, all subjects' standard space nonlinearly aligned images are merged into a single 4D image file called all_FA, created in a new subdirectory called stats. The mean of all FA images is created, called mean_FA, and this is then fed into the FA skeletonisation program to create mean_FA_skeleton.

Once the script has finished running, check that the mean FA image looks reasonable, and is well aligned with the MNI152 image:

cd stats
fslview $FSLDIR/data/standard/MNI152_T1_1mm mean_FA -l Red-Yellow -b .2,.6

As you move around in the image you should see that the mean FA image is indeed well aligned to standard space and corresponds to white matter in the MNI152 image. Now kill FSLView and restart it, loading the 4D aligned FA data underneath the skeleton:

fslview all_FA mean_FA_skeleton -l Green -b .2,.6

Now turn on the movie loop; you will see the mean FA skeleton on top of each different subject's aligned FA image. If all the processing so far has worked, the skeleton should look like the examples shown in the lecture. If the registration has worked well you should see that in general each subject's major tracts are reasonably well aligned to the relevant parts of the skeleton. If you set the skeleton threshold (in FSLView, the lower of the display range settings) much lower than 0.2, it will extend away towards extremes where there is too much cross-subject variability and where the nonlinear registration has not been able to give good alignments. In this dataset, with very few subjects, the mean FA image is quite noisy, so you probably want the threshold around 0.3.

4 - projecting all pre-aligned FA data onto the skeleton

The last TBSS script carries out the final steps necessary before you run the voxelwise cross-subject stats. It thresholds the mean FA skeleton image at the chosen threshold:

cd .. (if you're still in the stats directory)
tbss_4_prestats 0.3

This takes 4-5 minutes to run; read the rest of this section while it's running (though if you can't make much sense of it, don't worry - it's decribed in more detail in the paper!).

The thresholding creates a binary skeleton mask that defines the set of voxels used in all subsequent processing.

Next a "distance map" is created from the skeleton mask. This is used in the projection of each subject's FA onto the skeleton; when searching outwards from a skeleton voxel for the local tract centre, the search only continues while the distance map values keep increasing - this means that the search knows to stop when it has got more than halfway between the starting skeleton point and another separate part of the skeleton.

Finally, the script takes the 4D pre-aligned FA images in all_FA and, for each "timepoint" (subject ID), projects the FA data onto the mean FA skeleton. This results in a 4D image file containing the (projected) skeletonised FA data. It is this file that you will feed into voxelwise statistics in the next section.

Once the script has finished, cd into stats and have a look at all_FA_skeletonised in fslview - turn on movie mode to see the different timepoints of the skeletonised data.

5 - voxelwise statistics on the skeletonised FA data

The previous step resulted in the 4D skeletonised FA image. It is this that you now feed into voxelwise statistics, that, for example, tells you which FA skeleton voxels are significantly different between two groups of subjects.

The recommended way of doing the stats is to use the randomise tool. For more detail see the randomise manual. Before running randomise you will need to generate design matrix and contrast files (e.g., design.mat and design.con). We will use the Glm GUI to generate these. Note that the order of the entries (rows) in your design matrix must match the alphabetical order of your original FA images, as that determines the order of the aligned FA images in the final 4D file all_FA_skeletonised. In this case the order was: 3 young subjects and then 3 older subjects. Start the GUI:

cd stats
Glm

Change Timeseries design to Higher-level / non-timeseries design. Change the number of inputs to 6 (you may have to press Return after typing in 6) and then use the wizard to setup the two-group unpaired t-test. Reduce the number of contrasts to 2 (we're not interested in the group means on their own). Finally, save the design as filename design, and in the terminal use more to look at the design.mat and design.con files.

You are now ready to run the stats. Because this reduced dataset only contains 6 subjects, only 20 distinct permutations are possible, so by default randomise will run just these 20. Again, because this tiny dataset has so few subjects the raw t-stat values will not be very significant - so let's try setting a cluster-forming t threshold of 1.5 (in "real" datasets we would normally recommend using the --T2 option for TFCE instead of cluster-based thresholding):

randomise -i all_FA_skeletonised -o tbss -m mean_FA_skeleton_mask -d design.mat -t design.con -c 1.5 -V

Contrast 1 gives the young>older test. The raw unthresholded tstat image is tbss_tstat1 and the corresponding (p-values corrected for multiple comparisons) cluster image is tbss_clustere_corrp_tstat1.

Thresholding clusters at 0.949 (corresponding to thresholding the p-values at 0.05, because randomise outputs p-values as 1-p for convenience of display - so that higher values are more significant) shows a few areas of reduced FA in the older subjects. The following shows unthresholded t-stats in red-yellow and thresholded clusters in blue:

fslview $FSLDIR/data/standard/MNI152_T1_1mm mean_FA_skeleton -l Green -b .3,.7 tbss_tstat1 \
   -l Red-Yellow -b 1.5,3 tbss_clustere_corrp_tstat1 -l Blue-Lightblue -b 0.949,1

Look at the histogram of tstat1; it is clearly shifted to the right, suggesting a global decrease in FA with aging.

If you prefer to display "fattened" results (that are not just the one-voxel-thick skeleton, but "fill out" into the tracts as seen in mean FA image), whilst still masking with the mean FA image, there is a simple script which will do this for you, then run:

tbss_fill tbss_clustere_corrp_tstat1 0.949 mean_FA tbss_clustere_corrp_tstat1_filled

and then load this output into FSLView and choose an appropriate colourmap. This script smooths the thresholded input image (by a typical tract width according to the local tract structure), multiples by the mean FA image (in order to constrain the width) and then adds the original input back in so that the skeleton can still be seen within the thickened output.