Tabulated ReHo values don't match CIFTI outputs

[fjc20]

Summary

The tabulated (.tsv) outputs for ReHo data appear dramatically different, and sometimes missing entirely, compared to the actual CIFTI output.

Additional details

• xcp_d version: v0.7.3
• Docker version: 26.0.1, build d260a54

What were you trying to do?

Compute ReHo values for resting state fMRI data analyze regional values using the Glasser-atlas parcellated data. Use CIFTI outputs to confirm.

What did you expect to happen?

CIFTI and parcellated data should show similar spatial patterns.

What actually happened?

Ran various resting state fMRI files through the xcp-d pipeline and compared the tabulated regional ReHo values for the Glasser atlas (visualized with ggseg) to the actual CIFTIs (visualized in connectome workbench as well as the ciftiTools wb library). The vertex-wise mean values across subjects for the CIFTIs look very reasonable-- highly symmetric, consistent with what we've seen in volumetric analyses, and reasonably similar to previously published reports. However, the tabulated data is missing entirely from ~16 ROIs, and when visualizing with ggseg produces (apparently) randomly organized, asymmetric spatial maps. Extracting ROI level data manually from the CIFTI files produces outputs which do not match the tabulated files.

Reproducing the bug

Standard XCP-d processing using fMRIprep outputs.

[tsalo]

Can you share the code you used? Namely, (1) the fMRIPrep command you used, (2) the XCP-D command you used, and (3) any code you used to find and characterize the problem.

[WillForan]

Here are the two scripts for fmriprep and xcpd fmriprep+xcpd_luna7t-bash.zip and the example xcp-d and wb_command output

The scripts run

fmriprep

docker  run --rm \
     -v /Volumes:/Volumes  \
     -v "$FSlic:$FSlic" \
     nipreps/fmriprep:23.2.1 \
      --skip_bids_validation \
      `#--bids-filter-file "$filter"` \
      --fs-license-file "$FSlic" \
      --fs-subjects-dir "$FS" \
      --task-id "rest" \
      "$SINGLEBIDS" "$PREPOUT" participant \
      --participant-label "${ld8/_*/}" \
      --ignore slicetiming  \
      --skip-bids-validation \
      --cifti-output

(input rest nifti is 3depi, so we disable slicetiming)

xcp_d

docker run --rm `#-it` \
      -v "$ses_dataset_dir":/fmriprep:ro \
      -v "$FSlic:$FSlic:ro" \
      -v /Volumes/Hera/tmp/wkdir:/work:rw \
      -v "$out_dir":/out:rw \
      -v "${FS_dir}":/freesurfer:ro \
      -v "/opt/ni_tools/templateflow:/opt/ni_tools/templateflow" \
      -e "TEMPLATEFLOW_HOME=/opt/ni_tools/templateflow" \
      pennlinc/xcp_d:latest  \
      /fmriprep /out participant \
      --participant_label "$sub" \
      --fs-license-file "$FSlic" \
      --cifti --despike \
      --head_radius 40 -w /work \
      --smoothing 4 \
      -p 36P \
      --fd-thresh 0

---

Characterization

to compare xcp output with expected @valeriejill ran

wb_command -cifti-parcellate \
    sub-11676_task-rest_run-01_space-fsLR_den-91k_stat-reho_boldmap.dscalar.nii \
    ~/Software/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/atlas_dlabel_files/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii \
    COLUMN \
    sub-11676_task-rest_run-01_space-fsLR_den-91k_stat-reho_glasser.pscalar.nii

which I (maybe incorrectly) expanded to

# xcp-d outputs
reho_cifti=xcpd/ses-1/sub-11676/func/sub-11676_task-rest_run-01_space-fsLR_den-91k_stat-reho_boldmap.dscalar.nii
reho_glasser=xcpd/ses-1/sub-11676/func/sub-11676_task-rest_run-01_space-fsLR_seg-Glasser_stat-reho_bold.tsv

# atlas from xcp-d python package
glasser=/opt/ni_tools/python/virtualenv-xcpd/lib/python3.11/site-packages/xcp_d/data/atlases/tpl-fsLR_atlas-Glasser_den-32k_dseg.dlabel.nii

# prefix for temporary wb_command output files to compare 
prefix=/tmp/glasser-11676

