Running XCP-D

Warning

XCP-D may not run correctly on M1 chips.

Execution and Input Formats

The XCP-D workflow takes fMRIPRep, NiBabies, DCAN and HCP outputs in the form of BIDS derivatives. In these examples, we use an fmriprep output directory.

The outputs are required to include at least anatomical and functional outputs with at least one preprocessed BOLD image. Additionally, each of theseshould be in directories that can be parsed by the BIDS online validator (even if it is not BIDS valid - we do not require BIDS valid directories). The directories must also include a valid dataset_description.json.

The exact command to run in xcp_d depends on the installation method and data that needs to be processed. We start first with the bare-metal Manually Prepared Environment (Python 3.8+) installation, as the command line is simpler. xcp_d can be executed on the command line, processesing fMRIPrep outputs, using the following command-line structure, for example:

xcp_d <fmriprep_dir> <outputdir> --cifti --despike  --head_radius 40 -w /wkdir --smoothing 6

However, we strongly recommend using Container Technologies. Here, the command-line will be composed of a preamble to configure the container execution, followed by the xcp_d command-line options as if you were running it on a bare-metal installation.

Command-Line Arguments

xcp_d postprocessing workflow of fMRI data

usage: xcp_d [-h] [--version]
             [--participant_label PARTICIPANT_LABEL [PARTICIPANT_LABEL ...]]
             [-t TASK_ID] [-m] [-s] [--nthreads NTHREADS]
             [--omp-nthreads OMP_NTHREADS] [--mem_gb MEM_GB]
             [--use-plugin USE_PLUGIN] [-v]
             [--input-type {fmriprep,dcan,hpc,nibabies}]
             [--smoothing SMOOTHING] [--despike]
             [--nuissance-regressors {27P,36P,24P,acompcor,aroma,acompcor_gsr,aroma_gsr,custom} | -p {27P,36P,24P,acompcor,aroma,acompcor_gsr,aroma_gsr,custom}]
             [-c CUSTOM_CONFOUNDS] [-d DUMMYTIME]
             [--disable-bandpass-filter | --bandpass_filter BANDPASS_FILTER]
             [--lower-bpf LOWER_BPF] [--upper-bpf UPPER_BPF]
             [--bpf-order BPF_ORDER] [--motion-filter-type {lp,notch}]
             [--band-stop-min BPM] [--band-stop-max BPM]
             [--motion-filter-order MOTION_FILTER_ORDER] [-r HEAD_RADIUS]
             [-f FD_THRESH] [-w WORK_DIR] [--clean-workdir]
             [--resource-monitor] [--notrack] [--warp-surfaces-native2std]
             fmri_dir output_dir analysis_level

Positional Arguments

fmri_dir

the root folder of a preprocessed fMRI output .

output_dir

the output path for xcp_d

analysis_level

the analysis level for xcp_d, must be specified as “participant”.

Named Arguments

--version

show program’s version number and exit

Options for filtering BIDS queries

--participant_label, --participant-label

a space delimited list of participant identifiers or a single identifier (the sub- prefix can be removed)

-t, --task-id

select a specific task to be selected for the postprocessing

-m, --combineruns

this option combines all runs into one file

Options for cifti processing

-s, --cifti

postprocess cifti instead of nifti this is set default for dcan and hcp

Options to for resource management

--nthreads

maximum number of threads across all processes

--omp-nthreads

maximum number of threads per-process

--mem_gb, --mem_gb

upper bound memory limit for xcp_d processes

--use-plugin

nipype plugin configuration file. for more information see https://nipype.readthedocs.io/en/0.11.0/users/plugins.html

-v, --verbose

increases log verbosity for each occurence, debug level is -vvv

Input flags

--input-type

Possible choices: fmriprep, dcan, hpc, nibabies

The pipeline used to generate the preprocessed derivatives. The default pipeline is ‘fmriprep’. The ‘dcan’, ‘hcp’, and ‘nibabies’ pipelines are also supported. ‘nibabies’ assumes the same structure as ‘fmriprep’.

Parameters for postprocessing

--smoothing

