-
Notifications
You must be signed in to change notification settings - Fork 13
/
Copy pathproc_Flanker.sh
491 lines (396 loc) · 20.7 KB
/
proc_Flanker.sh
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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
#!/bin/tcsh -xef
echo "auto-generated by afni_proc.py, Fri Sep 28 14:08:08 2018"
echo "(version 6.14, May 25, 2018)"
echo "execution started: `date`"
# execute via :
# tcsh -xef proc.s01 |& tee output.proc.s01
# =========================== 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 3 May 2018
if ( $status ) then
echo "** this script requires newer AFNI binaries (than 3 May 2018)"
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 = s01
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 2`)
# create results and stimuli directories
mkdir $output_dir
mkdir $output_dir/stimuli
# copy stim files into stimulus directory
cp \
$PWD/{$subj}/func/congruent.1D \
$PWD/{$subj}/func/incongruent.1D \
$output_dir/stimuli
# copy anatomy to results dir
3dcopy \
$PWD/{$subj}/anat/{$subj}_T1w.nii.gz \
$output_dir/{$subj}_T1w
# ============================ auto block: tcat ============================
# apply 3dTcat to copy input dsets to results dir,
# while removing the first 0 TRs
3dTcat -prefix $output_dir/pb00.$subj.r01.tcat \
$PWD/{$subj}/func/{$subj}_task-flanker_run-1_bold.nii.gz'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r02.tcat \
$PWD/{$subj}/func/{$subj}_task-flanker_run-2_bold.nii.gz'[0..$]'
# and make note of repetitions (TRs) per run
set tr_counts = ( 146 146 )
# -------------------------------------------------------
# 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 2 -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
# get run number and TR index for minimum outlier volume
set minindex = `3dTstat -argmin -prefix - outcount_rall.1D\'`
set ovals = ( `1d_tool.py -set_run_lengths $tr_counts \
-index_to_run_tr $minindex` )
# save run and TR indices for extraction of vr_base_min_outlier
set minoutrun = $ovals[1]
set minouttr = $ovals[2]
echo "min outlier: run $minoutrun, TR $minouttr" | tee out.min_outlier.txt
# ================================= 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_min_outlier \
pb01.$subj.r$minoutrun.tshift+orig"[$minouttr]"
# ================================= align ==================================
# for e2a: compute anat alignment transformation to EPI registration base
# (new anat will be intermediate, stripped, {$subj}_T1w_ns+orig)
align_epi_anat.py -anat2epi -anat {$subj}_T1w+orig \
-save_skullstrip -suffix _al_junk \
-epi vr_base_min_outlier+orig -epi_base 0 \
-epi_strip 3dSkullStrip \
-giant_move \
-volreg off -tshift off
# ================================== tlrc ==================================
# warp anatomy to standard space
@auto_tlrc -base MNI_avg152T1+tlrc -input {$subj}_T1w_ns+orig -init_xform AUTO_CENTER -no_ss
# store forward transformation matrix in a text file
cat_matvec {$subj}_T1w_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 {$subj}_T1w_ns+tlrc.HEAD ) then
echo "** missing +tlrc warp dataset: {$subj}_T1w_ns+tlrc.HEAD"
exit
endif
# register and warp
foreach run ( $runs )
# register each volume to the base image
3dvolreg -verbose -zpad 1 -base vr_base_min_outlier+orig \
-1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run \
-cubic \
-1Dmatrix_save mat.r$run.vr.aff12.1D \
pb01.$subj.r$run.tshift+orig
# create an all-1 dataset to mask the extents of the warp
3dcalc -overwrite -a pb01.$subj.r$run.tshift+orig -expr 1 \
-prefix rm.epi.all1
# catenate volreg/epi2anat/tlrc xforms
cat_matvec -ONELINE \
{$subj}_T1w_ns+tlrc::WARP_DATA -I \
{$subj}_T1w_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 {$subj}_T1w_ns+tlrc \
-input pb01.$subj.r$run.tshift+orig \
-1Dmatrix_apply mat.r$run.warp.aff12.1D \
-mast_dxyz 3 \
-prefix rm.epi.nomask.r$run
# warp the all-1 dataset for extents masking
3dAllineate -base {$subj}_T1w_ns+tlrc \
-input rm.epi.all1+orig \
-1Dmatrix_apply mat.r$run.warp.aff12.1D \
-mast_dxyz 3 -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)
3dMean -datum short -prefix rm.epi.mean rm.epi.min.r*.HEAD
3dcalc -a rm.epi.mean+tlrc -expr 'step(a-0.999)' -prefix 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 pb02.$subj.r$run.volreg
end
# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE \
{$subj}_T1w_ns+tlrc::WARP_DATA -I \
{$subj}_T1w_al_junk_mat.aff12.1D -I > mat.basewarp.aff12.1D
3dAllineate -base {$subj}_T1w_ns+tlrc \
-input vr_base_min_outlier+orig \
-1Dmatrix_apply mat.basewarp.aff12.1D \
-mast_dxyz 3 \
-prefix final_epi_vr_base_min_outlier
# create an anat_final dataset, aligned with stats
3dcopy {$subj}_T1w_ns+tlrc anat_final.$subj
# record final registration costs
3dAllineate -base final_epi_vr_base_min_outlier+tlrc -allcostX \
-input anat_final.$subj+tlrc |& tee out.allcostX.txt
# -----------------------------------------
# warp anat follower datasets (affine)
3dAllineate -source {$subj}_T1w+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 4.0 -doall -prefix pb03.$subj.r$run.blur \
pb02.$subj.r$run.volreg+tlrc
end
# ================================== mask ==================================
# create 'full_mask' dataset (union mask)
foreach run ( $runs )
3dAutomask -prefix rm.mask_r$run pb03.$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 {$subj}_T1w_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 tighter EPI mask by intersecting with anat mask
3dmask_tool -input full_mask.$subj+tlrc mask_anat.$subj+tlrc \
-inter -prefix mask_epi_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 /Users/${USER}/abin/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
# ================================= 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+tlrc
3dcalc -a pb03.$subj.r$run.blur+tlrc -b rm.mean_r$run+tlrc \
-c mask_epi_extents+tlrc \
-expr 'c * 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 2 \
-demean -write motion_demean.1D
# compute motion parameter derivatives (just to have)
1d_tool.py -infile dfile_rall.1D -set_nruns 2 \
-derivative -demean -write motion_deriv.1D
# convert motion parameters for per-run regression
1d_tool.py -infile motion_demean.1D -set_nruns 2 \
-split_into_pad_runs mot_demean
# create censor file motion_${subj}_censor.1D, for censoring motion
1d_tool.py -infile dfile_rall.1D -set_nruns 2 \
-show_censor_count -censor_prev_TR \
-censor_motion 0.3 motion_${subj}
# note TRs that were not censored
set ktrs = `1d_tool.py -infile motion_${subj}_censor.1D \
-show_trs_uncensored encoded`
# ------------------------------
# run the regression analysis
3dDeconvolve -input pb04.$subj.r*.scale+tlrc.HEAD \
-censor motion_${subj}_censor.1D \
-mask mask_group+tlrc \
-polort 2 \
-num_stimts 14 \
-stim_times 1 stimuli/congruent.1D 'BLOCK(2,1)' \
-stim_label 1 congruent \
-stim_times 2 stimuli/incongruent.1D 'BLOCK(2,1)' \
-stim_label 2 incongruent \
-stim_file 3 mot_demean.r01.1D'[0]' -stim_base 3 -stim_label 3 roll_01 \
-stim_file 4 mot_demean.r01.1D'[1]' -stim_base 4 -stim_label 4 pitch_01 \
-stim_file 5 mot_demean.r01.1D'[2]' -stim_base 5 -stim_label 5 yaw_01 \
-stim_file 6 mot_demean.r01.1D'[3]' -stim_base 6 -stim_label 6 dS_01 \
-stim_file 7 mot_demean.r01.1D'[4]' -stim_base 7 -stim_label 7 dL_01 \
-stim_file 8 mot_demean.r01.1D'[5]' -stim_base 8 -stim_label 8 dP_01 \
-stim_file 9 mot_demean.r02.1D'[0]' -stim_base 9 -stim_label 9 roll_02 \
-stim_file 10 mot_demean.r02.1D'[1]' -stim_base 10 -stim_label 10 \
pitch_02 \
-stim_file 11 mot_demean.r02.1D'[2]' -stim_base 11 -stim_label 11 yaw_02 \
-stim_file 12 mot_demean.r02.1D'[3]' -stim_base 12 -stim_label 12 dS_02 \
-stim_file 13 mot_demean.r02.1D'[4]' -stim_base 13 -stim_label 13 dL_02 \
-stim_file 14 mot_demean.r02.1D'[5]' -stim_base 14 -stim_label 14 dP_02 \
-jobs 8 \
-gltsym 'SYM: incongruent -congruent' \
-glt_label 1 incongruent-congruent \
-gltsym 'SYM: congruent -incongruent' \
-glt_label 2 congruent-incongruent \
-fout -tout -x1D X.xmat.1D -xjpeg X.jpg \
-x1D_uncensored X.nocensor.xmat.1D \
-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
# -- execute the 3dREMLfit script, written by 3dDeconvolve --
tcsh -x stats.REML_cmd
# create an all_runs dataset to match the fitts, errts, etc.
3dTcat -prefix all_runs.$subj pb04.$subj.r*.scale+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}+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}+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
# create ideal files for fixed response stim types
1dcat X.nocensor.xmat.1D'[6]' > ideal_congruent.1D
1dcat X.nocensor.xmat.1D'[7]' > ideal_incongruent.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.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
# create directory for ACF curve files
mkdir files_ACF
# -- estimate blur for each run in epits --
touch blur.epits.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 \
-ACF files_ACF/out.3dFWHMx.ACF.epits.r$run.1D \
all_runs.$subj+tlrc"[$trs]" >> blur.epits.1D
end
# compute average FWHM blur (from every other row) and append
set blurs = ( `3dTstat -mean -prefix - blur.epits.1D'{0..$(2)}'\'` )
echo average epits FWHM blurs: $blurs
echo "$blurs # epits FWHM blur estimates" >> blur_est.$subj.1D
# compute average ACF blur (from every other row) and append
set blurs = ( `3dTstat -mean -prefix - blur.epits.1D'{1..$(2)}'\'` )
echo average epits ACF blurs: $blurs
echo "$blurs # epits ACF blur estimates" >> blur_est.$subj.1D
# -- 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 \
-ACF files_ACF/out.3dFWHMx.ACF.errts.r$run.1D \
errts.${subj}+tlrc"[$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
# ================== 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.3 -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 ..
mv {$subj}.results {$subj}
echo "execution finished: `date`"
# ==========================================================================
# script generated by the command:
#
# afni_proc.py -subj_id s01 -script proc.s01 -scr_overwrite -blocks tshift \
# align tlrc volreg blur mask scale regress -copy_anat \
# $PWD/{$subj}/anat/{$subj}_T1w.nii.gz \
# -dsets \
# $PWD/{$subj}/func/{$subj}_task-flanker_run-1_bold.nii.gz \
# $PWD/{$subj}/func/{$subj}_task-flanker_run-2_bold.nii.gz \
# -tcat_remove_first_trs 0 -align_opts_aea -giant_move -tlrc_base \
# MNI_avg152T1+tlrc -volreg_align_to MIN_OUTLIER -volreg_align_e2a \
# -volreg_tlrc_warp -blur_size 4.0 -regress_stim_times \
# $PWD/{$subj}/func/congruent.1D \
# $PWD/{$subj}/func/incongruent.1D \
# -regress_stim_labels congruent incongruent -regress_basis 'BLOCK(2,1)' \
# -regress_censor_motion 0.3 -regress_motion_per_run -regress_opts_3dD \
# -jobs 8 -regress_make_ideal_sum sum_ideal.1D -regress_est_blur_epits \
# -regress_est_blur_errts -regress_run_clustsim no