Problem with QSIPREP output space and calculating scalars (FA, MD, etc...) with Freesurfer and MRtrix

Software Support
1

Hi there!

I’ve been using QSIPREP to pre-process diffusion MRIs. Here’s the code I used:

singularity run --cleanenv \
    -B ${HOME}/EcriPark_Code:/code,${EcriPark}:/data,/scratch/lcorcos/EcriPark_QSIPREP/TESTV9:/out,/home/lcorcos/freesurfer/license.txt,/scratch/lcorcos/Temp_QSIPREP/TMP:/tmp \
    --nv ${HOME}/qsiprep-0.18.1.sif /data/ \
    /out/ participant --participant_label ${SUB} \
    --skip_bids_validation \
    --nthreads 24 \
    --omp_nthreads 12 \
    --fs_license_file /home/lcorcos/freesurfer/license.txt \
    --work_dir /tmp/ \
    --output_resolution 1 \
    --eddy_config /code/eddy_params.json \
    --anatomical-template MNI152NLin2009cAsym \
    --verbose \
    --anat_modality T1w \
    --b0_threshold 50 \
    --dwi_denoise_window 5 \
    --denoise_method dwidenoise \
    --unringing_method mrdegibbs \
    --distortion_group_merge average \
    --bo_motion_corr_to iterative \
    --hmc_transform Affine \
    --hmc_model eddy \
    --skull_strip_template OASIS \
    --skull_strip_fixed_seed \
    --write_graph

Version:

qsiprep 0.18.1
freesurfer-linux-centos8_x86_64-7.4.1-20230613-7eb8460
singularity 1.1.6-1.el7

Hardware:

Dell PowerEdge C4140 (32 cores) Intel Xeon CPU 5218
380 Go RAM
NVIDIA Tesla V100

The main objective is to use this data to calculate scalar maps (FA, MD, RD, etc…) and to do inter-group comparison. To do this, I imagined putting the data in the same space (MNI) and applying an atlas to it.

The problem is that the data processed by QSIPREP doesn’t seem to be in the MNI space. (I downloaded the ICBM 152 Nonlinear atlases version 2009 here: Templates and atlases — The Princeton Handbook for Reproducible Neuroimaging) and it didn’t fit with the data.
Moreover, each of my subjects didn’t fit between them, and I have the impression that there was no recalibration… (see photos)

Also, here’s the QSIPREP output structure:

├── anat
│   ├── sub-CTR01_desc-aseg_dseg.nii.gz
│   ├── sub-CTR01_desc-brain_mask.nii.gz
│   ├── sub-CTR01_desc-preproc_T1w.nii.gz
│   ├── sub-CTR01_dseg.nii.gz
│   ├── sub-CTR01_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5
│   ├── sub-CTR01_from-T1wACPC_to-T1wNative_mode-image_xfm.mat
│   ├── sub-CTR01_from-T1wNative_to-T1wACPC_mode-image_xfm.mat
│   └── sub-CTR01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5
├── figures
│   ├── sub-CTR01_seg_brainmask.svg
│   ├── sub-CTR01_ses-01_carpetplot.svg
│   ├── sub-CTR01_ses-01_coreg.svg
│   ├── sub-CTR01_ses-01_desc-initial_b0ref.svg
│   ├── sub-CTR01_ses-01_desc-sdc_b0.svg
│   ├── sub-CTR01_ses-01_dir-AP_dwi_denoise_ses_01_dir_AP_dwi_wf_denoising.svg
│   ├── sub-CTR01_ses-01_dir-AP_dwi_denoise_ses_01_dir_AP_dwi_wf_unringing.svg
│   ├── sub-CTR01_ses-01_dir-PA_dwi_denoise_ses_01_dir_PA_dwi_wf_denoising.svg
│   ├── sub-CTR01_ses-01_dir-PA_dwi_denoise_ses_01_dir_PA_dwi_wf_unringing.svg
│   ├── sub-CTR01_ses-01_sampling_scheme.gif
│   └── sub-CTR01_t1_2_mni.svg
└── ses-01
    ├── anat
    │   └── sub-CTR01_ses-01_from-orig_to-T1w_mode-image_xfm.txt
    └── dwi
        ├── sub-CTR01_ses-01_confounds.tsv
        ├── sub-CTR01_ses-01_desc-ImageQC_dwi.csv
        ├── sub-CTR01_ses-01_dwiqc.json
        ├── sub-CTR01_ses-01_space-T1w_desc-brain_mask.nii.gz
        ├── sub-CTR01_ses-01_space-T1w_desc-eddy_cnr.nii.gz
        ├── sub-CTR01_ses-01_space-T1w_desc-preproc_dwi.b
        ├── sub-CTR01_ses-01_space-T1w_desc-preproc_dwi.bval
        ├── sub-CTR01_ses-01_space-T1w_desc-preproc_dwi.bvec
        ├── sub-CTR01_ses-01_space-T1w_desc-preproc_dwi.nii.gz
        ├── sub-CTR01_ses-01_space-T1w_desc-preproc_dwi.txt
        └── sub-CTR01_ses-01_space-T1w_dwiref.nii.gz

