I want to use XCP-D’s aroma or aroma_gsr pipeline to post-process minimally preprocessed ICA-AROMA outputs from fmriprep to generate functional connectivity matrix and fALFF from aparc+aseg regions of Desikan Killiany atlas in FreeSurfer. Since ALFF are generated from bandpass filtered time series, what would be the best way to get fALFF from there on? Should I re-run XCP-D with no filtering and get power spectrum of post-processed aroma/aroma_gsr’ed time series without any filtering and use that to get fALFF?
Also, is it possible to run multiple pipelines at once, as in can I run aroma and aroma_gsr at the same time with -p flag? I tried doing it with -p aroma -p aroma_gsr flags to my command but aroma outputs were overwritten by aroma_gsr so not sure if it is possible?
In the document they suggest to use custom parcellation for Nifti processing pipeline using NiftiLabelsMasker, however it does not cover any steps or guidelines. Is there any documentation to be able to incorporate this with XCP-D?
Running XCP-D with filtering, then without filtering, and using that to calculate fALFF sounds like a reasonable way to do it.
I don’t know much about fALFF, but FWIW, based on a discussion I had with Dr. Valerie Sydnor, it seems like fALFF is correlated with motion, per Satterthwaite et al. (2013).
No, XCP-D can’t use multiple denoising strategies in the same run. I would recommend running XCP-D twice; once with each nuisance group and a unique output folder.
Nilearn has examples that show how to use the NiftiLabelsMasker (e.g., this one). Unfortunately, XCP-D does a bit of extra processing to deal with missing voxels/vertices, so using NiftiLabelsMasker directly will produce slightly different results. I can’t think of a particularly easy way to do this, but I can draft something up and post it in GitHub - PennLINC/xcp_d-examples: Examples of advanced use-cases for xcp_d. later today.
Running XCP-D with filtering, then without filtering, and using that to calculate fALFF sounds like a reasonable way to do it.
I don’t know much about fALFF, but FWIW, based on a discussion I had with Dr. Valerie Sydnor, it seems like fALFF is correlated with motion, per Satterthwaite et al. (2013).
No, XCP-D can’t use multiple denoising strategies in the same run. I would recommend running XCP-D twice; once with each nuisance group and a unique output folder.
Thank you for your prompt response and these pointers. We will be using both ALFF and fALFF to measure cognitive outcomes post stroke. Hoping ica-aroma and additional regressors from aroma and aroma_gsr pipeline will be enough to deal with motion and capture neuronal signal than noise.
Nilearn has examples that show how to use the NiftiLabelsMasker (e.g., this one). Unfortunately, XCP-D does a bit of extra processing to deal with missing voxels/vertices, so using NiftiLabelsMasker directly will produce slightly different results. I can’t think of a particularly easy way to do this, but I can draft something up and post it in GitHub - PennLINC/xcp_d-examples: Examples of advanced use-cases for xcp_d. later today.
There’s one cell for dealing with NIfTIs and one for CIFTIs. The NIfTI code requires ANTS, while the CIFTI code requires Connectome Workbench, so it’s best to copy over the code to a script and run it through the XCP-D Docker image.
No worries, and thank you so much. Let me give this a try and will let you know if I come across any hiccups or success. Before these steps, are we required to have our custom atlas in a certain orientation or something?
traits.trait_errors.TraitError: Cannot set the undefined 'temporal_mask' attribute of a '_NiftiConnectInputSpec' object.
Following are my steps, am I doing anything wrong? I am skipping steps to warp my atlas to BOLD space as I ran freesurfer via fmriprep pipeline which already aligns freesurfer’s aparc+aseg segmentations to my supplied output space in the fmriprep command (MNI152NLin6Asym) and resamples it to the BOLD grid.
import os
from xcp_d.interfaces.ants import ApplyTransforms
from xcp_d.interfaces.connectivity import NiftiConnect
from xcp_d.interfaces.nilearn import IndexImage
from xcp_d.utils.utils import get_std2bold_xfms
xcpd_dir = "/xcp_d/sub-fcs65c/ses-20091124/func"
fmriprep_dir = "/fmriprep_output/sub-fcs65c/ses-20091124/func"
mask = os.path.join(fmriprep_dir, "sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-brain_mask.nii.gz")
temporal_mask = os.path.join(xcpd_dir, "sub-fcs65c_ses-20091124_task-rest_run-002_outliers.tsv")
min_coverage = 0.5
correlate = True
denoise_bold_file = os.path.join(xcpd_dir, "sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-denoised_bold.nii.gz")
warped_aparcaseg_atlas = os.path.join(fmriprep_dir, "sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-aparcaseg_dseg.nii.gz")
aparcaseg_atlas_labels = "/fmriprep_output/desc-aparcaseg_dseg.tsv"
interface = NiftiConnect(
filtered_file=denoise_bold_file,
mask=mask,
temporal_mask=temporal_mask,
atlas=warped_aparcaseg_atlas ,
atlas_labels=aparcaseg_atlas_labels ,
min_coverage=min_coverage,
correlate=True,
)
traits.trait_errors.TraitError: Cannot set the undefined 'temporal_mask' attribute of a '_NiftiConnectInputSpec' object.
I was able to run xcp_d with aroma pipeline both with and without bandpass filtering (0.008 - 0.09 Hz). Would following AFNI commands be appropriate to use to compute fALFF from bandpassed ALFF and non-filtered aroma denoised time series?
# Get standard deviation of non-filtered denoised time series - this serves as the denominator for fALFF computation
3dTstat -stdev -mask /fmriprep_output/sub-fcs65c/ses-20091124/func/sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-brain_mask.nii.gz \
-prefix /xcp_d_aroma/xcp_d/sub-fcs65c/ses-20091124/func/stdev_sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-denoised_bold.nii.gz \
/xcp_d_aroma/xcp_d/sub-fcs65c/ses-20091124/func/sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-denoised_bold.nii.gz
# Compute fALFF from ALFF
3dcalc -a /fmriprep_output/sub-fcs65c/ses-20091124/func/sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-brain_mask.nii.gz \
-b /xcp_d_aroma_bpf/xcp_d/sub-fcs65c/ses-20091124/func/sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_alff.nii.gz \
-c /xcp_d_aroma/xcp_d/sub-fcs65c/ses-20091124/func/stdev_sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_desc-denoised_bold.nii.gz \
-expr '(1.0*bool(a))*((1.0*b)/(1.0*c))' \
-float -prefix /xcp_d_aroma_bpf/xcp_d/sub-fcs65c/ses-20091124/func/sub-fcs65c_ses-20091124_task-rest_run-002_space-MNI152NLin6Asym_res-2_falff.nii.gz
From init_alff_wf workflow of this source file, it seems that high-motion outlier volumes are already removed before filtering. My understanding is that this file is *space-MNI152NLin6Asym_res-2_desc-denoised_bold.nii.gz and I would not have to further postprocess this denoised file in terms of removing any outlier volumes or anything else to compute the denominator for fALFF computation?
Dear @tsalo I wanted to confirm if the above computation for fALFF made sense? Do I need to take any precautionary step to compute the denominator from denoised time series?
I’ve never calculated fALFF myself, but at a glance it looks like you’re essentially dividing the ALFF map by the standard deviation map from the unfiltered, denoised data. I don’t think the standard deviation is what you need. I think you need to essentially calculate ALFF on your unfiltered, denoised data. Perhaps the code XCP-D uses could be helpful for that: https://github.com/PennLINC/xcp_d/blob/e7d9bab99bcef4dc446cff348884137757ab6039/xcp_d/utils/restingstate.py#L94
Thank you for your response. Yes, it makes sense what you suggest and that’s the approach I am taking. I was initially trying to compute fALFF w/o having to recompute ALFF to eliminate unknowns and redundancy. I essentially added following two lines in compute_alff.py function of this script right after ALFF is computed to get fALFF from denoised func timeseries (*_desc-denoised_bold.nii.gz) generated from pipeline with bandpass filter enabled and censoring disabled --lower-bpf .008 --upper-bpf .09 --fd-thresh 0. I had to use *_desc-denoised_bold.nii.gz as an input file, which is the output of docker command with --lower-bpf .008 --upper-bpf .09 --fd-thresh 0 flags in it for the ALFF to match.
I am planning to either run this as a standalone script to just get fALFF in addition to ALFF or fork github repository and add these lines and make few downstream changes to reflect the change (if that is okay?)
# falff is sum(alff)/sum(full power spectrum range). Re-using the ff_alff lower bound so the denominator is power ([f_alff[0]-inf])
falff[ii] = np.sum(power_spec_density_sqrt[ff_alff[0]: ff_alff[1]]) / np.sum(power_spec_density_sqrt[ff_alff[0]:])
falff = np.reshape(falff, [len(falff), 1])
Thank you once again! If I may make a small request (knowing that fALFF is also used a lot in addition to other resting state derivatives) and if it is not too much of a hassle it might be best to add it in the next release.
Hello,
I am also trying to compute fALFF, but I am running XCP_D with singularity. So, in an attempt to recreate the code to compute fALFF, my current plan is to:
use the denoised image from a run of xcp_d without filtering as input (which would be the denoised, unfiltered data that @tsalo suggests, right?) to what the function read_ndata() does (in order to get the image in numpy.ndarray format)
a. since it’s a nifti, it will look to apply a mask, is the mask used here from fmriprep?
replicate compute_alff() using the numpy.ndarray output from above, manually enter TR, and use the same lowpass and hipass filter values from when I computed alff (0.01 - 0.08)
1] The way I replicated the pipeline so that everything is consistent up until ALFF estimation is that I used the pipeline with filter’s enabled and no threshold for frame displacement. My docker command for xcp_d is as follow for reference.
The mask that is used is from the output of fmriprep, the one that ends with *desc-brain_mask.nii.gz inside the func folder
2-3] I had to append the code for computation of fALFF within the routine of compute_alff() itself as it uses power_spec_density_sqrt in its computation which is used in ALFF estimation as well. I ran this routine outside of xcp_d pipeline as integrating it within would be a major task within a small time frame.
Hi Sneha, thank you for that response. I have two questions, 1. if you are using desc-denoised_bold.nii.gz generated from the pipeline with the bandpass filter enabled as input, isn’t this filtered denoised data? For fALFF, don’t you want the unfiltered data as a starting point since you are looking at the ratio of ALFF to the full frequency band? 2. Could you expand on what you mean by “I ran this routine outside of xcp_d pipeline”?
1] So I used bpf along with *_desc-denoised_bold.nii.gz as I needed ALFF computation as well which required it. However if you see the way fALFF is computed it does not directly take ALFF as its input rather it uses power_spec_density_sqrt function.
2] I basically created a new script (lets say compute_falff.py) and replicated everything that is within xcp_d’s compute_alff.py script in it. Right after ALFF is computed I added following 2 lines. I just commented out everything related to ALFF computation in it to avoid any redundancy (or you can delete those lines entirely), since ALFF was already computed by the docker routine I invoked. Once xcp_d is completed I then explicitly called this compute_falff.py script to get fALFF. You will need to pass bold file, bold mask, atlas files, minimum coverage and possibly few other parameters as well as an input to this script. I further used NiftiConnect workflow to extract regional fALFF from a warped atlas.