https://doi.org/10.5281/zenodo.18728812

TBSS Tutorial

TBSS Tutorial#

TBSS stands for Tract-Based Spatial Statistics, a suite of tools for analyzing diffusion data. TBSS comes installed as part of the FSL software package. TBSS uses a tensor-fitting method to generate different measures of diffusion, such as fractional anisotropy (FA) and mean diffusivity (MD). Once these measurements are created, you can then extract them using ROI tools like you would for fMRI data.

The package also includes tools for correctying distortions in the diffusion data; in particular, the commands topup and eddy will remove distortions caused by eddy currents and magnetic field inhomogeneities.

Authors: Andrew Jahn & Monika Doerig

Date: 20 June 2025

License:

Note: If this notebook uses neuroimaging tools from Neurocontainers, those tools retain their original licenses. Please see Neurodesk citation guidelines for details.

Citation and Resources:#

Tools included in this workflow#

TBSS:

MRtrix3:

  • Tournier, J.-D., Smith, R., Raffelt, D., Tabbara, R., Dhollander, T., Pietsch, M., Christiaens, D., Jeurissen, B., Yeh, C.-H., & Connelly, A. (2019). MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 202, 116137. https://doi.org/10.1016/j.neuroimage.2019.116137

Educational resources#

Andy’s Brain Book:

Dataset#

Opensource Data from OpenNeuro:

! Note:

Before beginning, you may want to review the basics of diffusion imaging, which can be found here. The chapter will review the principles of diffusion, how they are used to generate diffusion-weighted images, and the advantages and disadvantages of fitting tensors to diffusion data.

Goals of This Tutorial#

This tutorial will show you how to analyze a sample dataset from start to finish using TBSS. We will be closely following the steps provided on the FSL TBSS website; the goal is to provide you with more experience in how to use this package, and use figures and videos to broaden your understanding of how it works and how it can be applied in other scenarios.

Import Python Packages#

%%capture
! pip install nibabel nilearn numpy

import os
import subprocess
import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
import tempfile
import nilearn
from nilearn import plotting, image
from IPython.display import display, clear_output, Image
import time
from ipyniivue import NiiVue

Load Neuroimaging Software (MRtrix and TBSS)#

import module
await module.load('mrtrix3/3.0.4')
await module.load('fsl/6.0.3') # Need to use this version in order for TBSS commands to work
await module.list()

['mrtrix3/3.0.4', 'fsl/6.0.3']

Try typing one of the commands from the library, such as fslinfo. If FSL has been installed correctly, you should see the help page printed by default when no arguments are passed to the command:

!fslinfo

Usage: /opt/fsl-6.0.3/bin/fslinfo <filename>

If this works without any errors, you are ready to begin analyzing the dataset.

Downloading the Dataset#

For this example, we will be preprocessing a dataset from openneuro.org called BTC preop. It includes data from patients with gliomas, patients with meningiomas, and a group of control subjects.

To download the data of one subject, excute the next cell:

PATTERN = "sub-CON02"
!datalad install https://github.com/OpenNeuroDatasets/ds001226.git

!cd ds001226 && datalad get $PATTERN

#

If you want to download more subjects, here is an example how to get 4 controls and 4 patients:

PATTERN_CONTROLS="sub-CON0[1-4]"
PATTERN_PATIENTS="sub-PAT0[1-4]"

! datalad install https://github.com/OpenNeuroDatasets/ds001226.git

! cd ds001226 && datalad get $PATTERN_CONTROLS $PATTERN_PATIENTS

Looking at the Data#

In order to view the data, we will navigate to the folder sub-CON02/ses-preop/dwi, which contains the diffusion data. We import the os package to change directories, and also to list the contents of that directory:

# Navigate to the diffusion data
os.chdir("ds001226/sub-CON02/ses-preop/dwi")

os.listdir()

['sub-02_AP.bvec',
 'AP.nii.gz',
 'AP_Cor.nii.gz',
 'dti_MD.nii.gz',
 'dti_L3.nii.gz',
 'AP_brain_mask.nii.gz',
 'sub-02_PA.bval',
 'sub-CON02_ses-preop_acq-AP_dwi.nii.gz',
 'dti_V3.nii.gz',
 'dti_L1.nii.gz',
 'tbss',
 'PA.nii.gz',
 'acq_param.txt',
 'AP_1stVol.nii.gz',
 'dti_L2.nii.gz',
 'sub-02_PA.bvec',
 'dti_S0.nii.gz',
 'AP_brain.nii.gz',
 'AP_PA.nii.gz',
 'sub-CON02_ses-preop_acq-PA_dwi.json',
 'sub-02_AP.bval',
 'dti_MO.nii.gz',
 'sub-CON02_ses-preop_acq-PA_dwi.nii.gz',
 'dti_V1.nii.gz',
 'dti_V2.nii.gz',
 'AP_PA.topup_log',
 'AP_PA_topup_movpar.txt',
 'dti_FA.nii.gz',
 'sub-CON02_ses-preop_acq-AP_dwi.json',
 'AP_PA_topup_fieldcoef.nii.gz']
