afni_ci_test_data/AFNI_data6/FT_analysis/s11.proc.FT

#!/bin/tcsh -xef

echo "auto-generated by afni_proc.py, Mon Aug 22 15:47:09 2016"
echo "(version 5.01, August 22, 2016)"
echo "execution started: `date`"

# execute via : 
#   tcsh -xef proc.FT |& tee output.proc.FT

# =========================== 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 1 Dec 2015
if ( $status ) then
    echo "** this script requires newer AFNI binaries (than 1 Dec 2015)"
    echo "   (consider: @update.afni.binaries -defaults)"
    exit
endif

# the user may specify a single subject to run with
if ( $#argv > 0 ) then
    set subj = $argv[1]
else
    set subj = FT
endif

# assign output directory name
set output_dir = $subj.results

# 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 3`)

# create results and stimuli directories
mkdir $output_dir
mkdir $output_dir/stimuli

# copy stim files into stimulus directory
cp FT/AV1_vis.txt FT/AV2_aud.txt $output_dir/stimuli

# copy anatomy to results dir
3dcopy FT/FT_anat+orig $output_dir/FT_anat

# ============================ auto block: tcat ============================
# apply 3dTcat to copy input dsets to results dir, while
# removing the first 2 TRs
3dTcat -prefix $output_dir/pb00.$subj.r01.tcat FT/FT_epi_r1+orig'[2..$]'
3dTcat -prefix $output_dir/pb00.$subj.r02.tcat FT/FT_epi_r2+orig'[2..$]'
3dTcat -prefix $output_dir/pb00.$subj.r03.tcat FT/FT_epi_r3+orig'[2..$]'

# and make note of repetitions (TRs) per run
set tr_counts = ( 150 150 150 )

# -------------------------------------------------------
# enter the results directory (can begin processing data)
cd $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 3 -legendre                     \
                pb00.$subj.r$run.tcat+orig > outcount.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

# ================================= tshift =================================
# time shift data so all slice timing is the same 
foreach run ( $runs )
    3dTshift -tzero 0 -quintic -prefix pb01.$subj.r$run.tshift \
             pb00.$subj.r$run.tcat+orig
end

# --------------------------------
# extract volreg registration base
3dbucket -prefix vr_base pb01.$subj.r01.tshift+orig"[2]"

# ================================= volreg =================================
# align each dset to base volume
foreach run ( $runs )
    # register each volume to the base
    3dvolreg -verbose -zpad 1 -base vr_base+orig                    \
             -1Dfile dfile.r$run.1D -prefix pb02.$subj.r$run.volreg \
             -cubic                                                 \
             pb01.$subj.r$run.tshift+orig
end

# make a single file of registration params
cat dfile.r*.1D > dfile_rall.1D

# compute motion magnitude time series: the Euclidean norm
# (sqrt(sum squares)) of the motion parameter derivatives
1d_tool.py -infile dfile_rall.1D -set_nruns 3                       \
           -derivative  -collapse_cols euclidean_norm               \
           -write motion_${subj}_enorm.1D

# create an anat_final dataset, aligned with stats
3dcopy FT_anat+orig anat_final.$subj

# ================================== blur ==================================
# blur each volume of each run
foreach run ( $runs )
    3dmerge -1blur_fwhm 4.0 -doall -prefix pb03.$subj.r$run.blur \
            pb02.$subj.r$run.volreg+orig
end

# ================================== mask ==================================
# create 'full_mask' dataset (union mask)
foreach run ( $runs )
    3dAutomask -dilate 1 -prefix rm.mask_r$run pb03.$subj.r$run.blur+orig
end

# create union of inputs, output type is byte
3dmask_tool -inputs rm.mask_r*+orig.HEAD -union -prefix full_mask.$subj

# ================================= scale ==================================
# scale each voxel time series to have a mean of 100
# (be sure no negatives creep in)
# (subject to a range of [0,200])
foreach run ( $runs )
    3dTstat -prefix rm.mean_r$run pb03.$subj.r$run.blur+orig
    3dcalc -a pb03.$subj.r$run.blur+orig -b rm.mean_r$run+orig \
           -expr 'min(200, a/b*100)*step(a)*step(b)'           \
           -prefix pb04.$subj.r$run.scale
end

# ================================ regress =================================

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

# compute motion parameter derivatives (just to have)
1d_tool.py -infile dfile_rall.1D -set_nruns 3                            \
           -derivative -demean -write motion_deriv.1D

# ------------------------------
# run the regression analysis
3dDeconvolve -input pb04.$subj.r*.scale+orig.HEAD                        \
    -polort 3                                                            \
    -num_stimts 8                                                        \
    -stim_times 1 stimuli/AV1_vis.txt 'BLOCK(20,1)'                      \
    -stim_label 1 Vrel                                                   \
    -stim_times 2 stimuli/AV2_aud.txt 'BLOCK(20,1)'                      \
    -stim_label 2 Arel                                                   \
    -stim_file 3 motion_demean.1D'[0]' -stim_base 3 -stim_label 3 roll   \
    -stim_file 4 motion_demean.1D'[1]' -stim_base 4 -stim_label 4 pitch  \
    -stim_file 5 motion_demean.1D'[2]' -stim_base 5 -stim_label 5 yaw    \
    -stim_file 6 motion_demean.1D'[3]' -stim_base 6 -stim_label 6 dS     \
    -stim_file 7 motion_demean.1D'[4]' -stim_base 7 -stim_label 7 dL     \
    -stim_file 8 motion_demean.1D'[5]' -stim_base 8 -stim_label 8 dP     \
    -gltsym 'SYM: Vrel -Arel'                                            \
    -glt_label 1 V-A                                                     \
    -fout -tout -x1D X.xmat.1D -xjpeg X.jpg                              \
    -fitts fitts.$subj                                                   \
    -errts errts.${subj}                                                 \
    -bucket stats.$subj


# 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*.scale+orig.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+orig
3dTstat -stdev -prefix rm.noise.all errts.${subj}+orig
3dcalc -a rm.signal.all+orig                                             \
       -b rm.noise.all+orig                                              \
       -c full_mask.$subj+orig                                           \
       -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}+orig
