rsDenoise

A python-based, customisable and extensible pipeline for resting-state fMRI data.
Authors: Paola Galdi ([email protected]) and Julien Dubois ([email protected]).
The software is provided as is without warranty of any kind.

This configurable pipeline for the processing of resting-state fMRI was originally designed to work with data from the Human Connectome Project (and hence the HCP directory structure). The original code is still available here. The pipeline was first updated to work with data that are converted to an HCP-like format, using ciftify, and later extended to work with data minimally preprocessed with fmriprep.

The pipeline works with both volumetric (NIFTI) and surface data (either CIFTI or GIFTI).

All functions needed for processing are grouped in three helpers file:

• HCP_helpers.py contains the set of functions to work with the HCP directory structure
• fmriprepciftify_helpers.py is meant to work with data previously processed with fmriprep and freesurfer, converted to HCP-like structure with ciftify
• fmriprep_helpers.py is designed to work directly with fmriprep outputs, with or without freesurfer outputs

The IPython notebooks provide a few examples of use (more can be found here: HCP_MRI-behavior) A minimal script to launch the pipeline is provided in runPipeline_example.py. A command line wrapper is provided in rsDenoise.py. Example output for a HCP subject here.

References

• Dubois J, Galdi P, Han Y, Paul LK, Adolphs R. Resting-state functional brain connectivity best predicts the personality dimension of openness to experience. Personality neuroscience. 2018 Jul;1:e6.
• Dubois J, Galdi P, Paul LK, Adolphs R. A distributed brain network predicts general intelligence from resting-state human neuroimaging data. Philosophical Transactions of the Royal Society B: Biological Sciences. 2018 Sep 26;373(1756):20170284.
• Kliemann D, Adolphs R, Armstrong T, Galdi P, Kahn DA, Rusch T, Enkavi AZ, Liang D, Lograsso S, Zhu W, Yu R. Caltech Conte Center, a multimodal data resource for exploring social cognition and decision-making. Scientific Data. 2022 Mar 31;9(1):138.

Overview

Prerequisites

The code depends on a set of python libraries, including: scipy, pandas, matplotlib, statsmodels, numpy, nibabel, nilearn, nipype, seaborn, scikit_learn.

Preprocessing Pipeline

Preprocessing is performed by executing a sequence of steps or operations, each of which corresponds to a function in the helpers file. There are 7 available operations: Voxel Normalization, Detrending, Motion Regression, Scrubbing, Tissue Regression, Global Signal Regression and Temporal Filtering. The corresponding functions are built using the following template:

def OperationName(niiImg, flavor, masks, imgInfo):

• niiImg it is a list of two elements [data,volData]. If the file to be processed is a Nifti file, the variable data contains Nifti data and the variable volData is None. If the input file is a Cifti file, the variable data contains Cifti data and the variable volData contains volumetric data. Image data can be loaded with function load_img().
• flavor is a list containing the Operation parameters (detailed in the following section).
• masks is the output of the function makeTissueMasks(), a list of four elements corresponding to 1) whole brain mask, 2) white matter mask, 3) cerebrospinal fluid mask and 4) gray matter mask.
• imgInfo is a list of 6 elements corresponding to 1) no. of rows in a slice, 2) no. of columns in a slice, 3) no. of slices, 4) no. of time points, 5) the affine matrix and 6) repetion time. These data can be retrieved with function load_img().

A preprocessing pipeline can be fully described with the following data structure:

[
    ['Operation1',  1, [param1]],
    ['Operation2',  2, [param1, param2, param3]],
    ['Operation3',  3, [param1, param2]],
    ['Operation4',  4, [param1]],
    ['Operation5',  5, [param1, param2]],
    ['Operation6',  6, [param1, param2, param3]],
    ['Operation7',  7, [param]]
]

It is a list of structures. Each structure is a list of 3 elements:

1. The operation name.
2. The operation order.
3. The list of parameters for the specific operation.

