Checking Handedness in SPM?

ajschadler · August 7, 2021, 3:37am

I have a language task, Verbal Paired Associates, in which I am interested in lateralization of activation, primarily in language regions and the MTL. Data is preprocessed in fmriprep then smooth/scaled and deconvolved in AFNI.

I have been comparing two methods:

Afni

3dmerge --> 3dROIstats
This is performed on each hemispherical ROI. Voxel counts are then used to calculate an LI.

SPM

Using the LI-toolbox, default settings, using the same ROIs as in 1.

Unfortunately, both methods are giving opposing answers: one says a subject is left lateralized for a given ROI, the other says right lateralized!

One of my colleagues suggested this could be due to an orientation or flipping issue. The orientation of the fmri images is:

$ 3dinfo -handedness -orient sub-01_task-VPA_run-1_space-MNI152NLin2009cAsym_desc-preproc.nii.gz
R       LPI

$ fslinfo -getorient
NEUROLOGICAL

When using the LI toolbox, the output does mention that the “handedness” of my data differs from “standard conventions” (or something like that).

My main question is: how can I confirm whether SPM or AFNI is causing some form of orientation flip?

My only other thought, if it is not an orientation flip, it may be a thresholding issue. I.e. since I based my AFNI LI results on a single threshold, whereas the LI toolbox produces LI curves from a bootstrapped approach, the LI toolbox may produce more “robust results”.

ptaylor · August 7, 2021, 1:27pm

Howdy-

I think it would first help to describe what your Steps #1 and #2 are… What are the explicit commands in each case? (NB: I only know the AFNI side of this to be able to comment on.)

I doubt AFNI is flipping your data…

Actually left-right flipping is something we work to detect while processing with afni_proc.py—relative LR flips between EPI and anatomicals are a surprisingly common phenomenon, and worth checking out–see here for more details:
Glen DR, Taylor PA, Buchsbaum BR, Cox RW, Reynolds RC (2020). Beware (Surprisingly Common) Left-Right Flips in Your MRI Data: An Efficient and Robust Method to Check MRI Dataset Consistency Using AFNI. Front. Neuroinformatics 14. doi.org/10.3389/fninf.2020.00018
Beware (Surprisingly Common) Left-Right Flips in Your MRI Data: An Efficient and Robust Method to Check MRI Dataset Consistency Using AFNI

Another couple factors here are: there is a difference between dataset orientation (how the rows, cols and slices are stored on the disk) and the data coordinate orientation (do “left” side of the brain coords have more negative or more positive values?). Data can be stored in any orientation, and still show the brain with the same coords. Re. data orientation—another confusing factor is that different software and vendors use different conventions.

I believe that SPM will always use LPI for output; in AFNI, we tend to use RAI as default, but you can often specify what you would like to use.

Have you looked at your data to visually check which ROIs seem larger on a given side? In the AFNI GUI, for example, you can run clusterizing to see the cluster sizes in vivo (including thresholding, if you wish).

–pt

ajschadler · August 18, 2021, 5:29pm

Hello,

For specific steps and commands, data was processed through fMRIprep, using default settings. Data is outputted to both anatomical and mni152nlin2009casym template space.
After, images are smoothed using 3dBLURtoFWHM, to a 4.0 mm kernel, and scaled (using 3dTstat, 3dcalc, etc.). Data is then regressed using 3dDeconvolve

For 1.

3dmerge \
  -dxyz=1 \
  -1clust 1 20 \ # Corresponding to athr of ~ 0.05 @ pthr of 0.01
  -1dindex $dindex \
  -1tindex $tindex \
  -1thresh 2.59 \ # Rounded. T-val corresponds to pthr of ~ 0.01
  -1noneg \
  -prefix $tmpMap
  stats.$task.$sub+tlrc

for mask in {l,r}h.MTL_rs+tlrc.HEAD ; do
  3dROIstats \
    -mask $mask \
    -nomeanout \
    -nzvoxels \
    $tmpMap > $outfile
done

LI is then calculated using counts from $outfile.

For 2.:
Relevant T-maps from stats.$task.$sub+tlrc are converted to .nii; same for mask files. Mask files are additionally combined into a bilateral mask

3dAFNItoNIfTI \
  -input stats.$task.$sub+tlrc'[N_GLT#0_Tstat]' \
  -prefix 'afni_TN.nii'

for mask in {l,r}h.MTL_rs+tlrc ; do
  3dAFNItoNIfTI \
    -input $mask
done

3dcalc \
  -a 'lh_MTL_rs.nii' \
  -b rh_MTL_rs.nii \
  -expr 'a+b' \
  -prefix 'MTL_rs.nii'

Both the bilateral mask and relevant T-map are given to the LI() function from the LI SPM toolbox, which by default uses a bootstrap approach to construct LI curves.

ptaylor · August 20, 2021, 7:09pm

Hi-

Thanks for specifying that.

Some comments on “For 1.”:

I would recommend using 3dClusterize rather than 3dmerge—the latter is very old, and I find the syntax quite opaque. You can also input both a “data” volume (e.g., a beta weight or coefficient) as well as an accompanying stat (e.g., t-stat), for making maps of both simultaneously. Here is an example of the former from the help file:

 3dClusterize                  \
    -inset stats.FT+tlrc.      \
    -ithr 2                    \
    -idat 1                    \
    -mask mask_group+tlrc.     \
    -NN 1                      \
    -bisided -3.313 3.313      \
    -clust_nvox 157            \
    -pref_map ClusterMap       \
    -pref_dat ClusterEffEst

