13.2.8. Taylor et al. (2018). FMRI processing with AFNI: Some comments and corrections …

Introduction

Here we present commands used in the following paper:

  • Taylor PA, Chen G, Glen DR, Rajendra JK, Reynolds RC, Cox RW (2018). FMRI processing with AFNI: Some comments and corrections on ‘Exploring the Impact of Analysis Software on Task fMRI Results’. bioRxiv 308643; doi:10.1101/308643

Abstract: A recent study posted on bioRxiv by Bowring, Maumet and Nichols aimed to compare results of FMRI data that had been processed with three commonly used software packages (AFNI, FSL and SPM). Their stated purpose was to use “default” settings of each software’s pipeline for task-based FMRI, and then to quantify overlaps in final clustering results and to measure similarity/dissimilarity in the final outcomes of packages. While in theory the setup sounds simple (implement each package’s defaults and compare results), practical realities make this difficult. For example, different softwares would recommend different spatial resolutions of the final data, but for the sake of comparisons, the same value must be used across all. Moreover, we would say that AFNI does not have an explicit default pipeline available: a wide diversity of datasets and study designs are acquired across the neuroimaging community, often requiring bespoke tailoring of basic processing rather than a “one-size-fits-all” pipeline. However, we do have strong recommendations for certain steps, and we are also aware that the choice of a given step might place requirements on other processing steps. Given the very clear reporting of the AFNI pipeline used in Bowring et al. paper, we take this opportunity to comment on some of these aspects of processing with AFNI here, clarifying a few mistakes therein and also offering recommendations. We provide point-by-point considerations of using AFNI’s processing pipeline design tool at the individual level, afni_proc.py, along with supplementary programs; while specifically discussed in the context of the present usage, many of these choices may serve as useful starting points for broader processing. It is our intention/hope that the user should examine data quality at every step, and we demonstrate how this is facilitated in AFNI, as well.

Study keywords: task-block, EPI, MPRAGE, human, control, adult, MNI space, nonlinear align

Main programs: afni_proc.py, 3dmask_tool, recon-all (FS), @SUMA_Make_Spec_FS

Download scripts

To download, either:

  • ... click the link(s) in the following table (perhaps Rightclick -> “Save Link As…”):

    s.bmn_subject_level_02_ap.tcsh

    BMN-AFNI processing with afni_proc.py

    s.nimh_subject_level_01_qwarp.tcsh

    NIMH-AFNI processing: skullstripping and alignment to standard space (via @SSwarper)

    s.nimh_subject_level_02_ap.tcsh

    NIMH-AFNI processing with afni_proc.py

    s.nimh_group_level_02_mema.tcsh

    NIMH-AFNI group level processing: 3dMEMA, 3dClustSim, clusterizing

  • ... or copy+paste into a terminal:

    curl -O https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/codex/fmri/media/2018_TaylorEtal/s.bmn_subject_level_02_ap.tcsh
    curl -O https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/codex/fmri/media/2018_TaylorEtal/s.nimh_subject_level_01_qwarp.tcsh
    curl -O https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/codex/fmri/media/2018_TaylorEtal/s.nimh_subject_level_02_ap.tcsh
    curl -O https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/codex/fmri/media/2018_TaylorEtal/s.nimh_group_level_02_mema.tcsh
    

View scripts

These scripts describe different approaches for processing FMRI data with AFNI. Please read the comments at the tops of the scripts carefully, as well as the bioRxiv papers associated with each, in order to understand the steps.

These scripts were run on the NIH’s Biowulf HPC cluster, hence some of the variables like $SLURM_CPUS_PER_TASK, $SLURM_JOBID, etc.

s.bmn_subject_level_02_ap.tcsh

Comment: Several aspects of this particular script’s command would not be recommended for processing (e.g., the nonlinear alignment program; the selection of volreg reference, and more). Please see the main paper for further details.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#!/usr/bin/env tcsh