Operations are executed following the operation order specified by the user. Note that in case of regression operations, a single regression step combining multiple regressors (e.g., tissue regressors and global signal regressor) can be requested by assigning the same order to all regression steps.

There are several pipelines already specified helpers file:

config.operationDict = {
    'A': [ #Finn et al. 2015
        ['VoxelNormalization',      1, ['zscore']],
        ['Detrending',              2, ['legendre', 3, 'WMCSF']],
        ['TissueRegression',        3, ['WMCSF', 'GM']],
        ['MotionRegression',        4, ['R dR']],
        ['TemporalFiltering',       5, ['Gaussian', 1]],
        ['Detrending',              6, ['legendre', 3 ,'GM']],
        ['GlobalSignalRegression',  7, ['GS']]
        ],
    'B': [ #Satterthwaite et al. 2013 (Ciric7)
        ['VoxelNormalization',      1, ['demean']],
        ['Detrending',              2, ['poly', 2, 'wholebrain']],
        ['TemporalFiltering',       3, ['Butter', 0.01, 0.08]],
        ['MotionRegression',        4, ['R dR R^2 dR^2']],
        ['TissueRegression',        4, ['WMCSF+dt+sq', 'wholebrain']],
        ['GlobalSignalRegression',  4, ['GS+dt+sq']],
        ['Scrubbing',               4, ['RMS', 0.25]]
        ],
    'C': [ #Siegel et al. 2016 (SiegelB)
        ['VoxelNormalization',      1, ['demean']],
        ['Detrending',              2, ['poly', 1, 'wholebrain']],
        ['TissueRegression',        3, ['CompCor', 5, 'WMCSF', 'wholebrain']],
        ['TissueRegression',        3, ['GM', 'wholebrain']], 
        ['GlobalSignalRegression',  3, ['GS']],
        ['MotionRegression',        3, ['censoring']],
        ['Scrubbing',               3, ['FD+DVARS', 0.25, 5]], 
        ['TemporalFiltering',       4, ['Butter', 0.009, 0.08]]
        ],
        ...
    }

Pipeline Operations

Custom operations can be provided by users, using the template detailed above. New functions need to be added to the Hooks data structure in the helpers file:

Hooks={
   'MotionRegression'       : MotionRegression,
   'Scrubbing'              : Scrubbing,
   'TissueRegression'       : TissueRegression,
   'Detrending'             : Detrending,
   'TemporalFiltering'      : TemporalFiltering, 
   'GlobalSignalRegression' : GlobalSignalRegression, 
   'VoxelNormalization'     : VoxelNormalization,
   'newFunction'            : newFunction,
   }

In the following we detail those included in the helpers file.

Voxel Normalization

• zscore: convert each voxel’s time course to z-score (remove mean, divide by standard deviation).
• demean: substract the mean from each voxel's time series.
• pcSigCh: convert each voxel’s time course to % signal change (remove mean, divide by mean, multiply by 100)

Example: ['VoxelNormalization', 1, ['zscore']]

Detrending

• poly: polynomial regressors up to specified order.
  1. Specify polynomial order.
  2. Specify tissue, one of 'WMCSF' or 'GM'.
• legendre: Legendre polynomials up to specified order.
  1. Specify polynomial order
  2. Specify tissue, one of 'WMCSF' or 'GM'.

Example: ['Detrending', 2, ['poly', 3, 'GM']]

Motion Regression

Note: R = [X Y Z pitch yaw roll]
Note: if temporal filtering has already been performed, the motion regressors are filtered too.
Note: if scrubbing has been requested, a regressor is added for each volume to be censored; the censoring option performs only scrubbing.

