From 48145b411919e68dc5b79478468a422b4003e7dd Mon Sep 17 00:00:00 2001 From: Tor Wager Date: Sun, 8 May 2022 14:15:24 -0400 Subject: [PATCH] Updates to Meta_cluster_tools barplots and Meta_results_batch_script Some compat with new Matlab versions (2020b). CanlabCore also updated. --- densityUtility3/Meta_cluster_tools.m | 44 ++- .../Meta_results_batch_script.asv | 338 ++++++++++++++++++ .../Meta_results_batch_script.m | 221 +++++++----- 3 files changed, 496 insertions(+), 107 deletions(-) create mode 100644 densityUtility3/plotting_functions/Meta_results_batch_script.asv diff --git a/densityUtility3/Meta_cluster_tools.m b/densityUtility3/Meta_cluster_tools.m index bc8d3bf..8e258cb 100644 --- a/densityUtility3/Meta_cluster_tools.m +++ b/densityUtility3/Meta_cluster_tools.m @@ -40,11 +40,11 @@ % ------------------------------------------------------ % count studies by condition and plot [optional] % -% [prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition',dat,Xi,w,doplot,[xnames],[seriesnames], [colors]) +% [prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition',dat,Xi,w,doplot,[xnames],[seriesnames], [colors]) % -% [prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition',studybyset,MC_Setup.Xi,MC_Setup.wts,1) +% [prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition',studybyset,MC_Setup.Xi,MC_Setup.wts,1) % -% [prop_by_condition,se,num_by_condition,n] = ... +% [prop_by_condition,se,num_by_condition,n, table_obj] = ... % Meta_cluster_tools('count_by_condition',studybyroi,MC_Setup.Xi,MC_Setup.wts,1, ... % {'Right' 'Left'},MC_Setup.connames(1:5),{[1 0 0] [0 1 0] [1 0 1] [1 1 0] [0 0 1]}); % @@ -52,7 +52,7 @@ % nms = SOMResults.alltasknms(9:13) % w = SOMResults.w; % colors = {[1 0 0] [0 1 0] [1 0 1] [1 1 0] [0 0 1]}; -% [prop_by_condition,se,num_by_condition,n] = ... +% [prop_by_condition,se,num_by_condition,n, table_obj] = ... % Meta_cluster_tools('count_by_condition',studybyroi,Xi,w,1, ... % {'Right' 'Left'},nms,colors); % @@ -65,7 +65,7 @@ % R = Meta_Chisq_new('compute',MC_Setup,'mask',mask); % R = Meta_Chisq_new('write',R); % [studybyroi,studybyset] = Meta_cluster_tools('getdata',cl,R.dat,R.volInfo); -% [prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition',studybyset,R.Xi,R.w,1); +% [prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition',studybyset,R.Xi,R.w,1); % ------------------------------------------------------ switch meth @@ -92,14 +92,14 @@ end case 'count_by_condition' - %[prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition',dat,Xi,w,doplot,[xnames],[seriesnames], [colors]) - %[prop_by_condition,se,num_by_condition,n] = count_by_condition(dat,Xi,w,doplot) + %[prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition',dat,Xi,w,doplot,[regionnames],[seriesnames], [colors]) + %[prop_by_condition,se,num_by_condition,n, table_obj] = count_by_condition(dat,Xi,w,doplot) if length(varargin) < 4, varargin{4} = 0; end if length(varargin) < 5, varargin{5} = []; end if length(varargin) < 6, varargin{6} = []; end if length(varargin) < 7, varargin{7} = []; end %colors - [varargout{1},varargout{2},varargout{3},varargout{4}] = count_by_condition(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},varargin{6},varargin{7}); + [varargout{1},varargout{2},varargout{3},varargout{4}, varargout{5}] = count_by_condition(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5},varargin{6},varargin{7}); case 'activation_table' @@ -177,7 +177,7 @@ [studybyroi,studybyset] = Meta_cluster_tools('getdata',cl,MC_Setup.unweighted_study_data,MC_Setup.volInfo); disp('Counting studies by condition') -[prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition',studybyroi,MC_Setup.Xi,MC_Setup.wts,doplot); +[prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition',studybyroi,MC_Setup.Xi,MC_Setup.wts,doplot); % get field names for conditions fnames = MC_Setup.Xinms; @@ -255,7 +255,7 @@ % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- -function [prop_by_condition,se,num_by_condition,n] = count_by_condition(dat,Xi,w,doplot,varargin) +function [prop_by_condition,se,num_by_condition,n, table_obj] = count_by_condition(dat,Xi,w,doplot,varargin) if nargin < 4, doplot = 0; end @@ -273,12 +273,14 @@ se = ( (prop_by_condition .* (1-prop_by_condition) ) ./ n ).^.5; if doplot - xnames = []; - seriesnames = []; - mycolors = []; - if length(varargin) > 0, xnames = varargin{1}; end - if length(varargin) > 1, seriesnames = varargin{2}; end - if length(varargin) > 2, mycolors = varargin{3}; end + + xnames = repmat({'Region'}, 1, size(prop_by_condition, 1)); + seriesnames = repmat({'Cond'}, 1, size(prop_by_condition, 2)); + mycolors = scn_standard_colors(size(prop_by_condition, 2)); + + if length(varargin) > 0 && ~isempty(varargin{1}), xnames = varargin{1}; end + if length(varargin) > 1 && ~isempty(varargin{2}), seriesnames = varargin{2}; end + if length(varargin) > 2 && ~isempty(varargin{3}), mycolors = varargin{3}; end create_figure('barplot'); fprintf('Sample size is %3.0f\n', size(Xi, 1)); @@ -286,8 +288,18 @@ barplot_grouped(prop_by_condition, se, xnames, seriesnames, 'inputmeans', 'colors', mycolors); ylabel('Proportion of studies activating'); + end +% Table of contrast data +table_obj = table(); + +propdat = mat2cell(prop_by_condition, size(prop_by_condition, 1), ones(1, length(seriesnames))); +sedat = mat2cell(se, size(prop_by_condition, 1), ones(1, length(seriesnames))); + +senames = cellfun(@(x) ['SE_' x], seriesnames, 'UniformOutput', false); +table_obj = table(xnames', propdat{:}, sedat{:}, 'Variablenames', [{'Region'} seriesnames senames]); + end diff --git a/densityUtility3/plotting_functions/Meta_results_batch_script.asv b/densityUtility3/plotting_functions/Meta_results_batch_script.asv new file mode 100644 index 0000000..cf4d0a1 --- /dev/null +++ b/densityUtility3/plotting_functions/Meta_results_batch_script.asv @@ -0,0 +1,338 @@ +%% CANlab MKDA meta-analysis results report +% "What's in the box": A series of tables and plots for a standard +% coordinate-based meta-analysis created with Meta_Activation_FWE in the +% CANlab MKDA meta-analysis toolbox. +% +% The script Meta_results_batch_script is meant to be run stand-alone or +% used in publish_meta_analysis_report to create an HTML report with +% results and plots. + +% Display helper functions: Called by later scripts +% -------------------------------------------------------- + +dashes = '----------------------------------------------'; +printstr = @(dashes) disp(dashes); +printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes); + +%% Results orthviews and tables with multi-threshold and subclusters +% Multi-threshold: Show results with voxel-wise FWER and 'stringent' 'medium' and +% 'lenient' cluster extent-corrected results in different colors. +% +% Tables show numbers for contiguous clusters and arrows for sub-peaks +% within each contiguous cluster using SPM software's local maximum-finding +% algorithm. + +close all + +cl = Meta_Activation_FWE('results', 1); + +drawnow, snapnow + +% Preparing regions so that higher thresholds are included in region maps +% for lower thresholds (this is not done in the code above). + +h = fmri_data('Activation_FWE_height.img', 'noverbose'); +e1 = fmri_data('Activation_FWE_extent_stringent.img', 'noverbose'); +e2 = fmri_data('Activation_FWE_extent_medium.img', 'noverbose'); +e3 = fmri_data('Activation_FWE_extent_lenient.img', 'noverbose'); + +mydat = h; +r_height = region(mydat); + +mydat.dat = mydat.dat + e1.dat; +r_stringent = region(mydat); + +mydat.dat = mydat.dat + e2.dat; +r_medium = region(mydat); + +mydat.dat = mydat.dat + e3.dat; +r_lenient = region(mydat); + +% Set up to extract +% Study proportion data for voxel-wise FWER-corrected regions +% First, we will load the MC_Info.mat file, which contains lots of +% information about the analysis. The variable of main interest in this + +% file is MC_Setup, which contains the contrast indicator maps. +% This is stored in MC_Setup.unweighted_study_data + +load MC_Info + +% MC_Setup = +% unweighted_study_data: [231202×18 double] % Voxels x contrast indicator maps (1/0) +% volInfo: [1×1 struct] % Info for mapping back into brain space +% n: [2 3 2 2 1 2 1 2 4 2 2 3 3 2 2 8 2 1] % Coordinates per study +% wts: [18×1 double] % Contrast weights +% r: 5 % Radius in voxels +% cl: {1×18 cell} % contiguous blobs for each contrast, for Monte Carlo + +% Set up contrasts + +hascontrasts = isfield(MC_Setup, 'contrasts'); + +if hascontrasts + % These do not have the same (strange) property as the main activation + % images of excluding voxels at a more stringent threshold + + posfilenames{1} = [MC_Setup.connames{1} '_Pos_FWE_height.img']; + negfilenames{1} = [MC_Setup.connames{1} '_Neg_FWE_height.img']; + + posfilenames{2} = [MC_Setup.connames{1} '_Pos_FWE_extent_stringent.img']; + negfilenames{2} = [MC_Setup.connames{1} '_Neg_FWE_extent_stringent.img']; + + posfilenames{3} = [MC_Setup.connames{1} '_Pos_FWE_extent_medium.img']; + negfilenames{3} = [MC_Setup.connames{1} '_Neg_FWE_extent_medium.img']; + + posfilenames{4} = [MC_Setup.connames{1} '_Pos_FWE_extent.img']; + negfilenames{4} = [MC_Setup.connames{1} '_Neg_FWE_extent.img']; + + for i = 1:4 + r_pos_contrast{i} = region(posfilenames{i}, 'noverbose'); + r_neg_contrast{i} = region(negfilenames{i}, 'noverbose'); + end + +end + + + +%% Voxel-wise FWER: Table and slice montage with autolabeled regions + +printhdr('Voxel-wise FWER'); + +% Displays montages of the region object and a table of results. +% Table creation adds auto-labeled region names to the .shorttitle field of the region object. +% Extracted 1/0 data for activation by each contrast are saved in the region_object .dat field. + +r_height = extract_data_and_show_region_results(r_height, MC_Setup); + + +%% Stringent .001 cluster-extent corrected: Table and slice montage with autolabeled regions + +printhdr('Voxel-wise FWER + stringent primary threshold (.001) cluster-extent'); + +r_stringent = extract_data_and_show_region_results(r_stringent, MC_Setup); + + +%% Medium .01 cluster-extent corrected: Table and slice montage with autolabeled regions +% Excluding voxels already reported at more stringent/localizable thresholds + +printhdr('Voxel-wise FWER + medium primary threshold (.01) cluster-extent'); + +r_medium = extract_data_and_show_region_results(r_medium, MC_Setup); + +%% Lenient .05 cluster-extent corrected: Table and slice montage with autolabeled regions +% Excluding voxels already reported at more stringent/localizable thresholds + +printhdr('Voxel-wise FWER + lenient primary threshold (.05) cluster-extent'); + +r_lenient = extract_data_and_show_region_results(r_lenient, MC_Setup); + + +close all % save memory, etc. + + +% Save region_objects.mat file with regions and extracted data + +save region_objects r_height r_stringent r_medium r_lenient + + +%% All thresholds: Montage and individual points +% +% % Reload masked activation map and display slice montage +% % The code below uses the CANlab object-oriented tools, in +% % CANlab_Core_Tools repository. There are many more options for +% % visualization and data analysis + +disp('Activation_FWE_all.img') + +% load saved results mask (union of all thresholds) +img = fmri_data('Activation_FWE_all.img', 'noverbose'); % Create an fmri_data-class object + +create_figure('slice montage'); axis off +o2 = montage(img); % Display montage, and return an fmridisplay object called o2 + +drawnow, snapnow + +% Print transparent blobs and activation points on slice montage + +% load database with coordinates (also, % fix for bad xyz) +load SETUP DB +if ~isfield(DB, 'xyz') || size(DB.xyz, 2) ~= 3, DB.xyz = [DB.x DB.y DB.z]; end + +o2 = removeblobs(o2); +o2 = addblobs(o2, r_lenient, 'trans', 'transvalue', .6); +o2 = addpoints(o2, DB.xyz, 'Color', 'none', 'MarkerFaceColor', [1 1 0], 'Marker', 'o', 'MarkerSize', 4); + +drawnow, snapnow + +o2 = removeblobs(o2); + +drawnow, snapnow + +close all + +%% Contrasts + + +if hascontrasts + + clear poseffectname negeffectname + poseffectname{1} = sprintf('Contrast %s: Voxel FWER Pos Effects', MC_Setup.connames{1}); + negeffectname{1} = sprintf('Contrast %s: Voxel FWER Neg Effects', MC_Setup.connames{1}); + + poseffectname{2} = sprintf('Contrast %s: Voxel FWER + stringent cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{2} = sprintf('Contrast %s: Voxel FWER + stringent cluster-ext Neg Effects', MC_Setup.connames{1}); + + poseffectname{3} = sprintf('Contrast %s: Voxel FWER + medium cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{3} = sprintf('Contrast %s: Voxel FWER + medium cluster-ext Neg Effects', MC_Setup.connames{1}); + + poseffectname{4} = sprintf('Contrast %s: Voxel FWER + lenient cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{4} = sprintf('Contrast %s: Voxel FWER + lenient cluster-ext Neg Effects', MC_Setup.connames{1}); + + for i = 1:4 + + printhdr(poseffectname{i}); + + r_pos_contrast{i} = extract_data_and_show_region_results(r_pos_contrast{i}, MC_Setup); + + if ~isempty(r_pos_contrast{i}), r_pos_contrast{i}(1).title = poseffectname{i}; end + + % + + printhdr(negeffectname{i}); + + r_neg_contrast{i} = extract_data_and_show_region_results(r_neg_contrast{i}, MC_Setup); + + if ~isempty(r_neg_contrast{i}), r_neg_contrast{i}(1).title = negeffectname{i}; end + + end % loop + + + % Save region_objects.mat file with regions and extracted data + + save region_objects -append r_pos_contrast r_neg_contrast + + +end % if contrasts + +%% Surface renderings: Points on surfaces + +% No longer needed +% Create a cutaway surface with blobs +% surface(img, 'cutaway', 'ycut_mm', -30); +% drawnow, snapnow + +% Plot points as spheres on a canonical surface: + +plot_points_on_surface2(DB.xyz, {[1 1 0]}); + +drawnow, snapnow + +close all + + +%% Neurochemistry: PET binding map similarity + +% [pet, petlabels] = load_image_set('receptorbinding'); +% +% stats = image_similarity_plot(h, 'mapset', pet, 'cosine_similarity', 'bicolor', 'networknames', pet.metadata_table.target, 'plotstyle', 'polar', 'average'); +% +% title('PET tracer profile'); + + + + +% SUBFUNCTIONS +% --------------------------------------------------------------------------- + +function r_height = extract_data_and_show_region_results(r_height, MC_Setup) +% Displays montages of the region object and a table of results. +% Table creation adds auto-labeled region names to the .shorttitle field of the region object. +% Extracted 1/0 data for activation by each contrast are saved in the region_object .dat field. + +if ~isempty(r_height) + + % warning off, r = cluster2region(cl{1}); warning on + + montage(r_height, 'colormap'); + create_figure('surface'); axis off; surface(r_height, 'cutaway'); + drawnow, snapnow + + [rpos, rneg] = table(r_height); + montage([rpos rneg], 'regioncenters', 'colormap'); + + drawnow, snapnow + + r_height = [rpos rneg]; % save with region names, also for data extraction below + + r_height = extract_data_and_plot(r_height, MC_Setup); % extract data from regions + +else + + disp('No significant results.'); + +end + +end % subfunction + + + +function region_obj = extract_data_and_plot(region_obj, MC_Setup) +% This extracts data from a set of regions and makes a barplot of +% activation proportions. Extracted 1/0 data for activation by each +% contrast are saved in the region_object .dat field. + +indicator_maps = MC_Setup.unweighted_study_data; + +if ~isempty(region_obj) + + % Extract indicator for contrasts that activate within 10 mm (5 vox) of + % significant voxels region object r + studybyroi = Meta_cluster_tools('getdata', region_obj, indicator_maps, MC_Setup.volInfo); + + for i = 1:length(region_obj) + region_obj(i).dat = studybyroi(:, i); + end + + % studybyroi: contrasts x regions, values 1/0 for whether each contrast activates the region + % studybyset: contrasts x 1, values 1/0 for whether each contrast activates any region in the set + + create_figure('Activation_proportions') + + [n, k] = size(studybyroi); + prop_by_condition = sum(studybyroi) ./ n; % Proportion of contrasts activating each ROI + se = ( (prop_by_condition .* (1-prop_by_condition) ) ./ n ).^.5; % Standard Error based on binomial distribution + + han = bar(prop_by_condition); + ehan = errorbar(prop_by_condition, se); + + labels = format_strings_for_legend({region_obj.shorttitle}); + + set(han, 'EdgeColor', [0 0 .5], 'FaceColor', [.3 .3 .6]); + set(ehan, 'Color', [0 0 .5], 'LineStyle', 'none', 'LineWidth', 3); + set(gca, 'XTick', 1:k, 'XLim', [.5 k+.5], 'XTickLabel', labels); + set(gca, 'XTickLabelRotation', 70); + ylabel('Proportion of studies activating'); + + % Conditions for contrast plots, if entered + if isfield(MC_Setup, 'Xi') + + colors = seaborn_colors(size(MC_Setup.Xi, 2)); + regionnames = cat(2, {region_obj.shorttitle}); + + [prop_by_condition,se,num_by_condition,n] = Meta_cluster_tools('count_by_condition', studybyroi, MC_Setup.Xi, MC_Setup.wts, 1, regionnames, MC_Setup.Xinms, colors); + + % Table of contrast data + propdat = mat2cell(prop_by_condition, size(prop_by_condition, 1), ones(1, length(MC_Setup.Xinms))); + sedat = mat2cell(se, size(prop_by_condition, 1), ones(1, length(MC_Setup.Xinms))); + + senames = cellfun(@(x) ['SE_' x], MC_Setup.Xinms, 'UniformOutput', false); + t = table(regionnpropdat{:}, sedat{:}, 'Variablenames', [MC_Setup.Xinms senames], 'Rownames', regionnames); + + end % Conditions/contrasts plot + + drawnow, snapnow + +end % not empty + +end % function diff --git a/densityUtility3/plotting_functions/Meta_results_batch_script.m b/densityUtility3/plotting_functions/Meta_results_batch_script.m index eff4772..f41bc22 100644 --- a/densityUtility3/plotting_functions/Meta_results_batch_script.m +++ b/densityUtility3/plotting_functions/Meta_results_batch_script.m @@ -66,62 +66,51 @@ % r: 5 % Radius in voxels % cl: {1×18 cell} % contiguous blobs for each contrast, for Monte Carlo +% Set up contrasts -%% Voxel-wise FWER: Table and slice montage with autolabeled regions - -printhdr('Voxel-wise FWER'); +hascontrasts = isfield(MC_Setup, 'contrasts'); -if ~isempty(r_height) +if hascontrasts + % These do not have the same (strange) property as the main activation + % images of excluding voxels at a more stringent threshold - % warning off, r = cluster2region(cl{1}); warning on + posfilenames{1} = [MC_Setup.connames{1} '_Pos_FWE_height.img']; + negfilenames{1} = [MC_Setup.connames{1} '_Neg_FWE_height.img']; - montage(r_height, 'colormap'); - create_figure('surface'); axis off; surface(r_height, 'cutaway'); - drawnow, snapnow - - [rpos, rneg] = table(r_height); - montage([rpos rneg], 'regioncenters', 'colormap'); - - drawnow, snapnow - - r_height = [rpos rneg]; % save with region names, also for data extraction below + posfilenames{2} = [MC_Setup.connames{1} '_Pos_FWE_extent_stringent.img']; + negfilenames{2} = [MC_Setup.connames{1} '_Neg_FWE_extent_stringent.img']; - r_height = extract_data_and_plot(r_height, MC_Setup); % extract data from regions + posfilenames{3} = [MC_Setup.connames{1} '_Pos_FWE_extent_medium.img']; + negfilenames{3} = [MC_Setup.connames{1} '_Neg_FWE_extent_medium.img']; -else + posfilenames{4} = [MC_Setup.connames{1} '_Pos_FWE_extent.img']; + negfilenames{4} = [MC_Setup.connames{1} '_Neg_FWE_extent.img']; - disp('No significant results.'); + for i = 1:4 + r_pos_contrast{i} = region(posfilenames{i}, 'noverbose'); + r_neg_contrast{i} = region(negfilenames{i}, 'noverbose'); + end end + +%% Voxel-wise FWER: Table and slice montage with autolabeled regions + +printhdr('Voxel-wise FWER'); + +% Displays montages of the region object and a table of results. +% Table creation adds auto-labeled region names to the .shorttitle field of the region object. +% Extracted 1/0 data for activation by each contrast are saved in the region_object .dat field. + +r_height = extract_data_and_show_region_results(r_height, MC_Setup); + + %% Stringent .001 cluster-extent corrected: Table and slice montage with autolabeled regions printhdr('Voxel-wise FWER + stringent primary threshold (.001) cluster-extent'); -if ~isempty(r_stringent) - - %r = cluster2region(cl{2}); - - montage(r_stringent, 'colormap'); - create_figure('surface'); axis off; surface(r_stringent, 'cutaway'); - - drawnow, snapnow - - [rpos, rneg] = table(r_stringent); - montage([rpos rneg], 'regioncenters', 'colormap'); - - drawnow, snapnow - - r_stringent = [rpos rneg]; % save with region names, also for data extraction below - - r_stringent = extract_data_and_plot(r_stringent, MC_Setup); % extract data from regions - -else - - disp('No significant results.'); - -end +r_stringent = extract_data_and_show_region_results(r_stringent, MC_Setup); %% Medium .01 cluster-extent corrected: Table and slice montage with autolabeled regions @@ -129,59 +118,15 @@ printhdr('Voxel-wise FWER + medium primary threshold (.01) cluster-extent'); -if ~isempty(r_medium) - - %warning off, r = cluster2region(cl{3}); warning on - - montage(r_medium, 'colormap'); - create_figure('surface'); axis off; surface(r_medium, 'cutaway'); - - drawnow, snapnow - - [rpos, rneg] = table(r_medium); - montage([rpos rneg], 'regioncenters', 'colormap'); - - drawnow, snapnow - - r_medium = [rpos rneg]; % save with region names, also for data extraction below - - r_medium = extract_data_and_plot(r_medium, MC_Setup); % extract data from regions - -else - - disp('No significant results.'); - -end +r_medium = extract_data_and_show_region_results(r_medium, MC_Setup); %% Lenient .05 cluster-extent corrected: Table and slice montage with autolabeled regions % Excluding voxels already reported at more stringent/localizable thresholds printhdr('Voxel-wise FWER + lenient primary threshold (.05) cluster-extent'); -if ~isempty(r_lenient) - - %warning off, r = cluster2region(cl{4}); warning on - - montage(r_lenient, 'colormap'); - create_figure('surface'); axis off; surface(r_lenient, 'cutaway'); - - drawnow, snapnow - - [rpos, rneg] = table(r_lenient); - montage([rpos rneg], 'regioncenters', 'colormap'); - - drawnow, snapnow - - r_lenient = [rpos rneg]; % save with region names, also for data extraction below - - r_lenient = extract_data_and_plot(r_lenient, MC_Setup); % extract data from regions - - -else - - disp('No significant results.'); - -end +r_lenient = extract_data_and_show_region_results(r_lenient, MC_Setup); + close all % save memory, etc. @@ -210,10 +155,10 @@ % Print transparent blobs and activation points on slice montage -% load database with coordinates +% load database with coordinates (also, % fix for bad xyz) load SETUP DB -if ~isfield(DB, 'xyz'), DB.xyz = [DB.x DB.y DB.z]; end - +if ~isfield(DB, 'xyz') || size(DB.xyz, 2) ~= 3, DB.xyz = [DB.x DB.y DB.z]; end + o2 = removeblobs(o2); o2 = addblobs(o2, r_lenient, 'trans', 'transvalue', .6); o2 = addpoints(o2, DB.xyz, 'Color', 'none', 'MarkerFaceColor', [1 1 0], 'Marker', 'o', 'MarkerSize', 4); @@ -226,7 +171,49 @@ close all +%% Contrasts + + +if hascontrasts + + clear poseffectname negeffectname + poseffectname{1} = sprintf('Contrast %s: Voxel FWER Pos Effects', MC_Setup.connames{1}); + negeffectname{1} = sprintf('Contrast %s: Voxel FWER Neg Effects', MC_Setup.connames{1}); + + poseffectname{2} = sprintf('Contrast %s: Voxel FWER + stringent cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{2} = sprintf('Contrast %s: Voxel FWER + stringent cluster-ext Neg Effects', MC_Setup.connames{1}); + + poseffectname{3} = sprintf('Contrast %s: Voxel FWER + medium cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{3} = sprintf('Contrast %s: Voxel FWER + medium cluster-ext Neg Effects', MC_Setup.connames{1}); + + poseffectname{4} = sprintf('Contrast %s: Voxel FWER + lenient cluster-ext Pos Effects', MC_Setup.connames{1}); + negeffectname{4} = sprintf('Contrast %s: Voxel FWER + lenient cluster-ext Neg Effects', MC_Setup.connames{1}); + + for i = 1:4 + + printhdr(poseffectname{i}); + + r_pos_contrast{i} = extract_data_and_show_region_results(r_pos_contrast{i}, MC_Setup); + + if ~isempty(r_pos_contrast{i}), r_pos_contrast{i}(1).title = poseffectname{i}; end + + % + + printhdr(negeffectname{i}); + + r_neg_contrast{i} = extract_data_and_show_region_results(r_neg_contrast{i}, MC_Setup); + + if ~isempty(r_neg_contrast{i}), r_neg_contrast{i}(1).title = negeffectname{i}; end + + end % loop + + + % Save region_objects.mat file with regions and extracted data + + save region_objects -append r_pos_contrast r_neg_contrast + +end % if contrasts %% Surface renderings: Points on surfaces @@ -258,8 +245,42 @@ % SUBFUNCTIONS % --------------------------------------------------------------------------- +function r_height = extract_data_and_show_region_results(r_height, MC_Setup) +% Displays montages of the region object and a table of results. +% Table creation adds auto-labeled region names to the .shorttitle field of the region object. +% Extracted 1/0 data for activation by each contrast are saved in the region_object .dat field. + +if ~isempty(r_height) + + % warning off, r = cluster2region(cl{1}); warning on + + montage(r_height, 'colormap'); + create_figure('surface'); axis off; surface(r_height, 'cutaway'); + drawnow, snapnow + + [rpos, rneg] = table(r_height); + montage([rpos rneg], 'regioncenters', 'colormap'); + + drawnow, snapnow + + r_height = [rpos rneg]; % save with region names, also for data extraction below + + r_height = extract_data_and_plot(r_height, MC_Setup); % extract data from regions + +else + + disp('No significant results.'); + +end + +end % subfunction + + function region_obj = extract_data_and_plot(region_obj, MC_Setup) +% This extracts data from a set of regions and makes a barplot of +% activation proportions. Extracted 1/0 data for activation by each +% contrast are saved in the region_object .dat field. indicator_maps = MC_Setup.unweighted_study_data; @@ -293,6 +314,24 @@ set(gca, 'XTickLabelRotation', 70); ylabel('Proportion of studies activating'); + % Conditions for contrast plots, if entered + if isfield(MC_Setup, 'Xi') + + colors = seaborn_colors(size(MC_Setup.Xi, 2)); + regionnames = cat(2, {region_obj.shorttitle}); + + [prop_by_condition,se,num_by_condition,n, table_obj] = Meta_cluster_tools('count_by_condition', studybyroi, MC_Setup.Xi, MC_Setup.wts, 1, regionnames, MC_Setup.Xinms, colors); + + % Table of contrast data - now included in Meta_cluster_tools + disp(table_obj) + + disp('Attaching table obj to region(1).custom_info1') + region_obj(1).custom_info1 = table_obj; + + end % Conditions/contrasts plot + + drawnow, snapnow + end % not empty end % function