clusters of significant brain activation, max z score, MNI coordinates

[Franziska Motka]

Hi everyone, 

I wrote the following code to test the contrast [smoking-neutral] for each voxel against zero with TFCE correction:

%%% Hypothesis 1: For which voxels is the contrast [smoking-neutral] sig. different from zero? --> t-test (two-sided)
p_values_1 = NaN(size(ds_allsubs_pre.samples, 2), 1); % create an array with # rows = # voxel
[~, p_values_1,~,Stats_1] = ttest(ds_allsubs_pre.samples);

%%% Threshold-free cluster enhancement (TFCE)
% Find neighbors of each voxel
nh = cosmo_cluster_neighborhood(ds_allsubs_pre,'progress', false); %'progress', false deactivates code progress information
% Set options for monte carlo based clustering statistic
opt = struct(); % creates empty structure to save options for the TFCE analysis
opt.cluster_stat ='tfce'; % determine the TFCE method to calculate the sig. of clusters
opt.niter = 10000;  % Use 100 when testing the code and 10000 for results % # iterations for the Monte-Carlo-Simulation, to estimate the distribution of the cluster statistic for the null hypothesis                                  
opt.h0_mean = 0; % determine the mean of the null distribution for the Monte-Carlo-Simulation; =0 -> null hypothesis, that there are no systematic effects    
% Apply cluster-based correction
ds_allsubs_pre.sa.chunks = (1:98)'; %98 participants
ds_allsubs_pre.sa.targets = repmat(1, 98, 1); %one condition
z_ds=cosmo_montecarlo_cluster_stat(ds_allsubs_pre,nh,opt); % function performs the Monte-Carlo-Simulation, to determine the sig. of activity clusters and stores the results in a new data structure z_ds which contains the calcuated Z-value for each voxel
% Map it to a nifti file to visualize

cosmo_map2fmri(z_ds, 'Hypothesis1_TFCEcorrected.nii')

So i have a nii file which contains a z score for each voxel. Now I would like to get clusters of significant (abs(z score) > 1.96) brain activation, the max z score within each cluster and the MNI coordinates for the voxel with the max z score. Can anyone help with that? Here is what I tried so far, but this can not be correct since all z scores are zero and I receive way to many clusters (> 10.000).

%%% Identify clusters with CosmoMVPA
% Load the TFCE corrected NIfTI file
z_ds = cosmo_fmri_dataset('Hypothesis1_TFCEcorrected.nii');
% Define neighborhood
nh = cosmo_cluster_neighborhood(z_ds, 'progress', false);
% Define options for clusterizing
opt = struct();
opt.cluster_stat = 'maximal'; % Find maximal value in each cluster
% Apply clusterizing
nh_mat=cosmo_convert_neighborhood(nh,'matrix');
clusters_ds = cosmo_clusterize(z_ds.samples, nh_mat);
% Initialize array to store maximal values in each cluster
max_values_per_cluster = NaN(1, numel(clusters_ds));
% Iterate over each cluster
for i = 1:numel(clusters_ds)
    % Extract indices of voxels belonging to the current cluster
    voxel_indices = find(clusters_ds{i});
    % Extract Z-scores for voxels in the current cluster
    cluster_z_scores = z_ds.samples(voxel_indices);
    % Find the maximum Z-score in the cluster
    max_z_score = max(cluster_z_scores);
    % Store the maximum Z-score for the current cluster
    max_values_per_cluster(i) = max_z_score;
end
% Display the maximal values in each cluster
disp('Maximal values in each cluster:');
disp(max_values_per_cluster);

Thank you for any help!

Best, 

Franziska

[Nick Oosterhof]

Greetings,

> On 6 Apr 2024, at 13:33, Franziska Motka <franzis...@gmail.com> wrote:
>
> I wrote the following code to test the contrast [smoking-neutral] for each voxel against zero with TFCE correction:

> [snip]

>
> So i have a nii file which contains a z score for each voxel. Now I would like to get clusters of significant (abs(z score) > 1.96) brain activation, the max z score within each cluster and the MNI coordinates for the voxel with the max z score. Can anyone help with that? Here is what I tried so far, but this can not be correct since all z scores are zero and I receive way to many clusters (> 10.000).
>
> %%% Identify clusters with CosmoMVPA
> % Load the TFCE corrected NIfTI file
> z_ds = cosmo_fmri_dataset('Hypothesis1_TFCEcorrected.nii');
> % Define neighborhood
> nh = cosmo_cluster_neighborhood(z_ds, 'progress', false);
> % Define options for clusterizing
> opt = struct();
> opt.cluster_stat = 'maximal'; % Find maximal value in each cluster
> % Apply clusterizing
> nh_mat=cosmo_convert_neighborhood(nh,'matrix');
> clusters_ds = cosmo_clusterize(z_ds.samples, nh_mat);

Just a remark: cluster_ds is not a dataset as used in CoSMoMVPA; it is a cell. I find the _ds suffix a bit confusing.

> % Initialize array to store maximal values in each cluster
> max_values_per_cluster = NaN(1, numel(clusters_ds));
> % Iterate over each cluster
> for i = 1:numel(clusters_ds)
> % Extract indices of voxels belonging to the current cluster
> voxel_indices = find(clusters_ds{i});

Why are you using “find”? The vectors in cluster_ds should already contain the indices themselves.

> % Extract Z-scores for voxels in the current cluster
> cluster_z_scores = z_ds.samples(voxel_indices);
> % Find the maximum Z-score in the cluster
> max_z_score = max(cluster_z_scores);

Just a remark: max can return second output with the index of the maximum value.

> % Store the maximum Z-score for the current cluster
> max_values_per_cluster(i) = max_z_score;
> end
> % Display the maximal values in each cluster
> disp('Maximal values in each cluster:');
> disp(max_values_per_cluster);

You may want to look at cosmo_find_local_extrema. The documentation also shows how to get voxel maxima in world coordinates

Just a small ‘warning’ for finding maxima from Monte Carlo stats: with strong effects it can happen that a whole group of connected voxels shows the same z_max which is only limited by the number of iterations. In those cases the “max” function may return a voxel on the edge of the voxel group, rather than the center. Computing a center of mass of such a voxel group will be more representative when reporting world coordinates.

[Franziska Motka]

Thank you for your answer and comments! I will try it and may come back with some questions.

Best, 
Franziska