Fast regression and mediation analysis of vertex or voxel MRI data with TFCE
The pre-print manuscript is available here as well as the supporting information.
19-07-2017
- version 1.2.2 is now available on pypipe (PIP).
- fixes a bug in the permutation testing of surfaces with unequal sizes
5-07-2017
- version 1.2.1 is now available on pypipe (PIP).
- A number of small fixes have been made.
- Release of an alpha version of multi-surface, multi-modality regression (mmr) analysis using the tmi image. The scripts can be found here.
- tm_multimodality_multisurface_regression.py is still in the early stages but it already has many new features including:
- In the same set of subjects, cortical thickness, surface area, neurite density, and fMRI contrast maps could be packaged in the a tmi file using create_tmi.py. Statistical analysis (and randomization) across all modalities can be performed using tm_multimodality_multisurface_regression.py
- Study-wide FWE correction with TFCE. i.e., 1-p(FWE<0.05) images are produced that are corrected via a standardization algorithm and maximum TFCE values among all surfaces.
- The results can be exported as ply files and viewed using Blender or meshlab. Alternatively, the tmi result files can be converted to mgh or nifti.
- Many surfaces can be analyzed at once. For example, analysis among all subcortical structures can be performed with TFCE and corrected for all surfaces (or even subcortical, cortical surfaces can be combined together in a single analysis).
- Due to the flexibility of the tmi image format, many more features will be released in the new feature. For instance, the tmi format can include both voxel and vertex data in the same space. This allows voxel-based (such as TBSS) analyses and vertex-based (such cortical thickness) to be performed concurrently with TFCE and FWE correction among all data analyzed.
- TLDR, why do volumetric analyses when surface area and cortical thickness can be analyzed as the same time? Bonferroni correction among different imaging modalities is not necessary using tm_multimodality_multisurface_regression.py and the tmi image format.
21-06-2017
- version 1.2 is now available on pypipe (PIP)
- added the ability to read and write binarized PLY files using 'tm_tools convert-surface'
- a preview of the new tfce_mediation imaging (tmi) for format group analyses will be released before OHBM 2017. It is based on the the PLY file format. Some advantages include a fast, space efficient, and expandable neuroimaging that harmonizes voxel and vertex images as well as surfaces all with a file history function, but more on that later.
16-05-2017
- added an option to use geodesic FWHM smoothing in tfce_mediation step0-vertex. Geodesic smoothing requires a distances list that can either be created using fwhm_compute_distances_parallel.py or downloaded from tm_addons. Currently, only the midthickness surface is supported, but other surfaces will be supported soon.
13-05-2017
- convert-surface now supports painting freesurfer label files to .ply files as well as outputing a legend. An example of the lh aparc atlas on the midthickness surface is: left view, right view, legend, and the lh Glasser et al. (2016) HCP-MMC1 atlas: left view, right view, and a ridiculously long legend
- Support for loading 4D minc volumes with tfce_mediation step0-voxel. Right now the outputs will still be nifti (which could be converted to mnc e.g. by using nii2mnc).
05-05-2017
- convert-surface now supports exporting the Polygon File Format (PLY). Statistic files (*.mgh) can now be painted on the vertices with specifying a threshold range (e.g 0.95 1 for pFWE corrected images) with an either red-yellow or blue-lightblue color scheme (as well as all matplotlib color schemes). So now you can 3D print your surface based results!
- A few examples of the high quality (3840x2160) render of painted results of the ?h.mean.area.00.mgh from 0.1mm to 1.0mm using the 'jet' colormap right, rostral, and superior. The surfaces were rendered using Blender and exported *.ply objects.
- A sample blender scene for making figures can be downloaded from here. n.b., when importing a new *ply object, the scale will need to be reduced.
- wiki is slowly becoming alive
03-05-2017
- Added surface conversion tool (tm_tools convert-surface) for converting gifti / MNI objects / freesurfer .srf files to freesurfer .srf files for analysis with TFCE_mediation or waveform objects / Stl objects for importing the surface to 3D rendering software such as Blender or Meshlab. Among other things, this means that there is now full support for statistical analysis of neuroimages processed with CIVET. Here is an example of the midthickness surface LH and RH.
02-05-2017
- Added fwhm_compute_distances_parallel to the misc_scripts directory. It is necessary for geodesic FWHM smoothing.
- Full support for CIFTI-2 and GIFTI in the very near future.
- Version 1.1 is now available on pip
21-04-2017
- Geodesic FWHM smoothing on the vertex images. It is performed on the midthickness surface. The runtime is approximately 4 minutes per image, but it can probably be faster. No more fudge factors!
- FWHM smoothing on voxel images.
- Too many additions to tm_maths than can be listed. Some examples are surface-based ICA artefact removal, curve fitting for calculation of the max TFCE null distribution of FWE correction, pca compression smoothing, many clustering algorithms, and powerful normalization tools. Check the help function.
- Initial support for CIVET vertex-based analyses (see the issues section)
- Version 1.1 will be pushed to pip in the near future.
TFCE_mediation should work with any UNIX based system. It has been tested on Ubuntu 14.04 and 16.04, Arch Linux, and OSX
Required
- Freesurfer (for vertex-based analyses)
Recommended
Recommended for parallelization:
or
PIP Install (recommended) The python dependencies will automatically be installed.
sudo pip install tfce-mediation
OR
Install from source (requires the following python packages: numpy, scipy, matplotlib, nibabel, cython, scikit-learn, joblib):
git clone https://github.com/trislett/TFCE_mediation.git
cd TFCE_mediation
sudo python setup.py install
-
Run standard Freesurfer (for vertex-based analyses) recon-all preprocessing of T1-weighted images.
-
Spherical registration
tfce_mediation step0-vertex -i subjectslist.csv area -p
Explanation: Perform a spherical surface registration for each subject to the fsaverage resampling either area or thickness.
Inputs:
- subjectslist.csv (a text file contain one column of subject IDs. Note, the IDs must match those in the $SUBJECTS_DIR folder.)
- area (the surface to include in the analysis. For basic use, this should be area or thickness)
Outputs:
- ?h.all.area.00.mgh (all subjects to the fsaverage template using spherical registration)
- ?h.all.area.03B.mgh (the above file after 3mm FWHM smoothing of the surface)
- Optional: Box-Cox transformation
tm_tools vertex-box-cox-transform -i lh.all.area.00.mgh 8
tm_tools vertex-box-cox-transform -i rh.all.area.00.mgh 8
# replace the ?h.all.area.03B.mgh with ?h.all.area.03B.boxcox.mgh
for i in lh rh; do
mv ${i}.all.area.03B.mgh ${i}.all.area.03B.backup.mgh;
mv ${i}.all.area.03B.boxcox.mgh ${i}.all.area.03B.mgh;
done
Explanation: It has been suggested that surface area follows roughly a lognormal distribution (Winkler, et al., 2012); therefore, vertex-box-cox-transform normalizes the unsmoothed images using a power transformation (Box-Cox) transformation.
Inputs:
- ?h.all.area.00.mgh (The unsmoothed concatenated surface area image for subjects included)
- 8 (The number of processors to use)
Outputs:
- ?h.all.area.00.boxcox.mgh (Box-Cox transformed image)
- ?h.all.area.03B.boxcox.mgh (Box-Cox transformed image after 3mm FWHM smoothing)
- Optional: orthonormalizing the regressors
tm_tools regressor-tools -i predictors.csv covariates.csv -o -s
Explanation: For the two-step multiple regression and mediation analyses using TFCE_mediation, it is recommended to scale (or orthonormalization) the regressors. The input file(s) should be dummy coded, and comma deliminated. The program returns either the orthogonalization of the input file(s) or it returns the residuals from a least squares regression to remove the effect of covariates from variable. In this example, we using the orthonormalization option (-o -s).
Inputs:
- predictors.csv (dummy-coded regressors of interest)
- covariates.csv (dummy-coded regressors of no interest)
Outputs:
- predictors_orthogonized.csv (orthogonormalized regressors of interest)
- covariates_orthogonized.csv (orthogonormalized regressors of no interest)
- Multiple Regression
tfce_mediation step1-vertex-regress -i predictors_orthogonized.csv covariates_orthogonized.csv -s area
Explanation: A two-step multiple regression is performed, and the resulting T-statistic images then undergo TFCE.
Inputs:
- predictors_orthogonized.csv (orthogonormalized regressors of interest)
- covariates_orthogonized.csv (orthogonormalized regressors of no interest)
- area (the surface of interest)
Outputs:
- python_temp_area/*.npy (Memory mapped numpy objects used later for permutation testing)
- output_area/tstat_area_?h_con?.mgh (t-statistic image or the respective hemispehere and contrast)
- output_area/negtstat_area_?h_con?.mgh (negative t-statistic image or the respective hemispehere and contrast)
- output_area/tstat_area_?h_con?_TFCE.mgh (TFCE transformed t-statistic image)
- output_area/negtstat_area_?h_con?_TFCE.mgh (TFCE transformed negative t-statistic image)
- output_area/max_TFCE_contrast_values.csv (The max TFCE value for the positive and negative contrast(s))
- Permuation Testing (Randomization)
tfce_mediation step2-randomise-parallel --vertex area -n 10000 -p 8
Explanation: Permutation testing using parallel processing. This script is a wrapper to make n=200 permution chunks, and parallelizing processing of the vertex-regress-randomise scipt for each chunk (in this case). i.e., each chunk to a different processor.
Input:
- area (surface of interest)
- 10000 (total number of permutations)
- 8 (number of cores to used with GNU parallel)
Output:
- output/perm_Tstat_area/perm_tstat_con?_TFCE_maxVertex.csv (the maximum TFCE value among all vertices of the entire cortex for each permutation. It is used to correct for family-wise error)
- Apply family-wise error rate correction
cd output
for i in lh rh; do
tfce_mediation vertex-calculate-fwep -i tstat_con1_area_${i}_TFCE.mgh perm_Tstat_area/perm_tstat_con1_TFCE_maxVertex.csv
tfce_mediation vertex-calculate-fwep -i negtstat_con1_area_${i}_TFCE.mgh perm_Tstat_area/perm_tstat_con1_TFCE_maxVertex.csv
done
Explanation: Calculate 1-P(FWE) vertex image from max TFCE values from randomisation. Note, the max TFCE values is determined from both the left and right hemisphere; therefore, there is only one max TFCE file per contrast. It is appropiate to use it correct the positive and negative t-statistic values from each hemisphere.
Input:
- tstat_area_?h_con?_TFCE.mgh (TFCE transformed T-statistic surface image)
- perm_tstat_con?_TFCE_maxVertex.csv (List with maximum TFCE values)
Output:
- tstat_area_?h_con?_TFCE_FWEcorrP.mgh (1-P(FWE) corrected image)
- View results
tm_tools vertex-freeview-quick -i tstat_area_lh_con1_TFCE_FWEcorrP.mgh tstat_area_rh_con1_TFCE_FWEcorrP.mgh
-
Preprocess a 4D volume. For instance using TBSS.
-
Load 4D volume into tfce_mediation
tfce_mediation step0-voxel -i 4dvolume.nii.gz mask.nii.gz
Explanation: Load the 4D volume and its binarized mask into tfce_mediation
Inputs:
- 4Dvolume.nii.gz (3D NiFti volumes for all subjects. e.g. for TBSS output use all_FA_skeletonised.nii.gz)
- mask.nii.gz (binarized brain makes. e.g. for TBSS output use mean_FA_skeleton_mask.nii.gz)
Outputs:
- python_temp/*.npy (Numpy memory mapped objects used for all voxel analyses/actions)
- Remove the effect of covariates on regressors for mediation analysis
tm_tools regressor-tools -i predictor_variable.csv covariates.csv -r -s
tm_tools regressor-tools -i dependent_variable.csv covariates.csv -r -s
Explanation: For mediation analyses using TFCE_mediation, the effect of covariates should be removed from all regressors of interest. ‘regressor-tools -r’ returns the residuals after regressing out the effect of the covariates. In the next step (voxel-step1-????), the covariates will regressed out from the 4D_image.
Inputs:
- ?????_variable.csv (dummy-coded regressors of interest)
- covariates.csv (dummy-coded regressors of no interest)
Outputs:
- ?????_resids.csv (residuals of the variable of interest)
- covariates_std_dm.csv (unit variance and demean covariates file)
- Mediation analysis
tfce_mediation step1-voxel-mediation -i predictor_resids.csv dependent_resids.csv -c covariates_std_dm.csv -m M
Explanation: Mediation analyses using TFCE_mediation in which the ‘predictor_variable’ is the independent variable, the 4d_image is the mediator, and the ‘dependent_variable’ is the dependent variable.
Inputs:
- ?????_resids.csv (residuals of the variable of interest)
- covariates_std_dm.csv (unit variance and demean covariates file)
Outputs:
- output_med_M/SobelZ_M.nii.gz (the untransformed Sobel Z-statistics of the indirect effect)
- output_med_M/SobelZ_M_TFCE.nii.gz (the TFCE transformed Z-statistics)
- output_med_M/maxTFCE_contrast_value.csv (the maximum TFCE values of the statistics image)
Note, different mediation models can be performed.
- Permuation Testing (Randomization)
tfce_mediation step2-randomise-parallel --voxel -n 10000 -m M -p 8
Explanation: Permutation testing using parallel processing. This script is a wrapper to make n=200 permution chunks, and parallelizing processing of the voxel-mediation-randomise scipt for each chunk (in this case). i.e., each chunk to a different processor.
Input:
- M (mediation type)
- 10000 (total number of permutations)
- 8 (number of cores to used with GNU parallel)
Output:
- output_med_M/perm_SobelZ/perm_Zstat_TFCE_maxVoxel.csv (the maximum TFCE value among all voxels of the entire cortex for each permutation. It is used to correct for family-wise error)
- Apply family-wise error rate correction
cd output_med_M
tfce_mediation voxel-calculate-fwep -i SobelZ_M_TFCE.nii.gz perm_SobelZ/perm_Zstat_TFCE_maxVoxel.csv
Explanation: Calculate 1-P(FWE) voxel image from max TFCE values from randomisation.
Input:
- SobelZ_M_TFCE.nii.gz (TFCE transformed Z-statistic voxel image)
- perm_Zstat_TFCE_maxVoxel.csv (List with maximum TFCE values)
Output:
- SobelZ_M_TFCE_FWEcorrP.nii.gz (1-P(FWE) corrected image)
- View results
Use your favorite viewer. e.g., fslview, mricron, MRIcroGL, etc.