Adding RESTORE to main.nf

, ,
1

I have been tinkering with main.nf, so far successfully, to add steps such as Gibbs correction to the pipeline. Now I am trying to implement RESTORE robust tensor estimation. On paper, this seems to be an easy task, as the tool used to estimate the tensor (scil_compute_dti_metrics.py) has an option to run RESTORE built into it.

I tried to implement this by adding --method restore\ to the DTI_Metrics process, but I am thrown an error. Tracing back to the log file for the command that invokes scil_compute_dti_metrics, I see the following output message:

image
image2866×998 455 KB

Please let me know if there is a way to address this error. Also, should I be concerned about the b-vector normalization warning towards the beginning?

I am running Tractoflow 2.1.1 with the appropriate 2.1.1 singularity container on HPC.

Edit: The data have up to 30 gradients at b=1000 (depending on how much is thrown out during QC) and 4 b0 scans. The data have been denoised, de-Gibbed, and eddy/motion corrected. I have also tagged DIPY now, since SCILPY calls the RESTORE tensor fitting from DIPY.

Thanks,
Steven

2

Hi @Steven

What do you modified exactly ? Can you send me your code ?

Best

Guillaume

3

From main.nf

process DTI_Metrics {
cpus 3

input:
set sid, file(dwi), file(bval), file(bvec), file(b0_mask)\
    from dwi_and_grad_for_dti_metrics

output:
file "${sid}__ad.nii.gz"
file "${sid}__evecs.nii.gz"
file "${sid}__evecs_v1.nii.gz"
file "${sid}__evecs_v2.nii.gz"
file "${sid}__evecs_v3.nii.gz"
file "${sid}__evals.nii.gz"
file "${sid}__evals_e1.nii.gz"
file "${sid}__evals_e2.nii.gz"
file "${sid}__evals_e3.nii.gz"
file "${sid}__fa.nii.gz"
file "${sid}__ga.nii.gz"
file "${sid}__rgb.nii.gz"
file "${sid}__md.nii.gz"
file "${sid}__mode.nii.gz"
file "${sid}__norm.nii.gz"
file "${sid}__rd.nii.gz"
file "${sid}__tensor.nii.gz"
file "${sid}__nonphysical.nii.gz"
file "${sid}__pulsation_std_dwi.nii.gz"
file "${sid}__residual.nii.gz"
file "${sid}__residual_iqr_residuals.npy"
file "${sid}__residual_mean_residuals.npy"
file "${sid}__residual_q1_residuals.npy"
file "${sid}__residual_q3_residuals.npy"
file "${sid}__residual_residuals_stats.png"
file "${sid}__residual_std_residuals.npy"
set sid, "${sid}__fa.nii.gz" into\
    fa_for_reg

script:
"""
export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
export OMP_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
scil_compute_dti_metrics.py $dwi $bval $bvec --mask $b0_mask\
    --ad ${sid}__ad.nii.gz --evecs ${sid}__evecs.nii.gz\
    --evals ${sid}__evals.nii.gz --fa ${sid}__fa.nii.gz\
    --ga ${sid}__ga.nii.gz --rgb ${sid}__rgb.nii.gz\
    --md ${sid}__md.nii.gz --mode ${sid}__mode.nii.gz\
    --norm ${sid}__norm.nii.gz --rd ${sid}__rd.nii.gz\
    --tensor ${sid}__tensor.nii.gz\
    --non-physical ${sid}__nonphysical.nii.gz\
    --pulsation ${sid}__pulsation.nii.gz\
    --residual ${sid}__residual.nii.gz\
--method restore\
    -f
"""

}

Originally, the argument line with --method is not included. I have tried other DTI tensor models, such as --method NLLS, and it has run fine.

4

Hi! Yes: this looks like an error that is happening in DIPY. It’s something that does come up every once in a while… For example, see: https://github.com/dipy/dipy/issues/2062. This could also arise when the sigma noise estimate used in RESTORE is too small. @Guillaumeth: how is sigma estimated in scilpy? [EDITED: “tractoflow”=> “scilpy”]

5

The relevant lines of code are:

import dipy.denoise.noise_estimate as ne
import nibabel as nib
from dipy.reconst.dti import (TensorModel, color_fa, fractional_anisotropy,
                              geodesic_anisotropy, mean_diffusivity,
                              axial_diffusivity, norm,
                              radial_diffusivity, lower_triangular)
def _get_min_nonzero_signal(data):
    return np.min(data[data > 0])


img = nib.load(PREPROCESSED_DWI)
data = img.get_fdata(dtype=np.float32)
sigma = ne.estimate_sigma(data)
tenmodel = TensorModel(gtab, fit_method='restore', sigma=sigma, min_signal=_get_min_nonzero_signal(data))
6

Dear @Steven,

Don’t hesitate to suggest PR for tractoflow if a specific parameter is not coded.

If you feel this thread is solved can you “tag” it as solved ? It helps sort the different issues.

Thank you in advance
Arnaud