Qsiprep : how to toggle distortion correction properly with my dataset?

quentinduche · November 4, 2022, 3:06pm

Hi everyone,

I am new to both qsiprep and the BIDS format but I managed to make my dataset compatible with qsiprep.

I am using qsiprep-0.16.1.
Data were BIDS-validated with bids-validator 1.9.9.

However, I encountered problems on 66/69 subjects on my dataset. They all underwent a raw_rpe_concat node crash (more info below, in point ).

I have a dataset of 69 subjects where each participant underwent

a T1-weighted MPRAGE scan
diffusion weighted images with opposite phase encoding directions (PED)

A DW volume was obtained at b-value 0 in the PA PED, shape = (148, 148, 90)
67 DW volumes were obtained at different b-values in the AP PED, shape = (148, 148, 90, 67)

DW images are severely distorted and must be corrected for that using the PEPOLAR TOPUP technique.
I ran qsiprep in two different configurations for my dataset.

Configuration
I used the following BIDS-valid dataset :

├── anat
│   ├── sub-VS006_T1w.json
│   └── sub-VS006_T1w.nii.gz
├── dwi
│   ├── sub-VS006_acq-cusp66b0_dir-pa_sbref.bval
│   ├── sub-VS006_acq-cusp66b0_dir-pa_sbref.bvec
│   ├── sub-VS006_acq-cusp66b0_dir-pa_sbref.json
│   ├── sub-VS006_acq-cusp66b0_dir-pa_sbref.nii.gz (3D file, shape = (148,148,90))
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bval
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bvec
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.json
│   └── sub-VS006_acq-cusp66b3000_dir-ap_dwi.nii.gz  (4D file, shape = (148,148,90,67))

I had individual HTML reports for all subjects
Pre-processed DWI looked good except for…
distortion correction was not performed

Configuration
I used the following BIDS-valid dataset.
I adjusted the dataset by doing a tiny trick to try to toggle the PEPOLAR TOPUP distortion correction step.
My 3D PA single volume (shape = (148, 148, 90)) was converted into a 4D volume (shape = (148, 148, 90, 2)) by concatenating twice the same single PA volume.

├── anat
│   ├── sub-VS006_T1w.json
│   └── sub-VS006_T1w.nii.gz
├── dwi
│   ├── sub-VS006_acq-cusp66b0_dir-pa_dwi.bval
│   ├── sub-VS006_acq-cusp66b0_dir-pa_dwi.bvec
│   ├── sub-VS006_acq-cusp66b0_dir-pa_dwi.json
│   ├── sub-VS006_acq-cusp66b0_dir-pa_dwi.nii.gz (4D file, shape = (148,148,90,2))
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bval
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bvec
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.json
│   └── sub-VS006_acq-cusp66b3000_dir-ap_dwi.nii.gz (4D, file, shape = (148,148,90,67))

This indeed toggled the distortion correction but ended in a different system output.
I did not have individual HTML reports (only one subject succeeded in creating the html report)
However, pre-processed DWI looked good
DWI were unambiguously corrected for distortions
I had 66 raw_rpe_concat crashes (example crash file below)
For these 66 subjects two files were systematically missing compared to the outputs obtained for configuration sub-XXX_desc-ImageQC_dwi.csv and sub-XXX_desc-SliceQC_dwi.json

Content of an example crash file :

Node: qsiprep_wf.single_subject_VS006_wf.dwi_preproc_wf.pre_hmc_wf.raw_rpe_concat
Working directory: /work/qsiprep_wf/single_subject_VS006_wf/dwi_preproc_wf/pre_hmc_wf/raw_rpe_concat

Node inputs:

compress = True
dtype = f4
header_source = <undefined>
in_files = <undefined>
is_dwi = True

Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node
    result["result"] = node.run(updatehash=updatehash)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 527, in run
    result = self._run_interface(execute=True)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 645, in _run_interface
    return self._run_command(execute)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 771, in _run_command
    raise NodeExecutionError(msg)
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node raw_rpe_concat.