It differs from what’s marked in the doc (Preprocessing — qsiprep version documentation)

Is this due to out-of-date documentation or have I set the command options incorrectly?

Next, I use MRTRIX’s dwi2tensor and tensor2metric commands to obtain scalar maps. Is this the right choice?

Finally, do I need to run Freesurfer’s recon-all command to obtain a parcelation of the GM and WM? Or should I do a registration of an existing atlas?

Many thanks for your help.

MNI_and_sub1
MNI_and_sub11544×1005 90 KB

sub1_and_sub2
sub1_and_sub21544×1005 76.6 KB

2

Hi @MagicLudo,

I would calculate all scalars in native space. Thankfully that is what QSIPrep does. You might be interested in the multishell_scalarfest recon spec.

All DWI derivatives are in native space.

That is one of many correct choices (FSL, DIPY and other softwares have their own tensor fitting algorithm, and probably perform comparably well).

You can register an atlas to native space. Thankfully QSIPrep has also output the MNI-to-T1 transform, which you can invert to bring atlases into native space. And if you want to parcellate WM, you’d probably be better off using a tractography algorithm (e.g., pyafq_tractometry, mrtrix_multishell_msmt_pyafq_tractometry or dsi_studio_autotrack).

Best,
Steven

3

The results are in in AC-PC alignment, so not aligned to the raw BIDS data. The space-T1w results are in alignment with the acpc-aligned T1w.

4

I would calculate all scalars in native space. Thankfully that is what QSIPrep does. You might be interested in the multishell_scalarfest recon spec .

Thanks for the tip! I’ve been running the command for dozens of hours already and it’s generating scalar files for me.
However, do you know where I can find some documentation on what it does exactly? I understand that it uses a lot of DSI Studio’s calculation processes, but I’m having trouble understanding all the output metrics… Also, I can’t find the DSI Studio doc that tells me the output units. I have the impression that there’s a difference, because when I compare the mrtrix FA and the one generated (DTI and DKI), the images appear identical but don’t have the same values at all!

And if you want to parcellate WM, you’d probably be better off using a tractography algorithm (e.g., pyafq_tractometry , mrtrix_multishell_msmt_pyafq_tractometry or dsi_studio_autotrack ).

Thank you very much :slight_smile:
I’ll start with the GM before trying more complicated stuff haha

5

Looking through the recon spec will tell you a bit about the workflow that is happening and some of the parameters used. Some other resources:

When you’re ready to do this, you can chain output specs together (basically copying and pasting nodes of recon specs into their own file) and have one recon spec that does everything you want to do.

Best,
Steven

6

Hi @Steven

I tried to configure my rebuild configuration file, oriented on the different metrics according to the DTI, GQI and DKI models.

Here’s the json:

{
  "name": "multishell_scalarfest_perso",
  "space" : "T1w",
  "atlases": [ ],
  "anatomical": [ ],
  "nodes": [
    {
      "name": "dsistudio_gqi",
      "software": "DSI Studio",
      "action": "reconstruction",
      "input": "qsiprep",
      "output_suffix": "gqi",
      "parameters": {"method": "gqi"}
    },
    {
      "name": "dsistudio_dti",
      "software": "DSI Studio",
      "action": "reconstruction",
      "input": "qsiprep",
      "output_suffix": "dti",
      "parameters": {"method": "dti"}
    },
    {
      "name": "scalar_export_gqi",
      "software": "DSI Studio",
      "action": "export",
      "input": "dsistudio_gqi",
      "output_suffix": "gqiscalar"
    },
    {
      "name": "scalar_export_dti",
      "software": "DSI Studio",
      "action": "export",
      "input": "dsistudio_dti",
      "output_suffix": "dtiscalar"
    },
    {
      "name": "dki_recon",
      "software": "Dipy",
      "action": "DKI_reconstruction",
      "input": "qsiprep",
      "output_suffix": "DKI",
      "parameters": {
        "write_mif": false,
        "write_fibgz": false
      }
    }
  ]
}

In the end, the scalars reconstructed using DTI have the same values as those reconstructed using GQI… But it seems to me that the values should be different, don’t you think?

7

Hi @MagicLudo,

No, I do not think they should be different. GQI and DTI are completely separate models (or rather, GQI is model-free and DTI is based off the tensor model). Your DTI metrics (FA, MD, etc) are not going to come from a GQI reconstruction. As stated in the DSI Studio documentation, DTI metrics are output alongside GQI metrics, but the fitting process that goes into both sets of metrics should be performed independently.

Where I would expect a difference would be between DTI metrics output from DIPY DTI vs DIPY DKI. As shown in Veraart et al., 2011, DTI metrics become more stable when fit in the context of DKI models (as DKI models are actually composed of the DTI term and the kurtosis term; modeling the kurtosis gives a better estimate of the tensor).

Best,
Steven