resting state output different between AFNI versions


I seem to be having an issue with different versions of AFNI, where I get very different resting-state output between 2016 versions and 2019/2020. Using the same data and scripts, I get resting state results that make sense when using AFNI 2016 but the results appear random when using AFNI 2019. I have verified that it is not due to one subject or the hardware on one computer. It seems that the issue arises during alignment; output is the same until reaching the alignment routine.

Below is the script in question (; this file is called for each subject by a separate script). I tried to paste the output also but got an error that my message is too long.

Please help!


#!/bin/tcsh -xef

echo “auto-generated by, Wed Jan 15 19:33:11 2020”
echo “(version 4.72, June 30, 2016)”
echo “execution started: date

execute via :

tcsh -xef proc.NA8831 |& tee output.proc.NA8831

=========================== auto block: setup ============================

script setup

take note of the AFNI version

afni -ver

check that the current AFNI version is recent enough

afni_history -check_date 28 Oct 2015
if ( $status ) then
echo “** this script requires newer AFNI binaries (than 28 Oct 2015)”
echo " (consider: @update.afni.binaries -defaults)"

accept arguments from

set subj = $argv[1]
set nifti = $argv[2]
set anat = $argv[3]
set afnidir = $argv[4]

the user may specify a single subject to run with

#if ( $#argv > 0 ) then

set subj = $argv[1]


set subj = NA8831


assign output directory name

set output_dir = $subj.results

assign subject directory name

set subjdir=pwd

echo “subj = $subj”
echo “nifti = $nifti”
echo “anat = $anat”
echo “afnidir = $afnidir”
echo “subjdir = $subjdir”

verify that the results directory does not yet exist

if ( -d $output_dir ) then
echo output dir “$subj.results” already exists

set list of runs

set runs = (count -digits 2 1 1)

create results directory

mkdir $subjdir/$output_dir

copy anatomy to results dir

3dcopy $subjdir/$anat $subjdir/$output_dir/anat

============================ auto block: tcat ============================

apply 3dTcat to copy input dsets to results dir, while

removing the first 4 TRs

3dTcat -prefix $subjdir/$output_dir/pb00.$subj.r01.tcat $subjdir/$nifti’[4…$]’

and make note of repetitions (TRs) per run

set tr_counts = ( 266 )


enter the results directory (can begin processing data)

cd $subjdir/$output_dir

========================== auto block: outcount ==========================

data check: compute outlier fraction for each volume

touch out.pre_ss_warn.txt
foreach run ( $runs )
3dToutcount -automask -fraction -polort 4 -legendre
pb00.$subj.r$run.tcat+orig > outcount.r$run.1D

# censor outlier TRs per run, ignoring the first 0 TRs
# - censor when more than 0.25 of automask voxels are outliers
# - step() defines which TRs to remove via censoring
1deval -a outcount.r$run.1D -expr "1-step(a-0.25)" > rm.out.cen.r$run.1D

# outliers at TR 0 might suggest pre-steady state TRs
if ( `1deval -a outcount.r$run.1D"{0}" -expr "step(a-0.4)"` ) then
    echo "** TR #0 outliers: possible pre-steady state TRs in run $run" \
        >> out.pre_ss_warn.txt


catenate outlier counts into a single time series

cat outcount.r*.1D > outcount_rall.1D

catenate outlier censor files into a single time series

cat rm.out.cen.r*.1D > outcount_${subj}_censor.1D

================================ despike =================================

apply 3dDespike to each run

foreach run ( $runs )
3dDespike -NEW -nomask -prefix pb01.$subj.r$run.despike

================================= tshift =================================

time shift data so all slice timing is the same

foreach run ( $runs )
3dTshift -tzero 0 -quintic -prefix pb02.$subj.r$run.tshift


extract volreg registration base

3dbucket -prefix vr_base pb02.$subj.r01.tshift+orig"[0]"

================================= align ==================================

for e2a: compute anat alignment transformation to EPI registration base

(new anat will be intermediate, stripped, anat_ns+orig) -anat2epi -anat anat+orig
-save_skullstrip -suffix _al_junk
-epi vr_base+orig -epi_base 0
-epi_strip 3dAutomask
-volreg off -tshift off

================================== tlrc ==================================

warp anatomy to standard space

