Friday, July 1, 2022

fMRIPrep for a robust preprocessing pipeline

When dealing with neuroimaging studies, researchers often create their own preprocessing workflows that suit the individual labs and datasets. The existence of a variety of MRI analytical tools does not help with the standardization of the analysis. Recently, a group at Stanford University led by Russ Poldrack et al. came up with a robust and analysis-agnostic tool for robust and reproducible preprocessing called fMRIPrep. The software tool integrates existing software packages such as FSL, ANTs, AFNI, Nypipe, and SPM. It also works on the BIDS dataset, another useful framework to standardize neuroimaging data collection following the spirit of Open Science.

Preprocessing anatomical images

  1. Correction of intensity non-uniformity with N4BiasFieldCorrection (ANTs).
  2. Then skull-stripped the image with antsBrainExtraction.sh script (ANTs).
  3. Run recon-all (FreeSurfer) to reconstruct 2D cortical surface from the 3D T1-weighted image.
  4. Cortical grey matter segmentation is done using FreeSurfer.
  5. Do non-linear spatial normalization to the standard template* using antsRegistration (ANTs).
  6. Brain tissues: CSF, white matter/WM, and grey matter/GM are extracted using FAST (FSL).
* standard space follows the ICMB 152 nonlinear asymmetrical template (MNI, v2009c). You can also choose the slightly older version, MNI152 asymmetrical template 6th generation. 

Preprocessing functional images

  1. For each BOLD run found in the folder, we first find its reference volume, which is usually the 3D volume at midpoint (half-way through the run); and it should be skull-stripped also.
  2. Head-motion parameters are estimated using MCFLIRT (FSL).
  3. Slice timing correction is performed using 3dTshift (AFNI) if slice-time information is available.
  4. Fieldmap distortion correction is then applied for SDC around the air space.
  5. Co-registration into the subject's T1 image is done using FLIRT with BBS and 6 d.o.f (FSL).
  6. If surface reconstruction is selected, then bbresgiter (FreeSurfer) is applied for co-registration. The authors reported that this method yields the best results owing to the accuracy of GM/WM surfaces driving the process.
  7. A series of transforms of BOLD image up to the co-registration can be customized using antsApplyTransforms and Lanczos interpolation.
  8. BOLD images can be normalized to the MNI template of a different resolution: 1mm or 2mm, and to a different template version.
  9. Lastly, if ICA-AROMA* is enabled, then the pipeline will perform ICA analysis and do some denoising. At this stage, preprocessed data up to Step-7 has to be "smoothened" first before feeding it into FSL-MELODIC. 

* By default, fmriPrep will produce a non-aggresive denoised and preprocessed BOLD image in the standard space at the end of ICA-AROMA. The authors suggest that non-aggresive method of denoising is recommended as it removes less signal components that share variance with the nuisance regressors.  


Extraction of nuisance regressors

By default, there is no temporal denoising done in fMRIPrep but it produces a vast variety of nuisance regressors. This allows researchers to be more flexible in their denoising strategies.

  1. Component-based correction (CompCor) provides physiological noise regressors.
  2. PCA on high-pass filtered data is performed to yield tCompCor (top 5% of the variable voxels in the subcortical masks without any GM components) and aCompCor (intersection between mask and CSF/WM).
  3. Three nuisance regressors from CSF, WM, and whole brain signals.
  4. Framewise displacements and their first derivatives to signify head movements are also obtained (DVARS).
  5. Components produced by ICA-AROMA are also included.

Wednesday, March 30, 2022

Linear models: Getting the syntax right!

When conducting statistical analysis and mixed-effects modelling, I often either got confused or forgot how to write the syntax properly. I thought this entry will be useful. 

Simple linear models
In simple linear models in R using the lm library, we can write the formula as DV ~ IV (DV is a dependent variable). The command employs a few arithmetic symbols. A '+ sign' indicates more than one main effect or predictor (independent variable, IV) with no interaction. A '* sign' provides a short form of main effects with interaction. Number '1' with a '+ sign' refers to an intercept. There is no need to explicitly write this, as it is always implied and estimated by default. Sometimes, you want to omit this by simply writing '0'.


When we want to include fixed/random factors, we then use, e.g. the lmer library, and the syntax is slightly different.

DV ~ 1 + IV1 * IV2        DV ~ IV | grouping                 
  1. The dependent variable (DV) is the response variable to be predicted. The independent variable (IV) is the one whose effects we will assess. By default, IV is treated as a fixed factor.
  2. A vertical bar denotes a grouping factor or random factor. It separates expressions for design matrices from the grouping factor. For example: to fit a predictor for each random factor, you can write 1 + A|S, which means the same thing as (A|S).
  3. A '/ sign' indicates nesting. So (school/class) means classes are nested within school.
  4. Fixed factors can be included without any grouping. You can have additional random factors without any fixed factor (an intercept-only model). In each random/fixed factor, ask yourself whether you allow the intercept to change, the slope, or both!


    Hierarchical Linear Modeling
    Begin with the first and simplest model which is the intercept-only model. It has no predictor or independent variable (so, no slope!). Suppose we have one random factor, usually the participant factor. The model is also called the unconditional model.

    Y ~ (1|S)            ; intercept-only model 

    The next model is the random intercept model. We assign a new variable or predictor X as Level 1, which acts as a fixed factor. Here, the intercepts for different subjects will vary but not the slopes. Note that the intercept for X can be omitted as the function understands it to be 1+X. 

    Y ~ X + (1|S)        ; random intercept model

    Now, we can add complexity by allowing different subjects to have different intercepts and slopes. This is the random intercept + slope model. There are two possibilities: the variations of intercept and slope can be independent or correlated!

    Y ~ X + (1|S) + (0 + X|S)    ; independent
    Y ~ X + (X|S)                ; correlated

     Source: here