! Note:

To make the rest of the example easier to read as well, use the mv command to rename the .bval and .bvec files: 

! mv sub-CON02_ses-preop_acq-AP_dwi.bvec sub-02_AP.bvec
! mv sub-CON02_ses-preop_acq-AP_dwi.bval sub-02_AP.bval
! mv sub-CON02_ses-preop_acq-PA_dwi.bvec sub-02_PA.bvec
! mv sub-CON02_ses-preop_acq-PA_dwi.bval sub-02_PA.bval

mv: cannot stat 'sub-CON02_ses-preop_acq-AP_dwi.bvec': No such file or directory
mv: cannot stat 'sub-CON02_ses-preop_acq-AP_dwi.bval': No such file or directory
mv: cannot stat 'sub-CON02_ses-preop_acq-PA_dwi.bvec': No such file or directory
mv: cannot stat 'sub-CON02_ses-preop_acq-PA_dwi.bval': No such file or directory

! fslinfo sub-CON02_ses-preop_acq-AP_dwi.nii.gz

data_type	INT16
dim1		96
dim2		96
dim3		60
dim4		102
datatype	4
pixdim1		2.500000
pixdim2		2.500000
pixdim3		2.500000
pixdim4		8.700000
cal_max		0.000000
cal_min		0.000000
file_type	NIFTI-1+

Note that, since this is a 4-dimensional dataset, the last dimension is time; in other words, this file contains 102 volumes, each one with dimensions of 96 x 96 x 60 voxels. The last dimension of the Voxel size field - which in this case has a value of 8.7 - indicates the time it took to acquire each volume, the repetition time, or TR.

Bvals and Bvecs#

The other files we need to check are the bvals and bvecs files. Briefly, the bvals contain a single number per volume that indicates how large of a diffusion gradient was applied to the data; and the bvecs file contains a triplet of numbers per volume that shows in what directions the gradients were applied. In general, volumes with larger b-values will be more sensitive to diffusion changes, but the images will also be more susceptible to motion and physiological artifacts.

The most important check is to ensure that the number of bvals and the number of bvecs are the same as the number of volumes in the dataset. For example, we can find the number of volumes in the sub-CON02_ses-preop_acq-AP_dwi.nii.gz dataset by typing:

! fslinfo sub-CON02_ses-preop_acq-AP_dwi.nii.gz | awk 'FNR == 5 {print $2}'

102

Which returns a value of 102, the number in the 4th field of the dimensions header that corresponds to the number of time-points, or volumes, in the dataset. We then compare this with the number of bvals and bvecs by using awk to count the number of columns in each text file:

! awk '{print NF; exit}' sub-02_AP.bvec
! awk '{print NF; exit}' sub-02_AP.bval

102
102
! Note:

If the number of volumes in your dataset and the number of bvals and bvecs do not match, you should check with your scan technician about the discrepancy; the files may not have been properly uploaded to the server, or maybe the diffusion-weighted image wasn’t acquired correctly.

Visualization#

In order to view the image, we will use a couple of different packages, and you can determine which one you prefer. The first one we will use it matplotlib, which was loaded earlier in this notebook. In the code below, we will temporarily convert the data to NIFTI format and plot 3 volumes with Matplotlib.

dwi_file = 'sub-CON02_ses-preop_acq-AP_dwi.nii.gz'  

# Load the NIfTI file
nii_image = nib.load(dwi_file)
data = nii_image.get_fdata()

# Define volume indices for the three volumes you want to display
volume_indices = [0, 2, 3]  
slice_index = data.shape[2] // 2  # Choose the middle slice for visualization

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for i, vol_idx in enumerate(volume_indices):
    # Extract the middle slice from the selected volume and rotate it 90°
    rotated_slice = np.rot90(data[:, :, slice_index, vol_idx]) 
    axes[i].imshow(rotated_slice, cmap="gray")
    axes[i].set_title(f"Volume {vol_idx} - Middle Slice")
    axes[i].axis("off")

plt.suptitle("Middle Slices of Selected Volumes")
plt.show()
../../_images/046fa8509ac780428641da4b246e3639b9e91a46400c0f174356fa7594e679b3.png

Another way to view these images is through nilearn, which contains a plotting library. The code below will import the necessary packages, and then extract and plot the first volume of the time-series. This opens an interactive plot upon which you can click and drag your mouse cursor. Note in particular that the frontal pole is slightly deformed and squished inwards; we will fix that by using topup.

dwi_image = 'sub-CON02_ses-preop_acq-AP_dwi.nii.gz'
print(image.load_img(dwi_image).shape)
first_dwi_volume = image.index_img(dwi_image,0) # Load first volume of DWI time-series
plotting.view_img(first_dwi_volume, bg_img=False, threshold=None)

(96, 96, 60, 102)

/tmp/ipykernel_3623/3924413228.py:4: UserWarning: Casting data from int32 to float32
  plotting.view_img(first_dwi_volume, bg_img=False, threshold=None)