Error in ```makeTissueMasks``` Function: data shape issue

[littlehegemon]

I am encountering an error when running my script using the makeTissueMasks function from fmriprepciftify_helpers.py. The error seems to be related to the number of values being unpacked from the data shape. Below is a snippet of my code and the error message I receive:

config.pipelineName = 'A'
config.subject = '100206'
config.fmriRun = 'rfMRI_REST1_LR'
config.DATADIR = '<HCP path>'
config.isCifti = True
config.fmriFile = "<HCP path>/100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii"

print('Step 0 : Building WM, CSF and GM masks...')
masks = makeTissueMasks(overwrite=True)

masks = makeTissueMasks(overwrite=True)
  File "fmriprepciftify_helpers.py", line 634, in makeTissueMasks
    nRows, nCols, nSlices = tmpWM.header.get_data_shape()
ValueError: too many values to unpack (expected 3)```

Upon reviewing [the code](https://github.com/adolphslab/rsDenoise/blob/c28a0298af0eb6933a2a62067299b294ceb9e8b5/HCP_helpers.py#L467), it appears the issue might be related to the fsl.FLIRT function, as seen in the lines below:

```flirt_wmparc = fsl.FLIRT(in_file=wmparcFilein, out_file=wmparcFileout,
                         reference=fmriFile, apply_xfm=True,
                         in_matrix_file=eyeMat, out_matrix_file=wmparcMat, interp='nearestneighbour')
ref = nib.load(wmparcFileout)
img = nib.Nifti1Image(WMmask.reshape(ref.shape).astype('<f4'), ref.affine)
nib.save(img, WMmaskFileout)
tmpWM = nib.load(WMmaskFileout)

The shape of the image appears to be determined by fsl.FLIRT. I suspect the error might be due to differences in FSL versions, as the helper script was written approximately seven years ago.

I executed the flirt command directly in bash using FSL version 6.0.5.2, and it completed successfully:

flirt -in <HCP path>100206/MNINonLinear/wmparc.nii.gz -ref <HCP path>100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii -out test -applyxfm -init eye.mat -omat wmparc_flirt.mat -interp nearestneighbour

However, when I tried using FSL 5.0.10 (released 25th April 2017) and FSL 5.0.11 (released 23rd March 2018). the flirt command failed, likely due to incompatibility with the CIFTI format:

** ERROR (nifti_image_read): cannot create nifti image from header '<HCP path>/100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii'
** ERROR: nifti_image_open(<HCP path>/100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries): bad header info
ERROR: failed to open file <HCP path>/100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries
Image Exception : #22 :: ERROR: Could not open image <HCP path>/100206/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries
terminate called after throwing an instance of 'RBD_COMMON::BaseException'
Aborted (core dumped)

Could you please help me identify if the error is indeed due to version differences in FSL or if there is another issue at play? Any guidance on how to resolve this error would be greatly appreciated.

Thank you for your assistance!

[littlehegemon]

It seems flirt in different FSL versions are in the same version.

$ /usr/local/fsl/bin/flirt  -version
FLIRT version 6.0
$ /usr/local/fsl5.0.11/bin/flirt -version
FLIRT version 6.0

[paola-g]

Hello littlehegemon, and thanks for the interest in our pipeline! As you observed, the fmriprepciftify_helpers.py has not been updated in many years, because we moved to an implementation that starts directly from fmriprep outputs bypassing ciftify (fmriprep_helpers.py).

However if you are using HCP data, we have a separate helpers file that should work natively with cifti data (HCP_helpers.py). Is this your case?

[littlehegemon]

Hi Paola,

Thank you for your reply.

After several days of testing, I have successfully replicated our lab's previous results. I would like to document some points for future reference.

1. Versions

   Use the latest numpy. Although rsDenoise recommend to use numpy==1.12.1, when I use that version, fit = np.linalg.lstsq(X, data.T, rcond=None)[0] would show error. This is because in older versions, the rcond argument in np.linalg.lstsq() cannot be set to None. In more recent versions, the default value of rcond is None, which works fine.

   Python 3.10 works and 3.12 would not be compatible with a package.

   scipy==1.9.3 works and the latest version 1.14.1 would show module 'scipy.signal' has no attribute 'gaussian'.

