Scripting p-value threshold + binarize on GLT sub-brick for batch processing

AFNI Message Board
1

Howdy AFNI experts,

I am running a functional localizer GLM in each subject's T1w space using 3dDeconvolve. My design has three conditions (Animals, Foods, Tools) with a BLOCK(10,1) HRF across 3 concatenated runs, and I am using -tout to get t-statistics alongside my betas. Here is a simplified version of my call:

3dDeconvolve \
    -input  sub-${subj}_task-localizer_allrun_space-T1w_SI.nii.gz \
    -mask   sub-${subj}_space-T1w_GM_probseg_resampled.nii.gz \
    -concat '1D: 0 134 268' \
    -polort A \
    -num_stimts 3 \
    -stim_times 1 times_Animals.1D 'BLOCK(10,1)' -stim_label 1 Animals \
    -stim_times 2 times_Foods.1D   'BLOCK(10,1)' -stim_label 2 Foods \
    -stim_times 3 times_Tools.1D   'BLOCK(10,1)' -stim_label 3 Tools \
    -num_glt 7 \
    -gltsym 'SYM: Foods -0.5*Animals -0.5*Tools' -glt_label 3 Foods_vs_Others \
    -tout \
    -bucket sub-${subj}_T1w_SI_func_localizer.nii.gz

What I am trying to do:

For each subject, I want to take a specific contrast from the output bucket — for example Foods > Others (Foods_vs_Others) — threshold it at p < 0.001 (uncorrected), binarize the surviving voxels into a binary mask, and then intersect that binary mask with a subject-specific anatomical VTC ROI (also in T1w space). The goal is to get a subject-specific functional VTC mask showing only voxels that are both (a) within the anatomical VTC boundary and (b) significantly activated for that contrast.

I know this can be done interactively in the AFNI GUI — overlaying the stat map, setting the p-value threshold in the slider, and using the Write/Save function to extract the thresholded map. But I need to automate this across ~100 participants, so I am looking for the command-line equivalent.

My specific questions:

  1. What is the correct way to extract and threshold a specific GLT t-stat sub-brick from the bucket at p < 0.001? I was thinking 3dcalc or 3dmaskdump but I am not sure of the right approach.
  2. Once thresholded, what is the recommended way to binarize the surviving voxels into a 0/1 mask? Is fslmaths -bin appropriate here, or is there a preferred AFNI-native approach?
  3. Is there a way to have AFNI automatically determine the correct t-threshold for a given p-value, rather than hardcoding a t-value?

Any guidance would be much appreciated. Thank you!

Thanks,
Sahithyan

2

Hi, Sahityan-

AFNI programs general recognize subbrick selection, for one or more volumes (and the same for surfaces and 1D datasets) using squarebrackets. So, selecting the [10]th volume of a dset (counting from 0, of course!) would be done with: DSET"[10]". The need for quotes around the square bracket is a bit of shell stuff. You can specify intervals , like DSET"[10..15]" or even mixtures: DSET"[5,8,19..23,26..35]", etc.

Outputs of AFNI stats programs have string labels per volume, so you can specify with those. For example, a Bootcamp example stats dset from running 3dDeconvolve via afni_proc.py has this for selection: stats.sub-000.affine+tlrc.HEAD"[vis#0_Tstat]". To see what labels are used, run any of:

3dinfo DSET
3dinfo -verb DSET
3dinfo -label DSET

AFNI stats dsets also tend to have stat-to-p mapping internally, because they store the degree of freedom information for each stat. You have to also specify the sidedness of the test, which in 99.9% of cases should be 2-sided to match the standard hypotheses used, as per that paper.

The AFNI GUI makes direct use of this in various ways, having a p-value below the color/threshold bar, so you can see what a statistical threshold corresponds to. You can also right click there or the "Thr" button above the cbar, and have an entry menu to set the threshold value directly, or using a p-value.

Some programs make use of this directly, like 3dClusterize, and you can specify thresholds via p-values directly (as well as the stat value).

Additionally, you can use p2dsetstat to figure out the stat for a given p, like here:

p2dsetstat                                       \
    -inset stats.sub01+tlrc"[2]"                 \
    -pval 0.001                                  \
    -bisided

Adding the -quiet option makes it even more useful for scripting.

Consider 3dClusterize, if you want to binarize and find contiguous islands as a map of ROIs. To simply threshold and output a binarize mask, consider 3dcalc, which is an omnibus calculation program that is a Swiss Army Knife of handiness. For example, this puts together a couple things (using tcsh syntax for the variable, but you can change that accordingly):

set thr = `p2dsetstat                                                        \
              -inset       stats.sub-000.affine+tlrc.HEAD"[vis#0_Tstat]"     \
              -pval        0.001                                             \
              -bisided                                                       \
              -quiet`

3dcalc                                                                       \
    -a          stats.sub-000.affine+tlrc.HEAD"[vis#0_Tstat]"                \
    -expr       "ispositive(a-${thr})-isnegative(a+${thr})"                  \
    -prefix     mask_of_stat_regions.nii.gz

You can even use the stat mapping to

BUT please note that thresholding is often the wrong thing to focus on in a strict sense, and transparent thresholding is generally a much better, interpretable, meaningful and reproducible way to go. Please see Highlight, Don't Hide, Go Figure, Overcoming the Curse of Dimensionality, More than the Tip of the Iceberg and the Splendour of Unthresholded Maps, for compelling reasons why.

Importantly, you can also script making images with transparent thresholding, including specifying a p-value that gets translated in the above ways, using @chauffeur_afni. Many of the APQC HTML images are made this way, for example. Here is a tutorial showing how.

--pt

3

Hi @ptaylor

Thank you so much for the detailed response!