Preprocessing for atlas based connectivity matrices

1

Hi,

I just started using fmriprep and found it a great tool to preprocess resting state fMRI data. I am interested in using it to compute atlas-based connectivity matrices (conmats). Would the approach suggested below be adequate?

  1. Run fmriprep with T1w as output space (since conmats are computed in subject space):
    fmriprep-docker /raw_data /preproc participant --nthreads 10 --output-space T1w --use-aroma --work-dir /tmp_data --write-graph --config /raw_data/nipype.cfg

  2. Do non-agressive denoising, as suggested in fmriprep documentation:
    fsl_regfilt -i sub-<subject_label>_task-<task_id>_bold_space-T1w_preproc.nii.gz -f $(cat sub-<subject_label>_task-<task_id>_bold_AROMAnoiseICs.csv) -d sub-<subject_label>_task-<task_id>_bold_MELODICmix.tsv -o sub-<subject_label>_task-<task_id>_bold_space-<space>_AromaNonAggressiveDenoised.nii.gz

  3. Highpass / lowpass temporal filters + scrubbing (any particular default suggestions?)

  4. Mean of voxels in each ROI of the aparc+aseg atlas obtained in subject native space (/preproc//freesurfer/sub-<subject_label>/mri/aparc+aseg.mgz). Is there any transform to apply to that file to match sub-<subject_label>_task-<task_id>_bold_space-<space>_AromaNonAggressiveDenoised.nii.gz obtained in the previous step?

  5. (Partial-) Correlation of the resulting mean ROI signals to obtain the connectivity matrix.

As an alternative, would you suggest to do the calculations in surface+volume space (GII format using the same aparc+aseg atlas)? In this case, how would it change the above workflow and associated commands?

Any help would be greatly appreciated.

2

With regard to the FMRIPREP-related portions of this pipeline:

  1. ICA-AROMA is calculated in the MNI template space, which we do not calculate (to save time and disk space) if no MNI outputs are requested, so you’re going to need to use --output-space T1w template (output spaces are separated by a space, not a comma).

  2. That looks reasonable. I haven’t tried to denoise results in the subject native space, so I can’t attest that that works. Please let us know if you need to make any modifications, and we’ll updated the documentation.

  3. I don’t have experience here, so I’ll leave commentary to others.

  4. While we work to keep FreeSurfer and FMRIPREP outputs aligned, they aren’t always sampled exactly the same. For that you’re going to need to apply the FreeSurfer -> T1w transform. It’s not currently an output of fmriprep (though it probably ought to be), but you can create it with:

tkregister2 --noedit --mov /preproc/freesurfer/sub-<subject_label>/mri/T1.mgz \
    --targ /preproc/fmriprep/sub-<subject_label>/anat/sub-<subject_label>_T1w_preproc.nii.gz \
    --regheader --ltaout fs_to_t1w.lta

And then apply it:

mri_vol2vol --mov /preproc/freesurfer/sub-<subject_label>/mri/aparc+aseg.mgz \
    --targ /preproc/fmriprep/sub-<subject_label>/anat/sub-<subject_label>_T1w_preproc.nii.gz \
    --lta fs_to_t1w.lta --o aparc+aseg.nii.gz --nearest

Hopefully others can advise you more on the connectivity analysis-specific portions of your question.

1 Like
3

RE:
3. There was a good discussion of this on the fsl listserv, and would direct you to some papers. In brief, the predominant view is 0.01-0.08 Hz is the range of interest, but it depends on your TR. Additionally, see Beckman’s comment: https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind1203&L=fsl&D=0&P=231075

2 Likes
4

Dear @effigies, thank you very much for this answer, this is very useful. I will apply the steps including the tkregister2 and mri_vol2vol commands as you suggested. After some research, I found the pipeline suggestion below, and it would be great to know how to adapt fmriprep in order to implement it. Please let me know if a new post or an “issue” on github would be more appropriate.

Many thanks also @jdkent. As mentioned above, after a few reads, it seems the pipeline below is a good choice in general, and it would be great to know how to use fmriprep in order to implement it.

A manuscript uploaded a couple months ago on bioRxiv evaluated 16 different pipelines for noise reduction, and as a result recommended the following pipeline:
1. removal of the first four volumes of each acquisition
2. slicetime correction (implemented in SPM8)
3. despiking (using AFNI’s 3dDespike)
4. realignment (two-pass realignment: all volumes to the first volume, then to the mean volume, using SPM8)
5. co-registration of EPI data to the native, cropped, high-resolution structural image (via rigid-body registration using ANTs)
6. registration to MNI (application of the nonlinear transform T1 --> MNI to the T1-coregistered EPI data, using ANTs)
7. linear detrending of the spatially normalized BOLD time series (using rest_detrend function from REST toolbox)
8. intensity normalisation of the EPI data to mode 1000 units
9. spatial smoothing (with an 8mm FWHM kernel)
10. regression of ICA-AROMA noise ICs as well as averaged WM and CSF signals
11. bandpass filtering between 0.008 and 0.08 Hz
12. compute tDOF-loss to use as a covariate in group comparisons

It seems that the steps in bold (1, 3, 7, 8, 9, 10, 11) are not fully integrated in fmriprep (please correct me if I am wrong). In particular:

  • Step 1 would not be difficult to do manually (just remove the first 4 volumes before using fmriprep)
  • Step 3 and 7 seem not to be an option in fmriprep
  • Step 8 does not seem to be part of fmriprep either, and I am not sure of its significance/usefulness
  • For Step 9, it is not clear if and when smoothing is performed in fmriprep
  • Step 10 seems to require some modifications to be used in fmriprep:
    • I could see the mean signal of WM but not of CSF in the confounds file (sub-<subject>_task-rest_bold_confounds.tsv) although the CSF component seems to be computed in the confounds workflow
    • A regression would need to be performed using only the WM, CSF and ICA-AROMA columns of the sub-<subject>_task-rest_bold_confounds.tsv file. Is there a staightforward command for performing this step?
  • Since Step 11 and step 12 are last, they could also easily be performed independently of fmriprep, especially considering they are arguably not part of fMRI preprocessing (and Step 12 seems to be obtainable from the mriqc sister tool)

Would there be any suggestion on how to approach implementing this reference pipeline as of now and before fmriprep being modified in the future? (I’d be happy to contribute to creating this future fmriprep version)

5

Thank you for bringing this topic up. While designing FMRIPREP we thought hard on how to make it most usable for different types of analyses and any feedback is highly appreciated.

As you know there are many, many papers providing different recommendations on cleaning up data. Instead of implementing one default approach (or even a small set) we opted for generating outputs that would allow easy application of a wide range of different approaches. We are not married to this approach and happy to hear what do you think.

So to implement the closest version of this particular preprocessing you would have to do the following:

  1. Run FMRIPREP with --use-aroma option
  2. Take the *bold_space-MNI152NLin2009cAsym_preproc.nii.gz outputs.
  3. Either discard 4 first volumes, or as many volumes as the FMRIPREP steady non-state detector recommends as indicated by the number of NonSteadyStateOutlier. Alternatively you can use these columns in the denoising design matrix (see below).
  4. Smooth the data with 8mm (or whatever is appropriate for your data). Alternatively you can opt for not smoothing the 4D data and smoothing only the first level statistical maps prior to doing group level.
  5. Create a denoising design matrix that would include the following columns from the _confounds.tsv file: AROMAAggrComp* (those are the noise components from ICA AROMA), WhiteMatter, and aCompCor01 (first component of the anatomical compcor is a good approximation of the mean signal in WM and CSF). In addition to achieve the linear detrending you should add a column of ascending numbers (1,2,3,4 etc).
  6. Fit this model to the data and save the residuals - this is your denoised data.
  7. Save the number of columns for each run - this is your tDOF-loss you will use on group level.
  8. Perform band pass filtering (remember to do it after nuisance regression!).

In this solution I skipped intensity normalization since I don’t think it influences temporal noise and is more of convenience scaling operation.

One important thing that is missing in my recipe is despiking. As far as I know this cannot be easily done with nuissance regression since 3dDespike changes each voxel differently (so we cannot generate a single column for all voxels to include in the _confounds.tsv file). We could add 3dDespike as an optional (turned on by default or not) step in FMRIPREP. Would that be desirable? Do you know of any work recommending against using 3dDespike?

PS Most outputs of FMRIPREP are not smoothed. the only exception is variant-smoothAROMAnonaggr which is the non-aggressively ICA AROMA denoised data which is smoothed using 6mm kernel (as in the original AROMA paper).
PPS The reason why we output only WM but not CSF is that for participants with small ventricles we cannot reliably estimate CSF ROIs without running into partial volume issues. I would love to hear suggestions how to work around this.

3 Likes
6

@michael please check out https://github.com/poldracklab/fmriprep/issues/699

You’ll find a reference to a recent paper talking about ordering of head-motion estimation and corrections with time-interpolations (slice timing correction and despiking).

7

Glad to find this discussion here. Recently I worked on a analysis pipeline using FMRIPREP. And here’s what I did in preprocess after FMRIPREP for later modelling.

  1. Grab *bold_space-MNI152NLin2009cAsym_preproc.nii.gz outputs and normalized its median intensity to 10000 (follow FSL’s convention) .
  2. Use SUSAN smoothed the data with x mm FWHM.
  3. Optional Applied non-aggressive ICA-AROMA denoise as FMRIPREP document mentioned.
  4. Optional Use fsl_glm to regress out nuisance like WM from _confounds.tsv and taken the residuals images.
  5. Applied temporal filtering use fslmath -bptf (I think if you applied a highpass filter, it will taken care of the linear detrend.)
  6. Add back the mean image derived from 2 or 3 to make the filtered image a firm dataset.

I have several question regrading to this pipeline.

  1. As you mentioned,

Create a denoising design matrix that would include the following columns from the _confounds.tsv file: AROMAAggrComp* (those are the noise components from ICA AROMA), WhiteMatter, and aCompCor01 (first component of the anatomical compcor is a good approximation of the mean signal in WM and CSF). In addition to achieve the linear detrending you should add a column of ascending numbers (1,2,3,4 etc).

should I combined 3 and 4 into one step use fsl_glm? Is that an aggressive way to remove nuisance variables?

  1. I’m aware that the in FMRIPREP, ICA-AROMA was calculated on smoothed data with 6mm FWHM. Could I apply the derived regressors on smoothed data with different FWHM (or unsmoothed data for MVPA analysis)?

  2. Could I use these ICA-AROMA regressors on data in T1w space? Or some modification are needed?

Another question not quite related to this thread. If it’s inappropriate here I could open a new thread later.

I obtained several 7T MRI data as a benchmark for later study. For structure image, I got a MP2RAGE image with 0.7 mm resolution. After some tweaks, I still couldn’t get a good skullstrip result from FMRIPREP due to the image contains noise voxels around the brain with similar intensity. Is it possible to modifies anat workflow to handle MP2RAGE image?
Also, the partial FOV bold image corregistration was failed in this dataset (maybe due to problematic anat image). Should I scan a whole brain reference EPI image with same parameters (_sbref ?) to help registration? Could FMRIPREP support it?

Thanks.

8

Combining these steps does not make the denoising more or less aggressive, but will reduce the number of intermediate files you create.

In theory I don’t see anything wrong with that, but I have not looked into this issue too much.

Again - I don’t see anything wrong with this theoretically.

We would love to add support for MP2RAGE into FMRIPREP. What do you normally use for skullstripping?

We currently do not use sbref images, but we have seen successful coregistration with other FOV datasets. I would focus on issues with anat first.

Some more general questions:

  1. Why are you doing the median 10000 normalization? Would you prefer FMRIPREP to output such data by default?
  2. It seems that everyone is doing some sort of smoothing. Would it make sense to add an option in FMRIPREP to output smoothed data?
9

I used fsl_regfilt to remove ICA components. As far as I know, it will include all regressor in MELODICmix.tsv and subtract the effect of noise regressors you given from input data. So I could plug nuisance regressors into this file and specify its column number in fsl_regfilt to do a one step denoise. Is that correct?

For MP2RAGE image, I tried BET and AFNI 3dSkullStrip both on origin image. None of them given good result. Then I followed this note https://www.evernote.com/shard/s262/sh/11bd5998-9c8a-41a1-a452-c1c4e4876a4f/f44341b130b7ded5 create UNIxINV2 image (with/without N4BiasFieldCorrection) and do skullstrip on it. AFNI works ok but some gyrus in parietal cortex was missing. Bet still works poorly. For now I use AFNI with manually adjusted brainmask.

It just follow the FSL’s FEAT convention to make the preprocessed data looks like my old data. I have no idea why FSL do that.

I would love this feature. For task fMRI activation analysis, smoothing is necessary. And I have a general sense that applies smoothing on ‘raw’ data is better than smoothing a first-level/participant-level statistical map. For MVPA/RSA analysis, people usually used unsmoothed data (Although some paper says a little smooth actually have a little benefit for classification).
So I’m thinking of an --smooth-fwhm argument that could take a list fwhm for smoothing. For example, if I specify --smooth-fwhm 0.0 6.0, the output will contain a smoothed data for activation analysis and a unsmoothed data for MVPA analysis (also with ica denoised version). Neat!

10

Sorry - I misread your original message. To do non-agressive denoising two separate regressions are required so those steps cannot be combined. If you wanted to do aggressive denoising all confounds could be combined in one design matrix. Mind that FMRIPREP already outputs the non-aggressively denoised variant of the data in MNI space so you don’t have to run it yourself.

PS I’ll be on the lookout for good MP2RAGE skullstripping methods (CBS Tools might be worth checking out).

11

I only do it manually when I use a different smooth FWHM other than 6mm. That’s why I think a smooth argument is super helpful.

I will check out the CBS tools. Thanks for the suggestions.

12

Thanks a lot for these suggestions @ChrisGorgolewski. I am slightly confused by the subsequent answers to @Zhifang. I would like to implement your suggestions in order to obtained denoised data in T1 space at subject level (I have resting state data for which I would like to get the associated connectivity matrices).

  1. Run FMRIPREP with --use-aroma option

In this case, I assume I would need to run --use-aroma --output-space T1w template as suggested earlier by @effigies

  1. Take the *bold_space-MNI152NLin2009cAsym_preproc.nii.gz outputs.
  2. Either discard 4 first volumes, or as many volumes as the FMRIPREP steady non-state detector recommends as indicated by the number of NonSteadyStateOutlier. Alternatively you can use these columns in the denoising design matrix (see below).

I couldn’t find the NonSteadyStateOutlier information. I looked into the sub-<subject>_task-rest_bold_confounds.tsv file but could only find the columns WhiteMatter GlobalSignal stdDVARS non-stdDVARS vx-wisestdDVARS FramewiseDisplacement tCompCor0* aCompCor* X Y Z RotX RotY RotZ AROMAAggrComp0*

  1. Smooth the data with 8mm (or whatever is appropriate for your data). Alternatively you can opt for not smoothing the 4D data and smoothing only the first level statistical maps prior to doing group level.

From the subsequent posts related to the questions of @Zhifang it seems that a 6mm smoothing kernel is applied before the ICA-AROMA calculations. Are these smoothed data only temporary and discarded? That is, the 8mm smoothing kernel is to be applied to sub-<subject_label>_task-<task_id>_bold_space-T1w_preproc.nii.gz which is unsmoothed, is that correct?

  1. Create a denoising design matrix that would include the following columns from the _confounds.tsv file: AROMAAggrComp* (those are the noise components from ICA AROMA), WhiteMatter, and aCompCor01 (first component of the anatomical compcor is a good approximation of the mean signal in WM and CSF). In addition to achieve the linear detrending you should add a column of ascending numbers (1,2,3,4 etc).

I have a few questions regarding this point:

  • Is the associated command as below (confounds file transposed to be row-wise for fsl_regfit ):
    fsl_regfilt -i sub-<subject_label>_task-<task_id>_bold_space-T1w_preproc_smooth8mm.nii.gz -f <col nums for WhiteMatter, aCompCor01 and AROMAAggrComp*> -d sub-<subject_label>_task-<task_id>_bold_confounds_transpose.tsv -o sub-<subject_label>_task-<task_id>_bold_space-<space>_T1w_preproc_smooth8mm_Denoised.nii.gz
  • Linear detrending is done in the reference pipeline I suggested before the ICA-AROMA calculations, I assume this is not currently possible with fMRIprep? The fallback solution you proposed is doing it at the noise regression step by adding a row “1, 2, …, size of 4th dim” to sub-<subject_label>_task-<task_id>_bold_confounds_transpose.tsv if I understood correctly. I am not sure of the exact differences it creates compared to pre-AROMA calculations detrending.
  • Regarding WM and CSF regressors, from a statistical point of view, wouldn’t the correlation between WhiteMatter and aCompCor01 be a problem for the stability of the GLM (since aCompCor01 also includes WhiteMatter)? Maybe applying the solution proposed below would be helpful?

PPS The reason why we output only WM but not CSF is that for participants with small ventricles we cannot reliably estimate CSF ROIs without running into partial volume issues. I would love to hear suggestions how to work around this.

What about using a robust CSF estimation (e.g. tissue probability maps with 99% threshold, eroded CSF mask) and ignoring CSF signal if the number of voxels is below what you have in mind when you mentioned “small ventricles”? Could ignoring the signal be implemented as n/a values in the confounds file (i saw such n/a values in the first row of these files for the columns stdDVARS and non-stdDVARS )? This would allow to get a distinct CSF regressor column in the confounds file.

13

Thanks for the useful reference @oesteban. This point was actually addressed in the manuscript proposing the pipeline I referenced, where they tested despiking both before and after denoising, and did not see any significant difference in any of their pipeline performance metrics.

Considering it would be difficult at the moment to introduce despiking early on in the pipeline with fmriprep, it is then probably better (easier) to do it at the end.

14

Those columns are only added if the non-steady state was detected. If it wasn’t (probably because your scanner acquired dummy scans to deal with it before starting the sequence) no such columns will be added to the _confounds.tsv file. In such case I don’t see why would you need to discard any volumes, but it would not hurt (unless your sequence is very short).

The 6mm smoothed data is used to create the _variant-smoothAROMAnonaggr_preproc.nii.gz. The other _preproc.nii.gz file is unsmoothed.

Possibly - you would have to check if fsl_regfilt is not going to get stuck on headers (if it does it would be worth pinging FSL people if this could be fixed- selecting columns by names rather than numbers seems less error prone)

I don’t expect a huge difference. If you don’t detrent prior to ICA one of the components will most likely pick on the trend anyway.

Yes many of those regressors will be colinear, but since we only care about the residuals of this model this is not a problem. Colinearity will lead to variance in parameter estimates, but will not negatively impact the residuals.

Good idea. https://github.com/poldracklab/fmriprep/issues/702

15

Awesome, thanks @ChrisGorgolewski. Just two small remaining issues (let me know if I should write this somewhere else):

  1. I do not have any visualization output of the ICA-AROMA analysis, as described and shown in the fmriprep documentation

  2. You mentioned that:

The 6mm smoothed data is used to create the _variant-smoothAROMAnonaggr_preproc.nii.gz. The other _preproc.nii.gz file is unsmoothed.

But I do not have the smoothed file. In the func directory I have the following files:

sub-<subject>_task-rest_bold_AROMAnoiseICs.csv
sub-<subject>_task-rest_bold_confounds.tsv
sub-<subject>_task-rest_bold_MELODICmix.tsv
sub-<subject>_task-rest_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz
sub-<subject>_task-rest_bold_space-MNI152NLin2009cAsym_preproc.nii.gz
sub-<subject>_task-rest_bold_space-T1w_brainmask.nii.gz
sub-<subject>_task-rest_bold_space-T1w_preproc.nii.gz

And in the anat directory, I have (excluding surfaces):

sub-<subject>_T1w_brainmask.nii.gz
sub-<subject>_T1w_class-CSF_probtissue.nii.gz
sub-<subject>_T1w_class-GM_probtissue.nii.gz
sub-<subject>_T1w_class-WM_probtissue.nii.gz
sub-<subject>_T1w_dtissue.nii.gz
sub-<subject>_T1w_preproc.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_brainmask.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_class-CSF_probtissue.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_class-GM_probtissue.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_class-WM_probtissue.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_dtissue.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_preproc.nii.gz
sub-<subject>_T1w_space-MNI152NLin2009cAsym_warp.h5
16

This looks like a bug. Could you open an issue on github sharing an HTML report and the command line invocation you used?

17

Hi @michael. Sorry if I was unable to find it, but what version of fmriprep are you using? The outputs you’re missing were both added in 1.0.0-rc1 (a.k.a. 0.6.1), so if you’re using 0.6.0 or earlier, the issue is that you need to upgrade.

If not, I would encourage you to upgrade to the latest, and re-run. This should result in better error reporting, and make it easier to address this issue.

It would also be good to open a new issue on GitHub (or just a new thread in NeuroStars) to keep this thread from getting further off-topic.

18

Thank you @effigies, I now have to complement the pipeline discussed above with additional steps on the fmriprep outputs I already have, but for the next task I will certainly upgrade fmriprep to the latest version as you suggested (and if required open a new issue on github at that point).

19

RE: median 10000 normalization

I haven’t looked at the code to see the reason why, but FSL’s FILM_GLS does not like negative values or even values close to zero I believe, the performance of that algorithm depends on high(ish) voxel intensities, so making the median normalization makes it so that no voxels are in danger of being too close to zero or negative.

20

@ChrisGorgolewski Did you mean indeed aCompCor01 or aCompCor00? Also I have on average 30 +/- 10 AROMAAggrComp components per subject. Does that seem reasonable?