Integrating FreeSurfer and PALM
Contents
Introduction
fspalm is a FreeSurfer script written to aid in the use of Permutation Analysis of Linear Models (PALM) for the correction of multiple comparisons on the surface and in the volume. PALM, written by Anderson Winkler, offers a huge number of analysis options for permutation. See the PALM Users Guide for more details. The purpose of fspalm is to make it relatively easy to interface with PALM, though it is still up to the user to make most of the decisions. Those using PALM should cite the appropriate articles, eg:
Winkler AM, et al. Non-Parametric Combination and related permutation tests for neuroimaging. Hum Brain Mapp. 2016 Apr;37(4):1486-511.
Winkler AM, et al, . Multi-level block permutation. Neuroimage. 2015;123:253-68.
Winkler AM, et al, . Faster permutation inference in brain imaging. Neuroimage. 2016 Jun 7;141:502-516
Setting Up
Before using fspalm, download PALM from fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM and add it to your matlab path. Note: if you get errors, you man need a newer version of matlab (2018+ should work). You must also have these two FreeSurfer paths in your matlab path (if not there already): $FREESURFER_HOME/matlab and $FREESURFER_HOME/fsfast/toolbox. You will also need to get fspalm and a two other files from ftp://surfer.nmr.mgh.harvard.edu/pub/dist/freesurfer/6.0.0-patch/PALM.
To use fspalm, first run mri_glmfit to get a glm folder. This folder will contain the contrasts as specified with the --C option to mri_glmfit as well as all the other mri_glmfit output. You can use fspalm with both volume- and surface-based analyses; it will figure everything out. The configuration and output is meant to mimic the FreeSurfer mri_glmfit-sim program, which can also be used to perform permutation-based correction for multiple comparisons (though without all the options of PALM).
Use Cases
There are two use cases: (1) You just run fspalm to convert FreeSurfer data into something that PALM can understand and to create a matlab script to run PALM; you then edit the matlab script (or create a new one) to run the exact PALM analysis desired. This is done by adding the --monly flag to the script below. (2) You can use fspalm to run the permutation and generate cluster tables as shown below. You will not be able to exploit the full PALM functionality in this use case, but you might not need to.
Permutation can be tricky to apply. The user is responsible for making sure that their design is appropriate for permutation and for using the proper permutation settings.
Example fspalm command
fspalm --glmdir glmdir --cft 1.3 --twotail --name palm-twotail-1.3 --iters 1000 --2spaces --cwp .05 --glmdir is the output folder of mri_glmfit --cft 1.3 is the voxel-wise cluster forming threshold (CFT) as -log10(p). for cft = 1.3, p=.05 --name mypalm-twotail-1.3 : create a folder called mypalm-twotail-1.3 in glmdir for PALM output --iters 1000 means to run 1000 permutation iterations --2spaces : bonferroni correct clusterwise pvalues for two-spaces (eg, left and right hemispheres) --cwp .05 : cluster-wise p-value, filter out all clusters that have p>.05 (after any bonferroni correction)
This command will do several things. First it will convert the design and contrast matrices to design.mat and design.con files expected by PALM. It will create a matlab file called run_palm.m where it will create a palm command line. It will then run 1000 permutations and save the output with a base name of fsp (eg, fsp_clustere_tstat_fwep_c2.mgz). It will then create summary files as well as annotations (surface) or segmentations (volumes). For example, if there were a contrast called g1-g2, then it would create g1-g2.dpv.clustertable.summary with a list of the surviving clusters, their size, the location of the peak voxel in MNI305 space, and the segmentation/annotation that the peak voxel was located in. The clusterwise p-value is located in the "Max" column of the g1-g2.clustertable.summary. This file also has location information, but the values from the pv cluster table will be better since it uses the maximum from the uncorrected data as the location. There will also be a file called g1-g2.y.ocn.dat. This is a text file with mean input values for each subject and cluster (ie, each row is a subject, each column is a cluster). There will be a filed called g1-g2.ocn.{mgz,nii} that will be a segmentation (ie, the value of the voxel is the cluster number). For surface-based analyses, there will also be a file called g1-g2.ocn.annot which will be an annotation of the clusters.
Setting the voxel-wise Cluster-Forming Threshold (CFT)
The CFT is in log10 units, so if you want p<0.01, you would specify --cft 2 because 2= -log10(.01). Likewise for p<0.05, --cft 1.3; p<.001, --cft 3. The CFT is somewhat arbitrary, but if you set it too high, you may not have any voxels above threshold. If you set it too low, then there may be many clusters likely by chance making your clusters hard to find. In both cases, you are losing sensitivity (ie, increased false negatives). Permutation should assure that your specificity (false positive rate) is always a the level set by the --cwp argument.
One-tailed or two?
The "tailedness" of the analysis refers to whether you expect the contrast to be positive or negative or you have no expectation. If you have an expectation, then you should use a one-tailed test, otherwise use a two-tailed test. The one-tailed will be more powerful but you must have an a priori reason for using a one-tailed test. When selecting a two-tailed test, the CFT is increased by 0.3 behind the scenes because uncorrected, voxel-wise p-values in PALM are always one-tailed. Increasing the (log10) CFT by 0.3 is the same halving the p-value threshold. Eg, if you wanted a voxel-wise threshold of 0.05, the CFT would be -log10(.05)=1.3. For a two-tailed test, this would be changed to 1.6 which corresponds to .025. So you would specify --cft 1.3 --twotail and fspalm will automatically use 1.6. In some ways, this added complication is perhaps unnecessary because the CFT is mostly arbitrary anyway. But this is the way that mri_glmfit-sim behaves and we wanted to be consistent. One implication of this is that the cluster sizes will be different between the one- and two-tailed tests.
2spaces or 3spaces
When performing analyses in FreeSurfer, the hemispheres are analyzed separately. However, when performing correction for multiple comparisons, one must take into account the entire search space (eg, a cluster of a given size is twice as likely under the null when looking over two hemispheres as opposed to one). To take this into account, fspalm must be told how many search spaces there are. If you only have a single search space (eg, this is a pure volume analysis or a surface-based interhemispheric analysis), then do NOT specify --2 or --3spaces since the default is one space. If you have a surface-based analysis where you are going to be examining both left and right hemisphere (but not subcortical structures), then use --2spaces. If you are going to be looking at both surface- and volume-based results (eg, in an fMRI or PET analysis), then use --3spaces.
One-sample Group Mean (OSGM)
Those using an OSGM analysis (ie, the design matrix is a column of ones) should add --pargs="-ise" to the fspalm command line to perform sign flipping. This should actually be done for any contrast that is equivalent to OSGM. Eg, if you have two groups and a contrast of [+1 +1].
Other options
--octave : use octave instead of matlab. octave is basically a freely available equivalent of matlab. You must download and install octave on you system
--pargs SUBARGS : you can use this feature to add PALM-specific arguments to the palm command line in matlab. Note that if the argument has a dash ("-"), you must include it in quotes (eg, --pargs="-ise")
--centroid : report cluster centroid rather than max (surfaces only)