@auto_tlrc -base MNI_avg152T1+tlrc -input anat_ns+orig -no_ss

store forward transformation matrix in a text file

cat_matvec anat_ns+tlrc::WARP_DATA -I > warp.anat.Xat.1D

================================= volreg =================================

align each dset to base volume, align to anat, warp to tlrc space

verify that we have a +tlrc warp dataset

if ( ! -f anat_ns+tlrc.HEAD ) then
echo “** missing +tlrc warp dataset: anat_ns+tlrc.HEAD”

register and warp

foreach run ( $runs )
# register each volume to the base
3dvolreg -verbose -zpad 1 -base vr_base+orig
-1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run
-1Dmatrix_save mat.r$run.vr.aff12.1D

# create an all-1 dataset to mask the extents of the warp
3dcalc -overwrite -a pb02.$subj.r$run.tshift+orig -expr 1   \
       -prefix rm.epi.all1

# catenate volreg/epi2anat/tlrc xforms
cat_matvec -ONELINE                                         \
           anat_ns+tlrc::WARP_DATA -I                   \
           anat_al_junk_mat.aff12.1D -I                 \
           mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

# apply catenated xform: volreg/epi2anat/tlrc
3dAllineate -base anat_ns+tlrc                          \
            -input pb02.$subj.r$run.tshift+orig             \
            -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
            -mast_dxyz 3.5                                  \
            -prefix rm.epi.nomask.r$run

# warp the all-1 dataset for extents masking 
3dAllineate -base anat_ns+tlrc                          \
            -input rm.epi.all1+orig                         \
            -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
            -mast_dxyz 3.5 -final NN -quiet                 \
            -prefix rm.epi.1.r$run

# make an extents intersection mask of this run
3dTstat -min -prefix rm.epi.min.r$run rm.epi.1.r$run+tlrc


make a single file of registration params

cat dfile.r*.1D > dfile_rall.1D


create the extents mask: mask_epi_extents+tlrc

(this is a mask of voxels that have valid data at every TR)

(only 1 run, so just use 3dcopy to keep naming straight)

3dcopy rm.epi.min.r01+tlrc mask_epi_extents

and apply the extents mask to the EPI data

(delete any time series with missing data)

foreach run ( $runs )
3dcalc -a rm.epi.nomask.r$run+tlrc -b mask_epi_extents+tlrc
-expr ‘a*b’ -prefix pb03.$subj.r$run.volreg

warp the volreg base EPI dataset to make a final version

cat_matvec -ONELINE
anat_ns+tlrc::WARP_DATA -I
anat_al_junk_mat.aff12.1D -I > mat.basewarp.aff12.1D

3dAllineate -base anat_ns+tlrc
-input vr_base+orig
-1Dmatrix_apply mat.basewarp.aff12.1D
-mast_dxyz 3.5
-prefix final_epi_vr_base

create an anat_final dataset, aligned with stats

3dcopy anat_ns+tlrc anat_final.$subj


warp anat follower datasets (affine)

3dAllineate -source anat+orig
-master anat_final.$subj+tlrc
-final wsinc5 -1Dmatrix_apply warp.anat.Xat.1D
-prefix anat_w_skull_warped

================================== blur ==================================

blur each volume of each run

foreach run ( $runs )
3dmerge -1blur_fwhm 6.0 -doall -prefix rm.pb04.$subj.r$run.blur

# and apply extents mask, since no scale block
3dcalc -a rm.pb04.$subj.r$run.blur+tlrc -b mask_epi_extents+tlrc \
       -expr 'a*b' -prefix pb04.$subj.r$run.blur


================================== mask ==================================

create ‘full_mask’ dataset (union mask)

foreach run ( $runs )
3dAutomask -dilate 1 -prefix rm.mask_r$run pb04.$subj.r$run.blur+tlrc

create union of inputs, output type is byte

3dmask_tool -inputs rm.mask_r*+tlrc.HEAD -union -prefix full_mask.$subj

---- create subject anatomy mask, mask_anat.$subj+tlrc ----

(resampled from tlrc anat)

3dresample -master full_mask.$subj+tlrc -input anat_ns+tlrc
-prefix rm.resam.anat

convert to binary anat mask; fill gaps and holes

3dmask_tool -dilate_input 5 -5 -fill_holes -input rm.resam.anat+tlrc
-prefix mask_anat.$subj