2. How to run without qsub:

   To bypass the use of qsub, simply replace all instances of "qsub" with "bash" in HCP_helpers.py. I also replaced all occurrences of "wb_command" with the full path on my machine.

   Here’s the code I used:

   from HCP_helpers import *
   import glob
   
   config.pipelineName = 'A'
   config.Operations   = config.operationDict[config.pipelineName]
   config.DATADIR = '<HCP data path>'
   config.isCifti = True
   config.ext = '.dtseries.nii'
   config.useFIX = True
   config.movementRegressorsFile = 'Movement_Regressors_dt.txt'
   config.outDir = "<out path>"
   pattern = "<HCP data path>/100206/MNINonLinear/Results/rfMRI_REST*/rfMRI_REST*_Atlas_hp2000_clean.dtseries.nii"
   paths = glob.glob(pattern)
   
   for path in paths:
      path_parts = path.split("/")
      print(path_parts[-5], path_parts[-2])
      config.subject = path_parts[-5]
      config.fmriRun = path_parts[-2]
      config.fmriFile = path
      runPipelinePar()

   You can change "100206" to "*" if you wish to run the script for multiple subjects.

3. Difference between rsDenoise and HCP_MRI-behavior:

   I noticed a significant difference in results when using HCP_helpers.py from the rsDenoise repository versus the one in HCP_MRI-behavior.

   Below is a comparison for subject 100206, rfMRI_REST1_LR using the same code (above), but with different versions of HCP_helpers.py from each repository. The X axis is the 90000 vertex. Each point represent the correlation of 1200 timepoints of activation.

   

4. Necessary files to download:

   To download the required files from Amazon S3 to the current working directory:

   [ ! -f id_list.txt ] && aws s3 ls s3://hcp-openaccess/HCP_1200/ | awk '{print $2}' | sed 's#/##' > id_list.txt
   
   cat id_list.txt | parallel -j20 "/aws s3 sync s3://hcp-openaccess/HCP_1200/{}/ ./HCP_1200/{}/ --exclude '*'  \
   --include 'MNINonLinear/ribbon.nii.gz' --include 'MNINonLinear/wmparc.nii.gz' \
   --include 'MNINonLinear/Results/rfMRI_REST1_LR/Movement_Regressors_dt.txt' --include 'MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR.nii.gz' --include 'MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_hp2000_clean.dtseries.nii'  \
   --include 'MNINonLinear/Results/rfMRI_REST1_RL/Movement_Regressors_dt.txt' --include 'MNINonLinear/Results/rfMRI_REST1_RL/rfMRI_REST1_RL.nii.gz' --include 'MNINonLinear/Results/rfMRI_REST1_RL/rfMRI_REST1_RL_Atlas_hp2000_clean.dtseries.nii'  \
   --include 'MNINonLinear/Results/rfMRI_REST2_LR/Movement_Regressors_dt.txt' --include 'MNINonLinear/Results/rfMRI_REST2_LR/rfMRI_REST2_LR.nii.gz' --include 'MNINonLinear/Results/rfMRI_REST2_LR/rfMRI_REST2_LR_Atlas_hp2000_clean.dtseries.nii' \
   --include 'MNINonLinear/Results/rfMRI_REST2_RL/Movement_Regressors_dt.txt' --include 'MNINonLinear/Results/rfMRI_REST2_RL/rfMRI_REST2_RL.nii.gz' --include 'MNINonLinear/Results/rfMRI_REST2_RL/rfMRI_REST2_RL_Atlas_hp2000_clean.dtseries.nii'"

   You can adjust -j20 to set the desired number of concurrent threads.

Best regards.

[paola-g]

Many thanks for reporting these results. I have to confess I haven't updated the requirements in the README while moving to new versions of the underlying libraries so this is on me. Personally, I would always go with the latest libraries and adjust the rsDenoise code rather than the other way around.

About avoiding qsub, there is a way to switch off that feature by setting config.queue = False. Also, as long as wb_command is in your PATH it should be able to pick it up, but your workaround is also reasonable.

I will need to look into the difference between rsDenoise and HCP_MRI-behaviour. In the original repository, we kept a release of the original code we used in the publications to allow reproducibility, but since then most of the development has happened in this repository. However, we mostly focused on the fmriprep_helpers rather than the HCP_helpers because we have been working with datasets other than HCP, so I need to study the commit history to explain the discrepancies.