3dmaskave -quiet -mask full_mask.$subj+orig rm.errts.unit+orig           \
          > 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+orig -b gmean.errts.unit.1D -expr 'a*b' -prefix rm.DP
3dTstat -sum -prefix corr_brain rm.DP+orig

# create ideal files for fixed response stim types
1dcat X.xmat.1D'[12]' > ideal_Vrel.1D
1dcat X.xmat.1D'[13]' > ideal_Arel.1D

# --------------------------------------------------------
# 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.xmat.1D -show_indices_interest`
3dTstat -sum -prefix sum_ideal.1D X.xmat.1D"[$reg_cols]"

# also, create a stimulus-only X-matrix, for easy review
1dcat X.xmat.1D"[$reg_cols]" > X.stim.xmat.1D

# ============================ blur estimation =============================
# compute blur estimates
touch blur_est.$subj.1D   # start with empty file

# create directory for ACF curve files
mkdir files_ACF

# -- 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+orig                          \
            -ACF files_ACF/out.3dFWHMx.ACF.errts.r$run.1D                \
            errts.${subj}+orig"[$trs]" >> blur.errts.1D
end

# compute average FWHM blur (from every other row) and append
set blurs = ( `3dTstat -mean -prefix - blur.errts.1D'{0..$(2)}'\'` )
echo average errts FWHM blurs: $blurs
echo "$blurs   # errts FWHM blur estimates" >> blur_est.$subj.1D

# compute average ACF blur (from every other row) and append
set blurs = ( `3dTstat -mean -prefix - blur.errts.1D'{1..$(2)}'\'` )
echo average errts ACF blurs: $blurs
echo "$blurs   # errts ACF blur estimates" >> blur_est.$subj.1D


# add 3dClustSim results as attributes to any stats dset
mkdir files_ClustSim

# run Monte Carlo simulations using method 'ACF'
set params = ( `grep ACF blur_est.$subj.1D | tail -n 1` )
3dClustSim -both -mask full_mask.$subj+orig -acf $params[1-3]            \
           -cmd 3dClustSim.ACF.cmd -prefix files_ClustSim/ClustSim.ACF

# run 3drefit to attach 3dClustSim results to stats
set cmd = ( `cat 3dClustSim.ACF.cmd` )
$cmd stats.$subj+orig


# ================== 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 -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:
#
# afni_proc.py -subj_id FT -dsets FT/FT_epi_r1+orig.HEAD                     \
#     FT/FT_epi_r2+orig.HEAD FT/FT_epi_r3+orig.HEAD -copy_anat               \
#     FT/FT_anat+orig -tcat_remove_first_trs 2 -regress_stim_times           \
#     FT/AV1_vis.txt FT/AV2_aud.txt -regress_stim_labels Vrel Arel           \
#     -regress_basis 'BLOCK(20,1)' -regress_est_blur_errts -regress_opts_3dD \
#     -gltsym 'SYM: Vrel -Arel' -glt_label 1 V-A