Hello,
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 (MacRestproc.sh; 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!
Mike
–MacRestproc.sh–
#!/bin/tcsh -xef
echo “auto-generated by afni_proc.py, 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)"
exit
endif
accept arguments from rsRegress.sh
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]
#else
set subj = NA8831
#endif
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
exit
endif
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
endif
end
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
pb00.$subj.r$run.tcat+orig
end
================================= tshift =================================
time shift data so all slice timing is the same
foreach run ( $runs )
3dTshift -tzero 0 -quintic -prefix pb02.$subj.r$run.tshift
pb01.$subj.r$run.despike+orig
end
--------------------------------
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)
align_epi_anat.py -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”
exit
endif
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
-cubic
-1Dmatrix_save mat.r$run.vr.aff12.1D
pb02.$subj.r$run.tshift+orig
# 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
end
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
end
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
pb03.$subj.r$run.volreg+tlrc
# 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
end
================================== mask ==================================
create ‘full_mask’ dataset (union mask)
foreach run ( $runs )
3dAutomask -dilate 1 -prefix rm.mask_r$run pb04.$subj.r$run.blur+tlrc
end
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 ./rm.resam.group
-input $afnidir/MNI_avg152T1+tlrc
convert to binary group mask; fill gaps and holes
3dmask_tool -dilate_input 5 -5 -fill_holes -input rm.resam.group+tlrc
-prefix mask_group
================================ regress =================================
compute de-meaned motion parameters (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_nruns 1
-demean -write motion_demean.1D
compute motion parameter derivatives (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_nruns 1
-derivative -demean -write motion_deriv.1D
create censor file motion_${subj}_censor.1D, for censoring motion
1d_tool.py -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${subj}_combined_2.1D
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 = 1d_tool.py -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}
-x1D_stop
-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)’
exit
endif
display any large pairwise correlations from the X-matrix
1d_tool.py -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 1d_tool.py to get list of regressor colums)
set reg_cols = 1d_tool.py -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 = 1d_tool.py -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
end
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
gen_epi_review.py -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)
gen_ss_review_scripts.py -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”