smoothing the postprocessed output (fwhm)

--despike

despike the nifti/cifti before postprocessing

--nuissance-regressors

Possible choices: 27P, 36P, 24P, acompcor, aroma, acompcor_gsr, aroma_gsr, custom

Nuisance parameters to be selected, other options include 24P and 36P acompcor and aroma. See Ciric et. al (2007) for more information about regression strategies. This parameter is deprecated and will be removed in version 0.3.0. Please use -p or --nuisance-regressors.

-p, --nuisance-regressors

Possible choices: 27P, 36P, 24P, acompcor, aroma, acompcor_gsr, aroma_gsr, custom

Nuisance parameters to be selected. See Ciric et. al (2007).

-c, --custom_confounds

Custom confound to be added to nuisance regressors.

-d, --dummytime

first volume in seconds to be removed or skipped before postprocessing

Filtering parameters and default value

--disable-bandpass-filter, --disable_bandpass_filter

Disable bandpass filtering.

--bandpass_filter

Whether to Butterworth bandpass filter the data or not. This parameter is deprecated and will be removed in version 0.3.0. Bandpass filtering is performed by default, and if you wish to disable it, please use –disable-bandpass-filter`.

--lower-bpf

lower cut-off frequency (Hz) for the butterworth bandpass filter

--upper-bpf

upper cut-off frequency (Hz) for the butterworth bandpass filter

--bpf-order

number of filter coefficients for butterworth bandpass filter

--motion-filter-type

Possible choices: lp, notch

Type of band-stop filter to use for removing respiratory artifact from motion regressors. If not set, no filter will be applied.

If the filter type is set to “notch”, then both band-stop-min and band-stop-max must be defined. If the filter type is set to “lp”, then only band-stop-min must be defined.

--band-stop-min

Lower frequency for the band-stop motion filter, in breaths-per-minute (bpm). Motion filtering is only performed if motion-filter-type is not None. If used with the “lp” motion-filter-type, this parameter essentially corresponds to a low-pass filter (the maximum allowed frequency in the filtered data). This parameter is used in conjunction with motion-filter-order and band-stop-max.

Recommended values, based on participant age

Age Range (years)

Recommended Value (bpm)

< 1

30

1 - 2

25

2 - 6

20

6 - 12

15

12 - 18

12

19 - 65

12

65 - 80

12

> 80

10

When motion-filter-type is set to “lp” (low-pass filter), another commonly-used value for this parameter is 6 BPM (equivalent to 0.1 Hertz), based on Gratton et al. (2020).

--band-stop-max

Upper frequency for the band-stop motion filter, in breaths-per-minute (bpm). Motion filtering is only performed if motion-filter-type is not None. This parameter is only used if motion-filter-type is set to “notch”. This parameter is used in conjunction with motion-filter-order and band-stop-min.

Recommended values, based on participant age

Age Range (years)

Recommended Value (bpm)

< 1

60

1 - 2

50

2 - 6

35

6 - 12

25

12 - 18

20

19 - 65

18

65 - 80

28

> 80

30
--motion-filter-order

number of filter coeffecients for the band-stop filter

Censoring and scrubbing options

-r, --head_radius

head radius for computing FD, default is 50mm, 35mm is recommended for baby

-f, --fd-thresh

framewise displacement threshold for censoring, default is 0.2mm

Other options

-w, --work_dir

path where intermediate results should be stored

--clean-workdir

Clears working directory of contents. Use of this flag is notrecommended when running concurrent processes of xcp_d.

--resource-monitor

enable Nipype’s resource monitoring to keep track of memory and CPU usage

--notrack

Opt-out of sending tracking information

Experimental options

--warp-surfaces-native2std

If used, a workflow will be run to warp native-space (fsnative) reconstructed cortical surfaces (surf.gii files) produced by Freesurfer into standard (fsLR) space. These surface files are primarily used for visual quality assessment. By default, this workflow is disabled.

The surface files that are generated by the workflow

Filename

Description

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_pial.surf.gii

The gray matter / pial matter border.

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_smoothwm.surf.gii

The smoothed gray matter / white matter border for the cortex.

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_midthickness.surf.gii

The midpoints between wm and pial surfaces. This is derived from the FreeSurfer graymid (mris_expand with distance=0.5 applied to the WM surfs).

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_inflated.surf.gii

An inflation of the midthickness surface (useful for visualization). This file is only created if the input type is “hcp” or “dcan”.

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_desc-hcp_midthickness.surf.gii

The midpoints between wm and pial surfaces. This is created by averaging the coordinates from the wm and pial surfaces.

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_desc-hcp_inflated.surf.gii

An inflation of the midthickness surface (useful for visualization). This is derived from the HCP midthickness file. This file is only created if the input type is “fmriprep” or “nibabies”.

<source_entities>_space-fsLR_den-32k_hemi-<L|R>_desc-hcp_vinflated.surf.gii

A very-inflated midthicknesss surface (also for visualization). This is derived from the HCP midthickness file. This file is only created if the input type is “fmriprep” or “nibabies”.

see https://xcp-d.readthedocs.io/en/latest/generalworkflow.html

Running xcp_d via Docker containers

If you are running xcp_d locally, we recommend Docker. See Container Technologies for installation instructions.

In order to run Docker smoothly, it is best to prevent permissions issues associated with the root file system. Running Docker as user on the host will ensure the ownership of files written during the container execution.

A Docker container can be created using the following command:

docker run --rm -it \
    -v /fmriprepdata:/in/ \
    -v /tmp/wkdir:/wkdir/ \
    -v /tmp:/scrth/ \
    -v /tmp/xcpd_ciftiF/:/out. \
    pennlinc/xcp_d:latest \
    /in/ /out/ participant \
    --cifti --despike --head_radius 40 -w wkdir --smoothing 6

Running xcp_d via Singularity containers

If you are computing on an HPC, we recommend using Singularity. See Container Technologies for installation instructions.

If the data to be preprocessed is also on the HPC or a personal computer, you are ready to run xcp_d.

singularity run --cleanenv xcp_d.simg \
    path/to/data/fmri_dir  \
    path/to/output/dir \
    --participant-label label

Relevant aspects of the $HOME directory within the container

By default, Singularity will bind the user’s $HOME directory on the host into the /home/$USER directory (or equivalent) in the container. Most of the time, it will also redefine the $HOME environment variable and update it to point to the corresponding mount point in /home/$USER. However, these defaults can be overwritten in your system. It is recommended that you check your settings with your system’s administrator. If your Singularity installation allows it, you can work around the $HOME specification, combining the bind mounts argument (-B) with the home overwrite argument (--home) as follows:

singularity run -B $HOME:/home/xcp \
    --home /home/xcp \
    --cleanenv xcp_d.simg \
    <xcp_d arguments>

Therefore, once a user specifies the container options and the image to be run, the command line options are the same as the bare-metal installation.

Custom Confounds

XCP-D can implement custom confound regression (i.e., denoising). Here, you can supply your confounds, and optionally add these to a confound strategy already supported in XCP-D.

Task Regression

Here we document how to regress task block effects as well as the 36 parameter model confounds.

However, this method is still under development.

Regression of task effects from the BOLD timeseries is performed in 3 steps:

  1. Create a task event timing array

  2. Convolve task events with a gamma-shaped hemodynamic response function (HRF)

  3. Regress out the effects of task via a general linear model implemented with xcp_d

Create a task event array

First, for each condition (i.e., each separate contrast) in your task, create an Nx2 array where N is equal to the number of measurements (volumes) in your task fMRI run. Values in the first array column should increase sequentially by the length of the TR, with the first index = 0. Values in the second array column should equal either 0 or 1; each volume during which the condition/contrast was being tested should = 1, all others should = 0.

For example, for an fMRI task run with 210 measurements and a 3 second TR during which happy faces (events) were presented for 5.5 seconds at time = 36, 54, 90 seconds etc.

[  0.   0.]
[  3.   0.]
[  6.   0.]
[  9.   0.]
[ 12.   0.]
[ 15.   0.]
[ 18.   0.]
[ 21.   0.]
[ 24.   0.]
[ 27.   0.]
[ 30.   0.]
[ 33.   0.]
[ 36.   1.]
[ 39.   1.]
[ 42.   1.]
[ 45.   0.]
[ 48.   0.]
[ 51.   0.]
[ 54.   1.]
[ 57.   1.]
[ 60.   1.]
[ 63.   0.]
[ 66.   0.]
[ 69.   0.]
[ 72.   0.]
[ 75.   0.]
[ 78.   0.]
[ 81.   0.]
[ 84.   0.]
[ 87.   0.]
[ 90.   1.]
[ 93.   1.]
[ 96.   1.]
[ 99.   0.]

Convolve task events with the HRF

Next, the BOLD response to each event is modeled by convolving the task events with a canonical HRF. This can be done by first defining the HRF and then applying it to your task events array with numpy.convolve().

import numpy as np
from scipy.stats import gamma


def hrf(times):
   """Return values for HRF at given times."""
   # Gamma pdf for the peak
   peak_values = gamma.pdf(times, 6)
   # Gamma pdf for the undershoot
   undershoot_values = gamma.pdf(times, 12)
   # Combine them
   values = peak_values - 0.35 * undershoot_values
   # Scale max to 0.6
   return values / np.max(values) * 0.6

# Compute HRF with the signal
hrf_times = np.arange(0, 35, TR)  # TR = repetition time, in seconds
hrf_signal = hrf(hrf_times)
N = len(hrf_signal)-1
tt=np.convolve(taskevents[:,1], hrf_signal)  # taskevents = the array created in the prior step
realt=tt[:-N]  # realt = the output we need!

The code block above contains the following user-defined variables:

  • TR: a variable equal to the repetition time

  • taskevents: the Nx2 array created in the prior step

The code block above produces the numpy array realt, which must be saved to a file named ${subid}_${sesid}_task-${taskid}_desc-custom_timeseries.tsv. This tsv file will be used in the next step of xcp_d.

If you have multiple conditions/contrasts per task, steps 1 and 2 must be repeated for each such that you generate one taskevents Nx2 array per condition, and one corresponding realt numpy array. The realt outputs must all be combined into one space-delimited ${subid}_${sesid}_task-${taskname}_desc-custom_timeseries.tsv file. A task with 5 conditions (e.g., happy, angry, sad, fearful, and neutral faces) will have 5 columns in the custom .tsv file. Multiple realt outputs can be combined by modifying the example code below.

import pandas as pd

# Create an empty task array to save realt outputs to
taskarray = np.empty(shape=(measurements, 0))  # measurements = the number of fMRI volumes

# Create a taskevents file for each condition and convolve with the HRF, using the code above
## code to compute realt

# Write a combined custom.tsv file
taskarray = np.column_stack((taskarray, realt))
df = pd.DataFrame(taskarray)
df.to_csv(
   "{0}_{1}_task-{2}_desc-custom_timeseries.tsv".format(subid, sesid, taskid),
   index=False,
   header=False,
   sep='\t',
)

The space-delimited *desc-custom_timeseries.tsv file for a 5 condition task may look like:

0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.3957422940438729
0.0 0.0 0.0 0.0 0.9957422940438729
0.0 0.0 0.0 0.0 1.1009022019820307
0.0 0.0 0.0 0.0 0.5979640661963432
0.0 0.0 0.0 0.0 0.31017195439257517
0.0 0.0 0.0 0.0 0.7722398821320118
0.0 0.0 0.0 0.0 0.9755486196351566
0.0 0.0 0.0 0.0 0.9499183578181378
0.0 0.0 0.0 0.0 0.8987971115721047
0.0 0.0 0.0 0.0 0.8750149365335346
0.0 0.0 0.0 0.0 0.47218635162456457
0.0 0.0 0.0 0.0 -0.1294234695774829
0.3957422940438729 0.0 0.0 0.0 -0.23488535934344593
0.9957422940438729 0.0 0.0 0.0 -0.12773843588350925
1.1009022019820307 0.0 0.0 0.0 -0.04421213464698274
0.5979640661963432 0.0 0.0 0.0 -0.011439970324577234
-0.08557033965129775 0.0 0.0 0.0 -0.0023929581477495315
-0.22350241191186113 0.0 0.0 0.0 -0.00042413222490445587
0.27038871169699874 0.0 0.0 0.0 -6.512750410688139e-05
0.9519542916217947 0.0 0.0 0.0 -8.104611114467545e-06
1.0895273591615604 0.0 0.0 0.3957422940438729 0.0
0.5955792126597081 0.0 0.0 0.9957422940438729 0.0
-0.0859944718762022 0.0 0.0 1.1009022019820307 0.0
-0.223567539415968 0.0 0.0 0.5979640661963432 0.0
-0.12536168695798863 0.0 0.0 -0.08557033965129775 0.0
-0.04378800242207828 0.0 0.0 -0.22350241191186113 0.0
-0.011374842820470351 0.3957422940438729 0.0 -0.12535358234687416 0.0
-0.002384853536635064 0.9957422940438729 0.0 -0.04378800242207828 0.0
-0.00042413222490445587 1.1009022019820307 0.3957422940438729 -0.011374842820470351 0.0
-6.512750410688139e-05 0.5979640661963432 0.9957422940438729 -0.002384853536635064 0.0
0.39573418943275845 -0.08557033965129775 1.1009022019820307 -0.00042413222490445587 0.0
0.9957422940438729 -0.22350241191186113 0.5979640661963432 -6.512750410688139e-05 0.0
1.1009022019820307 -0.12535358234687416 -0.08557033965129775 -8.104611114467545e-06 0.0
0.5979640661963432 -0.04378800242207828 -0.22350241191186113 0.0 0.0

Command Line XCP-D with Custom Confounds

Last, supply the ${subid}_${sesid}_task-${taskid}_desc-custom_timeseries.tsv file to xcp_d with -c option. -c should point to the directory where this file exists, rather than to the file itself; xcp_d will identify the correct file based on the subid, sesid, and taskid. You can simultaneously perform additional confound regression by including, for example, -p 36P to the call.

singularity run --cleanenv -B /my/project/directory:/mnt xcpabcd_latest.simg \
   /mnt/input/fmriprep \
   /mnt/output/directory \
   participant \
   --despike --lower-bpf 0.01 --upper-bpf 0.08 \
   --participant_label $subid -p 36P -f 10 \
   -t emotionid -c /mnt/taskarray_file_dir

Custom Parcellations

While XCP-D comes with many built in parcellations, we understand that many users will want to use custom parcellations. We suggest running XCP-D with the -cifti option (assuming you have cifti files), and then using the Human Connectome Project wb_command to generate the time series:

wb_command \
   -cifti-parcellate \
   {SUB}_ses-{SESSION}_task-{TASK}_run-{RUN}_space-fsLR_den-91k_desc-residual_bold.dtseries.nii \
   your_parcels.dlabel \
   {SUB}_ses-{SESSION}_task-{TASK}_run-{RUN}_space-fsLR_den-91k_desc-residual_bold.ptseries.nii

After this, if one wishes to have a connectivity matrix:

wb_command \
   -cifti-correlation \
   {SUB}_ses-{SESSION}_task-{TASK}_run-{RUN}_space-fsLR_den-91k_desc-residual_bold.ptseries.nii \
   {SUB}_ses-{SESSION}_task-{TASK}_run-{RUN}_space-fsLR_den-91k_desc-residual_bold.pconn.nii

More information can be found at the HCP documentation

Troubleshooting

Logs and crashfiles are outputted into the <output dir>/xcp_d/sub-<participant_label>/log directory. Information on how to customize and understand these files can be found on the nipype debugging page.

Support and communication. The documentation of this project is found here: https://xcp-d.readthedocs.io/.

All bugs, concerns and enhancement requests for this software can be submitted here: https://github.com/PennLINC/xcp_d/issues.

If you have a question about using xcp_d, please create a new topic on NeuroStars with the “xcp_d” tag. The xcp_d developers follow NeuroStars, and will be able to answer your question there.