Quality of our multi-tissue multi-shell pipeline?

1

Dear developers, dear users,

we have a acquired multi-shell dataset of 31 Subjects (b=1000,1750, 2500) 1.75mm isotropic multi-band 3 (TR= 3800 ms, TE= 104.8 ms, flip angle= 90 deg., Refocus flip angle= 180 deg., FoV= 224 mm, 72 slices, a-p phase encoding) with 159 directions (6 b=0 images interspersed) including an additional reverse phase encoding session in a clinical 3T scanner (Siemens Skyra, 64-ch head/neck coil) under optimised head restriction conditions (Pearltec crania audult) including a HCP-like T1 image (0.75mm voxel resolution isotropic).

We then employed the following preprocessing pipeline before probabilistic tractography and would kindly ask you to comment on it (strongly ask you to criticise it):

#!/bin/bash

# MRtrix 3.0_RC3 pipeline
# PD Dr. med. M. Duering & Prof. Dr. med. P. zu Eulenburg
# June 25th 2018



## Preprocessing

# Denoising
dwidenoise ${id}_ses-ap_dwi.nii.gz temp_denoised_ap.nii.gz 2>&1
dwidenoise ${id}_ses-pa_dwi.nii.gz temp_denoised_pa.nii.gz 2>&1
mrdegibbs temp_denoised_ap.nii.gz temp_diffusion_ap.nii.gz 2>&1
mrdegibbs temp_denoised_pa.nii.gz temp_diffusion_pa.nii.gz 2>&1

# TOPUP and EDDY
# Extract b0s
select_dwi_vols temp_diffusion_ap.nii.gz ${id}_dwi.bval temp_b0s_ap.nii.gz 0
select_dwi_vols temp_diffusion_pa.nii.gz ${id}_dwi.bval temp_b0s_pa.nii.gz 0
fslmerge -t temp_allb0.nii.gz temp_b0s_ap.nii.gz temp_b0s_pa.nii.gz
# Run TOPUP
topup --imain=temp_allb0.nii.gz --datain=${param} --config=b02b0.cnf --out=topup --iout=temp_b0_hifi.nii.gz
fslmaths temp_b0_hifi.nii.gz -Tmean temp_b0_hifi_avg.nii.gz
bet temp_b0_hifi_avg.nii.gz temp_b0_hifi_avg_brain -f 0.3 -R -m
fslmaths temp_b0_hifi_avg_brain_mask -fillh temp_b0_hifi_avg_brain_mask
# Run EDDY (5.0.11)
fslmerge -t temp_diffusion.nii.gz temp_diffusion_ap.nii.gz temp_diffusion_pa.nii.gz 
eddy_cpu -v --imain=temp_diffusion.nii.gz --mask=temp_b0_hifi_avg_brain_mask --acqp=${param} --index=${index} --bvecs=temp_diffusion.bvec --bvals=temp_diffusion.bval --topup=topup --out=temp_eddy --resamp=lsr --fep --data_is_shelled

# Bias correction
N4BiasFieldCorrection -d 3 -i temp_b0_corrected_avg.nii.gz -w temp_b0_hifi_avg_brain_mask -o [temp_b0_hifi_avg_nobias.nii.gz,${id}_dwi-bias.nii.gz] -b [150,3] -c [1000x1000,0.0]
fslmaths temp_eddy.nii.gz -div ${id}_dwi-bias.nii.gz ${id}_dwi-data.nii.gz


## MRtrix pipeline

# Convert to mif format
mrconvert ${id}_dwi-data.nii.gz dwi.mif -fslgrad ${id}_dwi.bvec ${id}_dwi.bval -datatype float32 -stride 0,0,0,1  2>&1

# Mask
dwiextract dwi.mif - -bzero | mrmath - mean meanb0.nii.gz -axis 3  2>&1
maskfilter ${id}_dwi-mask.nii.gz erode mask.mif -npass 1  2>&1

# Register T1 to DWI using ANTS
antsRegistrationSyN.sh -d 3 -f meanb0.nii.gz -m ${id}_T1w.nii.gz -o t1_to_dwi -t r -p f
antsApplyTransforms -d 3 -e 3 -i ${id}_T1w.nii.gz -r ${id}_T1w.nii.gz -t t1_to_dwi0GenericAffine.mat -o t1_to_dwi.nii.gz --float

# 5TT image
5ttgen fsl t1_to_dwi.nii.gz 5tt.mif  2>&1
5tt2vis 5tt.mif 5tt_vis.mif 2>&1