# ----------------------------------------------------------------------
# Script: s.bmn_subject_level_02_ap.tcsh
# Run on: openfmri/ds001_R2.0.4
# Date  : April, 2018
#
# Time series analysis for one subject, basically unmodified from the
# BMN afni_proc.py command.  (Does not require pre-running of any
# script.)
#
# Used for "BMN-AFNI" processing in:
#
#   Some comments and corrections on FMRI processing with AFNI in
#   "Exploring the Impact of Analysis Software on Task fMRI Results"
# 
#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
#   Richard C. Reynolds, Robert W. Cox.
#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
#   Bethesda, MD, USA.
#
# NOTE: This afni_proc.py command follows the Bowringer et al. (2018)
# command, described on bioRxiv.  Variable names have been adjusted
# for convenience of running on local computers, but not selected
# options.  See their bioRxiv draft and github page for reference.
#
# ----------------------------------------------------------------------
#
# To run for a single subject, for example:
#
#   tcsh -ef s.bmn_subject_level_02_ap.tcsh "01"
#
# After this, one can run the following for QC and checking processing
# steps:
#
#   @ss_review_basic
#   @ss_review_driver
#
# ... and after all subjects in the group have been processed, the
# following can be run to combine the single subject processing
# numbers into one table:
#
#   gen_ss_review_table.py                       \
#      -tablefile review_table_bmn.xls           \
#      -infiles data_proc_bmn/sub*/out.ss_review*
#
#########################################################################

# ===================== Set inputs and paths =========================

