1#!/bin/tcsh 2 3# Make a table of NAcc properties 4 5# NB: this script is just for *single* echo data at the moment 6 7### ver = 2.0 8# [PT: Mar 30, 2021] Rewrote to put in loops over ROIs and dsets: more 9# generalizable. Added mean volreg dset info, too. 10### ver = 2.1 11# [PT: Mar 30, 2021] Fixed a bug in the winf* calcs 12### ver = 2.2 13# [PT: Mar 30, 2021] add in Harv-Ox atlas ROIs for stats 14 15# ------------------------------------------------------------------------ 16# user specifies inputs and parameters 17 18set narg_good = 3 # just the num of args needed, for checking 19set subj = $1 # subj ID 20set dir_ap_ss = $2 # something/${subj}.results, or ${tempdir} (for SLURM) 21set dir_fs_ss = $3 # something/SUMA 22 23# ------------------------------------------------------------------------ 24# define other variables of interest: dsets, intermediate prefixes and 25# final outputs; note that we will jump to AP results directory for 26# most processing, so most names are local to there. 27 28set here = $PWD 29 30# these should be the follower dset names used in AP 31set foll_fs2000 = follow_ROI_FSall00+tlrc.HEAD 32set foll_fs2009 = follow_ROI_FSall09+tlrc.HEAD 33 34# this is the Harv-Ox atl, resampled to final EPI grid, with just some 35# particular ROIs present+labeled 36set hox_atl = /data/NIMH_Haskins/dylan/template_hox/HOX_rois_333.nii.gz 37 38# tmp prefix for intermediates; cleaned up at end 39set tpref = __rep_tmp 40 41# output name of report 42set nacc_rep = report_nacc.txt 43 44# ------------------------------------------------------------------------ 45 46set narg = ${#argv} 47if ( ${narg} == ${narg_good} ) then 48 echo "" 49 echo "++ Good number of input args (${narg}) for subj: ${subj}" 50 echo "" 51else 52 echo "** ERROR: Should have ${narg_good} args input, but have: ${narg}" 53 exit 1 54endif 55 56# ------------------------------------------------------------------------- 57# before jumping to AP results, finalize ROI dataset names and properties 58 59# the dsets output by @SUMA_Make_Spec_FS, i.e., ones with labels 60set orig_2000 = ( ${dir_fs_ss}/aparc+aseg_REN_all.nii* ) 61set orig_2009 = ( ${dir_fs_ss}/aparc.a2009s+aseg_REN_all.nii* ) 62 63# reattach labeltables, if they don't have them yet 64if ( ! `3dinfo -is_atlas_or_labeltable ${dir_ap_ss}/${foll_fs2000}` ) then 65 3drefit -copytables ${orig_2000} ${dir_ap_ss}/${foll_fs2000} 66 3drefit -cmap INT_CMAP ${dir_ap_ss}/${foll_fs2000} 67endif 68 69if ( ! `3dinfo -is_atlas_or_labeltable ${dir_ap_ss}/${foll_fs2009}` ) then 70 3drefit -copytables ${orig_2009} ${dir_ap_ss}/${foll_fs2009} 71 3drefit -cmap INT_CMAP ${dir_ap_ss}/${foll_fs2009} 72endif 73 74# ------------------------------------------------------------------------- 75# from here on, everything is local to the AP results directory, so 76# move there, for simpler naming 77 78cd ${dir_ap_ss} 79 80# ------------------------------------------------------------------------- 81# define all ROIs of interest. Each dset is on final EPI grid already. 82 83# corresponding dsets (ROIs within which to get stats) and labels 84 85set all_rois = ( ${foll_fs2000}"<Left-Accumbens-area>" \ 86 ${foll_fs2000}"<Right-Accumbens-area>" \ 87 ${hox_atl}"<lh_hox_nacc>" \ 88 ${hox_atl}"<rh_hox_nacc>" ) 89 90set lab_rois = ( lh_nacc_fs rh_nacc_fs lh_nacc_hox rh_nacc_hox ) 91 92# ------------------------------------------------------------------------- 93# finalize dsets from within which we will get stats 94 95# the specific volreg dset to grab. just the run one one at 96# present. Would also need to adjust for multi-echo data, too. 97set dset_vr_4d = ( pb*${subj}*.r01*volreg+*.HEAD ) 98 99echo "++ Make intermediate dset: mean across time of volreg data" 1003dTstat \ 101 -overwrite \ 102 -prefix ${tpref}_vr.nii.gz \ 103 "${dset_vr_4d}" 104 105# corresponding dsets (all from which we will get stats) and labels 106 107set all_vols = ( ${dir_ap_ss}/TSNR.vreg.r01.${subj}+tlrc.HEAD \ 108 ${dir_ap_ss}/TSNR.${subj}+tlrc.HEAD \ 109 ${tpref}_vr.nii.gz ) 110 111set lab_vols = ( tsnr_vr \ 112 tsnr_fin \ 113 mean_vr ) 114 115# ------------------------------------------------------------------------- 116# start report, and populate with values 117 118echo "" 119echo "Number of ROI : $#all_rois" 120echo "Number of volumes : $#all_vols" 121echo "" 122 123printf '' > ${nacc_rep} 124 125printf "%-30s : %10s \n" \ 126 "subject ID" "${subj}" \ 127 >> ${nacc_rep} 128 129printf "\n" >> ${nacc_rep} 130 131# Value of interest: size of each ROI 132 133foreach ii ( `seq 1 1 ${#all_rois}` ) 134 set rdset = "${all_rois[$ii]}" 135 set rlab = "${lab_rois[$ii]}" 136 137 set val = `3dROIstats -quiet \ 138 -nzvoxels \ 139 -nobriklab \ 140 -nomeanout \ 141 -mask "${rdset}" \ 142 "${rdset}"` 143 144 printf "%-30s : %10s \n" \ 145 "${rlab} nvox" "${val}" \ 146 >> ${nacc_rep} 147end 148 149printf "\n" >> ${nacc_rep} 150 151# Value of interest: stats and percentiles of data within each ROI. 152# Loop over each volume of interest 153 154foreach hh ( `seq 1 1 ${#all_vols}` ) 155 set vdset = "${all_vols[$hh]}" 156 set vlab = "${lab_vols[$hh]}" 157 158 foreach ii ( `seq 1 1 ${#all_rois}` ) 159 set rdset = "${all_rois[$ii]}" 160 set rlab = "${lab_rois[$ii]}" 161 162 # calculate all stats of interest, then reorder for presenting 163 set sss = `3dBrickStat \ 164 -slow \ 165 -perc_quiet \ 166 -percentile 0 25 100 \ 167 -mean \ 168 -stdev \ 169 -mask "${rdset}" \ 170 "${vdset}"` 171 set all_stat = ( $sss[6-7] $sss[1-5] ) 172 set lab_stat = ( mean std min p25 med p75 max ) 173 174 # write out each stat of interest into the report 175 foreach jj ( `seq 1 1 ${#all_stat}` ) 176 set stat = "${all_stat[$jj]}" 177 set slab = "${lab_stat[$jj]}" 178 179 # full label for table row 180 set flab = "${rlab} ${vlab} ${slab}" 181 182 printf "%-30s : %10.1f \n" \ 183 "${flab}" "${stat}" \ 184 >> ${nacc_rep} 185 end 186 187 # Extra 'stat' of interest: 5-95%ile window scaled by median 188 # this gives a dimensionless value of the range of values in 189 # the ROI. That is, large values have a wide range, small 190 # values have a small one. Also do this for a 25-75%ile 191 # window, in case that is more elucidating. We call these 192 # "window fractions" (winf) 193 194 set www = `3dBrickStat \ 195 -slow \ 196 -perc_quiet \ 197 -perclist 5 5 25 50 75 95 \ 198 -mask "${rdset}" \ 199 "${vdset}"` 200 201 set w1 = `echo "scale=4; (${www[5]}-${www[1]})/${www[3]}" | bc` 202 set w2 = `echo "scale=4; (${www[4]}-${www[2]})/${www[3]}" | bc` 203 204 set all_win = ( $w1 $w2 ) 205 set lab_win = ( winf90 winf50 ) 206 207 # write out each stat of interest into the report 208 foreach jj ( `seq 1 1 ${#all_win}` ) 209 set wval = "${all_win[$jj]}" 210 set wlab = "${lab_win[$jj]}" 211 212 # full label for table row 213 set flab = "${rlab} ${vlab} ${wlab}" 214 215 printf "%-30s : %13.4f \n" \ 216 "${flab}" "${wval}" \ 217 >> ${nacc_rep} 218 end 219 220 printf "\n" >> ${nacc_rep} 221 end 222end 223 224# ------------------------------------------------------------------------ 225# done with most work 226 227# clean up 228if ( 1 ) then 229 echo "++ Clean up intermediate file(s)" 230 \rm ${tpref}_* 231endif 232 233# jump back 234cd ${here} 235 236echo "" 237echo "++ Done with report. Check out file:" 238echo " ${dir_ap_ss}/${nacc_rep}" 239echo "" 240 241exit 0 242