# actually extract
wb_command -cifti-parcellate $reho_cifti $glasser COLUMN $prefix.pscalar.nii
wb_command -cifti-convert -to-text $prefix.pscalar.nii  $prefix.tsv

# quick compare
paste $prefix.tsv <(sed '1d;s/\t/\n/g' $reho_glasser)|head

I'm not sure if the roi names line up, but there are no n/a in the wb_command generated tsv but are in the xcp_d glasser reho output.

0.761217        0.5979426967424731
0.608849        0.7017811064265276
0.93336 0.6262983058283969
0.801778        0.6258578027772196
0.815781        0.6652791527515337
0.804304        0.6991992888317721
0.812782        0.58908518931518
0.671818        0.6957912081364289
0.69467 0.7073237402876328
0.789031        n/a

[tsalo]

Thank you for sharing your code and results. The n/a is a result of the minimum coverage threshold (--min-coverage). >50% of the vertices in that parcel are probably zeros or NaNs in the ReHo CIFTI file. I'm looking into the other differences now to see if they stem from the way XCP-D masks out zeros and NaNs or if they're the result of a bug on XCP-D's side.

[fjc20]

I don't think the n/a's are just coverage. For one thing, there are 13 ROIs that have NAs for all 244 of our sessions:

label n_NA 1 lh_L_23d 244 2 lh_L_24dv 244 3 lh_L_33pr 244 4 lh_L_46 244 5 lh_L_8Ad 244 6 lh_L_8Av 244 7 lh_L_ProS 244 8 lh_L_RSC 244 9 lh_L_SCEF 244 10 lh_L_p24pr 244 11 rh_R_6d 244 12 rh_R_FEF 244 13 rh_R_SCEF 244 The CIFTI maps don't have missingness near these ROIs. For comparison, I'm showing the mean ReHo values across all 244 sessions for the glasser atlases visualized with ggseg, compared to the vertex-wise mean of the CIFTI files themselves. Note not only the missing ROIs in the tabulated Glasser data, but also the lack of correspondence in other regions.

[tsalo]

Thanks. I'm trying to reproduce this locally, but I don't use R. Can you share the ggseg commands you used?

[tsalo]

Okay, I'm seeing something weird in the ReHo CIFTIs. Can you check if you're seeing the same problem in ALFF maps (assuming you kept them)?

[tsalo]

I've merged a PR (#1175) that should fix the problem, but I want to make sure the results look good before I make a new release.

[WillForan]

rerunning xcp_d now! Will update comment when it finishes.

source /opt/ni_tools/python/virtualenv-xcpd/bin/activate                                                                            
pip install git+https://github.com/pennlinc/xcp_d.git
pip show xcp_d|grep Version  #  Version: 0.7.4.dev10+gf61ac56

mv xcpd/ses-1/sub-11676{,.bad_reho_tsv}
xcp_d fmriprep/ses-1/ xcpd/ses-1/ participant --participant_label 11676 \
      --fs-license-file /opt/ni_tools/freesurfer/license.txt \
      --cifti --despike --head_radius 40 -w /Volumes/Hera/tmp/wkdir2/ \
      --smoothing 4  -p 36P  --fd-thresh 0

update: failing to run probably b/c I did a bad job cleaning up the previous run

240524-16:38:15,438 nipype.workflow ERROR:                                                                                            
         could not run node: xcp_d_0_7_wf.sub_11676_wf.load_atlases_wf.atlas_file_grabber                                             
240524-16:38:15,611 nipype.workflow CRITICAL:                                                                                         
         XCP-D failed: 20 raised. Re-raising first.                                                                                   
240524-16:38:15,619 nipype.utils WARNING:                                                                                             
         Previous output generated by version 0.7.3 found. 

[tsalo]

I just merged #1174, which further refactors the CIFTI parcellation workflow to make it more transparent by using wb_command for most of the parcellation and connectivity steps instead of a Python function I wrote. That PR also adds plots of (1) coverage, (2) dense ALFF, (3) parcellated ALFF, (4) dense ReHo, and (5) parcellated ReHo.

[valeriejill]

Wow, that was a major overhaul and so much work!! It's awesome to have the plots for dense and parcellated metrics added. Thank you so much @tsalo!!!!

[tsalo]

Thanks. Please reopen this if you see problems in 0.7.4 (which has #1175) or unstable (which has #1174 merged in).