NB: If the stats volume was created by AFNI, then it would have header information to convert from a specified p-value to an equivalent t-stat (see 3dClusterize help examples for more info on that). Note when converting between p- and statistics-values, one needs to have an appropriate sided hypothesis specified—in AFNI, our default is 2sided testing (because most researchers are looking for differences between things), while in other software the default tends to be 1sided (because of the kind of statistic output). So, this is something to be aware of in mixing software outputs and commands. Not correcting for having two, 1sided tests will artificially double FPR values; see more info here: A tail of two sides: Artificially doubled false positive rates in neuroimaging due to the sidedness choice with t‐tests.
If you want your positive and negative supra-threshold regions to not count as a single cluster (which I guess you wouldn’t), then using the “-bisided” option, as above, is good. I am not sure what the equivalent would be in other software.

Also note, when you run 3dClusterize, you will get a table per cluster that tells you how many voxels there are per cluster (3dROIstats can also do the same).

Another difference between software is what it means to be a neighbor of a voxel. Do they have to share a face, or a face or an edge, or a face or an edge or a node? AFNI, FSL and SPM have typically used different definitions as a default—in AFNI you can specify what you want, though, using “NN” (nearest neighbor value). See here:
https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/FATCAT/MakingROIs.html#getting-to-know-your-neighbors
I am not sure what method you used to calculate your clustersize threshold, but this might have been calculated using one “neighborhood” definition, and different software might have a different one.

Note that there has also been a lot of discussion in the field about what appropriate p-values should be used for voxelwise thresholding. Several papers (and journals, too) have pushed for using p=0.001 and not p=0.01—this is of course entirely up to you, but note that that might be a point of discussion with a reviewer.

On a scripting note: you are using “>” to redirect output in your for-loop. I think this will overwrite whatever was there previously, and it doesn’t look like $outfile changes value within that loop. So you might want to adjust that by having the $outfile change each time, or by using “>>” to redirect with appending.

Some comments on “For 2.”:

Converting BRIK/HEAD to NIFTI is fine, sure. One can use 3dAFNItoNIFTI or 3dcopy, etc. AFNI programs do read/write NIFTI files directly, in most cases, if specified (afni_proc.py uses all BRIK/HEAD, but copying output is pretty straightforward, as you have done).

The for-loop here is missing a “-prefix …” for the command, I think.

In the bilateral mask, if the hemispheric masks are each binary, then the L and R hemispheres have lost their identity in this sum; if you wanted to be able to tell them apart in the output, you could use:

3dcalc \
  -a 'lh_MTL_rs.nii' \
  -b rh_MTL_rs.nii \
  -expr 'a+2*b' \
  -prefix 'MTL_rs.nii'

for example. Of course, if they aren’t binary masks, this may already be accounted for. (NB: this assumes that the ‘a’ and ‘b’ datasets don’t overlap—you can check for that visually, if it might be a concern, or by seeing what the max value in the dset is after.)

In general—
You had initially asked about potential left-right flipping in each case. I don’t see how that could happen with the above commands. You can also verify by looking at the clusters and hemispheric maps in the GUI. You can even use the “Clusterize” button in the GUI, to see the clusterizing at work.

Sorry for the long and possibly extraneous reply!

–pt

ajschadler · August 20, 2021, 8:29pm

My apologies for not specifying that. Here is how t-values are calculated for the 3merge command:

# Please forgive the atrocious use of cut
df=$(3dinfo -VERB ${statmap} | grep -A 1 "'${contrast}_GLT#0_Tstat" | grep statpar | cut -d ';' -f 2 | cut -d ' ' -f 5)

This is fed to a python function like so

from scipy import stats
criT = stats.t.ppf(1 - (alpha/2), df)

Apologies for not clarifying that. There’s a conditional statement that checks if the output file exists, switching to >> if $outfile does.

Apologies, the -prefix is included in my original script.

Some additional clarification:
The conversion is only to use the SPM LI-toolbox, which requires a bilateral-mask.

It was super helpful. Thank you so much! After meeting with my colleagues, as well as following the suggesting of Ruff et al. 2008 (DOI: 10.3174/ajnr.A0841), we think it might be either 1. a thresholding issue, or 2. a flipping issue may be happening when entering the data into SPM.

I’m currently investigating 2. and would be happy to follow up on this thread when I find the answer.

ptaylor · August 20, 2021, 9:15pm

OK, glad some of that was useful, and I see that some of the scripting things were taken care of elsewhere—great.

I would still use 3dClusterize for simplicity, rather than 3dmerge. And if that stat volume is an AFNI-created one (which it appears to be), then you can use something like the following to have the internal header information convert the p-to-stat for you—note again the sidedness being directly stated, as well as the NN value:

 3dClusterize                  \
    -inset stats.FT+tlrc.      \
    -ithr "${contrast}_GLT#0_Tstat" \
    -idat 1                    \
    -mask MASK     \
    -NN 1                      \
    -bisided p=0.001           \
    -clust_nvox 157            \
    -pref_map ClusterMap       \
    -pref_dat ClusterEffEst

You can also calculate the p-to-any stat using AFNI’s p2dsetstat, which also seems simpler. The degrees of freedom are stored internally, so rather than use 3dinfo with grep and then a Python function, the value pops out directly. And you can run it with or without the “-quiet” option to for scripting convenience, e.g.:

p2dsetstat  \
    -inset stats.sub01+tlrc"[${contrast}_GLT#0_Tstat]"            \
    -pval 0.001                            \
    -2sided                                 \
    -quiet

Again, the sidedness issue between softwares should be noted in any conversion, as well as the neighborhood-ness. But some of that can be doublechecked by visualizing the data in a GUI, too.

–pt