• R: translational and rotational realignment parameters (R) are used as explanatory variables in motion regression.
• R dR: translational and rotational realignment parameters (R) and their temporal derivatives (dR) are used as explanatory variables in motion regression.
• R R^2: translational and rotational realignment parameters (R) and their quadratic terms (R^2) are used as explanatory variables in motion regression.
• R dR R^2 dR^2: realignment parameters (R) with their derivatives (dR), quadratic terms (R^2) and square of derivatives (dR^2) are used in motion regression.
• R R^2 R-1 R-1^2: realignment parameters at time t (R) and realignment parameters at time t-1 (R-1) with their square terms are used in motion regression.
• R R^2 R-1 R-1^2 R-2 R-2^2: realignment parameters at time t (R), t-1 (R-1) and t-2 (R-2) with their square terms are used in motion regression.
• censoring: for each volume tagged for scrubbing, a unit impulse function with a value of 1 at that time point and 0 elsewhere is included as a regressor.
• ICA-AROMA: ICA is used to identify noise components that are used as regressors.

Example: ['MotionRegression', 3, ['R dR']]

Scrubbing

Note: this step only flags the volumes to be censored, that are then regressed out in the MotionRegression step.
Note: uncensored segments of data lasting fewer than 5 contiguous volumes, are flagged for removal as well.

• FD
  1. Specify a threshold for framewise displacement (FD) in mm.
  2. Specify number of adjacent volumes to exclude (optional).
• DVARS
  1. Specify a threshold for variance of differentiated signal (DVARS).
  2. Specify number of adjacent volumes to exclude (optional).
• FD+DVARS
  1. Specify a threshold for framewise displacement (FD) in mm.
  2. Specify a threshold t s.t. volumes with a variance of differentiated signal (DVARS) greater than (100 + t)% of the run median DVARS are flagged for removal (as in Siegel et al, Cerebral Cortex, 2016).
  3. Specify number of adjacent volumes to exclude (optional).
• RMS
  1. Specify threshold for root mean square displacement in mm.
  2. Specify number of adjacent volumes to exclude (optional).

Example: ['Scrubbing', 4, ['FD+DVARS', 0.25, 5]]

Tissue Regression

• GM: the gray matter signal is added as a regressor.
  1. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain')
• WM: the white matter signal is added as a regressor.
  1. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain')
• WMCSF: white matter and cerebrospinal fluid signals are added as regressors.
  1. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain')
• WMCSF+dt: white matter and cerebrospinal fluid signals with their derivatives are added as regressors.
  1. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain')
• WMCSF+dt+sq: white matter and cerebrospinal fluid signals with their derivatives and quadratic terms are added as regressors.
  1. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain')
• CompCor: a PCA-based method (Behzadi et al., 2007) is used to derive N components from CSF and WM signals.
  1. Specify no. of components to compute for specified tissue mask (see following parameter).
  2. Specify if components should be computed using a single mask for white matter and cerebrospinal fluid ('WMCSF') or separatedly for white matter and cerebrospinal fluid ('WM+CSF')
  3. Specify if regression should be performed on gray matter signal ('GM') or whole brain signal ('wholebrain').

Example: ['TissueRegression', 5, ['CompCor', 5, 'WMCSF', 'wholebrain']]

Global Signal Regression

• GS: the global signal is added as a regressor.
• GS+dt+sq: the global signal with its derivatives are added as regressors.
• GS+dt+sq: the global signal with its derivatives and square term are added as regressors.

Example: ['GlobalSignalRegression', 6, ['GS+dt+sq']]

Temporal Filtering

Note: if scrubbing has been requested, censored volumes are replaced by linear interpolation.

• Butter: Butterworth band pass filtering.
  1. Specify high pass threshold.
  2. Specify low pass threshold.
• Gaussian: Low pass Gaussian smoothing.
  1. Specify standard deviation.
• DCT: The discrete cosine transform is used to compute regressors to perform temporal filtering.
  1. Specify high pass threshold.
  2. Specify low pass threshold.

Example: ['TemporalFiltering', 7, ['Butter', 0.009, 0.08]]

Advanced functionalities

• QC plots: FD trace and grayplot
• XML reports
• Basic support for sge
• Computation of functional connectivy matrices, given a parcellation, or seed-connectivity matrices, given a seed
• Prediction of individual scores from functional connectivity matrices