I wanted to share a code I picked up at the Afni Bootcamp! One can use fmriprep output with afni_proc to do the following analysis (My example: blur scaling and regressions). The code might throw an error when generating the QC page. But no worries, the main analysis will still wrap up just fine. The QC page is a bit tricky and might need to wait for an update in the future.
Best
#!/usr/bin/env tcsh
# set data directory
set top_dir = /media/test/afni_fmriprep/sub-LM027
# run afni_proc.py to create a single subject processing script
afni_proc.py -subj_id sub-LM027 \
-blocks mask blur scale regress \
-radial_correlate_blocks tcat regress \
-copy_anat $top_dir/sub-LM027_ses-T2_desc-preproc_T1w.nii.gz \
-dsets \
$top_dir/sub-LM027_ses-T2_task-Prototype1_space-T1w_desc-preproc_bold.nii.gz \
-blur_size 4.0 \
-regress_motion_file \
$top_dir/afni_sub-LM027_ses-T2_task-Prototype1_confounds.txt \
-regress_stim_times \
$top_dir/afni_sub-LM027_ses-T2_task-Prototype1_allTrials.txt \
# I used the following to transform fsl-type regressor into afni-type. \
# timing_tool.py -fsl_timing_files input -write_timing output \
-regress_stim_labels \
allTrials \
-regress_basis 'BLOCK(5,1)' \
-regress_censor_motion 1.0 \
-regress_censor_outliers 0.05 \
-regress_motion_per_run \
-regress_opts_3dD \
-jobs 1 \
-gltsym 'SYM: allTrials ' -glt_label 1 V-A \
-regress_compute_fitts \
-regress_make_ideal_sum sum_ideal.1D \
-regress_est_blur_epits \
-regress_est_blur_errts \
-regress_run_clustsim no \
-html_review_style pythonic \
-execute
Curious does the error with the QC report look like the one I described in my post? Trying to figure out if I need to solve that problem or if I can keep going with my analysis and just be sad I don't get a pretty QC output.
I'm on vacation at the moment, but will just comment briefly:
As the afni_proc.py proc script processes, a dictionary of useful file (inputs, intermediate files and outputs) is built up.
The programs that build the APQC HTML use that dictionary to know what to build. For example, knowing what dset/file defines the final space (either template, anat_final or volreg_epi reference vol).
The fmriprep output covers the first few processing blocks of a typical afni_proc.py run. But items from that are not part of the reference dictionary above.
This particular error this causes for the HTML build issue that is referenced in this post, is that the volume defining the final space isn't known for the reference dictionary.
The tricky thing is that it appears there could be multiple spaces to choose from at the point fmriprep finishes. So, the user would have to specify one, and the question of how best to do that remains a bit.
This is extremely useful! I know it is months later, but I have a few questions in case anyone is still able to respond to this thread:
I am wondering about the option -regress_motion_file: are you including all confounds exported by fmriprep (e.g. all translations, rotations, derivatives, GSR, AROMA, aCompCor, cosines, etc)? Is it a subset of parameters for denoising? Is it just the 6 translation/rotation parameters?
I see that the anatomical file is in T1w (i.e. native) space-- is it recommended to do the regression in native space and then warp the results to another space, e.g. MNI non-linear 2009c? Or could one simply use the already-warped data as an input?
This script has the block order as mask, blur, scale and regress. In another example on OSF, I found that in their script used blur, mask, scale and regress. Are there pros and cons to the mask/blur steps being switched around?
When using -regress_motion_file, one should pass a file with only the 6 rigid body motion parameters. This file is used (as one mechanism) to detect motion and to censor, and it is expected to have 6 columns.
Other parameters can be passed in a separate file using -regress_extra_ortvec (and -regress_extra_ortvec_labels), or in multiple files if you wanted to pass cosines separately from aCompCor, for example. But putting them together into one file is okay, too.
I expect the fmriprep outputs are in MNI space already, aren't they? The NIFTI datasets should specify this. Check using something like:
The order between mask and blur is not too important, as the mask is really just used for QC. But it seems preferable to do mask then blur. We had older examples that had blur first.
Thank you so much @rickr for your speedy and precise response! Good to know that the motion regressors only want 6 terms, and that mask then blur is preferred but it's not a big deal.
fMRIPrep can output data in multiple spaces if you ask it nicely. By default it uses the MNI152NLin6Asym, but we asked for both native space and also MNI152NLin2009cAsym as well. It's handily specified right in the output filenames using the "space-" key, e.g.: sub-02_ses-01_task-mytask_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz (normalized to template) sub-02_ses-01_task-mytask_space-T1w_desc-preproc_bold.nii.gz (aligned to subject's anatomical T1, i.e. native space)
For good measure, I did check your code snippet as well, and based on the information here, it looks consistent.
I have a couple of follow-up questions, if you'll permit me.
The fMRIPrep confounds file already has regressors that identify non-steady-state outliers (i.e. high signal TRs near the beginning of the scan when it is still warming up). I think AFNI favors discarding these using the -tcat_remove_first_trs option. Would it be okay to include the fMRIPrep regressors (there is one regressor per outlier scan) in -regress_extra_ortvec instead of dropping those initial scans, so that I don't have to edit all the other confounds to remove the first few values?
Similarly, fMRIPrep identifies motion outliers based on criteria specified when it was run (by default, FD > 0.5 mm or standardized DVARS > 1.5). Each motion outlier has its own regressor. Would it be acceptable to use these regressors in -regress_extra_ortvec instead of -regress_censor_motion and -regress_censor_outliers?
Note that with something like "tail -n +5 FILE", it will print out all lines from 5 to the end, for example. That is a scriptable way to remove the first 4 lines from a text file.
But far more importantly, the reason that we prefer to throw away the pre-steady state TRs is that they affect any temporal operations, which we might rather not happen. Some issues, off the top of my head:
outliers and the MIN_OUTLIER TR: The notion of what it means to be an outlier would be highly impacted by pre-SS volumes.
tshift : This will blur spikes, including pre-SS, and esp with -fourier interpolation.
combine/tedana : This has 2 aspects, one is that we might not want the pre-SS affecting combining echoes (via OC), but tedana (or any PCA/ICA-based regression) is going to have to fight with that pre-SS spike, and there is no reason to add such confusion to the software.
scaling : The pre-SS time will slightly affect the scaling, reducing it a bit.
other PC regression : Any computation of PCs will be heavily affected by those initial volumes.
QC, such as radcor : Correlations for QC, such as the helpful radial correlations would probably just end up with mostly red images, where the entire volume looks correlated.
other QC: Surely this would affect a lot of QC, such as variance lines, TSNR evaluations, censoring and motion averages and statistics.
Anyway, more examples came quickly to mind than I anticipated. But after this nice little exercise I will more strongly say, please omit the pre-SS volumes. :)
Oh, and regarding censoring, you can do it however you feel is best. There is a -regress_censor_extern file that you could pass, or use the fmriprep regressors. But I suggest comparing them to what afni_proc.py would use (for censoring motion and outliers). Well, I guess you cannot properly compare outliers with DVARS, since afni_proc.py will not get to compute outliers on unregistered data, but still, there are details you can think about.
Interesting, these are definitely some good examples to consider. I will probably drop these frames before the GLM because of the scaling issue and radcor usefulness. Also, now that I had an idea of what to look for, I did find some extra information:
Based on this neurostars thread, it looks like fMRIPrep ignores any detected dummy scans when computing things like slice time correction and aCompCor.
"CompCor will ignore the dummy volumes to avoid giving you regressors that capture non-steady-state effects that will be accounted for by the non_steady_state_XX regressors."
I also found this this neurostars thread specifically for tedana. Looks like dropping those volumes is the right way to go for now in that instance (but it's not part of what I'm trying to do).
Cheers!
(Edit: here is one more relevant post that might help anyone looking at this in the future.)
Thanks a lot for the links.
Maybe I missed it, but what I did not see in any of those posts is a good reason to keep the dummy scans in the data to be processed. Because of the differential signal, they should not benefit the regression, so it is not clear why they would be kept in any of the preprocessing, since they then need to be avoided for so many of the steps (processing/QC/PCA). The only place we might use a pre-SS volume is in using index #0 for anat/epi registration, if the epi has poor contrast. But as that is provided as a separate volume, it still does not come with the EPI time series.
I agree, I can't think of a lot of reasons to keep those volumes around. But I can see why fMRIPrep would be conservative and not remove data unilaterally, especially because it may cause misalignment with experiment timing files that people may have generated using another program. At least this way, you must remove the volumes yourself, and you're more likely to remember to adjust experiment timing files to match. That's my best guess.
I have checked in my own files generated by fmriprep version 20.2.7 LTS, and confirmed that most fmriprep confound regressors (aCompCor, AROMA, etc) have values of zero for the non-steady-state volumes, and the regressors are also normalized so that their mean is zero. This was not the case for movement parameters, such as x/y/z_rot, x/y/z_trans, their derivatives/powers, and dvars and framewise_displacement. These have values even during the non-steady-state scans, but they are also not normalized (they just have their raw values), so removing the first few would not mess with any normalization. Anyway, I'm now confident that the fmriprep regressors are set up so that those first frames can indeed be dropped.
And finally, one more question on the original code posted: if I am going to use ETAC later to determine the significance of clusters, I would want to completely eliminate the blur block and its options (e.g. -blur_size 4.0), correct?
Sure, it does sound like the early time points can just be dropped.
Blurring with ETAC depends on what results and images you want to show in a paper. Volumetric results would be shown at a single blur level (or none), yet in that case ETAC would be evaluated across multiple blur levels. That gives a bit of a question regarding what is to be shown. But ETAC can be used without multiple blur levels (omitting -ETAC_blur), in which case it is fine to apply whatever you prefer initially.
rick
The
National Institute of Mental Health (NIMH) is part of the National Institutes of
Health (NIH), a component of the U.S. Department of Health and Human
Services.