compute overlaps between anat and EPI masks

3dABoverlap -no_automask full_mask.$subj+tlrc mask_anat.$subj+tlrc
|& tee out.mask_ae_overlap.txt

note Dice coefficient of masks, as well

3ddot -dodice full_mask.$subj+tlrc mask_anat.$subj+tlrc
|& tee out.mask_ae_dice.txt

---- create group anatomy mask, mask_group+tlrc ----

(resampled from tlrc base anat, MNI_avg152T1+tlrc)

3dresample -master full_mask.$subj+tlrc -prefix ./
-input $afnidir/MNI_avg152T1+tlrc

convert to binary group mask; fill gaps and holes

3dmask_tool -dilate_input 5 -5 -fill_holes -input
-prefix mask_group

================================ regress =================================

compute de-meaned motion parameters (for use in regression) -infile dfile_rall.1D -set_nruns 1
-demean -write motion_demean.1D

compute motion parameter derivatives (for use in regression) -infile dfile_rall.1D -set_nruns 1
-derivative -demean -write motion_deriv.1D

create censor file motion_${subj}_censor.1D, for censoring motion -infile dfile_rall.1D -set_nruns 1
-show_censor_count -censor_prev_TR
-censor_motion 0.2 motion_${subj}

combine multiple censor files

1deval -a motion_${subj}censor.1D -b outcount${subj}censor.1D
-expr “a*b” > censor

create bandpass regressors (instead of using 3dBandpass, say)

1dBport -nodata 266 2.0 -band 0.009 0.08 -invert -nozero > bandpass_rall.1D

note TRs that were not censored

set ktrs = -infile censor_${subj}_combined_2.1D \ -show_trs_uncensored encoded


run the regression analysis

3dDeconvolve -input pb04.$subj.r*.blur+tlrc.HEAD
-censor censor_${subj}_combined_2.1D
-ortvec bandpass_rall.1D bandpass
-polort 4
-num_stimts 12
-stim_file 1 motion_demean.1D’[0]’ -stim_base 1 -stim_label 1 roll_01
-stim_file 2 motion_demean.1D’[1]’ -stim_base 2 -stim_label 2 pitch_01
-stim_file 3 motion_demean.1D’[2]’ -stim_base 3 -stim_label 3 yaw_01
-stim_file 4 motion_demean.1D’[3]’ -stim_base 4 -stim_label 4 dS_01
-stim_file 5 motion_demean.1D’[4]’ -stim_base 5 -stim_label 5 dL_01
-stim_file 6 motion_demean.1D’[5]’ -stim_base 6 -stim_label 6 dP_01
-stim_file 7 motion_deriv.1D’[0]’ -stim_base 7 -stim_label 7 roll_02
-stim_file 8 motion_deriv.1D’[1]’ -stim_base 8 -stim_label 8 pitch_02
-stim_file 9 motion_deriv.1D’[2]’ -stim_base 9 -stim_label 9 yaw_02
-stim_file 10 motion_deriv.1D’[3]’ -stim_base 10 -stim_label 10 dS_02
-stim_file 11 motion_deriv.1D’[4]’ -stim_base 11 -stim_label 11 dL_02
-stim_file 12 motion_deriv.1D’[5]’ -stim_base 12 -stim_label 12 dP_02
-fout -tout -x1D X.xmat.1D -xjpeg X.jpg
-x1D_uncensored X.nocensor.xmat.1D
-fitts fitts.$subj
-errts errts.${subj}
-bucket stats.$subj

– use 3dTproject to project out regression matrix –

3dTproject -polort 0 -input pb04.$subj.r*.blur+tlrc.HEAD
-censor censor_${subj}_combined_2.1D -cenmode ZERO
-ort X.nocensor.xmat.1D -prefix errts.${subj}.tproject

if 3dDeconvolve fails, terminate the script

if ( $status != 0 ) then
echo ‘---------------------------------------’
echo ‘** 3dDeconvolve error, failing…’
echo ’ (consider the file 3dDeconvolve.err)’

display any large pairwise correlations from the X-matrix -show_cormat_warnings -infile X.xmat.1D |& tee out.cormat_warn.txt

create an all_runs dataset to match the fitts, errts, etc.

3dTcat -prefix all_runs.$subj pb04.$subj.r*.blur+tlrc.HEAD


