#!/usr/bin/env bash
# tat2: time-averaged T2*
# glue for 3dcalc 3dTcat and 3dTstat (+ 3dNotes)
# to compute time average T2
# default values
t2_inv=0 # before volnorm (3dcalc, 3dTstat)
MASK_REL="subject_mask.nii.gz" # 3dROIstats: use to calculate volnorm.1D
MASK="" # "
MKMASK=0 # "
volnorm_opt="-nzmedian" # 3dcalc: run_tat2.nii.gz (vol normalized)
SCALE=1000 # "
vox_scale=1 # "
use_zscore=0 # "
use_ln=0 # " (added 20240719)
MAXVOLS=-1 # " truncate how many volumes are in run_tat2.nii.gz
censor_rel="" # 3dTcat: truncate input
MAXVOLSTOTAL=-1 # " truncate how many volumes are ultimately averaged.
IDX_SAMPLE_METHOD="first" # how to sample indexes (first, last, random) when MAXVOLS*
tnorm_opt="-nzmean" # 3dTstat: final output (time normalized)
OUTPUT=tat2star.nii.gz
TMPLOCATION=/tmp
CLEAN=1
TMPD="" # will be folder specific to single invocation of tat2. removed at the end
WRITE_JSON_LOG=1
# all the inputs (possibly globs) we want to work on
declare -a GLOBFILES
filedesc_whichvol=""
TAT2_VER="20240719-calc_ln+BUGFIX-no_vol"
# 20201116
# back to mean default
# add
# across-TR median (instead of mean)
# whole brain normalization to mean (instead of sum)
# and whole brain z-score normalization instead
##
# initially from
# https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0016093
# author correspondence
# > The normalization is within the brain for each time point.
# > Not all voxels in the volume, but only those that belong to the brain mask.
# > Then we normalize so that the sum of all brain voxels is some fixed number,
# > e.g., 10000. The number doesn't really matter.
# also see "Relative Concentration of Brain Iron (rcFe)"
# https://www.biorxiv.org/content/biorxiv/early/2019/03/16/579763.full.pdf
usage(){
cat <<-HEREDOC
USAGE
tat2 '/paths/to/*/t2.nii.gz' [ -mask_rel $MASK_REL | -mask_rel s/_bold.nii.gz/_bold_mask.nii.gz/ | -mask /the/only/mask.nii.gz | -mkmask] [-output $OUTPUT] [-scale $SCALE] [-censor_rel relative/censor.1D | -censor_rel s/epi/censor/] [-median_time|-mean_time] [-median_vol|-mean_vol|-no_vol] [-calc_zscore|-calc_ln] [-no_voxscale] [-inverse] [-noclean] [-verbose] [-tmp $TMPLOCATION]
SYNOPSIS
calculates time averaged T2* like:
1) 3dROIstats > volnorm.1D
# tune with -{median,mean,no}_vol, and -mask, -mkmask, or -mask_rel
2) 3dcalc -x input_1 -m volnorm_1.1D -expr "(x/m)*SCALE/numvox #rep. f.ea. input
# tune with -no_voxscale, -scale, -calc_zscore, -calc_ln
# numvox derived from mask (but output not masked)
3) 3dTstat allnormed_concated.nii.gz
# tune with -{median,mean}_time
final '$OUTPUT' will be in the cwd
OPTIONS (in order of relevance)
-censor_rel FILE
specify the motion censor file. single column no header file. nLines = nVols. 0=exclude, 1=keep.
either a filename or regexp
as a filename in input nifti directory (resolves to \$(dirname \$input)/\$censor_rel) OR
as a perl regexp run on \$input. like 's/(rest-\\d)_bold.nii.gz/\\1_cesnor.1d/'
despite switch's name, FILE also be an absolute path (start with '/')
-mask_rel volnorm mask is relative to the input nifit (default: subject_mask.nii.gz)
like censor_rel, if given 's/find/replace/' perl regexp
will replace 'find' in the input.nii.gz with 'replace'
-mask use mask to create volnorm.1D (-median_vol or -mean_vol)
default none, defer to mask_rel
-mkmask ensures the number of voxels used to normalize is the same
across runs by create a coverage mask from all inputs
-nzmedian/nzmean and '-count -non-zero' could maybe get away with no mask,
but option is not provided
-maxvols NVOL
limit each run to NVOL volumes. run after and only after censor_rel.
useful when groups (e.g. young, old) have disproportonate censored volumes
-maxvolstotal NVOL
in final combined timeseries of all runs, take only the first NVOL volumes
-sample_method {first,last,random}
how to sample indexes when using maxvols/maxvolstotal. take first n, last n, or random n
-median_time or -mean_time
set how voxels are collapsed across time (3dTstat)
time norm default: '$tnorm_opt'
(history: mean until 20201117)
-median_vol or -mean_vol or -no_vol
set value used to scale across volume (each TR, 3dROIstats)
vol norm default: '$volnorm_opt'
(history: median [20210302], mean [20201116], median [20201016]. no_vol added 20210921 but *x instead of *m. fixed 20240719)
-calc_zscore or -calc_log
for 3dcalc per volume voxelwise normalization:
instead of 'x/m'
-calc_zscore '(x-m)/s'
-calc_log '-ln(x/m)'
-no_voxscale disables scaling. (SCALE=1, not normalized by number of nzvoxels)
-inverse create 1/t2* and use as initial input
-noclean keeps tmp folder around, useful for debugging
-tmp DIR sets the directory to use for temporary calculations to DIR.
useful if forking enough jobs to overwhelm /tmp or using with -noclean
REFS
https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0016093
https://www.biorxiv.org/content/biorxiv/early/2019/03/16/579763.full.pdf
HEREDOC
exit 1
}
err() { echo -e "$@" >&2; exit 1; }
msg() { echo -e "# [$(date)] $@"; }
find_rel_file(){
local input="$1"; shift
local rel_str="$1"; shift
local rel_file="$(dirname "$input")/$rel_str"
# if given an absolute path. try that
[ ! -r "$rel_file" -a ${rel_str:0:1} == "/" ] &&
rel_file="$rel_str"
# try regexp pattern replace if rel_str looks like s/
[ ! -r "$rel_file" -a ${rel_str:0:2} == "s/" ] &&
rel_file="$(echo "$input" | perl -pe "$rel_str")"
# make sure we didn't just find the input file
if [ "$rel_file" == "$input" ]; then
err "ERR: relative file subsitition failed '$rel_file' matches inputfile! rel '$rel_str' for '$input'"
fi
# and make sure we actually found a file
if [ ! -r "$rel_file" ]; then
err "ERR: file '$rel_file' (from '$rel_str') DNE!"
fi
echo $rel_file
}
# pull out good index (censor=1 means keep). zero based
where1csv(){ perl -lne 'push @i, $.-1 if /^\s?1/; END{print join ",", @i;}' "$@"; }
firstn_csv() { perl -F, -slane 'print $n<=0?$_:(join ",", grep {!/^$/} @F[0..$n-1])' -- -n=$1; }
idx_shuffle() {
local method="$1"; shift
local cmd
case "$method" in
first) cmd=cat;;
last) cmd=tac;;
random) cmd=shuf;;
esac
tr , '\n' | $cmd |tr '\n' ,|sed s/,$//
}
collapse_seq_idx(){
# for very large volumes, the argument list is too large!
# > ERROR: (The longest filename allowed in AFNI is 5095 characters)
perl -F, -slane '
BEGIN{
$p=-2;
@IDX=();
$n="";
sub ndot{ push(@IDX, "..@_") if @_[0];}
}
for $i (@F){
if($i==$p+1){
$n=$i;
} else{
ndot $n; $n="";
}
push(@IDX, $i) unless $n;
$p=$i;
}
END{
ndot($n);
$_=join(",",@IDX);
s/,(\.\.)/\1/g;
print
}'
}
parse_args(){
# keep args around for 3dNotes history
allargs="$*"
# read in any arguments/paramaters
while [ $# -gt 0 ]; do
case $1 in
-mask) MASK="$2"; shift 2;;
-mask_rel) MASK_REL="$2"; shift 2;;
-output) OUTPUT="$2"; shift 2;;
-scale) SCALE="$2"; shift 2;;
-censor_rel) censor_rel="$2"; shift 2;;
-tmp) TMPLOCATION="$2"; shift 2;;
-maxvols) MAXVOLS="$2"; shift 2;;
-maxvolstotal) MAXVOLSTOTAL="$2"; shift 2;;
-sample_method) IDX_SAMPLE_METHOD=$2; shift 2;;
-mkmask) MKMASK=1; MASK_REL=""; shift;;
-median_vol) volnorm_opt="-nzmedian"; shift;;
-mean_vol) volnorm_opt="-nzmean"; shift;;
-no_vol) volnorm_opt="none"; shift;;
-calc_zscore|-zscore_vol) use_zscore=1; shift;;
-calc_ln) use_ln=1; shift;; # added 20240719
-mean_time) tnorm_opt="-nzmean"; shift;;
-median_time)tnorm_opt="-nzmedian"; shift;;
-inverse) t2_inv=1; shift;;
-no_voxscale)vox_scale=0; shift;;
-noclean) CLEAN=0; shift;;
-verbose) set -x; shift;;
-h*) usage;;
-*) echo "unknown option '$1'"; usage;;
*) GLOBFILES+=("$1"); shift;;
esac
done
return 0
}
args_are_sane(){
# help if nothing given
[ $# -eq 0 ] && usage
parse_args "$@"
! [[ "$IDX_SAMPLE_METHOD" =~ first|last|random ]] &&
err "-sample_method must be first, last, or random; not '$IDX_SAMPLE_METHOD'"
[ "$IDX_SAMPLE_METHOD" != first -a $MAXVOLSTOTAL = -1 -a $MAXVOLS = -1 ] &&
err "'-sample_method $IDX_SAMPLE_METHOD' doesn't makes sense without -maxvols*"
local no_vol=0
[[ $volnorm_opt == "none" ]] && no_vol=1
if [[ "$use_zscore$use_ln$no_vol" =~ 1.*1 ]]; then
err "cannot combine -calc_zscore or -calc_ln or -no_vol"
fi
if [[ $use_zscore$use_ln =~ 1 && $vox_scale -eq 1 ]]; then
err "must use -no_voxscale with -calc_zscore or -calc_ln"
fi
echo "#files/globs: ${#GLOBFILES[@]}"
# need to have at least one file to average
[ -z "${GLOBFILES[*]}" ] && usage
# how many files do we have
nfiles=$(find -L ${GLOBFILES[@]} -maxdepth 0 | wc -l)
[ $nfiles -eq 0 ] && err "no files match input GLOBFILESs: ${GLOBFILES[@]}"
[ $nfiles -eq 1 ] && echo "WARNING: only one file matches '${GLOBFILES[@]}'. expected all (>1) runs"
[ $nfiles -gt 10 ] && echo "WARNING: running on $nfiles epi files! Are you sure you don't want to run one visit at a time?"
[ -r "$OUTPUT" ] && echo "# have $OUTPUT; rm $OUTPUT # to redo" && exit 0
[ "$MAXVOLS" == "1" ] && echo "-maxvols cannot be 1. 3dcalc doesn't like applying a single value 1D file to a 4D dataset" && exit 1
return 0
}
update_with_censor(){
declare -g input filedesc_whichvol censor_file
filedesc_whichvol=""
readonly censor_rel
local idxs nkeep tat2inputfile
# update 'input' and 'runoutput'
# censor input
if [ -n "$censor_rel" ]; then
censor_file=$(find_rel_file "$input" "$censor_rel")
# if MAXVOLS is default "-1", firstn_cvs does nothing
idxs=$(where1csv $censor_file| idx_shuffle $IDX_SAMPLE_METHOD | firstn_csv $MAXVOLS)
nkeep=$(echo $idxs| tr ',' '\n' |wc -l)
filedesc_whichvol=_keep-${nkeep}
msg "censor: using $nkeep/$(3dinfo -nt $input) timepoints. (last vol: ${idxs/*,}) [$censor_file]"
if [ $MAXVOLS -gt 0 ]; then
[ $nkeep -ne $MAXVOLS ] &&
msg "WARNING:run$cnt: have $nkeep != maxvols ($MAXVOLS). censored more volumes than min required?!"
# nkeep next used for file names.
filedesc_whichvol=_lastidx-${idxs/*,/}
fi
# truncate input
# previous used 3dTcat to truncation. but 3dROIstats and 3dcalc understand '[$index]' notation
# so no need to make another (potentially large) file
tat2inputfile=$input
input="$input[$idxs]"
# > ** ERROR: (The longest filename allowed in AFNI is 5095 characters)
# NB will be weird (do nothing?) when IDX_SAMPLE_METHOD is not "first"
#
# should remove conditional (>5k) so if this caues a bug
# it will show up everywhere and will be tested (?)
[ ${#input} -ge 5095 ] &&
input="$tat2inputfile[$(echo $idxs | collapse_seq_idx)]"
fi
return 0
}
args_to_3dcalc_expr(){
readonly use_zscore use_ln vox_scale volnorm_opt numvox
## what does 3dcalc's expression look like
# default is x/m*$SCALE/$numvox
# 4 possibilities. product of
# '(x-m)/s' vs 'x/m' AND '$SCALE/$numvox' vs '1'
# zscore vs not * vox_scale vs not
calc_expr="(x/m)"
calc_scale="$SCALE/$numvox"
if [ $use_zscore -eq 1 ]; then
#sd_input="-s $vol_sd" # handled elsehwere to avoid sideffects here
calc_expr="(x-m)/s"
elif [ $use_ln -eq 1 ]; then
calc_expr="-1*log(x/m)"
fi
[ $vox_scale -eq 0 ] && calc_scale="1"
# 20210921 - BL measure effect of volume normalization
# 20240719 - BUGFIX
# 1. had 'x*1' but want to undo x/m. is now (x/m)*m
# 2. did not respect scaling
# now can do: x/m*m * scale/nvox (but thats a bad idea?)
[ $volnorm_opt == "none" ] && calc_scale="m*$calc_scale"
# normalize each voxel by the volume average or median
# scale by SCALE and number of good (within mask) voxels in the run
# sd_input is only non-empty when '-zscore'
# '$calc_expr*$calc_scale' defaults to (x/m)*$SCALE/$numvox
# s.t.
# v=(x/mean(x))*(1000/length(x))
# sum(v) = 1000
echo "$calc_expr*$calc_scale"
}
one_ta(){
# GLOBALS:
# MASK, MASK_REL, censor_rel, MAXVOLS,
# t2_inv, volnorm_opt, use_zscore, SCALE, vox_scale, all_numvox
#
# makes $runoutput and $volnorm_1D
# looks like ${cnt}_tat2.nii.gz or $tmpd/${cnt}_keep$n_tat2.nii.gz if censoring
local cnt="$1"; shift
local input="$1"; shift
local tmpd="$1"; shift
local ROIstats_cmd;
# 20231207 - getting bash quoting right is hard
! [[ $input =~ .nii.gz$|.nii$|.HEAD$ ]] &&
err "4D EPI input file '$input' must end in nii, nii.gz, or HEAD
Do you have a bad glob for eg -censor_rel? try putting 'echo' in front of your command to debug"
[ -n "$MASK" ] && mask="$MASK" || mask=$(find_rel_file "$input" "$MASK_REL")
[ ! -f $mask ] && err "mask '$mask' DNE; add -mkmask to create"
local filedesc_volnorm=_volnorm-${volnorm_opt//[- ]/}
# change $input and $filedesc_whichvol if $cesnor_rel exist
censor_file="" # updated version of '$censor_rel' if used
update_with_censor
# must match *_tat2.nii.gz to be picked up by 3dTstat at the end
runoutput="$tmpd/${cnt}${filedesc_whichvol}${filedesc_volnorm}_tat2.nii.gz"
# calc number of voxels
numvox=$(3dBrickStat -count -non-zero $mask|sed 's/[\t ]//g')
# verbose
msg "tat2: $input -> $runoutput"
volnorm_1D=$tmpd/${cnt}${filedesc_volnorm}.1D
vol_sd=$tmpd/${cnt}_sd.1D
# if we want 1/t2*, switch up input to that
if [ $t2_inv -eq 1 ]; then
inv_out="$tmpd/${cnt}_inv.nii.gz"
3dcalc -datum float -x $input -expr 1/x -prefix $inv_out
input="$inv_out"
fi
# calc sum of each volume. output is normalize value per TR
# as long as we want to do some volume normalization
if [ $volnorm_opt != "none" ]; then
ROIstats_cmd="3dROIstats -nomeanout $volnorm_opt -mask $mask -1Dformat '$input' > '$volnorm_1D';"
eval "$ROIstats_cmd"
else
ROIstats_cmd="# no vol norm cmd;"
# this is redundant with later calc_expr change. but here for constancy
perl -se 'print "1\n"x$nt' -- -nt=$(3dinfo -nt "$input") > "$volnorm_1D"
fi
# do the same for stddev (sigma) if zscoring
local calc_sd_input=""
if [ $use_zscore -eq 1 ]; then
local sd_cmd="3dROIstats -nomeanout -nzsigma -mask '$mask' -1Dformat '$input' > '$vol_sd';"
ROIstats_cmd="$ROIstats_cmd$sd_cmd"
eval "$sd_cmd"
calc_sd_input="-s $vol_sd"
fi
full_expr="$(args_to_3dcalc_expr)" # side-effect
calc_cmd="3dcalc \
-x $input \
-m $volnorm_1D \
$calc_sd_input \
-datum float -overwrite \
-expr '$full_expr'\
-prefix $runoutput;"
eval "$calc_cmd"
# for logging we store the total number of voxels used
# and the list of censor files
all_numvox="$all_numvox$numvox,"
all_censors="$all_censors\"$censor_file\","
all_roistats="$all_roistats$ROIstats_cmd"
all_calc_cmd="$all_calc_cmd$calc_cmd"
}
# make list of input arguments formated as json
json_list(){
local list=$(paste -sd, <<<$(printf '"%s"\n' "$@"))
echo -n "[${list}]"
}
_tat2(){
TMPD=$(mktemp -d $TMPLOCATION/tat2star_XXXX)
cnt=0
volcount=0
all_numvox="" # will be populated with each subjectmask (csv)
all_censors="" # appened list of all censor files used (csv)
all_roistats="" # all 3dROIStats volume normalization commands
all_calc_cmd="" # 3dcalc commands
local nvols_nocen=0
if [ $MKMASK -eq 1 ]; then
MASK="$TMPD/subj_mask.nii.gz"
coverage_mask $MASK "${GLOBFILES[@]}"
fi
for input in ${GLOBFILES[@]}; do
one_ta $cnt "$input" "$TMPD"
# now have e.g. 0_keep200_tat2.nii.gz
prev_cnt=$cnt
let ++cnt
# track total uncensored volume count
let nvols_nocen+=$(3dinfo -nt "$input")
# we can stop processing runs if we have enough
if [ $MAXVOLSTOTAL -gt 0 -a $IDX_SAMPLE_METHOD == "first" ]; then
this_run=$TMPD/${prev_cnt}_*tat2.nii.gz
n_vol_run=$(3dinfo -nt $this_run)
let volcount+=$n_vol_run
[ $volcount -lt $MAXVOLSTOTAL ] && continue
msg "have $volcount >= $MAXVOLSTOTAL and -sample_method $IDX_SAMPLE_METHOD. stopping at $cnt/$nfiles runs"
nfiles=$cnt
break
fi
done
actual_nfiles=$(ls $TMPD/*_tat2.nii.gz|wc -l)
[ $actual_nfiles -ne $nfiles ] && err "ERROR: create $actual_nfiles/$nfiles ($cnt iterations) in $TMPD"
msg "tat2star $OUTPUT"
# combine all the normalized runs if we need to
if [ $nfiles -gt 1 ]; then
3dTcat -overwrite -prefix $TMPD/tat2_all.nii.gz $TMPD/*_tat2.nii.gz
concated_file=$TMPD/tat2_all.nii.gz
else
concated_file=$TMPD/*_tat2.nii.gz
fi
# truncate them
# count is 1-indexed, indexing is 0-indexed
# what to do if too few. maybe should die
# instead use min and update MAXTOTAL to be current maxtotal
concat_nvol=$(3dinfo -nt $concated_file)
if [ $MAXVOLSTOTAL -gt 0 ]; then
use_vols0=$MAXVOLSTOTAL
concat_nvol0=$concat_nvol
# 0 index
let concat_nvol0--
let use_vols0--
if [ $concat_nvol0 -lt $use_vols0 ]; then
msg "WARNING: have fewer volumes than requested! $use_vols0 < $concat_nvol0; 0-idx"
MAXVOLSTOTAL=$concat_nvol
use_vols0=$concat_nvol0
fi
maxvolsubset="[0..$use_vols0]"
[ "$IDX_SAMPLE_METHOD" != first ] &&
maxvolsubset="[$(seq -s, 0 $concat_nvol0 |idx_shuffle $IDX_SAMPLE_METHOD | firstn_csv $MAXVOLSTOTAL)]"
msg "using $MAXVOLSTOTAL from $concat_nvol ($maxvolsubset)"
else
maxvolsubset=""
MAXVOLSTOTAL=$concat_nvol
fi
# time norm
local Tstat_cmd="3dTstat -prefix $OUTPUT -overwrite $tnorm_opt $concated_file$maxvolsubset"
eval "$Tstat_cmd"
### HISTORY/PROVENANCE
# 1D create file provenance is not tracked and tmp files will be remove.
# put final tat2 file. also in log.json
3dNotes -h "$all_roistats" "$OUTPUT"
all_nts="$(3dinfo -nt $TMPD/*_tat2.nii.gz|tr '\n' ',')"
[ $vox_scale -eq 0 ] && scale_msg="SCALE=NA" || scale_msg="SCALE=$SCALE"
# use gitver (lncdtool) when on path and tat2 tracked in git repo (e.g. in development)
# otherwise use hardcoded TAT2_VER (e.g. copied just tat2/no repo; or using from system package)
gitver="$(gitver "$0"||:)"
[[ -z "$gitver" || "$gitver" =~ :NA$ ]] && gitver=$TAT2_VER
note_msg="$0 $allargs # $(pwd); $scale_msg, nvoxes=$all_numvox, censors=$all_censors; nt=$all_nts; vol:$volnorm_opt, time:$tnorm_opt; $MAXVOLSTOTAL/$concat_nvol/$nvols_nocen used/not-cen/total nvols; $cnt/${#GLOBFILES[@]} files; $gitver"
3dNotes -h "$note_msg" "$OUTPUT"
# computer-friendly output for parsing when running many option permutations
[ "${WRITE_JSON_LOG:-1}" -eq 1 ] &&
cat <<- HEREDOC > "${OUTPUT/.nii.gz/.log.json}"
{
"cmd": "$0 $allargs",
"roistats_cmds": "${all_roistats}",
"volume_norm_cmds": "${all_calc_cmd}",
"collapse_cmd": "${Tstat_cmd}",
"expr": "$full_expr",
"volume_norm_3dROIstats": "${volnorm_opt//-/}",
"time_norm_3dcalc": "${tnorm_opt//-/}",
"ref_mask": "$mask",
"nvox": [${all_numvox::-1}],`# -1 to remove trailing comma`
"censor_files": [${all_censors::-1}],
"nt": [${all_nts::-1}],
"tat2_input": $(json_list $TMPD/*_tat2.nii.gz),
"raw_inputs": $(json_list ${GLOBFILES[@]}),
"n_raw_used": $cnt,
"output": "$OUTPUT",
"concat_nvol": $concat_nvol,
"nvolx_nocen": $nvols_nocen,
"scale": "$SCALE",
"git": "$gitver",
"tat2_ver": "$TAT2_VER",
"use_inverse": "$t2_inv",
"use_scale": "$vox_scale",
"maxvols": $MAXVOLSTOTAL
}
HEREDOC
[ $CLEAN -eq 1 ] && \
rm -r "$TMPD" ||\
echo "find temporary files in '$TMPD'; remove '-noclean' if not testing"
return 0
}
# run if not sourced
if ! [[ "$(caller)" != "0 "* ]]; then
set -eou pipefail
# 20210218WF - hygenic file usage -- tmp files after crash clog fileystem
trap '[ -n "$TMPD" -a -d "$TMPD" -a $CLEAN -eq 1 ] && [[ "$TMPD" =~ $TMPLOCATION ]] && rm -r "$TMPD"' EXIT
args_are_sane "$@"
_tat2
fi