Demo: Threshold-Free Cluster Enhancement (TFCE) on surface dataset
The data used here is available from http://cosmomvpa.org/datadb.zip
This example uses the following dataset: + 'digit' A participant made finger pressed with the index and middle finger of the right hand during 4 runs in an fMRI study. Each run was divided in 4 blocks with presses of each finger and analyzed with the GLM, resulting in 2*4*4=32 t-values
This example illustrates the use of Threshold-Free Cluster Enhancement with a permutation test to correct for multiple comparisons.
TFCE reference: Stephen M. Smith, Thomas E. Nichols, Threshold-free cluster enhancement: Addressing problems of smoothing, threshold dependence and localisation in cluster inference, NeuroImage, Volume 44, Issue 1, 1 January 2009, Pages 83-98.
- For CoSMoMVPA's copyright information and license terms, #
- see the COPYING file distributed with CoSMoMVPA. #
Contents
Check externals
cosmo_check_external({'surfing','afni'});
Set data paths
The function cosmo_config() returns a struct containing paths to tutorial data. (Alternatively the paths can be set manually without using cosmo_config.)
config=cosmo_config(); digit_study_path=fullfile(config.tutorial_data_path,'digit'); readme_fn=fullfile(digit_study_path,'README'); cosmo_type(readme_fn); output_path=config.output_data_path; % resolution parameter for input surfaces % 64 is for high-quality results; use 16 for fast execution ld=16; % reset citation list cosmo_check_external('-tic'); % load single surface intermediate_fn=fullfile(digit_study_path,... sprintf('ico%d_mh.intermediate_al.asc', ld)); [vertices,faces]=surfing_read(intermediate_fn);
Overview
--------
fMRI responses to a human participant pressing buttons with the index and middle finger.
Contents
--------
- glm_T_stats_allruns+orig.{BRIK,HEAD}:
t-statistics associated with finger presses.
There are 4 runs, each with 4 blocks.
Each block has two t-statistics, one for each of button presses for the
index and middle finger (in that order).
- epi+orig.{HEAD,BRIK}:
A single EPI image.
- icoXX_mh.YY.asc:
Surface anatomy meshes of two hemispheres with standard topology, generated
using AFNI SUMA's MapIcosahedron. Left and right hemisphere surfaces were
merged in the order left, right; the first half of the nodes and faces refer
to the left hemisphere, the second half to the right hemisphere.
XX={16,64} is the number of linear divisions of MapIcosahedron; surfaces
have (XX^2)*10+2 nodes and XX^2*20 faces in each hemisphere.
YY={pial_al,white_al} are the outer and inner surfaces around the grey matter
generated by FreeSurfer. YY=intermediate_al is the node-wise average of the
pial and white surfaces; this surface can be used for single-surface
analyses when surfaces are generated using Caret or BrainVoyager
YY=inflated_alCoMmedial are inflated surfaces with the center of mass
along the medial side of the two hemispheres. These surfaces are suitable
for visualization purposes.
Methods
-------
This dataset contains data from a fingerpress experiment where a participant
(28y right-handed male) pressed buttons on a box with the index and middle
finger of the right hand while fMRI volumes were acquired. Over four runs, in
total sixteen blocks of 72s finger presses were acquired; in each block the
participant pressed the index and middle finger during different `mini-blocks'
of 8 s each. Each block was preceded and followed by a 16s rest period.
Functional data was preprocessed in AFNI with despiking, time slice correction,
motion correction, and scaling to percent signal change by dividing the signal
for each volume by the mean of the run. No spatial smoothing or interpolation
was applied to the functional data, except for interpolation during motion
correction. The preprocessed data was analyzed with a general linear model
(GLM) with separate regressors for each finger and each block (and some
regressors of no interest) to obtain t-values for each finger in each block.
License
-------
The contents are made available by Nikolaas N. Oosterhof <nikolaas.oosterhof
|at| unitn.it> under the Creative Commons CC0 1.0 Universal Public Domain
Dedication ("CC0"). See the LICENSE file for details, or visit
http://creativecommons.org/publicdomain/zero/1.0/deed.en.
Contact
-------
Nikolaas N. Oosterhof <nikolaas.oosterhof |at| unitn.it>
Load functional data
data_path=digit_study_path; data_fn=fullfile(data_path,'glm_T_stats_perblock+orig'); targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2 chunks=floor(((1:32)-1)/4)+1; % half-run: 1 1 1 1 2 2 2 2 3 3 ... 8 8 vol_ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);
Map univariate response data to surface
% this measure averages the data near each node to get a surface dataset radius=0; surf_band_range=[-2 2]; % get voxel data within 2mm from surface surf_def={vertices,faces,[-2 2]}; nbrhood=cosmo_surficial_neighborhood(vol_ds,surf_def,'radius',radius); measure=@(x,opt) cosmo_structjoin('samples',mean(x.samples,2),'sa',x.sa); surf_ds=cosmo_searchlight(vol_ds,nbrhood,measure); fprintf('Univariate surface data:\n'); cosmo_disp(surf_ds);
Warning: found 415 / 5124 center nodes outside the volume, these will be
ignored.
Using 210084 / 210084 voxels in functional volume mask
+00:00:00 [####################] -00:00:00 r=0.00, 3.2 vox
+00:00:00 [####################] -00:00:00
Univariate surface data:
.a
.fdim
.labels
{ 'node_indices' }
.values
{ [ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124 }
.fa
.node_indices
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
.center_ids
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
.samples
[ NaN 0.11 0.722 ... 5.24 -0.0639 -0.632
NaN -0.253 -1.03 ... 7.41 0.846 -1.7
NaN 0.421 -0.71 ... -0.813 -0.349 0.72
: : : : : :
NaN 0.738 -0.143 ... 0.827 -1.21 -2.01
NaN -0.00235 -0.294 ... 0.19 -0.393 0.203
NaN 0.17 -0.765 ... 2.27 -0.393 -0.0536 ]@32x5124
.sa
.chunks
[ 1
1
1
:
8
8
8 ]@32x1
.labels
{ 'fi_i_R1_B01#0_Tstat'
'fi_m_R1_B01#0_Tstat'
'fi_i_R1_B02#0_Tstat'
:
'fi_m_R4_B03#0_Tstat'
'fi_i_R4_B04#0_Tstat'
'fi_m_R4_B04#0_Tstat' }@32x1
.stats
{ 'Ttest(168)'
'Ttest(168)'
'Ttest(168)'
:
'Ttest(168)'
'Ttest(168)'
'Ttest(168)' }@32x1
.targets
[ 1
2
1
:
2
1
2 ]@32x1
Average data in each chunk
for this example only consider the samples in the first condition (targets==1), and average the samples in each chunk
for group analysis: set chunks to (1:nsubj)', assuming each sample is data from a single participant
surf_ds=cosmo_slice(surf_ds,surf_ds.sa.targets==1); surf_ds=cosmo_average_samples(surf_ds); fn_surf_ds=fullfile(output_path, 'digit_target1.niml.dset'); % save to disc cosmo_map2surface(surf_ds, fn_surf_ds); fprintf('Input data saved to %s\n', fn_surf_ds);
Input data saved to /Users/nick/outputs/digit_target1.niml.dset
Run Threshold-Free Cluster Enhancement (TFCE)
% All data is prepared; surf_ds has 8 samples and 5124 nodes. We want to % see if there are clusters that show a significant difference from zero in % their response. Thus, .sa.targets is set to all ones (the same % condition), whereas .sa.chunks is set to (1:8)', indicating that all % samples are assumed to be independent. % % (While this is a within-subject analysis, exactly the same logic can be % applied to a group-level analysis) % define neighborhood for each feature % (cosmo_cluster_neighborhood can be used also for meeg or volumetric % fmri datasets) cluster_nbrhood=cosmo_cluster_neighborhood(surf_ds,... 'vertices',vertices,'faces',faces); fprintf('Cluster neighborhood:\n'); cosmo_disp(cluster_nbrhood); opt=struct(); % number of null iterations. for publication-quality, use >=1000; % 10000 is even better opt.niter=250; % in this case we run a one-sample test against a mean of 0, and it is % necessary to specify the mean under the null hypothesis % (when testing classification accuracies, h0_mean should be set to chance % level, assuming a balanced crossvalidation scheme was used) opt.h0_mean=0; % this example uses the data itself (with resampling) to obtain cluster % statistcs under the null hypothesis. This is (in this case) somewhat % conservative due to how the resampling is performed. % Alternatively, and for better estimates (at the cost of computational % cost), one can generate a set of (say, 50) datasets using permuted data % e.g. using cosmo_randomize_targets), put them in a cell and provide % them as the null argument. opt.null=[]; fprintf('Running multiple-comparison correction with these options:\n'); cosmo_disp(opt); % Run TFCE-based cluster correction for multiple comparisons. % The output has z-scores for each node indicating the probablity to find % the same, or higher, TFCE value under the null hypothesis tfce_ds=cosmo_montecarlo_cluster_stat(surf_ds,cluster_nbrhood,opt);
Cluster neighborhood:
.neighbors
{ [ 1 13 43 163 298 433 ]
[ 2 688 703 823 958 1.09e+03 ]
[ 3 972 987 1.11e+03 1.35e+03 1.48e+03 ]
:
[ 4.55e+03 4.55e+03 5.12e+03 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x7
[ 4.9e+03 4.9e+03 5.12e+03 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x7
[ 4.55e+03 4.55e+03 4.9e+03 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x7 }@5124x1
.fa
.radius
[ 5.74 9.59 5.81 ... 10.6 10.6 9.3 ]@1x5124
.node_indices
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
.sizes
[ 7.29 29.4 9.05 ... 27.4 33.1 31.3 ]@1x5124
.a
.fdim
.values
{ [ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124 }
.labels
{ 'node_indices' }
.origin
.a
.fdim
.labels
{ 'node_indices' }
.values
{ [ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124 }
.fa
.node_indices
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
.center_ids
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
Running multiple-comparison correction with these options:
.niter
[ 250 ]
.h0_mean
[ 0 ]
.null
[ ]
+00:00:09 [####################] -00:00:00 p = 0.004 / 0.004 [+/-0.032] (left/right)
Show results
fprintf('TFCE z-score dataset\n'); cosmo_disp(tfce_ds); nfeatures=size(tfce_ds.samples,2); percentiles=(1:nfeatures)/nfeatures*100; plot(percentiles,sort(tfce_ds.samples)) title('sorted TFCE z-scores'); xlabel('feature percentile'); ylabel('z-score'); nvertices=size(vertices,1); disp_opt=struct(); disp_opt.DataRange=[-2 2]; DispIVSurf(vertices,faces,1:nvertices,tfce_ds.samples',0,disp_opt); % store results fn_tfce_ds=fullfile(output_path, 'digit_target1_tfce.niml.dset'); cosmo_map2surface(tfce_ds, fn_tfce_ds); surf_fn=fullfile(output_path, 'digit_intermediate.asc'); surfing_write(surf_fn,vertices,faces); % show citation information cosmo_check_external('-cite');
TFCE z-score dataset
.samples
[ 0 0 0.745 ... 0 0 0 ]@1x5124
.sa
.stats
{ 'Zscore()' }
.a
.fdim
.labels
{ 'node_indices' }
.values
{ [ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124 }
.fa
.node_indices
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
.center_ids
[ 1 2 3 ... 5.12e+03 5.12e+03 5.12e+03 ]@1x5124
DispIVSurf verbose: z-buffer mode.
DispIVSurf verbose: Scaling ITvect to fit colormap
DispIVSurf verbose: Set hold to off
DispIVSurf verbose: Drawing Patches
DispIVSurf verbose: Applied interp for Shade.
DispIVSurf verbose: Set Edgecolor to black
DispIVSurf verbose: Set view to 2D
DispIVSurf verbose: Displaying colormap
DispIVSurf verbose: Seting Axes properties
DispIVSurf verbose: Done
If you use CoSMoMVPA and/or some other toolboxes for a publication, please cite:
Z. Saad, G. Chen. AFNI Matlab library. available online from https://github.com/afni/AFNI
N. N. Oosterhof, T. Wiestler, J. Diedrichsen (2011). A comparison of volume-based and surface-based multi-voxel pattern analysis. Neuroimage 56 (2), 593-600. Surfing toolbox available online from http://github.com/nno/surfing
Gabriel Peyre. Toolbox Fast Marching - A toolbox Fast Marching and level sets computations [https://www.ceremade.dauphine.fr/~peyre/matlab/fast-marching/content.html]. toolbox fast marching [included in surfing] available online from http://github.com/nno/surfing
N. N. Oosterhof, A. C. Connolly, J. V. Haxby (2016). CoSMoMVPA: multi-modal multivariate pattern analysis of neuroimaging data in Matlab / GNU Octave. Frontiers in Neuroinformatics, doi:10.3389/fninf.2016.00027.. CoSMoMVPA toolbox available online from http://cosmomvpa.org
The Mathworks, Natick, MA, United States. Matlab 9.1.0.441655 (R2016b) (September 7, 2016). available online from http://www.mathworks.com