create a temporal signal to noise ratio dataset

signal: if ‘scale’ block, mean should be 100

noise : compute standard deviation of errts

3dTstat -mean -prefix rm.signal.all all_runs.$subj+tlrc"[$ktrs]"
3dTstat -stdev -prefix rm.noise.all errts.${subj}.tproject+tlrc"[$ktrs]"
3dcalc -a rm.signal.all+tlrc
-b rm.noise.all+tlrc
-c full_mask.$subj+tlrc
-expr ‘c*a/b’ -prefix TSNR.$subj


compute and store GCOR (global correlation average)

(sum of squares of global mean of unit errts)

3dTnorm -norm2 -prefix rm.errts.unit errts.${subj}.tproject+tlrc
3dmaskave -quiet -mask full_mask.$subj+tlrc rm.errts.unit+tlrc
> gmean.errts.unit.1D
3dTstat -sos -prefix - gmean.errts.unit.1D' > out.gcor.1D
echo “-- GCOR = cat out.gcor.1D


compute correlation volume

(per voxel: average correlation across masked brain)

(now just dot product with average unit time series)

3dcalc -a rm.errts.unit+tlrc -b gmean.errts.unit.1D -expr ‘a*b’ -prefix rm.DP
3dTstat -sum -prefix corr_brain rm.DP+tlrc


compute sum of non-baseline regressors from the X-matrix

(use to get list of regressor colums)

set reg_cols = -infile X.nocensor.xmat.1D -show_indices_interest
3dTstat -sum -prefix sum_ideal.1D X.nocensor.xmat.1D"[$reg_cols]"

also, create a stimulus-only X-matrix, for easy review

1dcat X.nocensor.xmat.1D"[$reg_cols]" > X.stim.xmat.1D

============================ blur estimation =============================

compute blur estimates

touch blur_est.$subj.1D # start with empty file

– estimate blur for each run in errts –

touch blur.errts.1D

restrict to uncensored TRs, per run

foreach run ( $runs )
set trs = -infile X.xmat.1D -show_trs_uncensored encoded \ -show_trs_run $run
if ( $trs == “” ) continue
3dFWHMx -detrend -mask full_mask.$subj+tlrc
errts.${subj}.tproject+tlrc"[$trs]" >> blur.errts.1D

compute average blur and append

set blurs = ( cat blur.errts.1D )
echo average errts blurs: $blurs
echo “$blurs # errts blur estimates” >> blur_est.$subj.1D

================== auto block: generate review scripts ===================

generate a review script for the unprocessed EPI data -script @epi_review.$subj
-dsets pb00.$subj.r*.tcat+orig.HEAD

generate scripts to review single subject results

(try with defaults, but do not allow bad exit status) -mot_limit 0.2 -out_limit 0.25
-errts_dset errts.${subj}.tproject+tlrc.HEAD -exit0

========================== auto block: finalize ==========================

remove temporary files

\rm -f rm.*

if the basic subject review script is here, run it

(want this to be the last text output)

if ( -e @ss_review_basic ) ./@ss_review_basic |& tee out.ss_review.$subj.txt

return to parent directory

cd …

echo “execution finished: date


script generated by the command: -subj_id NA8831 -script proc.NA8831 -scr_overwrite -blocks \

despike tshift align tlrc volreg blur mask regress -copy_anat \

/home/mjm693/Desktop/NA8831a/anat8831.nii.gz -tcat_remove_first_trs 4 \

-dsets /home/mjm693/Desktop/NA8831a/rs8831.nii.gz -tlrc_base \

MNI_avg152T1+tlrc -volreg_align_to first -volreg_align_e2a \

-volreg_tlrc_warp -blur_size 6.0 -regress_censor_motion 0.2 \

-regress_censor_outliers 0.25 -regress_bandpass 0.009 0.08 \

-regress_apply_mot_types demean deriv -regress_est_blur_errts \

-regress_run_clustsim no


I think posting your command would be more illustrative. Listing the version numbers of what you are testing would also be useful.

Note that probably multiple things, perhaps even some defaults, have changed since 2016. If you post your command, it might be possible to get a better sense of that.

Also, you should state explicitly what you are doing to create your results-- t-tests, seedbased correlation maps, something else? Stating that “the results appear random [now]” is not enough to judge what undrelying differences might have changed.