# The only input argument is the subject number (zero-paddeed, as
# "01", "02", etc., if that is how the file naming conventions are)
if( $#argv < 1 )then
    echo "** ERROR: no command line argument?" 
    exit 1
endif

# Set subject ID from input
set ss         = $argv[1]
set subj       = sub-${ss}  
set grp        = "bmn"

# Set reference template for standard space (location found below)
set itemp      = MNI_avg152T1+tlrc.HEAD

# Make path/file structure variables for readability.

# Top level directory
set p_0        = /data/NIMH_SSCC
set path_main  = ${p_0}/openfmri/ds001_R2.0.4
# Inputs
set prep_onset = ${path_main}/data_basic/ONSETS
set path_func  = ${path_main}/data_basic/sub-${ss}/func
set path_anat  = ${path_main}/data_basic/sub-${ss}/anat
# Outputs
set path_ogrp  = ${path_main}/data_proc_${grp}
set path_out   = ${path_ogrp}/sub-${ss}

mkdir -p ${path_ogrp}

# =================== Environment variables ===========================

# NB: need this because of initial data's qform code
setenv AFNI_NIFTI_VIEW orig

# ================== Optional: Biowulf cluster =========================

# This section for possibly running on Biowulf cluster.

# Set thread count
if( $?SLURM_CPUS_PER_TASK )then
  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
endif

# Set temporary output directory; then requires using something like
# this on the swarm command line: --sbatch '--gres=lscratch:100'.
# These variables used again *after* afni_proc.py command, if Biowulfing.
if( $?SLURM_JOBID )then
  set tempdir = /lscratch/$SLURM_JOBID/${subj}
  set usetemp = 1
else
  set tempdir = $path_out
  set usetemp = 0
endif

# ================== Find standard space template ======================

# Selecting/finding reference template (for defining final standard
# space); the chosen template file is specified above

set tpath    = `@FindAfniDsetPath $itemp`
set basedset = $tpath/$itemp

if( "$tpath" == "" )then
    echo "** ERROR: can't find my template $itemp"
    exit 1
endif

if( ! -f $basedset )then
    echo "** ERROR: template $basedset is not present" 
    exit 1
endif

# ========================= Specify processing =========================

# [PT: July 10, 2019] AP command amended, to have ${tempdir} provided
#    to the '-out_dir ..' opt; allows actual use of scratch disk, if
#    asked for above. No change in processing results, just speed of
#    processing in some cases.  Thanks, Yoichi M!

# FINALLY: the afni_proc.py command itself, which makes a single
# subject processing script
afni_proc.py                                                                \
    -scr_overwrite                                                          \
    -subj_id sub${ss}                                                       \
    -script  proc_${grp}.sub${ss}                                           \
    -out_dir ${tempdir}                                                     \
    -scr_overwrite                                                          \
    -blocks tshift align tlrc volreg blur mask scale regress                \
    -copy_anat ${path_anat}/sub-${ss}_T1w.nii                               \
    -dsets                                                                  \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-01_bold.nii.gz \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-02_bold.nii.gz \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-03_bold.nii.gz \
    -tcat_remove_first_trs 2                                                \
    -align_opts_aea -cost mi -ginormous_move -check_flip                    \
    -anat_uniform_method unifize                                            \
    -tlrc_base ${basedset}                                                  \
    -volreg_warp_dxyz 2                                                     \
    -volreg_align_to third                                                  \
    -volreg_align_e2a                                                       \
    -volreg_tlrc_warp                                                       \
    -blur_size 5.0                                                          \
    -regress_stim_times                                                     \
       ${prep_onset}/sub-${ss}_combined_cash_demean_afni.1d                 \
       ${prep_onset}/sub-${ss}_combined_cash_RT_afni.1d                     \
       ${prep_onset}/sub-${ss}_combined_control_pumps_demean_afni.1d        \
       ${prep_onset}/sub-${ss}_combined_control_pumps_RT_afni.1d            \
       ${prep_onset}/sub-${ss}_combined_explode_demean_afni.1d              \
       ${prep_onset}/sub-${ss}_combined_pumps_demean_afni.1d                \
       ${prep_onset}/sub-${ss}_combined_pumps_RT_afni.1d                    \
    -regress_stim_labels                                                    \
       cash_demean cash_RT control_pumps_demean                             \
       control_pumps_RT explode_demean pumps_demean                         \
       pumps_RT                                                             \
    -regress_basis_multi                                                    \
       'BLOCK(0.772,1)' 'dmBLOCK' 'BLOCK(0.772,1)' 'dmBLOCK'                \
       'BLOCK(0.772,1)' 'BLOCK(0.772,1)' 'dmBLOCK'                          \
    -regress_stim_types                                                     \
       AM2 AM1 AM2 AM1 AM2 AM2 AM1                                          \
    -regress_censor_motion 0.3                                              \
    -regress_opts_3dD                                                       \
    -gltsym 'SYM: pumps_demean -control_pumps_demean'                       \
    -glt_label 1 pumps_demean_vs_ctrl_demean                                \
    -regress_make_ideal_sum sum_ideal.1D                                    \
    -regress_est_blur_epits                                                 \
    -regress_est_blur_errts                                                 \
    -execute

# =============================================================================

# Again, Biowulf-running considerations: if processing went fine and
# we were using a temporary directory, copy data back.
if( $usetemp && -d $tempdir )then
    echo "Copying data from $tempdir to $path_out"
    mkdir -p $path_out
    \cp -pr $tempdir/* $path_out/
endif

# If it worked, run some volreg snapshots and compress outputs:
# compare chosen EPI reference volume (from alignment to anatomical)
# to the final anatomical.
if( -d $path_out )then
    cd $path_out
    @snapshot_volreg                       \
        anat_final.${subj}+tlrc.HEAD       \
        final_epi_vr_base*+tlrc.HEAD       \
        ${subj}
    cd ..
endif

# =============================================================================

echo "++ DONE with afni_proc for subject: $subj"

time ; exit 0

Script: s.nimh_subject_level_01_qwarp.tcsh

Comment: @SSwarper syntax has changed a lot since these days.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#!/usr/bin/env tcsh

# ----------------------------------------------------------------------
# Script: s.nimh_subject_level_01_qwarp.tcsh
# Run on: openfmri/ds001_R2.0.4
# Date  : April, 2018
#
# This script runs @SSwarper for one subject, whose cognomen (ID) is
# the sole command line argument. The results are used in the analysis
# script.
#
# From @SSwarper's help: 
#   Script to skull-strip and warp to the MNI 2009 template.
#
# Used for "NIMH-AFNI" processing in:
#
#   Some comments and corrections on FMRI processing with AFNI in
#   "Exploring the Impact of Analysis Software on Task fMRI Results"
# 
#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
#   Richard C. Reynolds, Robert W. Cox.
#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
#   Bethesda, MD, USA.
#
# ----------------------------------------------------------------------
#
# To run for a single subject, for example:
#
#   tcsh -ef s.nimh_subject_level_01_qwarp.tcsh "01"
#
#########################################################################

# The only input argument is the subject number (zero-paddeed, as
# "01", "02", etc., if that is how the file naming conventions are)
if( $#argv < 1 )then
    echo "** ERROR: no command line argument?" 
    exit 1
endif

# Set subject ID 
set ss = $argv[1]

# Make path/file structure variables for readability.

# Top level directory
set p_0        = /data/NIMH_SSCC
set path_main  = ${p_0}/openfmri/ds001_R2.0.4
# Inputs
set prep_onset = ${path_main}/data_basic/ONSETS
set path_func  = ${path_main}/data_basic/sub-${ss}/func
set path_anat  = ${path_main}/data_basic/sub-${ss}/anat

# ========================================================================

# This section for possibly running on Biowulf cluster.

# Set thread count
if( $?SLURM_CPUS_PER_TASK )then
  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
endif

# ========================================================================

# switch to the subject’s anat directory
if( ! -e $path_anat/sub-${ss}_T1w.nii.gz ) then
    echo "** ERROR: can’t find file $path_anat/sub-${ss}_T1w.nii" 
    exit 1
endif

cd $path_anat

# Do the work
@SSwarper sub-${ss}_T1w.nii.gz sub-${ss}

# =============================================================================

echo "++ DONE with warping for subject: $subj"

time ; exit 0

Script: s.nimh_subject_level_02_ap.tcsh

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
#!/usr/bin/env tcsh

# ----------------------------------------------------------------------
# Script: s.nimh_subject_level_02_ap.tcsh
# Run on: openfmri/ds001_R2.0.4
# Date  : April, 2018
#
# Time series analysis for one subject, modified from BMN afni_proc.py
# command to adhere to current AFNI usage recommendations. The
# accompanying "s.nimh_subject_level_01_qwarp.tcsh" script must have
# been run previously to get the warped-to-MNI files for the current
# subject via @SSwarper.
#
# Used for "NIMH-AFNI" processing in:
#
#   Some comments and corrections on FMRI processing with AFNI in
#   "Exploring the Impact of Analysis Software on Task fMRI Results"
# 
#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
#   Richard C. Reynolds, Robert W. Cox.
#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
#   Bethesda, MD, USA.
#
# NOTE: This afni_proc.py command includes several features that would
# have been chosen in preference to those of BMN, as described in the
# above text.  *However*, it is not a complete set.  Further features
# should also likely be adjusted, such as in the GLT specification
# (see points A-5 and A-6 in the accompanying text) and upsampling to
# 2x2x2 mm^3 would also note be chosen.  These features were not
# changed *here*, for the purposes of matching the BMN specifications
# more closely.  In general, please see the descriptions in the above
# bioRxiv text for details on settings/options that were chosen.
#
# ----------------------------------------------------------------------
#
# To run for a single subject, for example:
#
#   tcsh -ef s.nimh_subject_level_02_ap.tcsh "01"
#
# After this, one can run the following for QC and checking processing
# steps:
#
#   @ss_review_basic
#   @ss_review_driver
#
# ... and after all subjects in the group have been processed, the
# following can be run to combine the single subject processing
# numbers into one table:
#
#   gen_ss_review_table.py                       \
#      -tablefile review_table_nimh.xls          \
#      -infiles data_proc_nimh/sub*/out.ss_review*
#
#########################################################################

# ===================== Set inputs and paths =========================

# The only input argument is the subject number (zero-paddeed, as
# "01", "02", etc., if that is how the file naming conventions are)
if( $#argv < 1 )then
    echo "** ERROR: no command line argument?" 
    exit 1
endif

# Set subject ID from input
set ss         = $argv[1]
set subj       = sub-${ss}  
set grp        = "nimh"

# Set reference template for standard space (location found below)
set itemp      = MNI152_2009_template.nii.gz

# Make path/file structure variables for readability.

# Top level directory
set p_0        = /data/NIMH_SSCC
set path_main  = ${p_0}/openfmri/ds001_R2.0.4
# Inputs
set prep_onset = ${path_main}/data_basic/ONSETS
set path_func  = ${path_main}/data_basic/sub-${ss}/func
set path_anat  = ${path_main}/data_basic/sub-${ss}/anat
# Outputs
set path_ogrp  = ${path_main}/data_proc_${grp}
set path_out   = ${path_ogrp}/sub-${ss}

mkdir -p ${path_ogrp}

# =================== Environment variables ===========================

# NB: need this because of initial data's qform code
setenv AFNI_NIFTI_VIEW orig

# ================== Optional: Biowulf cluster =========================

# This section for possibly running on Biowulf cluster.

# Set thread count
if( $?SLURM_CPUS_PER_TASK )then
  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
endif

# Set temporary output directory; then requires using something like
# this on the swarm command line: --sbatch '--gres=lscratch:100'.
# These variables used again *after* afni_proc.py command, if Biowulfing.
if( $?SLURM_JOBID )then
  set tempdir = /lscratch/$SLURM_JOBID/${subj}
  set usetemp = 1
else
  set tempdir = $path_out
  set usetemp = 0
endif

# ================== Find standard space template ======================

# Selecting/finding reference template (for defining final standard
# space); the chosen template file is specified above

set tpath    = `@FindAfniDsetPath $itemp`
set basedset = $tpath/$itemp

if( "$tpath" == "" )then
    echo "** ERROR: can't find my template $itemp"
    exit 1
endif

if( ! -f $basedset )then
    echo "** ERROR: template $basedset is not present" 
    exit 1
endif

# ========================= Specify processing =========================

# [PT: July 10, 2019] AP command amended, to have ${tempdir} provided
#    to the '-out_dir ..' opt; allows actual use of scratch disk, if
#    asked for above. No change in processing results, just speed of
#    processing in some cases.  Thanks, Yoichi M!

# FINALLY: the afni_proc.py command itself, which makes a single
# subject processing script
afni_proc.py                                                                \
    -scr_overwrite                                                          \
    -subj_id  sub${ss}                                                      \
    -script   proc_${grp}.sub${ss}                                          \
    -out_dir  ${tempdir}                                                    \
    -blocks tshift align tlrc volreg blur mask scale regress                \
    -copy_anat ${path_anat}/anatSS.${subj}.nii                              \
          -anat_has_skull no                                                \
    -dsets                                                                  \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-01_bold.nii.gz \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-02_bold.nii.gz \
       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-03_bold.nii.gz \
    -tcat_remove_first_trs 2                                                \
    -align_opts_aea -cost lpc+ZZ -ginormous_move -check_flip                \
    -volreg_warp_dxyz 2                                                     \
    -volreg_align_to MIN_OUTLIER                                            \
    -volreg_align_e2a                                                       \
    -volreg_tlrc_warp                                                       \
    -tlrc_base $basedset                                                    \
    -tlrc_NL_warp                                                           \
    -tlrc_NL_warped_dsets                                                   \
         $path_anat/anatQQ.${subj}.nii                                      \
         $path_anat/anatQQ.${subj}.aff12.1D                                 \
         $path_anat/anatQQ.${subj}_WARP.nii                                 \
    -blur_size 5.0                                                          \
    -regress_stim_times                                                     \
       ${prep_onset}/sub-${ss}_combined_cash_demean_afni.1d                 \
       ${prep_onset}/sub-${ss}_combined_cash_RT_afni.1d                     \
       ${prep_onset}/sub-${ss}_combined_control_pumps_demean_afni.1d        \
       ${prep_onset}/sub-${ss}_combined_control_pumps_RT_afni.1d            \
       ${prep_onset}/sub-${ss}_combined_explode_demean_afni.1d              \
       ${prep_onset}/sub-${ss}_combined_pumps_demean_afni.1d                \
       ${prep_onset}/sub-${ss}_combined_pumps_RT_afni.1d                    \
    -regress_stim_labels                                                    \
       cash_demean cash_RT control_pumps_demean                             \
       control_pumps_RT explode_demean pumps_demean                         \
       pumps_RT                                                             \
    -regress_basis_multi                                                    \
       'BLOCK(0.772,1)' 'dmBLOCK' 'BLOCK(0.772,1)' 'dmBLOCK'                \
       'BLOCK(0.772,1)' 'BLOCK(0.772,1)' 'dmBLOCK'                          \
    -regress_stim_types                                                     \
       AM2 AM1 AM2 AM1 AM2 AM2 AM1                                          \
    -regress_censor_motion 0.3                                              \
    -regress_motion_per_run                                                 \
    -regress_opts_3dD                                                       \
    -gltsym 'SYM: pumps_demean -control_pumps_demean'                       \
    -glt_label 1 pumps_demean_vs_ctrl_demean                                \
    -regress_make_ideal_sum sum_ideal.1D                                    \
    -regress_3dD_stop                                                       \
    -regress_reml_exec                                                      \
    -regress_est_blur_epits                                                 \
    -regress_est_blur_errts                                                 \
    -execute

# =============================================================================

# Again, Biowulf-running considerations: if processing went fine and
# we were using a temporary directory, copy data back.
if( $usetemp && -d $tempdir )then
    echo "Copying data from $tempdir to $path_out"
    mkdir -p $path_out
    \cp -pr $tempdir/* $path_out/
endif

# If it worked, run some volreg snapshots and compress outputs:
# compare chosen EPI reference volume (from alignment to anatomical)
# to the final anatomical.
if( -d $path_out )then
    cd $path_out
    @snapshot_volreg                       \
        anat_final.${subj}+tlrc.HEAD       \
        final_epi_vr_base*+tlrc.HEAD       \
        ${subj}
    cd ..
endif

# =============================================================================

echo "++ DONE with afni_proc for subject: $subj"

time ; exit 0

Script: s.nimh_group_level_02_mema.tcsh

Comment: Instead of using 3dclust here, the more modern and convenient 3dClusterize would now be recommended.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
#!/bin/tcsh

# ----------------------------------------------------------------------
# Script: s.nimh_group_level_02_mema.tcsh
# Run on: openfmri/ds001_R2.0.4
# Date  : April, 2018
#
# Take a group of subjects processed for task-based FMRI, run 3dMEMA
# on their final statistical results.  Clustering is also performed
# here, using 3dClustSim (with the ACF blur estimation).  Some
# individual volumes of t-statistics and effect estimate maps are also
# output.
#
# Used for "NIMH-AFNI" processing in:
#
#   Some comments and corrections on FMRI processing with AFNI in
#   "Exploring the Impact of Analysis Software on Task fMRI Results"
# 
#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
#   Richard C. Reynolds, Robert W. Cox.
#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
#   Bethesda, MD, USA.
#
# NOTE: This set of commands includes some features that would have
# been chosen in preference to those of BMN, as described in the above
# text (e.g., using the '-missing_data 0' option in 3dMEMA).
# *However*, there are still some things that should be altered, such
# as using paired one-sided testing-- the search for clusters here is
# two-sided, and so the two-sided test should be used for correctness.
#
# ----------------------------------------------------------------------
#
# To run for a single subject, for example:
#
#   tcsh -ef s.nimh_group_level_02_mema.tcsh
#
# ----------------------------------------------------------------------

# ================ Set up paths and many params ======================

set here        = $PWD
set p_0         = /data/NIMH_SSCC
set path_main   = ${p_0}/openfmri/ds001_R2.0.4
set path_proc   = ${path_main}/data_proc_nimh
set odir        = GROUP_LEVEL_nimh
set path_mema   = ${path_proc}/${odir}

# Running 3dMEMA
set script_mema = do_group_level_nimh.tcsh # generated command file name
set omema       = mema_results.nii.gz      # output effect+stats filename
set omema_mskd  = mema_results_mskd.nii.gz # masked, stats in brain

# Cluster parameters
set csim_NN     = "NN1"    # neighborhood; could be NN1, NN2, NN3
set csim_sided  = "1sided" # could be 1sided, 2sided or bisided
set csim_pthr   = 0.01     # voxelwise thr in 3dClustSim
set csim_alpha  = 0.05     # nominal FWE

# ================== Optional: Biowulf cluster =========================

# This section for possibly running on Biowulf cluster.

#### The following 2 lines could be uncommented to load these modules;
#### that would likely not be necessary outside of NIH's Biowulf cluster.
# source /etc/profile.d/modules.csh   
# module load afni R

# set thread count if we are running SLURM
if ( $?SLURM_CPUS_PER_TASK ) then
  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
endif

# Set temporary output directory; then requires using something like
# this on the swarm command line: --sbatch '--gres=lscratch:50'.
# These variables used again *after* afni_proc.py command, if running
# on Biowulf.
if( $?SLURM_JOBID )then
  set tempdir = /lscratch/$SLURM_JOBID/${odir}
  set usetemp = 1
else
  set tempdir = $path_mema
  set usetemp = 0
endif

mkdir -p $tempdir

# ====================== Voxelwise modeling ==========================

# Create and execute 3dMEMA command.  Note the inclusion of the
# "-options -missing_data 0" flags, to avoid some numerical badness of
# possible FOV edge effects, etc.
gen_group_command.py                                           \
    -command      3dMEMA                                       \
    -write_script ${tempdir}/${script_mema}                    \
    -prefix       ${tempdir}/${omema}                          \
    -dsets        ${path_proc}/sub*/stats.sub*_REML+tlrc.HEAD  \
    -set_labels   "controls"                                   \
    -subs_betas   'pumps_demean_vs_ctrl_demean#0_Coef'         \
    -subs_tstats  'pumps_demean_vs_ctrl_demean#0_Tstat'        \
    -options      -missing_data 0

tcsh -ef ${tempdir}/${script_mema}
echo "++ MEMAed!"

cd ${tempdir}
echo "++ Now in the group stat dir: $PWD"

# ====================== Make group mask ==========================

# Make a group level mask for reporting results.  This also determines
# the volume over which multiple comparisons corrections for voxelwise
# stats are done.
3dmask_tool                                             \
    -prefix mask.nii.gz                                 \
    -input `ls ${path_proc}/sub-*/mask_epi_anat.*.HEAD` \
    -frac 1.0

# ==================== Calc average smoothness =======================

# Get the line from each blur estimate file that has our desired
# parameters; we are running REML here and using the ACF smoothness
# estimation, so the gpattern variable reflects both.  For N subjects,
# the output *1D file contains N rows of the 3 ACF parameters.
set gpattern = 'err_reml ACF'
grep -h "$gpattern" ${path_proc}/sub*/blur_est*   \
    | cut -d\  -f1-3                              \
    > group_ACF_ests.1D

# Get the group average ACF parameters.
set blur_est = ( `3dTstat -mean -prefix - group_ACF_ests.1D\'` )
echo "++ The group average ACF params are: $blur_est"

# ==================== Cluster simulations =======================

# Simulations for FWE corrected cluster-size inference 
3dClustSim                                       \
    -both                                        \
    -mask   mask.nii.gz                          \
    -acf    $blur_est                            \
    -prefix ClustSim 

# Get the volume threshold value from the ClustSim results, based on
# user's choice of parameters as set above.  The "-verb 0" means that
# just the threshold number of voxels is returned, so we can save that
# as a variable.
set clust_thrvol = `1d_tool.py -verb 0                                \
                        -infile ClustSim.${csim_NN}_${csim_sided}.1D  \
                        -csim_pthr   $csim_pthr                       \
                        -csim_alpha "$csim_alpha"`

# Get the statistic value equivalent to the desired voxelwise p-value,
# for thresholding purposes.  Using the same p-value and sidedness
# that were selected in the ClustSim results.  This program also gets
# the number of degrees of freedom (DOF) from the header of the volume
# containing the statistic. The "-quiet" means that only the
# associated statistic value is returned, so we can save it to a
# variable.
set voxstat_thr = `p2dsetstat -quiet                    \
                    -pval $csim_pthr                    \
                    "-${csim_sided}"                    \
                    -inset ${omema}'[controls:t]'`

echo "++ The final cluster volume threshold is:  $clust_thrvol"
echo "++ The voxelwise stat value threshold is:  $voxstat_thr"

# ================== Make cluster maps =====================

# Apply WB mask to the statistics results, so we only find clusters
# that lie within the mask of interest.
3dcalc                              \
    -a ${omema}                     \
    -b mask.nii.gz                  \
    -expr 'a*b'                     \
    -prefix ${omema_mskd} 

# Following the previous paper, these commands extract the "positive"
# and "negative" cluster info, respectively, into separate files.
# Each applies clustering parameters above and makes: a *map* of
# cluster ROIs (regions numbered 1, 2, 3, ...); a map of the effect
# estimate (EE) values within those clusters; and a table report of
# the clusters.  
# Technical note: '-1tindex 1' means that the threshold is applied to
# the [1]st volume of the input dset, which here holds the statistic
# values; '-1dindex 0' means that data values saved in '-prefix ...'
# come from the [0]th volume of the input dset, which here holds the
# EE values.  
set csim_pref = "clust_tpos"
3dclust                                                        \
    -1Dformat -nosum -1tindex 1 -1dindex 0                     \
    -2thresh -1e+09 $voxstat_thr  -dxyz=1                      \
    -savemask ${csim_pref}_map.nii.gz                          \
    -prefix   ${csim_pref}_EE.nii.gz                           \
    1.01 ${clust_thrvol} ${omema_mskd} > ${csim_pref}_table.1D
# Also make a mask of t-stats (not necessary, but matching previous
# work)
3dcalc \
    -a ${csim_pref}_map.nii.gz                \
    -b ${omema}'[1]'                          \
    -expr "step(a)*b"                         \
    -prefix ${csim_pref}_t.nii.gz             \
    -float

set csim_pref = "clust_tneg"
3dclust                                                        \
    -1Dformat -nosum -1tindex 1 -1dindex 0                     \
    -2thresh -$voxstat_thr 1e+09 -dxyz=1                       \
    -savemask ${csim_pref}_map.nii.gz                          \
    -prefix   ${csim_pref}_EE.nii.gz                           \
    1.01 ${clust_thrvol} ${omema_mskd} > ${csim_pref}_table.1D
# Also make a mask of t-stats (not necessary, but matching previous
# work)
3dcalc \
    -a ${csim_pref}_map.nii.gz                \
    -b ${omema}'[1]'                          \
    -expr "step(a)*b"                         \
    -prefix ${csim_pref}_t.nii.gz             \
    -float

cd $here

# ===================================================================

# Again, Biowulf-running considerations: if processing went fine and
# we were using a temporary directory, copy data back.
if( $usetemp && -d $tempdir )then
    echo "Copying data from $tempdir to $path_mema"
    mkdir -p $path_mema
    \cp -pr $tempdir/* $path_mema/
endif

# ===================================================================

echo "\n++ DONE with group level stats + clustering!\n"

time ; exit 0