Taylor et al. (2018). FMRI processing with AFNI: Some comments and corrections on “Exploring the Impact of Analysis Software on Task fMRI Results

Taylor PA, Chen GC, 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”

See also

-> brief entry

Tag Label
FMRI paradigm: task-block
FMRI dset: EPI
Anatomical dset: MPRAGE
Subject population: human
Subject characteristic: control
Subject age: adult
Template space: MNI
Template align method: nonlinear
Tissue segmentation:
Tissue regression:
Comments:  

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
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.
  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
 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
  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
  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