14.2.10. 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#!/usr/bin/env tcsh
  2
  3# ----------------------------------------------------------------------
  4# Script: s.bmn_subject_level_02_ap.tcsh
  5# Run on: openfmri/ds001_R2.0.4
  6# Date  : April, 2018
  7#
  8# Time series analysis for one subject, basically unmodified from the
  9# BMN afni_proc.py command.  (Does not require pre-running of any
 10# script.)
 11#
 12# Used for "BMN-AFNI" processing in:
 13#
 14#   Some comments and corrections on FMRI processing with AFNI in
 15#   "Exploring the Impact of Analysis Software on Task fMRI Results"
 16# 
 17#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
 18#   Richard C. Reynolds, Robert W. Cox.
 19#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
 20#   Bethesda, MD, USA.
 21#
 22# NOTE: This afni_proc.py command follows the Bowringer et al. (2018)
 23# command, described on bioRxiv.  Variable names have been adjusted
 24# for convenience of running on local computers, but not selected
 25# options.  See their bioRxiv draft and github page for reference.
 26#
 27# ----------------------------------------------------------------------
 28#
 29# To run for a single subject, for example:
 30#
 31#   tcsh -ef s.bmn_subject_level_02_ap.tcsh "01"
 32#
 33# After this, one can run the following for QC and checking processing
 34# steps:
 35#
 36#   @ss_review_basic
 37#   @ss_review_driver
 38#
 39# ... and after all subjects in the group have been processed, the
 40# following can be run to combine the single subject processing
 41# numbers into one table:
 42#
 43#   gen_ss_review_table.py                       \
 44#      -tablefile review_table_bmn.xls           \
 45#      -infiles data_proc_bmn/sub*/out.ss_review*
 46#
 47#########################################################################
 48
 49# ===================== Set inputs and paths =========================
 50
 51# The only input argument is the subject number (zero-paddeed, as
 52# "01", "02", etc., if that is how the file naming conventions are)
 53if( $#argv < 1 )then
 54    echo "** ERROR: no command line argument?" 
 55    exit 1
 56endif
 57
 58# Set subject ID from input
 59set ss         = $argv[1]
 60set subj       = sub-${ss}  
 61set grp        = "bmn"
 62
 63# Set reference template for standard space (location found below)
 64set itemp      = MNI_avg152T1+tlrc.HEAD
 65
 66# Make path/file structure variables for readability.
 67
 68# Top level directory
 69set p_0        = /data/NIMH_SSCC
 70set path_main  = ${p_0}/openfmri/ds001_R2.0.4
 71# Inputs
 72set prep_onset = ${path_main}/data_basic/ONSETS
 73set path_func  = ${path_main}/data_basic/sub-${ss}/func
 74set path_anat  = ${path_main}/data_basic/sub-${ss}/anat
 75# Outputs
 76set path_ogrp  = ${path_main}/data_proc_${grp}
 77set path_out   = ${path_ogrp}/sub-${ss}
 78
 79mkdir -p ${path_ogrp}
 80
 81# =================== Environment variables ===========================
 82
 83# NB: need this because of initial data's qform code
 84setenv AFNI_NIFTI_VIEW orig
 85
 86# ================== Optional: Biowulf cluster =========================
 87
 88# This section for possibly running on Biowulf cluster.
 89
 90# Set thread count
 91if( $?SLURM_CPUS_PER_TASK )then
 92  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
 93endif
 94
 95# Set temporary output directory; then requires using something like
 96# this on the swarm command line: --sbatch '--gres=lscratch:100'.
 97# These variables used again *after* afni_proc.py command, if Biowulfing.
 98if( $?SLURM_JOBID )then
 99  set tempdir = /lscratch/$SLURM_JOBID/${subj}
100  set usetemp = 1
101else
102  set tempdir = $path_out
103  set usetemp = 0
104endif
105
106# ================== Find standard space template ======================
107
108# Selecting/finding reference template (for defining final standard
109# space); the chosen template file is specified above
110
111set tpath    = `@FindAfniDsetPath $itemp`
112set basedset = $tpath/$itemp
113
114if( "$tpath" == "" )then
115    echo "** ERROR: can't find my template $itemp"
116    exit 1
117endif
118
119if( ! -f $basedset )then
120    echo "** ERROR: template $basedset is not present" 
121    exit 1
122endif
123
124# ========================= Specify processing =========================
125
126# [PT: July 10, 2019] AP command amended, to have ${tempdir} provided
127#    to the '-out_dir ..' opt; allows actual use of scratch disk, if
128#    asked for above. No change in processing results, just speed of
129#    processing in some cases.  Thanks, Yoichi M!
130
131# FINALLY: the afni_proc.py command itself, which makes a single
132# subject processing script
133afni_proc.py                                                                \
134    -scr_overwrite                                                          \
135    -subj_id sub${ss}                                                       \
136    -script  proc_${grp}.sub${ss}                                           \
137    -out_dir ${tempdir}                                                     \
138    -scr_overwrite                                                          \
139    -blocks tshift align tlrc volreg blur mask scale regress                \
140    -copy_anat ${path_anat}/sub-${ss}_T1w.nii                               \
141    -dsets                                                                  \
142       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-01_bold.nii.gz \
143       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-02_bold.nii.gz \
144       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-03_bold.nii.gz \
145    -tcat_remove_first_trs 2                                                \
146    -align_opts_aea -cost mi -ginormous_move -check_flip                    \
147    -anat_uniform_method unifize                                            \
148    -tlrc_base ${basedset}                                                  \
149    -volreg_warp_dxyz 2                                                     \
150    -volreg_align_to third                                                  \
151    -volreg_align_e2a                                                       \
152    -volreg_tlrc_warp                                                       \
153    -blur_size 5.0                                                          \
154    -regress_stim_times                                                     \
155       ${prep_onset}/sub-${ss}_combined_cash_demean_afni.1d                 \
156       ${prep_onset}/sub-${ss}_combined_cash_RT_afni.1d                     \
157       ${prep_onset}/sub-${ss}_combined_control_pumps_demean_afni.1d        \
158       ${prep_onset}/sub-${ss}_combined_control_pumps_RT_afni.1d            \
159       ${prep_onset}/sub-${ss}_combined_explode_demean_afni.1d              \
160       ${prep_onset}/sub-${ss}_combined_pumps_demean_afni.1d                \
161       ${prep_onset}/sub-${ss}_combined_pumps_RT_afni.1d                    \
162    -regress_stim_labels                                                    \
163       cash_demean cash_RT control_pumps_demean                             \
164       control_pumps_RT explode_demean pumps_demean                         \
165       pumps_RT                                                             \
166    -regress_basis_multi                                                    \
167       'BLOCK(0.772,1)' 'dmBLOCK' 'BLOCK(0.772,1)' 'dmBLOCK'                \
168       'BLOCK(0.772,1)' 'BLOCK(0.772,1)' 'dmBLOCK'                          \
169    -regress_stim_types                                                     \
170       AM2 AM1 AM2 AM1 AM2 AM2 AM1                                          \
171    -regress_censor_motion 0.3                                              \
172    -regress_opts_3dD                                                       \
173    -gltsym 'SYM: pumps_demean -control_pumps_demean'                       \
174    -glt_label 1 pumps_demean_vs_ctrl_demean                                \
175    -regress_make_ideal_sum sum_ideal.1D                                    \
176    -regress_est_blur_epits                                                 \
177    -regress_est_blur_errts                                                 \
178    -execute
179
180# =============================================================================
181
182# Again, Biowulf-running considerations: if processing went fine and
183# we were using a temporary directory, copy data back.
184if( $usetemp && -d $tempdir )then
185    echo "Copying data from $tempdir to $path_out"
186    mkdir -p $path_out
187    \cp -pr $tempdir/* $path_out/
188endif
189
190# If it worked, run some volreg snapshots and compress outputs:
191# compare chosen EPI reference volume (from alignment to anatomical)
192# to the final anatomical.
193if( -d $path_out )then
194    cd $path_out
195    @snapshot_volreg                       \
196        anat_final.${subj}+tlrc.HEAD       \
197        final_epi_vr_base*+tlrc.HEAD       \
198        ${subj}
199    cd ..
200endif
201
202# =============================================================================
203
204echo "++ DONE with afni_proc for subject: $subj"
205
206time ; exit 0

Script: s.nimh_subject_level_01_qwarp.tcsh

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

 1#!/usr/bin/env tcsh
 2
 3# ----------------------------------------------------------------------
 4# Script: s.nimh_subject_level_01_qwarp.tcsh
 5# Run on: openfmri/ds001_R2.0.4
 6# Date  : April, 2018
 7#
 8# This script runs @SSwarper for one subject, whose cognomen (ID) is
 9# the sole command line argument. The results are used in the analysis
10# script.
11#
12# From @SSwarper's help: 
13#   Script to skull-strip and warp to the MNI 2009 template.
14#
15# Used for "NIMH-AFNI" processing in:
16#
17#   Some comments and corrections on FMRI processing with AFNI in
18#   "Exploring the Impact of Analysis Software on Task fMRI Results"
19# 
20#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
21#   Richard C. Reynolds, Robert W. Cox.
22#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
23#   Bethesda, MD, USA.
24#
25# ----------------------------------------------------------------------
26#
27# To run for a single subject, for example:
28#
29#   tcsh -ef s.nimh_subject_level_01_qwarp.tcsh "01"
30#
31#########################################################################
32
33# The only input argument is the subject number (zero-paddeed, as
34# "01", "02", etc., if that is how the file naming conventions are)
35if( $#argv < 1 )then
36    echo "** ERROR: no command line argument?" 
37    exit 1
38endif
39
40# Set subject ID 
41set ss = $argv[1]
42
43# Make path/file structure variables for readability.
44
45# Top level directory
46set p_0        = /data/NIMH_SSCC
47set path_main  = ${p_0}/openfmri/ds001_R2.0.4
48# Inputs
49set prep_onset = ${path_main}/data_basic/ONSETS
50set path_func  = ${path_main}/data_basic/sub-${ss}/func
51set path_anat  = ${path_main}/data_basic/sub-${ss}/anat
52
53# ========================================================================
54
55# This section for possibly running on Biowulf cluster.
56
57# Set thread count
58if( $?SLURM_CPUS_PER_TASK )then
59  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
60endif
61
62# ========================================================================
63
64# switch to the subject’s anat directory
65if( ! -e $path_anat/sub-${ss}_T1w.nii.gz ) then
66    echo "** ERROR: can’t find file $path_anat/sub-${ss}_T1w.nii" 
67    exit 1
68endif
69
70cd $path_anat
71
72# Do the work
73@SSwarper sub-${ss}_T1w.nii.gz sub-${ss}
74
75# =============================================================================
76
77echo "++ DONE with warping for subject: $subj"
78
79time ; exit 0

Script: s.nimh_subject_level_02_ap.tcsh

  1#!/usr/bin/env tcsh
  2
  3# ----------------------------------------------------------------------
  4# Script: s.nimh_subject_level_02_ap.tcsh
  5# Run on: openfmri/ds001_R2.0.4
  6# Date  : April, 2018
  7#
  8# Time series analysis for one subject, modified from BMN afni_proc.py
  9# command to adhere to current AFNI usage recommendations. The
 10# accompanying "s.nimh_subject_level_01_qwarp.tcsh" script must have
 11# been run previously to get the warped-to-MNI files for the current
 12# subject via @SSwarper.
 13#
 14# Used for "NIMH-AFNI" processing in:
 15#
 16#   Some comments and corrections on FMRI processing with AFNI in
 17#   "Exploring the Impact of Analysis Software on Task fMRI Results"
 18# 
 19#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
 20#   Richard C. Reynolds, Robert W. Cox.
 21#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
 22#   Bethesda, MD, USA.
 23#
 24# NOTE: This afni_proc.py command includes several features that would
 25# have been chosen in preference to those of BMN, as described in the
 26# above text.  *However*, it is not a complete set.  Further features
 27# should also likely be adjusted, such as in the GLT specification
 28# (see points A-5 and A-6 in the accompanying text) and upsampling to
 29# 2x2x2 mm^3 would also note be chosen.  These features were not
 30# changed *here*, for the purposes of matching the BMN specifications
 31# more closely.  In general, please see the descriptions in the above
 32# bioRxiv text for details on settings/options that were chosen.
 33#
 34# ----------------------------------------------------------------------
 35#
 36# To run for a single subject, for example:
 37#
 38#   tcsh -ef s.nimh_subject_level_02_ap.tcsh "01"
 39#
 40# After this, one can run the following for QC and checking processing
 41# steps:
 42#
 43#   @ss_review_basic
 44#   @ss_review_driver
 45#
 46# ... and after all subjects in the group have been processed, the
 47# following can be run to combine the single subject processing
 48# numbers into one table:
 49#
 50#   gen_ss_review_table.py                       \
 51#      -tablefile review_table_nimh.xls          \
 52#      -infiles data_proc_nimh/sub*/out.ss_review*
 53#
 54#########################################################################
 55
 56# ===================== Set inputs and paths =========================
 57
 58# The only input argument is the subject number (zero-paddeed, as
 59# "01", "02", etc., if that is how the file naming conventions are)
 60if( $#argv < 1 )then
 61    echo "** ERROR: no command line argument?" 
 62    exit 1
 63endif
 64
 65# Set subject ID from input
 66set ss         = $argv[1]
 67set subj       = sub-${ss}  
 68set grp        = "nimh"
 69
 70# Set reference template for standard space (location found below)
 71set itemp      = MNI152_2009_template.nii.gz
 72
 73# Make path/file structure variables for readability.
 74
 75# Top level directory
 76set p_0        = /data/NIMH_SSCC
 77set path_main  = ${p_0}/openfmri/ds001_R2.0.4
 78# Inputs
 79set prep_onset = ${path_main}/data_basic/ONSETS
 80set path_func  = ${path_main}/data_basic/sub-${ss}/func
 81set path_anat  = ${path_main}/data_basic/sub-${ss}/anat
 82# Outputs
 83set path_ogrp  = ${path_main}/data_proc_${grp}
 84set path_out   = ${path_ogrp}/sub-${ss}
 85
 86mkdir -p ${path_ogrp}
 87
 88# =================== Environment variables ===========================
 89
 90# NB: need this because of initial data's qform code
 91setenv AFNI_NIFTI_VIEW orig
 92
 93# ================== Optional: Biowulf cluster =========================
 94
 95# This section for possibly running on Biowulf cluster.
 96
 97# Set thread count
 98if( $?SLURM_CPUS_PER_TASK )then
 99  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
100endif
101
102# Set temporary output directory; then requires using something like
103# this on the swarm command line: --sbatch '--gres=lscratch:100'.
104# These variables used again *after* afni_proc.py command, if Biowulfing.
105if( $?SLURM_JOBID )then
106  set tempdir = /lscratch/$SLURM_JOBID/${subj}
107  set usetemp = 1
108else
109  set tempdir = $path_out
110  set usetemp = 0
111endif
112
113# ================== Find standard space template ======================
114
115# Selecting/finding reference template (for defining final standard
116# space); the chosen template file is specified above
117
118set tpath    = `@FindAfniDsetPath $itemp`
119set basedset = $tpath/$itemp
120
121if( "$tpath" == "" )then
122    echo "** ERROR: can't find my template $itemp"
123    exit 1
124endif
125
126if( ! -f $basedset )then
127    echo "** ERROR: template $basedset is not present" 
128    exit 1
129endif
130
131# ========================= Specify processing =========================
132
133# [PT: July 10, 2019] AP command amended, to have ${tempdir} provided
134#    to the '-out_dir ..' opt; allows actual use of scratch disk, if
135#    asked for above. No change in processing results, just speed of
136#    processing in some cases.  Thanks, Yoichi M!
137
138# FINALLY: the afni_proc.py command itself, which makes a single
139# subject processing script
140afni_proc.py                                                                \
141    -scr_overwrite                                                          \
142    -subj_id  sub${ss}                                                      \
143    -script   proc_${grp}.sub${ss}                                          \
144    -out_dir  ${tempdir}                                                    \
145    -blocks tshift align tlrc volreg blur mask scale regress                \
146    -copy_anat ${path_anat}/anatSS.${subj}.nii                              \
147          -anat_has_skull no                                                \
148    -dsets                                                                  \
149       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-01_bold.nii.gz \
150       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-02_bold.nii.gz \
151       ${path_func}/sub-${ss}_task-balloonanalogrisktask_run-03_bold.nii.gz \
152    -tcat_remove_first_trs 2                                                \
153    -align_opts_aea -cost lpc+ZZ -ginormous_move -check_flip                \
154    -volreg_warp_dxyz 2                                                     \
155    -volreg_align_to MIN_OUTLIER                                            \
156    -volreg_align_e2a                                                       \
157    -volreg_tlrc_warp                                                       \
158    -tlrc_base $basedset                                                    \
159    -tlrc_NL_warp                                                           \
160    -tlrc_NL_warped_dsets                                                   \
161         $path_anat/anatQQ.${subj}.nii                                      \
162         $path_anat/anatQQ.${subj}.aff12.1D                                 \
163         $path_anat/anatQQ.${subj}_WARP.nii                                 \
164    -blur_size 5.0                                                          \
165    -regress_stim_times                                                     \
166       ${prep_onset}/sub-${ss}_combined_cash_demean_afni.1d                 \
167       ${prep_onset}/sub-${ss}_combined_cash_RT_afni.1d                     \
168       ${prep_onset}/sub-${ss}_combined_control_pumps_demean_afni.1d        \
169       ${prep_onset}/sub-${ss}_combined_control_pumps_RT_afni.1d            \
170       ${prep_onset}/sub-${ss}_combined_explode_demean_afni.1d              \
171       ${prep_onset}/sub-${ss}_combined_pumps_demean_afni.1d                \
172       ${prep_onset}/sub-${ss}_combined_pumps_RT_afni.1d                    \
173    -regress_stim_labels                                                    \
174       cash_demean cash_RT control_pumps_demean                             \
175       control_pumps_RT explode_demean pumps_demean                         \
176       pumps_RT                                                             \
177    -regress_basis_multi                                                    \
178       'BLOCK(0.772,1)' 'dmBLOCK' 'BLOCK(0.772,1)' 'dmBLOCK'                \
179       'BLOCK(0.772,1)' 'BLOCK(0.772,1)' 'dmBLOCK'                          \
180    -regress_stim_types                                                     \
181       AM2 AM1 AM2 AM1 AM2 AM2 AM1                                          \
182    -regress_censor_motion 0.3                                              \
183    -regress_motion_per_run                                                 \
184    -regress_opts_3dD                                                       \
185    -gltsym 'SYM: pumps_demean -control_pumps_demean'                       \
186    -glt_label 1 pumps_demean_vs_ctrl_demean                                \
187    -regress_make_ideal_sum sum_ideal.1D                                    \
188    -regress_3dD_stop                                                       \
189    -regress_reml_exec                                                      \
190    -regress_est_blur_epits                                                 \
191    -regress_est_blur_errts                                                 \
192    -execute
193
194# =============================================================================
195
196# Again, Biowulf-running considerations: if processing went fine and
197# we were using a temporary directory, copy data back.
198if( $usetemp && -d $tempdir )then
199    echo "Copying data from $tempdir to $path_out"
200    mkdir -p $path_out
201    \cp -pr $tempdir/* $path_out/
202endif
203
204# If it worked, run some volreg snapshots and compress outputs:
205# compare chosen EPI reference volume (from alignment to anatomical)
206# to the final anatomical.
207if( -d $path_out )then
208    cd $path_out
209    @snapshot_volreg                       \
210        anat_final.${subj}+tlrc.HEAD       \
211        final_epi_vr_base*+tlrc.HEAD       \
212        ${subj}
213    cd ..
214endif
215
216# =============================================================================
217
218echo "++ DONE with afni_proc for subject: $subj"
219
220time ; 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#!/bin/tcsh
  2
  3# ----------------------------------------------------------------------
  4# Script: s.nimh_group_level_02_mema.tcsh
  5# Run on: openfmri/ds001_R2.0.4
  6# Date  : April, 2018
  7#
  8# Take a group of subjects processed for task-based FMRI, run 3dMEMA
  9# on their final statistical results.  Clustering is also performed
 10# here, using 3dClustSim (with the ACF blur estimation).  Some
 11# individual volumes of t-statistics and effect estimate maps are also
 12# output.
 13#
 14# Used for "NIMH-AFNI" processing in:
 15#
 16#   Some comments and corrections on FMRI processing with AFNI in
 17#   "Exploring the Impact of Analysis Software on Task fMRI Results"
 18# 
 19#   Paul A. Taylor, Gang Chen, Daniel R. Glen, Justin K. Rajendra,
 20#   Richard C. Reynolds, Robert W. Cox.
 21#   Scientific and Statistical Computing Core, NIMH/NIH/DHHS,
 22#   Bethesda, MD, USA.
 23#
 24# NOTE: This set of commands includes some features that would have
 25# been chosen in preference to those of BMN, as described in the above
 26# text (e.g., using the '-missing_data 0' option in 3dMEMA).
 27# *However*, there are still some things that should be altered, such
 28# as using paired one-sided testing-- the search for clusters here is
 29# two-sided, and so the two-sided test should be used for correctness.
 30#
 31# ----------------------------------------------------------------------
 32#
 33# To run for a single subject, for example:
 34#
 35#   tcsh -ef s.nimh_group_level_02_mema.tcsh
 36#
 37# ----------------------------------------------------------------------
 38
 39# ================ Set up paths and many params ======================
 40
 41set here        = $PWD
 42set p_0         = /data/NIMH_SSCC
 43set path_main   = ${p_0}/openfmri/ds001_R2.0.4
 44set path_proc   = ${path_main}/data_proc_nimh
 45set odir        = GROUP_LEVEL_nimh
 46set path_mema   = ${path_proc}/${odir}
 47
 48# Running 3dMEMA
 49set script_mema = do_group_level_nimh.tcsh # generated command file name
 50set omema       = mema_results.nii.gz      # output effect+stats filename
 51set omema_mskd  = mema_results_mskd.nii.gz # masked, stats in brain
 52
 53# Cluster parameters
 54set csim_NN     = "NN1"    # neighborhood; could be NN1, NN2, NN3
 55set csim_sided  = "1sided" # could be 1sided, 2sided or bisided
 56set csim_pthr   = 0.01     # voxelwise thr in 3dClustSim
 57set csim_alpha  = 0.05     # nominal FWE
 58
 59# ================== Optional: Biowulf cluster =========================
 60
 61# This section for possibly running on Biowulf cluster.
 62
 63#### The following 2 lines could be uncommented to load these modules;
 64#### that would likely not be necessary outside of NIH's Biowulf cluster.
 65# source /etc/profile.d/modules.csh   
 66# module load afni R
 67
 68# set thread count if we are running SLURM
 69if ( $?SLURM_CPUS_PER_TASK ) then
 70  setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
 71endif
 72
 73# Set temporary output directory; then requires using something like
 74# this on the swarm command line: --sbatch '--gres=lscratch:50'.
 75# These variables used again *after* afni_proc.py command, if running
 76# on Biowulf.
 77if( $?SLURM_JOBID )then
 78  set tempdir = /lscratch/$SLURM_JOBID/${odir}
 79  set usetemp = 1
 80else
 81  set tempdir = $path_mema
 82  set usetemp = 0
 83endif
 84
 85mkdir -p $tempdir
 86
 87# ====================== Voxelwise modeling ==========================
 88
 89# Create and execute 3dMEMA command.  Note the inclusion of the
 90# "-options -missing_data 0" flags, to avoid some numerical badness of
 91# possible FOV edge effects, etc.
 92gen_group_command.py                                           \
 93    -command      3dMEMA                                       \
 94    -write_script ${tempdir}/${script_mema}                    \
 95    -prefix       ${tempdir}/${omema}                          \
 96    -dsets        ${path_proc}/sub*/stats.sub*_REML+tlrc.HEAD  \
 97    -set_labels   "controls"                                   \
 98    -subs_betas   'pumps_demean_vs_ctrl_demean#0_Coef'         \
 99    -subs_tstats  'pumps_demean_vs_ctrl_demean#0_Tstat'        \
100    -options      -missing_data 0
101
102tcsh -ef ${tempdir}/${script_mema}
103echo "++ MEMAed!"
104
105cd ${tempdir}
106echo "++ Now in the group stat dir: $PWD"
107
108# ====================== Make group mask ==========================
109
110# Make a group level mask for reporting results.  This also determines
111# the volume over which multiple comparisons corrections for voxelwise
112# stats are done.
1133dmask_tool                                             \
114    -prefix mask.nii.gz                                 \
115    -input `ls ${path_proc}/sub-*/mask_epi_anat.*.HEAD` \
116    -frac 1.0
117
118# ==================== Calc average smoothness =======================
119
120# Get the line from each blur estimate file that has our desired
121# parameters; we are running REML here and using the ACF smoothness
122# estimation, so the gpattern variable reflects both.  For N subjects,
123# the output *1D file contains N rows of the 3 ACF parameters.
124set gpattern = 'err_reml ACF'
125grep -h "$gpattern" ${path_proc}/sub*/blur_est*   \
126    | cut -d\  -f1-3                              \
127    > group_ACF_ests.1D
128
129# Get the group average ACF parameters.
130set blur_est = ( `3dTstat -mean -prefix - group_ACF_ests.1D\'` )
131echo "++ The group average ACF params are: $blur_est"
132
133# ==================== Cluster simulations =======================
134
135# Simulations for FWE corrected cluster-size inference 
1363dClustSim                                       \
137    -both                                        \
138    -mask   mask.nii.gz                          \
139    -acf    $blur_est                            \
140    -prefix ClustSim 
141
142# Get the volume threshold value from the ClustSim results, based on
143# user's choice of parameters as set above.  The "-verb 0" means that
144# just the threshold number of voxels is returned, so we can save that
145# as a variable.
146set clust_thrvol = `1d_tool.py -verb 0                                \
147                        -infile ClustSim.${csim_NN}_${csim_sided}.1D  \
148                        -csim_pthr   $csim_pthr                       \
149                        -csim_alpha "$csim_alpha"`
150
151# Get the statistic value equivalent to the desired voxelwise p-value,
152# for thresholding purposes.  Using the same p-value and sidedness
153# that were selected in the ClustSim results.  This program also gets
154# the number of degrees of freedom (DOF) from the header of the volume
155# containing the statistic. The "-quiet" means that only the
156# associated statistic value is returned, so we can save it to a
157# variable.
158set voxstat_thr = `p2dsetstat -quiet                    \
159                    -pval $csim_pthr                    \
160                    "-${csim_sided}"                    \
161                    -inset ${omema}'[controls:t]'`
162
163echo "++ The final cluster volume threshold is:  $clust_thrvol"
164echo "++ The voxelwise stat value threshold is:  $voxstat_thr"
165
166# ================== Make cluster maps =====================
167
168# Apply WB mask to the statistics results, so we only find clusters
169# that lie within the mask of interest.
1703dcalc                              \
171    -a ${omema}                     \
172    -b mask.nii.gz                  \
173    -expr 'a*b'                     \
174    -prefix ${omema_mskd} 
175
176# Following the previous paper, these commands extract the "positive"
177# and "negative" cluster info, respectively, into separate files.
178# Each applies clustering parameters above and makes: a *map* of
179# cluster ROIs (regions numbered 1, 2, 3, ...); a map of the effect
180# estimate (EE) values within those clusters; and a table report of
181# the clusters.  
182# Technical note: '-1tindex 1' means that the threshold is applied to
183# the [1]st volume of the input dset, which here holds the statistic
184# values; '-1dindex 0' means that data values saved in '-prefix ...'
185# come from the [0]th volume of the input dset, which here holds the
186# EE values.  
187set csim_pref = "clust_tpos"
1883dclust                                                        \
189    -1Dformat -nosum -1tindex 1 -1dindex 0                     \
190    -2thresh -1e+09 $voxstat_thr  -dxyz=1                      \
191    -savemask ${csim_pref}_map.nii.gz                          \
192    -prefix   ${csim_pref}_EE.nii.gz                           \
193    1.01 ${clust_thrvol} ${omema_mskd} > ${csim_pref}_table.1D
194# Also make a mask of t-stats (not necessary, but matching previous
195# work)
1963dcalc \
197    -a ${csim_pref}_map.nii.gz                \
198    -b ${omema}'[1]'                          \
199    -expr "step(a)*b"                         \
200    -prefix ${csim_pref}_t.nii.gz             \
201    -float
202
203set csim_pref = "clust_tneg"
2043dclust                                                        \
205    -1Dformat -nosum -1tindex 1 -1dindex 0                     \
206    -2thresh -$voxstat_thr 1e+09 -dxyz=1                       \
207    -savemask ${csim_pref}_map.nii.gz                          \
208    -prefix   ${csim_pref}_EE.nii.gz                           \
209    1.01 ${clust_thrvol} ${omema_mskd} > ${csim_pref}_table.1D
210# Also make a mask of t-stats (not necessary, but matching previous
211# work)
2123dcalc \
213    -a ${csim_pref}_map.nii.gz                \
214    -b ${omema}'[1]'                          \
215    -expr "step(a)*b"                         \
216    -prefix ${csim_pref}_t.nii.gz             \
217    -float
218
219cd $here
220
221# ===================================================================
222
223# Again, Biowulf-running considerations: if processing went fine and
224# we were using a temporary directory, copy data back.
225if( $usetemp && -d $tempdir )then
226    echo "Copying data from $tempdir to $path_mema"
227    mkdir -p $path_mema
228    \cp -pr $tempdir/* $path_mema/
229endif
230
231# ===================================================================
232
233echo "\n++ DONE with group level stats + clustering!\n"
234
235time ; exit 0