# Response function
dwi2response dhollander dwi.mif wm_response.txt gm_response.txt csf_response.txt  2>&1

# FOD
dwi2fod msmt_csd -mask mask.mif dwi.mif wm_response.txt wmfod.mif gm_response.txt gm.mif csf_response.txt csf.mif  2>&1
mrconvert -coord 3 0 wmfod.mif - | mrcat csf.mif gm.mif - vf.mif  2>&1

# Multi-tissue informed log-domain intensity normalisation
mtnormalise wmfod.mif wmfod_norm.mif gm.mif gm_norm.mif csf.mif csf_norm.mif -mask mask.mif  2>&1

# Tractography
tckgen wmfod_norm.mif 10M.tck -act 5tt.mif -backtrack -crop_at_gmwmi -seed_dynamic wmfod_norm.mif -maxlength 250 -select 10M -cutoff 0.06

# SIFT2
tcksift2 10M.tck wmfod_norm.mif sift2_weights.txt -act 5TT.mif

# Connectome; nodes_in in .mif format
tck2connectome 10M.tck sift2_weights.txt [nodes_in] connectome.csv

thanks in advance for any feedback
Peter zu Eulenburg

2

Dear Peter,

Your pipeline looks very robust to me, I would modify a couple of things.

Denoising
dwidenoise {id}_ses-ap_dwi.nii.gz temp_denoised_ap.nii.gz 2>&1 dwidenoise {id}_ses-pa_dwi.nii.gz temp_denoised_pa.nii.gz 2>&1
mrdegibbs temp_denoised_ap.nii.gz temp_diffusion_ap.nii.gz 2>&1
mrdegibbs temp_denoised_pa.nii.gz temp_diffusion_pa.nii.gz 2>&1

As your data contains 159 volumes, maybe you want to add the flag -extent 7,7,7 to obtain better results.

Run EDDY (5.0.11)
fslmerge -t temp_diffusion.nii.gz temp_diffusion_ap.nii.gz temp_diffusion_pa.nii.gz
eddy_cpu -v --imain=temp_diffusion.nii.gz --mask=temp_b0_hifi_avg_brain_mask --acqp={param} --index={index} --bvecs=temp_diffusion.bvec --bvals=temp_diffusion.bval --topup=topup --out=temp_eddy --resamp=lsr --fep --data_is_shelled

The option --repol in eddy is also very interesting, you can take a look here.

Bias correction
N4BiasFieldCorrection -d 3 -i temp_b0_corrected_avg.nii.gz -w temp_b0_hifi_avg_brain_mask -o [temp_b0_hifi_avg_nobias.nii.gz,{id}_dwi-bias.nii.gz] -b [150,3] -c [1000x1000,0.0] fslmaths temp_eddy.nii.gz -div {id}_dwi-bias.nii.gz ${id}_dwi-data.nii.gz

Why don’t use the dwibiascorrect script? now provides updated parameters and also doesn’t introduce big scale differences between before and after bias correction.

Register T1 to DWI using ANTS
antsRegistrationSyN.sh -d 3 -f meanb0.nii.gz -m {id}_T1w.nii.gz -o t1_to_dwi -t r -p f antsApplyTransforms -d 3 -e 3 -i {id}_T1w.nii.gz -r ${id}_T1w.nii.gz -t t1_to_dwi0GenericAffine.mat -o t1_to_dwi.nii.gz --float

As the data has been previously corrected for all distortions, a rigid body transformation should be enough. In this scenario the bbr registration has shown to provide better results that a simply rigid body. In the forum this issue has been discussed a couple of times, here.

I hope all of this helps.

Best regards,

Manuel

3

Thanks a lot for the more than detailed feedback.

What exactly does the -extent 7,7,7 flag do?

Best regards
Peter

4

Hi,

In the manual it says:

The kernel size, default 5x5x5, can be chosen by the user (option: -extent). For maximal SNR gain we suggest to choose N>M for which M is typically the number of DW images in the data (single or multi-shell), where N is the number of kernel elements. However, in case of spatially varying noise, it might be beneficial to select smaller sliding kernels, e.g. N~M, to balance between precision, accuracy, and resolution of the noise map.

With 159 directions, the default kernel size is too small (125). Here there is a little bit of information about this, and how it works.

Best regards,

Manuel

5

Thanks again for the additional information. People told me this forum would be really helpful and fast.

Best wishes
Peter