Traceback:
	Traceback (most recent call last):
	  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/base/core.py", line 398, in run
	    runtime = self._run_interface(runtime)
	  File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/interfaces/nilearn.py", line 134, in _run_interface
	    new_nii = concat_imgs(self.inputs.in_files, dtype=self.inputs.dtype)
	  File "/usr/local/miniconda/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py", line 487, in concat_niimgs
	    for index, (size, niimg) in enumerate(zip(lengths, _iter_check_niimg(
	  File "/usr/local/miniconda/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py", line 160, in _iter_check_niimg
	    raise ValueError(
	ValueError: Field of view of image #1 is different from reference FOV.
	Reference affine:
	array([[-1.49992088e+00,  1.51489924e-06,  1.54058855e-02,
	         1.04946426e+02],
	       [-3.61095565e-03,  1.45818010e+00, -3.51706936e-01,
	        -8.49841766e+01],
	       [ 1.49767256e-02,  3.51725472e-01,  1.45810318e+00,
	        -8.69976807e+01],
	       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
	         1.00000000e+00]])
	Image affine:
	array([[-1.49992108e+00, -8.02594400e-08,  1.53924711e-02,
	         1.04946426e+02],
	       [-3.60936113e-03,  1.45818007e+00, -3.51706922e-01,
	        -8.49841766e+01],
	       [ 1.49633121e-02,  3.51725459e-01,  1.45810342e+00,
	        -8.69976807e+01],
	       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
	         1.00000000e+00]])
	Reference shape:
	(148, 148, 90)
	Image shape:
	(148, 148, 90, 67)

Here are my questions :

Can the trick (3D to 4D volume) be the reason of the crash ?
How should I organize my dataset to properly toggle the PEPOLAR TOPUP distortion correction with qsiprep ?
Can I use the pre-processed DWI file ?

Thank you very much for your help.

Cheers
Quentin

jsein · November 4, 2022, 4:34pm

Hi Quentin, welcome to Neurostars!

I think you should use the extension _dwi for your PA scans. And no need to concatenate your b0 PA scan to become a 4D image. QSIPREP will detect the b0 images in both sets (dir-ap and dir-pa) and use them automatically for topup.

Could you look at the affine of the AP and PA images to check if their affine differ as it is said in the error message (with fslhd for instance). Are your images AP and PA acquired with the same FOV? If so the affines should not differ.
In your case it looks like the affine are differing by a very small amount, which looks like what is reported here and recently corrected in FMRIPREP (relaxing the affine difference tolerance):

github.com/nipreps/fmriprep

concat_blips

opened 04:19PM - 07 Feb 22 UTC

closed 04:08PM - 25 Apr 22 UTC

juansanchezpena

bug

### What happened? I am unable to produce func results for a data that looks… standard. ### What command did you use? ```shell /usr/local/pipelines/subfmri-arg_latest -i /MRI_DATA/nyspi/test/rawdata -o /MRI_DATA/nyspi/test/derivatives -l /MRI_DATA/nyspi/test/derivatives/fmriprep/logs/gold.log -a "--participant-label test03gold --output-spaces T1w MNI152NLin2009cAsym:res-2 --cifti-output --output-layout legacy” ``` ### What version of fMRIPrep are you running? 21.0 ### How are you running fMRIPrep? Docker ### Is your data BIDS valid? Yes ### Are you reusing any previously computed results? No ### Please copy and paste any relevant log output. ```shell Node Name: fmriprep_wf.single_subject_test03gold_wf.fmap_preproc_wf.wf_auto_00000.concat_blips File: /output/fmriprep/sub-test03gold/log/20220204-140358_cfdc6d40-fc00-456a-9b71-4a44a81323ce/crash-20220204-140521-root-concat_blips-138212d8-935a-4164-a554-5ca62213c680.txt Working Directory: /tmp/work/fmriprep_wf/single_subject_test03gold_wf/fmap_preproc_wf/wf_auto_00000/concat_blips Inputs: allow_4D: True in_files: Traceback (most recent call last): File "/opt/conda/lib/python3.8/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node result["result"] = node.run(updatehash=updatehash) File "/opt/conda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 516, in run result = self._run_interface(execute=True) File "/opt/conda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 635, in _run_interface return self._run_command(execute) File "/opt/conda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 741, in _run_command result = self._interface.run(cwd=outdir) File "/opt/conda/lib/python3.8/site-packages/nipype/interfaces/base/core.py", line 428, in run runtime = self._run_interface(runtime) File "/opt/conda/lib/python3.8/site-packages/niworkflows/interfaces/nibabel.py", line 191, in _run_interface img_4d = nb.concat_images(nii_list) File "/opt/conda/lib/python3.8/site-packages/nibabel/funcs.py", line 141, in concat_images raise ValueError(f'Affine for image {i} does not match affine for first image') ValueError: Affine for image 2 does not match affine for first image ``` ### Additional information / screenshots RAWDATA: anat: sub-test03gold_ses-21117_acq-cubepromo_T2w.json sub-test03gold_ses-21117_acq-nopromo_T1w.json sub-test03gold_ses-21117_acq-cubepromo_T2w.nii.gz sub-test03gold_ses-21117_acq-nopromo_T1w.nii.gz fmap: sub-test03gold_ses-21117_dir-fpe_epi.json sub-test03gold_ses-21117_dir-rpe_epi.json sub-test03gold_ses-21117_dir-fpe_epi.nii.gz sub-test03gold_ses-21117_dir-rpe_epi.nii.gz func: sub-test03gold_ses-21117_task-rest_run-01_bold.json sub-test03gold_ses-21117_task-rest_run-01_bold.nii.gz This same raw data went through abcd

A workaround could be to use fslcpgeom to force the affine of all your images to be identical.

In your case, you also should use the --separate_all_dwis argument

mattcieslak · November 4, 2022, 5:23pm

In this case, where you only have one 3d volume with the opposite phase encoding direction, I’d put it in the fmap/ directory with the _epi suffix.

Steven · November 4, 2022, 6:29pm

Hi @quentinduche,

Since you are new to BIDS, I feel I should preemptively add a bit of detail to @mattcieslak’s answer. For your new “fmap” file (your relocated b0 image), you should add a field to your _epi.json file called IntendedFor. The argument for this field should map to your DWI image (relative to the inside of a subject folder). For example:

IntendedFor: ["dwi/sub-VSXXX_acq-cusp66b3000_dir-ap_dwi.nii.gz"]

That way, BIDS apps such as QSIPrep will know to associate that specific fieldmap to a particular DWI scan.

Best,
Steven

quentinduche · November 7, 2022, 11:04am

Hello Julien, @mattcieslak and @Steven and thank you very much for your feedback.

I moved the 3D volume into the fmap directory and added the "IntendedFor" key in the associated .json file and that worked like a charm .

Not a single raw_rpe_concat crash
I have HTML reports
Pre-processed DWI look good
DWI are unambiguously corrected for distortions (animation in html report is )

This is what my dataset looked like :

├── anat
│   ├── sub-VS006_T1w.json
│   └── sub-VS006_T1w.nii.gz
├── dwi
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bval
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.bvec
│   ├── sub-VS006_acq-cusp66b3000_dir-ap_dwi.json
│   └── sub-VS006_acq-cusp66b3000_dir-ap_dwi.nii.gz
├── fmap
│   ├── sub-VS006_acq-cusp66b0_dir-pa_epi.json (with "IntendedFor" : ["dwi/sub-VS006_acq-cusp66b3000_dir-ap_dwi.nii.gz"])
│   ├── sub-VS006_acq-cusp66b0_dir-pa_epi.nii.gz  (3D file, shape = (148,148,90))

I encountered other crashes but they are not related to my original question :

14 raw_gqi crashes (similar bug than reported here)
1 dwi_rpt crash

Node: qsiprep_wf.single_subject_VS056_wf.dwi_preproc_acq_cusp66b3000_dir_ap_wf.fmap_unwarp_report_wf.dwi_rpt
Working directory: /work/qsiprep_wf/single_subject_VS056_wf/dwi_preproc_acq_cusp66b3000_dir_ap_wf/fmap_unwarp_report_wf/dwi_rpt

Node inputs:

after = <undefined>
before = <undefined>
compress_report = auto
out_report = report.svg
wm_seg = <undefined>

Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node
    result["result"] = node.run(updatehash=updatehash)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 527, in run
    result = self._run_interface(execute=True)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 645, in _run_interface
    return self._run_command(execute)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 771, in _run_command
    raise NodeExecutionError(msg)
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node dwi_rpt.

Traceback:
        Traceback (most recent call last):
          File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/base/core.py", line 399, in run
            runtime = self._post_run_hook(runtime)
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/niworkflows/interfaces/registration.py", line 312, in _post_run_hook
            return super(SimpleBeforeAfterRPT, self)._post_run_hook(runtime)
          File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/mixins/reporting.py", line 50, in _post_run_hook
            self._generate_report()
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/niworkflows/interfaces/report_base.py", line 59, in _generate_report
            cuts = cuts_from_bbox(contour_nii, cuts=n_cuts)
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/niworkflows/viz/utils.py", line 205, in cuts_from_bbox
            start_coords = B.min(0)
          File "/usr/local/miniconda/lib/python3.8/site-packages/numpy/core/_methods.py", line 44, in _amin
            return umr_minimum(a, axis, None, out, keepdims, initial, where)
        ValueError: zero-size array to reduction operation minimum which has no identity

1 gather_inputs crash

ode: qsiprep_wf.single_subject_VS054_wf.dwi_preproc_acq_cusp66b3000_dir_ap_wf.hmc_sdc_wf.gather_inputs
Working directory: /work/qsiprep_wf/single_subject_VS054_wf/dwi_preproc_acq_cusp66b3000_dir_ap_wf/hmc_sdc_wf/gather_inputs

Node inputs:

b0_threshold = 100
bval_file = <undefined>
bvec_file = <undefined>
dwi_file = <undefined>
epi_fmaps = ['/bids_dataset/sub-VS054/fmap/sub-VS054_acq-cusp66b0_dir-pa_epi.nii.gz']
original_files = <undefined>
raw_image_sdc = True
topup_max_b0s_per_spec = 3
topup_requested = True

Traceback (most recent call last):
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/plugins/multiproc.py", line 67, in run_node
    result["result"] = node.run(updatehash=updatehash)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 527, in run
    result = self._run_interface(execute=True)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 645, in _run_interface
    return self._run_command(execute)
  File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py", line 771, in _run_command
    raise NodeExecutionError(msg)
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node gather_inputs.

Traceback:
        Traceback (most recent call last):
          File "/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/base/core.py", line 398, in run
            runtime = self._run_interface(runtime)
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/interfaces/eddy.py", line 76, in _run_interface
            get_best_b0_topup_inputs_from(
          File "/usr/local/miniconda/lib/python3.8/site-packages/qsiprep/interfaces/epi_fmap.py", line 212, in get_best_b0_topup_inputs_from
            raise Exception("Unable to run TOPUP: not enough distortion groups. "
        Exception: Unable to run TOPUP: not enough distortion groups. Check "IntendedFor" fields or consider using --ignore fieldmaps.

mattcieslak · November 10, 2022, 4:38pm

I think you may be running out of memory or using too many cpus here. What are your arguments for --nthreads and --omp-nthreads?

quentinduche · November 10, 2022, 6:28pm

I did not specify arguments for --nthreads and --omp-nthreads

mattcieslak · January 31, 2023, 1:19am

what are the PhaseEncodingDirections of your fmaps and dwi images?