From 9fd1dfba5d80dd38aaa53a9fd5a9b85c40989bf7 Mon Sep 17 00:00:00 2001 From: arno Date: Tue, 23 May 2023 16:45:31 -0700 Subject: [PATCH 01/25] snippet calculation --- plotconnectivity.m | 5 ++-- pop_roi_connect.m | 72 +++++++++++++++++++++++++-------------------- roi_lookupregions.m | 7 +++-- 3 files changed, 48 insertions(+), 36 deletions(-) diff --git a/plotconnectivity.m b/plotconnectivity.m index 1533e72..d7b14ed 100644 --- a/plotconnectivity.m +++ b/plotconnectivity.m @@ -80,7 +80,7 @@ if isempty(g.labels) if strcmpi(g.brainimg, 'on') disp('Cannot plot on brain with area labels') - g.brainimg = 'off'; + g.brainimg = 'off';g end end @@ -193,6 +193,7 @@ % rename labels % ------------- +if ~strcmpi(g.brainimg, 'off') % CTC added .... to branch around code if not plotting if isempty(g.labels) for iPnt = 1:length(anglesInit) g.labels{iPnt} = sprintf(' Area %d', iPnt); @@ -205,7 +206,7 @@ g.labels{iPnt} = [ ' ' g.labels{iPnt} ]; end end - +end % CTC added .... to branch around above code if not plotting warning off; if isempty(limits) limits(2) = max(array(:)); diff --git a/pop_roi_connect.m b/pop_roi_connect.m index b386718..f8a0987 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -28,6 +28,7 @@ % 'MIM' : Multivariate Interaction Measure for each ROI % 'MIC' : Maximized Imaginary Coherency for each ROI % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast comptuation). Default is 'off'. % 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' % 'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). If specified, the signal is zero padded accordingly. Default is 0 (means no padding. @@ -144,6 +145,7 @@ 'naccu' 'integer' { } 0; 'methods' 'cell' { } {}; 'snippet' 'string' { 'on', 'off' } 'off'; + 'firstsnippet' 'string' { 'on', 'off' } 'off'; 'nepochs' 'real' {} []; 'snip_length' 'integer' { } 60; 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; @@ -205,43 +207,49 @@ snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet nsnips = floor(EEG.trials/snip_eps); if nsnips < 1 - error('Snippet length cannot exceed data length.') - end - diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); - if diff ~= 0 - warning(strcat(int2str(diff), ' seconds are thrown away.')); - end + warning('Snippet length cannot exceed data length. Using the whole data length.') + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); + else + diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); + if diff ~= 0 + warning(strcat(int2str(diff), ' seconds are thrown away.')); + end - source_roi_data_save = EEG.roi.source_roi_data; - for isnip = 1:nsnips - roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets - EEG.roi.source_roi_data = single(roi_snip); - EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); % compute connectivity over one snippet - for fc = 1:n_conn_metrics - fc_name = options{2}{fc}; - fc_matrix = EEG.roi.(fc_name); - conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure + if strcmpi(g.firstsnippet, 'on') + nsnips = 1; end - end - % compute mean over connectivity of each snippet - for fc = 1:n_conn_metrics - fc_name = options{2}{fc}; - [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); - - conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell - mat = cell2mat(conn_cell); - reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); - reshaped = squeeze(permute(reshaped, [2,1,3,4])); - if strcmpi(g.fcsave_format, 'all_snips') - EEG.roi.(fc_name) = reshaped; - else - if nsnips > 1 - mean_conn = squeeze(mean(reshaped, 1)); + source_roi_data_save = EEG.roi.source_roi_data; + for isnip = 1:nsnips + roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets + EEG.roi.source_roi_data = single(roi_snip); + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); % compute connectivity over one snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + fc_matrix = EEG.roi.(fc_name); + conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure + end + end + + % compute mean over connectivity of each snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); + + conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell + mat = cell2mat(conn_cell); + reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); + reshaped = squeeze(permute(reshaped, [2,1,3,4])); + if strcmpi(g.fcsave_format, 'all_snips') + EEG.roi.(fc_name) = reshaped; else - mean_conn = reshaped; + if nsnips > 1 + mean_conn = squeeze(mean(reshaped, 1)); + else + mean_conn = reshaped; + end + EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end - EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end end else diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 56d5ab6..5f526a3 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -75,10 +75,13 @@ end indNet = []; +len = 0; for iNet = 1:length(allnetworknames) if contains(networkname, allnetworknames{iNet}) - indNet = iNet; - break + if length( allnetworknames{iNet} ) > len + indNet = iNet; + len = length( allnetworknames{iNet} ); + end end end From d767b29b5f3db258f848867e97c743b34e2e8397 Mon Sep 17 00:00:00 2001 From: fpellegrini <38397852+fpellegrini@users.noreply.github.com> Date: Mon, 19 Jun 2023 08:15:14 +0200 Subject: [PATCH 02/25] Update README.md Updated publication --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d6c4e75..418a8f1 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ ROIconnect is a freely available open-source plugin to [EEGLAB](https://github.com/sccn/eeglab) for EEG data analysis. It allows you to perform functional connectivity analysis between regions of interests (ROIs) on source level. The results can be visualized in 2-D and 3-D. ROIs are defined based on popular fMRI atlases, and source localization can be performed through LCMV beamforming or eLORETA. Connectivity analysis can be performed between all pairs of brain regions using Granger Causality, Time-reversed Granger Causality, Multivariate Interaction Measure, Maximized Imaginary Coherency, and other methods. This plugin is compatible with Fieldtrip, Brainstorm and NFT head models. 📚 Check out the following paper to learn about recommended methods and pipelines for connectivity experiments: -> Pellegrini, F., Delorme, A., Nikulin, V. & Haufe, S., 2022. Identifying best practices for detecting inter-regional functional connectivity from EEG. bioRxiv 2022.10.05.510753. https://doi.org/10.1101/2022.10.05.510753 +> Pellegrini, F., Delorme, A., Nikulin, V., & Haufe, S. (2023). Identifying good practices for detecting inter-regional linear functional connectivity from EEG. NeuroImage, 120218. [doi: 10.1016/j.neuroimage.2023.120218](https://doi.org/10.1016/j.neuroimage.2023.120218) ⚠️ Disclaimer: This plugin implements the best-practice pipeline that we identified for our studied setting. We believe it can be used in such environments without hesitation. Additional code to reproduce our experiments entirely is provided in a [separate repository](https://github.com/fpellegrini/FCsim). In the medium term, we intend to extend this plugin for other use cases, which will be backed up by respective validation studies. From d31304568562a405c8a75a40f78965c3ba083279 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen <48970646+nguyen-td@users.noreply.github.com> Date: Thu, 22 Jun 2023 20:00:01 +0200 Subject: [PATCH 03/25] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d6c4e75..24f57e6 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ # What is ROIconnect? -ROIconnect is a freely available open-source plugin to [EEGLAB](https://github.com/sccn/eeglab) for EEG data analysis. It allows you to perform functional connectivity analysis between regions of interests (ROIs) on source level. The results can be visualized in 2-D and 3-D. ROIs are defined based on popular fMRI atlases, and source localization can be performed through LCMV beamforming or eLORETA. Connectivity analysis can be performed between all pairs of brain regions using Granger Causality, Time-reversed Granger Causality, Multivariate Interaction Measure, Maximized Imaginary Coherency, and other methods. This plugin is compatible with Fieldtrip, Brainstorm and NFT head models. +ROIconnect is a freely available open-source plugin to [EEGLAB](https://github.com/sccn/eeglab) for EEG data analysis. It allows you to perform functional connectivity analysis between regions of interest (ROIs) on source level. The results can be visualized in 2-D and 3-D. ROIs are defined based on popular fMRI atlases, and source localization can be performed through LCMV beamforming or eLORETA. Connectivity analysis can be performed between all pairs of brain regions using Granger Causality, Time-reversed Granger Causality, Multivariate Interaction Measure, Maximized Imaginary Coherency, and other methods. This plugin is compatible with Fieldtrip, Brainstorm and NFT head models. 📚 Check out the following paper to learn about recommended methods and pipelines for connectivity experiments: > Pellegrini, F., Delorme, A., Nikulin, V. & Haufe, S., 2022. Identifying best practices for detecting inter-regional functional connectivity from EEG. bioRxiv 2022.10.05.510753. https://doi.org/10.1101/2022.10.05.510753 From 1e3f5dec548496aad294cdeca81b34ee23add950 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:24:46 +0200 Subject: [PATCH 04/25] added MIM statistics incl. plotting --- libs/pellegrini/shuffle_MIM.m | 2 + pop_roi_statsplot.m | 81 +++++++++++++++++++++++++++++++++++ roi_connstats.m | 4 +- 3 files changed, 84 insertions(+), 3 deletions(-) create mode 100644 pop_roi_statsplot.m diff --git a/libs/pellegrini/shuffle_MIM.m b/libs/pellegrini/shuffle_MIM.m index 132b6ce..6a2463b 100644 --- a/libs/pellegrini/shuffle_MIM.m +++ b/libs/pellegrini/shuffle_MIM.m @@ -8,6 +8,7 @@ ninds = length(inds); warning('One iteration takes about 90 seconds.') + fprintf('Generating null distribution using %d shuffles...\n', nshuf) fprintf('Progress of %d:', nshuf); for ishuf = 1:nshuf %one iteration takes ~90 sec on my local laptop if mod(ishuf, 10) == 0 @@ -52,6 +53,7 @@ MIM_s(:, :, :, ishuf) = get_connect_mat(conn.MIM, nroi, +1); % [MIM_s(:,:,ishuf), ~, ~, ~, ~, ~] = fp_unwrap_conn(conn,nroi,filt1,PCA_inds); end + fprintf('\n'); end function measure = get_connect_mat( measureOri, nROI, signVal) diff --git a/pop_roi_statsplot.m b/pop_roi_statsplot.m new file mode 100644 index 0000000..649b1da --- /dev/null +++ b/pop_roi_statsplot.m @@ -0,0 +1,81 @@ +% pop_roi_statsplot() - Generate p-values from FC null distributions and plots them. +% +% Inputs: +% EEG - EEGLAB dataset with ROI activity computed. +% +% Optional inputs: +% 'measure' - [cell] Cell of strings corresponding to methods. +% 'CS' : Cross spectrum +% 'aCOH' : Coherence +% 'cCOH' : (complex-valued) coherency +% 'iCOH' : absolute value of the imaginary part of coherency +% 'wPLI' : Weighted Phase Lag Index +% 'PDC' : Partial directed coherence +% 'TRPDC' : Time-reversed partial directed coherence +% 'DTF' : Directed transfer entropy +% 'TRDTF' : Time-reversed directed transfer entropy +% 'MIM' : Multivariate Interaction Measure for each ROI +% 'MIC' : Maximized Imaginary Coherency for each ROI +% 'freqrange' - [min max] frequency range in Hz. Default is to plot broadband power. + +function EEG = pop_roi_statsplot(EEG, varargin) + + if nargin < 2 + help roi_connstats; + return + end + + if ~isfield(EEG, 'roi') || ~isfield(EEG.roi, 'source_roi_data') + error('Cannot find ROI data - compute ROI data first'); + end + + % decode input parameters + % ----------------------- + g = finputcheck(varargin, { + 'measure' 'string' { } ''; + 'freqrange' 'real' { } []}, 'pop_roi_statsplot'); + if ischar(g), error(g); end + S = EEG.roi; + + if isempty(g.measure) + error('You must define a measure to plot'); + end + + % extract frequency indices + if ~isempty(g.freqrange) + try + frq_inds = find(S.freqs >= g.freqrange(1) & S.freqs < g.freqrange(2)); + title = sprintf('%1.1f-%1.1f Hz frequency band', g.freqrange(1), g.freqrange(2)); + catch + frq_inds = find(S.freqs == g.freqrange(1)); % if a single frequency is passed + title = sprintf('%1.1f Hz', g.freqrange(1)); + end + else + frq_inds = 1:length(S.freqs); + title = 'broadband'; + end + + % select frequency or frequency band + if length(frq_inds) > 1 + matrix = squeeze(mean(S.(g.measure)(frq_inds, :, :, :))); + else + matrix = squeeze(S.(g.measure)(frq_inds, :, :, :)); + end + + % generate p-values by comparing the true FC (first shuffle) to null distribution + netFC = squeeze(mean(matrix, 2)); + FC_pn = sum(netFC(:, 1) < netFC(:, 2:end), 2)./(size(matrix, 3) - 1); + + % use FDR-correction for multiple comparison's correction + [p_fdr, ~] = fdr(FC_pn, 0.05); + FC_pn(FC_pn > p_fdr) = 1; + + % plot + load cm17; + load cortex; + FC_pn(FC_pn==0) = 1 / (size(netMIM, 2) - 1); % 1 / nshuf + data = -log10(FC_pn); + allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3); + h = textsc(title, 'title'); + set(h, 'fontsize', 20); +end \ No newline at end of file diff --git a/roi_connstats.m b/roi_connstats.m index 436fb15..b98ac47 100644 --- a/roi_connstats.m +++ b/roi_connstats.m @@ -10,8 +10,6 @@ % 'aCOH' : Coherence % 'cCOH' : (complex-valued) coherency % 'iCOH' : absolute value of the imaginary part of coherency -% 'GC' : Granger Causality -% 'TRGC' : Time-reversed Granger Causality % 'wPLI' : Weighted Phase Lag Index % 'PDC' : Partial directed coherence % 'TRPDC' : Time-reversed partial directed coherence @@ -19,7 +17,7 @@ % 'TRDTF' : Time-reversed directed transfer entropy % 'MIM' : Multivariate Interaction Measure for each ROI % 'MIC' : Maximized Imaginary Coherency for each ROI -% nshuf' - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 101. +% nshuf - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 1001. % roi_selection - [cell array of integers] Cell array of ROI indices {1, 2, 3, ...} indicating for which regions (ROIs) connectivity should be computed. % Default is all (set to EEG.roi.nROI). From ea2408d5cb68c91cbe920880005948bce40ebe3e Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:46:17 +0200 Subject: [PATCH 05/25] replaced matlab version check, adressing #651 issue in sccn/eeglab --- pop_roi_connectplot.m | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 3499481..0b47ac8 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -763,10 +763,15 @@ function roi_plotcoloredlobes( EEG, matrix, titleStr, measure, hemisphere, group imagesc(matrix); colormap(cmap); end cb = colorbar; - tf = isMATLABReleaseOlderThan("R2022a"); - if tf +% tf = isMATLABReleaseOlderThan("R2022a"); +% if tf +% caxis([clim_min clim_max]) +% else +% clim([clim_min clim_max]) +% end + try caxis([clim_min clim_max]) - else + catch clim([clim_min clim_max]) end if isDKatlas == 1 From 257ed5f83be4ade0f41d722f547dd8ed1245bc06 Mon Sep 17 00:00:00 2001 From: arno Date: Mon, 17 Jul 2023 09:56:16 -0700 Subject: [PATCH 06/25] update warning message --- roi_lookupregions.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 5f526a3..d46ac09 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -86,7 +86,7 @@ end if isempty(indNet) - fprintf('Network %s not found\n', networkname); + fprintf('The brain region "%s" was not found (it is ok if it is not a brain region)\n', networkname); elseif ~isstruct(networkdefs) regions = roiTable.(allnetworknames{iNet}); else From bb95f1efdfdcc59ad8e056e43a35169677de326a Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:31 -1000 Subject: [PATCH 07/25] implement fix for LORETA coordinates --- roi_sourceplot.m | 166 ++++++++++++++++++++++++++++++----------------- 1 file changed, 108 insertions(+), 58 deletions(-) diff --git a/roi_sourceplot.m b/roi_sourceplot.m index e33c23c..0a7d4ae 100644 --- a/roi_sourceplot.m +++ b/roi_sourceplot.m @@ -15,6 +15,10 @@ % 'freqselect' - [real] frequency of interest or frequency range of interest. % Defaut is all frequencies. % +% Example: +% % run ROI activity using ROI connect to compute activity in ROI (and voxels) +% roi_sourceplot(EEG.roi.freqs, EEG.roi.source_voxel_power, EEG.roi.cortex ); +% % Author: Arnaud Delorme % Copyright (C) Arnaud Delorme, arnodelorme@gmail.com @@ -51,13 +55,17 @@ g = finputcheck(varargin, { ... 'freqrange' {'cell' 'real'} { [] [] } {}; 'saveasfile' 'string' { } ''; + 'mask' 'string' { } ''; 'noplot' 'string' { 'on' 'off' } 'off'; 'precomputed' 'struct' { } struct([]); 'limits' '' [] []; - 'slice' 'integer' [] [5 8 10 15 18]}, 'roi_sourceplot'); + 'slice' 'integer' [] []}, 'roi_sourceplot'); if ischar(g) error(g); end +if isempty(g.freqrange) + g.freqrange = [freqs(1) freqs(end)]; +end if ~iscell(g.freqrange) g.freqrange = { g.freqrange }; end @@ -68,23 +76,7 @@ alldata = []; for iFreq = 1:length(g.freqrange) - g.freqselect = g.freqrange{iFreq}; - if isempty(g.freqselect) - indFreq = 1:length(freqs); - elseif length(g.freqselect) == 1 - [~,indFreq] = min(abs(freqs-g.freqselect)); - indFreq(2) = indFreq; - elseif length(g.freqselect) == 2 - [~,indFreq1] = min(abs(freqs-g.freqselect(1))); - [~,indFreq2] = min(abs(freqs-g.freqselect(2))); - indFreq = indFreq1:indFreq2; - else - error('Frequency selection must be an array of 1 or 2 elements'); - end - if isempty(indFreq) - error('No frequency found'); - end - + % transform to volume if ischar(sourcemodel) sourceProjtmp = load('-mat', sourcemodel); @@ -94,62 +86,96 @@ if ~isfield(sourceProjtmp, 'Vertices') sourceProjtmp.Vertices = sourceProjtmp.pos; end - xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; - yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; - zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; - downscale = 5; - volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); - clear sourcemodelout; - sourcemodelout.dim = size(volMat); - [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); - sourcemodelout.pos = [r,c,v]; - sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); - sourcemodelout.unit = 'mm'; - %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); - sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid - allInds = zeros(size(sourceProjtmp.Vertices)); - allIndVolume = zeros(length(sourceProjtmp.Vertices),1); - for iVert = 1:length(sourceProjtmp.Vertices) - xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; - yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; - zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; - ind = sub2ind(size(volMat), xVal, yVal, zVal); - volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); - allIndVolume(iVert) = ind; - allInds(iVert,:) = [xVal yVal zVal]; - sourcemodelout.inside(ind) = true; + + stride = unique(abs(diff(sourceProjtmp.Vertices))); + if iscell(sourceact) + % select precomputed activity matrix + volMat = loreta2volume(sourceact{iFreq}, sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); + else + % select frequency range in activity matrix + g.freqselect = g.freqrange{iFreq}; + if isempty(g.freqselect) + indFreq = 1:length(freqs); + elseif length(g.freqselect) == 1 + [~,indFreq] = min(abs(freqs-g.freqselect)); + indFreq(2) = indFreq; + elseif length(g.freqselect) == 2 + [~,indFreq1] = min(abs(freqs-g.freqselect(1))); + [~,indFreq2] = min(abs(freqs-g.freqselect(2))); + indFreq = indFreq1:indFreq2; + else + error('Frequency selection must be an array of 1 or 2 elements'); + end + if isempty(indFreq) + error('No frequency found'); + end + volMat = loreta2volume(mean(sourceact(:,indFreq), 2), sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); end + if stride(2) ~= floor(stride(2)) + error('Non-integer stride') + end + +% xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; +% yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; +% zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; +% downscale = 10; +% volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); +% clear sourcemodelout; +% sourcemodelout.dim = size(volMat); +% [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); +% sourcemodelout.pos = [r,c,v]; +% sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); +% sourcemodelout.unit = 'mm'; +% %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); +% sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid +% allInds = zeros(size(sourceProjtmp.Vertices)); +% allIndVolume = zeros(length(sourceProjtmp.Vertices),1); +% for iVert = 1:length(sourceProjtmp.Vertices) +% xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; +% yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; +% zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; +% ind = sub2ind(size(volMat), xVal, yVal, zVal); +% volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); +% allIndVolume(iVert) = ind; +% allInds(iVert,:) = [xVal yVal zVal]; +% sourcemodelout.inside(ind) = true; +% end % put precomputed data in VolMat - for iSlice = 1:length(g.slice) - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + if ~isempty(g.slice) + tmpSlice = g.slice; + else + tmpSlice = ceil(linspace(3, size(volMat,3)-2, 5)); + end + for iSlice = 1:length(tmpSlice) + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); if isfield(g.precomputed, fieldVal) res = g.precomputed.(fieldVal); res(isnan(res)) = 0; res(:,:) = res(end:-1:1,:); res = res'; - volMat(:,:,g.slice(iSlice)) = res; + volMat(:,:,tmpSlice(iSlice)) = res; end end % plot - res = squeeze(volMat(:,:,g.slice)); + res = squeeze(volMat(:,:,tmpSlice)); if isempty(g.limits) - mi = min(res(:)); - mx = max(res(:)); + mi = min(res(:),[],'omitnan'); + mx = max(res(:),[],'omitnan'); fprintf('Loreta limits from data: %1.2f to %1.2f\n', mi, mx); else mi = 0; mx = g.limits(iFreq); end - cmap = colormap('jet'); - for iSlice = 1:length(g.slice) - res = squeeze(volMat(:,:,g.slice(iSlice))); + cmap = colormap('turbo'); + for iSlice = 1:length(tmpSlice) + res = squeeze(volMat(:,:,tmpSlice(iSlice))); res = res'; res(:,:) = res(end:-1:1,:); % save and retreive data - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); alldata.(fieldVal) = res; if strcmpi(g.noplot, 'off') @@ -157,15 +183,19 @@ for iPix1 = 1:size(res,1) for iPix2 = 1:size(res,2) if res(iPix1,iPix2) ~= 0 - ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; - ind = max(1, ind); - ind = min(size(cmap,1), ind); - resrgb(iPix1,iPix2,:) = cmap(ind,:); + if isnan(res(iPix1,iPix2)) + resrgb(iPix1,iPix2,:) = [0.9 0.9 0.9]; + else + ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; + ind = max(1, ind); + ind = min(size(cmap,1), ind); + resrgb(iPix1,iPix2,:) = cmap(ind,:); + end end end end - subplot(length(g.slice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); + subplot(length(tmpSlice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); imagesc(resrgb); axis equal; axis off; if iSlice == 1 if length(g.freqrange{iFreq}) == 2 @@ -216,4 +246,24 @@ cfg.location = [26 8 10]; cfg.latency = 10; cfg.slicepos = [8 9]; -ft_sourceplot(cfg, sourceProj); \ No newline at end of file +ft_sourceplot(cfg, sourceProj); + +% convert coordinates to volume +% ----------------------------- +function vol = loreta2volume(act, pos, stride) + +minx = min(pos(:,1)); maxx = max(pos(:,1)); +miny = min(pos(:,2)); maxy = max(pos(:,2)); +minz = min(pos(:,3)); maxz = max(pos(:,3)); + +vol = zeros((maxx-minx)/stride(1), (maxy-miny)/stride(2), (maxz-minz)/stride(3)); + +for iPos = 1:size(pos,1) + vol( (pos(iPos,1)-minx)/stride(1)+1, ... + (pos(iPos,2)-miny)/stride(2)+1, ... + (pos(iPos,3)-minz)/stride(3)+1) = act(iPos); + % in case using the sourcemodel with shrunk coordinates + % vol( floor((pos(iPos,1)-minx)/stride(1)+1), ... + % floor((pos(iPos,2)-miny)/stride(2)+1), ... + % floor((pos(iPos,3)-minz)/stride(3)+1)) = act(iPos); +end From b9c8f2fc3afea7720a09129bd197d3a6dae2ee58 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:48 -1000 Subject: [PATCH 08/25] now output voxel power --- roi_activity.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/roi_activity.m b/roi_activity.m index 05b4840..2c65838 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -332,7 +332,8 @@ tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels - + source_voxel_power = tmpWelch; + % fooof settings if strcmpi(g.fooof, 'on') f_range = g.fooof_frange; % freq range where 1/f should be fitted @@ -377,6 +378,7 @@ if strcmpi(g.exportvoxact, 'on') EEG.roi.source_voxel_data = source_voxel_data; % large (takes lots of RAM) end +EEG.roi.source_voxel_power = single(source_voxel_power); EEG.roi.source_roi_data = single(source_roi_data); EEG.roi.source_roi_power = source_roi_power; % used for plotting EEG.roi.source_roi_power_norm = source_roi_power_norm; % used for cross-sprectum From 8b5daa0653f809952baeffefb3c2a5e5f5523d1e Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:31 -1000 Subject: [PATCH 09/25] implement fix for LORETA coordinates --- roi_sourceplot.m | 166 ++++++++++++++++++++++++++++++----------------- 1 file changed, 108 insertions(+), 58 deletions(-) diff --git a/roi_sourceplot.m b/roi_sourceplot.m index e33c23c..0a7d4ae 100644 --- a/roi_sourceplot.m +++ b/roi_sourceplot.m @@ -15,6 +15,10 @@ % 'freqselect' - [real] frequency of interest or frequency range of interest. % Defaut is all frequencies. % +% Example: +% % run ROI activity using ROI connect to compute activity in ROI (and voxels) +% roi_sourceplot(EEG.roi.freqs, EEG.roi.source_voxel_power, EEG.roi.cortex ); +% % Author: Arnaud Delorme % Copyright (C) Arnaud Delorme, arnodelorme@gmail.com @@ -51,13 +55,17 @@ g = finputcheck(varargin, { ... 'freqrange' {'cell' 'real'} { [] [] } {}; 'saveasfile' 'string' { } ''; + 'mask' 'string' { } ''; 'noplot' 'string' { 'on' 'off' } 'off'; 'precomputed' 'struct' { } struct([]); 'limits' '' [] []; - 'slice' 'integer' [] [5 8 10 15 18]}, 'roi_sourceplot'); + 'slice' 'integer' [] []}, 'roi_sourceplot'); if ischar(g) error(g); end +if isempty(g.freqrange) + g.freqrange = [freqs(1) freqs(end)]; +end if ~iscell(g.freqrange) g.freqrange = { g.freqrange }; end @@ -68,23 +76,7 @@ alldata = []; for iFreq = 1:length(g.freqrange) - g.freqselect = g.freqrange{iFreq}; - if isempty(g.freqselect) - indFreq = 1:length(freqs); - elseif length(g.freqselect) == 1 - [~,indFreq] = min(abs(freqs-g.freqselect)); - indFreq(2) = indFreq; - elseif length(g.freqselect) == 2 - [~,indFreq1] = min(abs(freqs-g.freqselect(1))); - [~,indFreq2] = min(abs(freqs-g.freqselect(2))); - indFreq = indFreq1:indFreq2; - else - error('Frequency selection must be an array of 1 or 2 elements'); - end - if isempty(indFreq) - error('No frequency found'); - end - + % transform to volume if ischar(sourcemodel) sourceProjtmp = load('-mat', sourcemodel); @@ -94,62 +86,96 @@ if ~isfield(sourceProjtmp, 'Vertices') sourceProjtmp.Vertices = sourceProjtmp.pos; end - xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; - yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; - zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; - downscale = 5; - volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); - clear sourcemodelout; - sourcemodelout.dim = size(volMat); - [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); - sourcemodelout.pos = [r,c,v]; - sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); - sourcemodelout.unit = 'mm'; - %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); - sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid - allInds = zeros(size(sourceProjtmp.Vertices)); - allIndVolume = zeros(length(sourceProjtmp.Vertices),1); - for iVert = 1:length(sourceProjtmp.Vertices) - xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; - yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; - zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; - ind = sub2ind(size(volMat), xVal, yVal, zVal); - volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); - allIndVolume(iVert) = ind; - allInds(iVert,:) = [xVal yVal zVal]; - sourcemodelout.inside(ind) = true; + + stride = unique(abs(diff(sourceProjtmp.Vertices))); + if iscell(sourceact) + % select precomputed activity matrix + volMat = loreta2volume(sourceact{iFreq}, sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); + else + % select frequency range in activity matrix + g.freqselect = g.freqrange{iFreq}; + if isempty(g.freqselect) + indFreq = 1:length(freqs); + elseif length(g.freqselect) == 1 + [~,indFreq] = min(abs(freqs-g.freqselect)); + indFreq(2) = indFreq; + elseif length(g.freqselect) == 2 + [~,indFreq1] = min(abs(freqs-g.freqselect(1))); + [~,indFreq2] = min(abs(freqs-g.freqselect(2))); + indFreq = indFreq1:indFreq2; + else + error('Frequency selection must be an array of 1 or 2 elements'); + end + if isempty(indFreq) + error('No frequency found'); + end + volMat = loreta2volume(mean(sourceact(:,indFreq), 2), sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); end + if stride(2) ~= floor(stride(2)) + error('Non-integer stride') + end + +% xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; +% yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; +% zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; +% downscale = 10; +% volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); +% clear sourcemodelout; +% sourcemodelout.dim = size(volMat); +% [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); +% sourcemodelout.pos = [r,c,v]; +% sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); +% sourcemodelout.unit = 'mm'; +% %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); +% sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid +% allInds = zeros(size(sourceProjtmp.Vertices)); +% allIndVolume = zeros(length(sourceProjtmp.Vertices),1); +% for iVert = 1:length(sourceProjtmp.Vertices) +% xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; +% yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; +% zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; +% ind = sub2ind(size(volMat), xVal, yVal, zVal); +% volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); +% allIndVolume(iVert) = ind; +% allInds(iVert,:) = [xVal yVal zVal]; +% sourcemodelout.inside(ind) = true; +% end % put precomputed data in VolMat - for iSlice = 1:length(g.slice) - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + if ~isempty(g.slice) + tmpSlice = g.slice; + else + tmpSlice = ceil(linspace(3, size(volMat,3)-2, 5)); + end + for iSlice = 1:length(tmpSlice) + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); if isfield(g.precomputed, fieldVal) res = g.precomputed.(fieldVal); res(isnan(res)) = 0; res(:,:) = res(end:-1:1,:); res = res'; - volMat(:,:,g.slice(iSlice)) = res; + volMat(:,:,tmpSlice(iSlice)) = res; end end % plot - res = squeeze(volMat(:,:,g.slice)); + res = squeeze(volMat(:,:,tmpSlice)); if isempty(g.limits) - mi = min(res(:)); - mx = max(res(:)); + mi = min(res(:),[],'omitnan'); + mx = max(res(:),[],'omitnan'); fprintf('Loreta limits from data: %1.2f to %1.2f\n', mi, mx); else mi = 0; mx = g.limits(iFreq); end - cmap = colormap('jet'); - for iSlice = 1:length(g.slice) - res = squeeze(volMat(:,:,g.slice(iSlice))); + cmap = colormap('turbo'); + for iSlice = 1:length(tmpSlice) + res = squeeze(volMat(:,:,tmpSlice(iSlice))); res = res'; res(:,:) = res(end:-1:1,:); % save and retreive data - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); alldata.(fieldVal) = res; if strcmpi(g.noplot, 'off') @@ -157,15 +183,19 @@ for iPix1 = 1:size(res,1) for iPix2 = 1:size(res,2) if res(iPix1,iPix2) ~= 0 - ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; - ind = max(1, ind); - ind = min(size(cmap,1), ind); - resrgb(iPix1,iPix2,:) = cmap(ind,:); + if isnan(res(iPix1,iPix2)) + resrgb(iPix1,iPix2,:) = [0.9 0.9 0.9]; + else + ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; + ind = max(1, ind); + ind = min(size(cmap,1), ind); + resrgb(iPix1,iPix2,:) = cmap(ind,:); + end end end end - subplot(length(g.slice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); + subplot(length(tmpSlice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); imagesc(resrgb); axis equal; axis off; if iSlice == 1 if length(g.freqrange{iFreq}) == 2 @@ -216,4 +246,24 @@ cfg.location = [26 8 10]; cfg.latency = 10; cfg.slicepos = [8 9]; -ft_sourceplot(cfg, sourceProj); \ No newline at end of file +ft_sourceplot(cfg, sourceProj); + +% convert coordinates to volume +% ----------------------------- +function vol = loreta2volume(act, pos, stride) + +minx = min(pos(:,1)); maxx = max(pos(:,1)); +miny = min(pos(:,2)); maxy = max(pos(:,2)); +minz = min(pos(:,3)); maxz = max(pos(:,3)); + +vol = zeros((maxx-minx)/stride(1), (maxy-miny)/stride(2), (maxz-minz)/stride(3)); + +for iPos = 1:size(pos,1) + vol( (pos(iPos,1)-minx)/stride(1)+1, ... + (pos(iPos,2)-miny)/stride(2)+1, ... + (pos(iPos,3)-minz)/stride(3)+1) = act(iPos); + % in case using the sourcemodel with shrunk coordinates + % vol( floor((pos(iPos,1)-minx)/stride(1)+1), ... + % floor((pos(iPos,2)-miny)/stride(2)+1), ... + % floor((pos(iPos,3)-minz)/stride(3)+1)) = act(iPos); +end From 12d215f5357e30ffc02a9f0f15f63a5ef87d876b Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:48 -1000 Subject: [PATCH 10/25] now output voxel power --- roi_activity.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/roi_activity.m b/roi_activity.m index c43b0e2..01febf6 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -332,7 +332,8 @@ tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels - + source_voxel_power = tmpWelch; + % fooof settings if strcmpi(g.fooof, 'on') f_range = g.fooof_frange; % freq range where 1/f should be fitted @@ -377,6 +378,7 @@ if strcmpi(g.exportvoxact, 'on') EEG.roi.source_voxel_data = source_voxel_data; % large (takes lots of RAM) end +EEG.roi.source_voxel_power = single(source_voxel_power); EEG.roi.source_roi_data = single(source_roi_data); EEG.roi.source_roi_power = source_roi_power; % used for plotting EEG.roi.source_roi_power_norm = source_roi_power_norm; % used for cross-sprectum From f7221b690bc6a44a887c6a07b0f6bb005e86521c Mon Sep 17 00:00:00 2001 From: arno Date: Tue, 23 May 2023 16:45:31 -0700 Subject: [PATCH 11/25] snippet calculation --- plotconnectivity.m | 5 +++-- pop_roi_connect.m | 54 ++++++++++++++++++++++++++++++--------------- roi_lookupregions.m | 7 ++++-- 3 files changed, 44 insertions(+), 22 deletions(-) diff --git a/plotconnectivity.m b/plotconnectivity.m index 1533e72..d7b14ed 100644 --- a/plotconnectivity.m +++ b/plotconnectivity.m @@ -80,7 +80,7 @@ if isempty(g.labels) if strcmpi(g.brainimg, 'on') disp('Cannot plot on brain with area labels') - g.brainimg = 'off'; + g.brainimg = 'off';g end end @@ -193,6 +193,7 @@ % rename labels % ------------- +if ~strcmpi(g.brainimg, 'off') % CTC added .... to branch around code if not plotting if isempty(g.labels) for iPnt = 1:length(anglesInit) g.labels{iPnt} = sprintf(' Area %d', iPnt); @@ -205,7 +206,7 @@ g.labels{iPnt} = [ ' ' g.labels{iPnt} ]; end end - +end % CTC added .... to branch around above code if not plotting warning off; if isempty(limits) limits(2) = max(array(:)); diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 106a34e..b2fa1dc 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -31,6 +31,7 @@ % 'MIC' : Maximized Imaginary Coherency for each ROI % 'PAC' : Phase-amplitude coupling between ROIs % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast comptuation). Default is 'off'. % 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets % (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' @@ -164,7 +165,8 @@ 'naccu' 'integer' { } 0; 'methods' 'cell' { } { }; 'snippet' 'string' { 'on', 'off' } 'off'; - 'nepochs' 'real' {} [ ]; + 'firstsnippet' 'string' { 'on', 'off' } 'off'; + 'nepochs' 'real' {} []; 'snip_length' 'integer' { } 60; 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; 'freqresolution' 'integer' { } 0; @@ -245,26 +247,42 @@ fc_matrix = EEG.roi.(fc_name); conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure end - end - - % compute mean over connectivity of each snippet - for fc = 1:n_conn_metrics - fc_name = options{2}{fc}; - [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); - conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell - mat = cell2mat(conn_cell); - reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); - reshaped = squeeze(permute(reshaped, [2,1,3,4])); - if strcmpi(g.fcsave_format, 'all_snips') - EEG.roi.(fc_name) = reshaped; - else - if nsnips > 1 - mean_conn = squeeze(mean(reshaped, 1)); + if strcmpi(g.firstsnippet, 'on') + nsnips = 1; + end + + source_roi_data_save = EEG.roi.source_roi_data; + for isnip = 1:nsnips + roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets + EEG.roi.source_roi_data = single(roi_snip); + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); % compute connectivity over one snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + fc_matrix = EEG.roi.(fc_name); + conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure + end + end + + % compute mean over connectivity of each snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); + + conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell + mat = cell2mat(conn_cell); + reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); + reshaped = squeeze(permute(reshaped, [2,1,3,4])); + if strcmpi(g.fcsave_format, 'all_snips') + EEG.roi.(fc_name) = reshaped; else - mean_conn = reshaped; + if nsnips > 1 + mean_conn = squeeze(mean(reshaped, 1)); + else + mean_conn = reshaped; + end + EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end - EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end end else diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 56d5ab6..5f526a3 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -75,10 +75,13 @@ end indNet = []; +len = 0; for iNet = 1:length(allnetworknames) if contains(networkname, allnetworknames{iNet}) - indNet = iNet; - break + if length( allnetworknames{iNet} ) > len + indNet = iNet; + len = length( allnetworknames{iNet} ); + end end end From b84bc91582a0b5807bd3e958cfa2ca80c0c6b444 Mon Sep 17 00:00:00 2001 From: arno Date: Mon, 17 Jul 2023 09:56:16 -0700 Subject: [PATCH 12/25] update warning message --- roi_lookupregions.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 5f526a3..d46ac09 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -86,7 +86,7 @@ end if isempty(indNet) - fprintf('Network %s not found\n', networkname); + fprintf('The brain region "%s" was not found (it is ok if it is not a brain region)\n', networkname); elseif ~isstruct(networkdefs) regions = roiTable.(allnetworknames{iNet}); else From 2e733ce2e4219fffd2f6eb41934be9681e44a20e Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:31 -1000 Subject: [PATCH 13/25] implement fix for LORETA coordinates --- roi_sourceplot.m | 166 ++++++++++++++++++++++++++++++----------------- 1 file changed, 108 insertions(+), 58 deletions(-) diff --git a/roi_sourceplot.m b/roi_sourceplot.m index e33c23c..0a7d4ae 100644 --- a/roi_sourceplot.m +++ b/roi_sourceplot.m @@ -15,6 +15,10 @@ % 'freqselect' - [real] frequency of interest or frequency range of interest. % Defaut is all frequencies. % +% Example: +% % run ROI activity using ROI connect to compute activity in ROI (and voxels) +% roi_sourceplot(EEG.roi.freqs, EEG.roi.source_voxel_power, EEG.roi.cortex ); +% % Author: Arnaud Delorme % Copyright (C) Arnaud Delorme, arnodelorme@gmail.com @@ -51,13 +55,17 @@ g = finputcheck(varargin, { ... 'freqrange' {'cell' 'real'} { [] [] } {}; 'saveasfile' 'string' { } ''; + 'mask' 'string' { } ''; 'noplot' 'string' { 'on' 'off' } 'off'; 'precomputed' 'struct' { } struct([]); 'limits' '' [] []; - 'slice' 'integer' [] [5 8 10 15 18]}, 'roi_sourceplot'); + 'slice' 'integer' [] []}, 'roi_sourceplot'); if ischar(g) error(g); end +if isempty(g.freqrange) + g.freqrange = [freqs(1) freqs(end)]; +end if ~iscell(g.freqrange) g.freqrange = { g.freqrange }; end @@ -68,23 +76,7 @@ alldata = []; for iFreq = 1:length(g.freqrange) - g.freqselect = g.freqrange{iFreq}; - if isempty(g.freqselect) - indFreq = 1:length(freqs); - elseif length(g.freqselect) == 1 - [~,indFreq] = min(abs(freqs-g.freqselect)); - indFreq(2) = indFreq; - elseif length(g.freqselect) == 2 - [~,indFreq1] = min(abs(freqs-g.freqselect(1))); - [~,indFreq2] = min(abs(freqs-g.freqselect(2))); - indFreq = indFreq1:indFreq2; - else - error('Frequency selection must be an array of 1 or 2 elements'); - end - if isempty(indFreq) - error('No frequency found'); - end - + % transform to volume if ischar(sourcemodel) sourceProjtmp = load('-mat', sourcemodel); @@ -94,62 +86,96 @@ if ~isfield(sourceProjtmp, 'Vertices') sourceProjtmp.Vertices = sourceProjtmp.pos; end - xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; - yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; - zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; - downscale = 5; - volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); - clear sourcemodelout; - sourcemodelout.dim = size(volMat); - [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); - sourcemodelout.pos = [r,c,v]; - sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); - sourcemodelout.unit = 'mm'; - %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); - sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid - allInds = zeros(size(sourceProjtmp.Vertices)); - allIndVolume = zeros(length(sourceProjtmp.Vertices),1); - for iVert = 1:length(sourceProjtmp.Vertices) - xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; - yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; - zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; - ind = sub2ind(size(volMat), xVal, yVal, zVal); - volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); - allIndVolume(iVert) = ind; - allInds(iVert,:) = [xVal yVal zVal]; - sourcemodelout.inside(ind) = true; + + stride = unique(abs(diff(sourceProjtmp.Vertices))); + if iscell(sourceact) + % select precomputed activity matrix + volMat = loreta2volume(sourceact{iFreq}, sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); + else + % select frequency range in activity matrix + g.freqselect = g.freqrange{iFreq}; + if isempty(g.freqselect) + indFreq = 1:length(freqs); + elseif length(g.freqselect) == 1 + [~,indFreq] = min(abs(freqs-g.freqselect)); + indFreq(2) = indFreq; + elseif length(g.freqselect) == 2 + [~,indFreq1] = min(abs(freqs-g.freqselect(1))); + [~,indFreq2] = min(abs(freqs-g.freqselect(2))); + indFreq = indFreq1:indFreq2; + else + error('Frequency selection must be an array of 1 or 2 elements'); + end + if isempty(indFreq) + error('No frequency found'); + end + volMat = loreta2volume(mean(sourceact(:,indFreq), 2), sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]); end + if stride(2) ~= floor(stride(2)) + error('Non-integer stride') + end + +% xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ]; +% yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ]; +% zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ]; +% downscale = 10; +% volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1)); +% clear sourcemodelout; +% sourcemodelout.dim = size(volMat); +% [r,c,v] = ind2sub(size(volMat),find(volMat == 0)); +% sourcemodelout.pos = [r,c,v]; +% sourcemodelout.inside = zeros(size(sourcemodelout.pos,1),1); +% sourcemodelout.unit = 'mm'; +% %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]); +% sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid +% allInds = zeros(size(sourceProjtmp.Vertices)); +% allIndVolume = zeros(length(sourceProjtmp.Vertices),1); +% for iVert = 1:length(sourceProjtmp.Vertices) +% xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1; +% yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1; +% zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1; +% ind = sub2ind(size(volMat), xVal, yVal, zVal); +% volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2); +% allIndVolume(iVert) = ind; +% allInds(iVert,:) = [xVal yVal zVal]; +% sourcemodelout.inside(ind) = true; +% end % put precomputed data in VolMat - for iSlice = 1:length(g.slice) - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + if ~isempty(g.slice) + tmpSlice = g.slice; + else + tmpSlice = ceil(linspace(3, size(volMat,3)-2, 5)); + end + for iSlice = 1:length(tmpSlice) + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); if isfield(g.precomputed, fieldVal) res = g.precomputed.(fieldVal); res(isnan(res)) = 0; res(:,:) = res(end:-1:1,:); res = res'; - volMat(:,:,g.slice(iSlice)) = res; + volMat(:,:,tmpSlice(iSlice)) = res; end end % plot - res = squeeze(volMat(:,:,g.slice)); + res = squeeze(volMat(:,:,tmpSlice)); if isempty(g.limits) - mi = min(res(:)); - mx = max(res(:)); + mi = min(res(:),[],'omitnan'); + mx = max(res(:),[],'omitnan'); fprintf('Loreta limits from data: %1.2f to %1.2f\n', mi, mx); else mi = 0; mx = g.limits(iFreq); end - cmap = colormap('jet'); - for iSlice = 1:length(g.slice) - res = squeeze(volMat(:,:,g.slice(iSlice))); + cmap = colormap('turbo'); + for iSlice = 1:length(tmpSlice) + res = squeeze(volMat(:,:,tmpSlice(iSlice))); res = res'; res(:,:) = res(end:-1:1,:); % save and retreive data - fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice)); + fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice)); alldata.(fieldVal) = res; if strcmpi(g.noplot, 'off') @@ -157,15 +183,19 @@ for iPix1 = 1:size(res,1) for iPix2 = 1:size(res,2) if res(iPix1,iPix2) ~= 0 - ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; - ind = max(1, ind); - ind = min(size(cmap,1), ind); - resrgb(iPix1,iPix2,:) = cmap(ind,:); + if isnan(res(iPix1,iPix2)) + resrgb(iPix1,iPix2,:) = [0.9 0.9 0.9]; + else + ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1; + ind = max(1, ind); + ind = min(size(cmap,1), ind); + resrgb(iPix1,iPix2,:) = cmap(ind,:); + end end end end - subplot(length(g.slice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); + subplot(length(tmpSlice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1)); imagesc(resrgb); axis equal; axis off; if iSlice == 1 if length(g.freqrange{iFreq}) == 2 @@ -216,4 +246,24 @@ cfg.location = [26 8 10]; cfg.latency = 10; cfg.slicepos = [8 9]; -ft_sourceplot(cfg, sourceProj); \ No newline at end of file +ft_sourceplot(cfg, sourceProj); + +% convert coordinates to volume +% ----------------------------- +function vol = loreta2volume(act, pos, stride) + +minx = min(pos(:,1)); maxx = max(pos(:,1)); +miny = min(pos(:,2)); maxy = max(pos(:,2)); +minz = min(pos(:,3)); maxz = max(pos(:,3)); + +vol = zeros((maxx-minx)/stride(1), (maxy-miny)/stride(2), (maxz-minz)/stride(3)); + +for iPos = 1:size(pos,1) + vol( (pos(iPos,1)-minx)/stride(1)+1, ... + (pos(iPos,2)-miny)/stride(2)+1, ... + (pos(iPos,3)-minz)/stride(3)+1) = act(iPos); + % in case using the sourcemodel with shrunk coordinates + % vol( floor((pos(iPos,1)-minx)/stride(1)+1), ... + % floor((pos(iPos,2)-miny)/stride(2)+1), ... + % floor((pos(iPos,3)-minz)/stride(3)+1)) = act(iPos); +end From 6e17fa95297a9b5d5ba815ee023d2fd71ae8a1e9 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:41:48 -1000 Subject: [PATCH 14/25] now output voxel power --- roi_activity.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/roi_activity.m b/roi_activity.m index c43b0e2..01febf6 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -332,7 +332,8 @@ tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels - + source_voxel_power = tmpWelch; + % fooof settings if strcmpi(g.fooof, 'on') f_range = g.fooof_frange; % freq range where 1/f should be fitted @@ -377,6 +378,7 @@ if strcmpi(g.exportvoxact, 'on') EEG.roi.source_voxel_data = source_voxel_data; % large (takes lots of RAM) end +EEG.roi.source_voxel_power = single(source_voxel_power); EEG.roi.source_roi_data = single(source_roi_data); EEG.roi.source_roi_power = source_roi_power; % used for plotting EEG.roi.source_roi_power_norm = source_roi_power_norm; % used for cross-sprectum From 854d9ff1e0eb2b5a19d8f7397fdc25b6b313bc49 Mon Sep 17 00:00:00 2001 From: arno Date: Tue, 23 May 2023 16:45:31 -0700 Subject: [PATCH 15/25] snippet calculation --- plotconnectivity.m | 5 +++-- pop_roi_connect.m | 54 ++++++++++++++++++++++++++++++--------------- roi_lookupregions.m | 7 ++++-- 3 files changed, 44 insertions(+), 22 deletions(-) diff --git a/plotconnectivity.m b/plotconnectivity.m index 1533e72..d7b14ed 100644 --- a/plotconnectivity.m +++ b/plotconnectivity.m @@ -80,7 +80,7 @@ if isempty(g.labels) if strcmpi(g.brainimg, 'on') disp('Cannot plot on brain with area labels') - g.brainimg = 'off'; + g.brainimg = 'off';g end end @@ -193,6 +193,7 @@ % rename labels % ------------- +if ~strcmpi(g.brainimg, 'off') % CTC added .... to branch around code if not plotting if isempty(g.labels) for iPnt = 1:length(anglesInit) g.labels{iPnt} = sprintf(' Area %d', iPnt); @@ -205,7 +206,7 @@ g.labels{iPnt} = [ ' ' g.labels{iPnt} ]; end end - +end % CTC added .... to branch around above code if not plotting warning off; if isempty(limits) limits(2) = max(array(:)); diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 985210c..f14ca3a 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -31,6 +31,7 @@ % 'MIC' : Maximized Imaginary Coherency for each ROI % 'PAC' : Phase-amplitude coupling between ROIs % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast comptuation). Default is 'off'. % 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets % (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' @@ -170,7 +171,8 @@ 'naccu' 'integer' { } 0; 'methods' 'cell' { } { }; 'snippet' 'string' { 'on', 'off' } 'off'; - 'nepochs' 'real' {} [ ]; + 'firstsnippet' 'string' { 'on', 'off' } 'off'; + 'nepochs' 'real' {} []; 'snip_length' 'integer' { } 60; 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; 'freqresolution' 'integer' { } 0; @@ -254,26 +256,42 @@ fc_matrix = EEG.roi.(fc_name); conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure end - end - - % compute mean over connectivity of each snippet - for fc = 1:n_conn_metrics - fc_name = g.methods{fc}; - [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); - conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell - mat = cell2mat(conn_cell); - reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); - reshaped = squeeze(permute(reshaped, [2,1,3,4])); - if strcmpi(g.fcsave_format, 'all_snips') - EEG.roi.(fc_name) = reshaped; - else - if nsnips > 1 - mean_conn = squeeze(mean(reshaped, 1)); + if strcmpi(g.firstsnippet, 'on') + nsnips = 1; + end + + source_roi_data_save = EEG.roi.source_roi_data; + for isnip = 1:nsnips + roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets + EEG.roi.source_roi_data = single(roi_snip); + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); % compute connectivity over one snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + fc_matrix = EEG.roi.(fc_name); + conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure + end + end + + % compute mean over connectivity of each snippet + for fc = 1:n_conn_metrics + fc_name = options{2}{fc}; + [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc}); + + conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell + mat = cell2mat(conn_cell); + reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim); + reshaped = squeeze(permute(reshaped, [2,1,3,4])); + if strcmpi(g.fcsave_format, 'all_snips') + EEG.roi.(fc_name) = reshaped; else - mean_conn = reshaped; + if nsnips > 1 + mean_conn = squeeze(mean(reshaped, 1)); + else + mean_conn = reshaped; + end + EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end - EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end end elseif strcmpi(g.conn_stats, 'on') diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 56d5ab6..5f526a3 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -75,10 +75,13 @@ end indNet = []; +len = 0; for iNet = 1:length(allnetworknames) if contains(networkname, allnetworknames{iNet}) - indNet = iNet; - break + if length( allnetworknames{iNet} ) > len + indNet = iNet; + len = length( allnetworknames{iNet} ); + end end end From 2b66a22f5c2e311eaf93aba906a3ab71da31b54b Mon Sep 17 00:00:00 2001 From: arno Date: Mon, 17 Jul 2023 09:56:16 -0700 Subject: [PATCH 16/25] update warning message --- roi_lookupregions.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_lookupregions.m b/roi_lookupregions.m index 5f526a3..d46ac09 100644 --- a/roi_lookupregions.m +++ b/roi_lookupregions.m @@ -86,7 +86,7 @@ end if isempty(indNet) - fprintf('Network %s not found\n', networkname); + fprintf('The brain region "%s" was not found (it is ok if it is not a brain region)\n', networkname); elseif ~isstruct(networkdefs) regions = roiTable.(allnetworknames{iNet}); else From a5da10156faebe8daa1eb55483984779065cd357 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:59:05 -1000 Subject: [PATCH 17/25] remnove typo --- plotconnectivity.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plotconnectivity.m b/plotconnectivity.m index d7b14ed..924d4a0 100644 --- a/plotconnectivity.m +++ b/plotconnectivity.m @@ -80,7 +80,7 @@ if isempty(g.labels) if strcmpi(g.brainimg, 'on') disp('Cannot plot on brain with area labels') - g.brainimg = 'off';g + g.brainimg = 'off'; end end From d31bb4a101273eac8525bfce5dc1a6a2f4078b04 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 12:00:21 -1000 Subject: [PATCH 18/25] remnove typo --- pop_roi_connect.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index f14ca3a..a7d6fec 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -31,7 +31,7 @@ % 'MIC' : Maximized Imaginary Coherency for each ROI % 'PAC' : Phase-amplitude coupling between ROIs % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. -% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast comptuation). Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast computation). Default is 'off'. % 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets % (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' From e3f7ef1fecd68ec79be853fa0ca7e2dc49c27f82 Mon Sep 17 00:00:00 2001 From: arno Date: Tue, 23 May 2023 16:45:31 -0700 Subject: [PATCH 19/25] snippet calculation --- pop_roi_connect.m | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index f14ca3a..580f35d 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -239,22 +239,12 @@ snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet nsnips = floor(EEG.trials/snip_eps); if nsnips < 1 - error('Snippet length cannot exceed data length.') - end - diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); - if diff ~= 0 - warning(strcat(int2str(diff), ' seconds are thrown away.')); - end - - source_roi_data_save = EEG.roi.source_roi_data; - for isnip = 1:nsnips - roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets - EEG.roi.source_roi_data = single(roi_snip); - EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection); % compute connectivity over one snippet - for fc = 1:n_conn_metrics - fc_name = g.methods{fc}; - fc_matrix = EEG.roi.(fc_name); - conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure + warning('Snippet length cannot exceed data length. Using the whole data length.') + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); + else + diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); + if diff ~= 0 + warning(strcat(int2str(diff), ' seconds are thrown away.')); end if strcmpi(g.firstsnippet, 'on') From 66b34baf0b5300fd655ea9692a9fc1768561c961 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 11:59:05 -1000 Subject: [PATCH 20/25] remnove typo --- plotconnectivity.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plotconnectivity.m b/plotconnectivity.m index d7b14ed..924d4a0 100644 --- a/plotconnectivity.m +++ b/plotconnectivity.m @@ -80,7 +80,7 @@ if isempty(g.labels) if strcmpi(g.brainimg, 'on') disp('Cannot plot on brain with area labels') - g.brainimg = 'off';g + g.brainimg = 'off'; end end From 90634f4bde7be79ee8cc943259270755e8719cda Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Tue, 18 Jul 2023 12:00:21 -1000 Subject: [PATCH 21/25] remnove typo --- pop_roi_connect.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 580f35d..9bdbfd9 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -31,7 +31,7 @@ % 'MIC' : Maximized Imaginary Coherency for each ROI % 'PAC' : Phase-amplitude coupling between ROIs % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. -% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast comptuation). Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast computation). Default is 'off'. % 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets % (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' From 318224df395acb701d48534a10f74591ae971b20 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Thu, 20 Jul 2023 16:02:57 -1000 Subject: [PATCH 22/25] avoiding to change parent dataset and passing on necessary values --- pop_roi_activity.m | 23 ++- pop_roi_connect.m | 284 ++++++++++++++++++++++++++++- roi_activity.m | 12 +- roi_connect.m | 8 +- roi_pac.m | 2 +- roi_sourceplot.m | 2 +- test_pipes/pipeline_connectivity.m | 2 +- 7 files changed, 312 insertions(+), 21 deletions(-) diff --git a/pop_roi_activity.m b/pop_roi_activity.m index 8255c37..ae38514 100644 --- a/pop_roi_activity.m +++ b/pop_roi_activity.m @@ -17,7 +17,12 @@ % 'resample' - [integer] resample to the desired sampling rate. Default % is 100. Adujst the model order accordingly. ROIconnect % has only be tested with 100 Hz sampling rates. -% 'fooof' - ['on'|'off'] enable FOOOF analysis (this method can be used to parameterize neural power spectra and is described here: https://fooof-tools.github.io/fooof/). Default is 'off'. +% 'effectchanges' - ['on'|'off'] apply resampling and epoching to the +% dataset given as input. Otherwise, only update +% EEG.roi structure. Default is 'off'. +% 'fooof' - ['on'|'off'] enable FOOOF analysis (this method can be +% used to parameterize neural power spectra and is described here: +% https://fooof-tools.github.io/fooof/). Default is 'off'. % 'fooof_frange' - [ ] FOOOF fitting range. Default is [1 30] like in the MATLAB example: % https://github.com/fooof-tools/fooof_mat/blob/main/examples/fooof_example_one_spectrum.m. % 'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). If @@ -244,6 +249,7 @@ 'atlas' 'string' {} ''; 'resample' 'real' {} 100; 'regepochs' 'string' { 'on' 'off'} 'off'; % ignored + 'effectchanges' 'string' { 'on' 'off'} 'off'; 'nPCA' 'real' {} 3; 'epochlen' 'real' {} 2; 'fooof' 'string' { 'on' 'off'} 'off'; @@ -251,12 +257,13 @@ 'fooof_frange' '' {} [1 30]}, 'pop_roi_activity', 'ignore'); if ischar(g), error(g); end +EEGOUT = EEG; if ~isempty(g.resample) && ~isequal(EEG.srate, g.resample) - EEG = pop_resample(EEG, g.resample); + EEGOUT = pop_resample(EEGOUT, g.resample); end -if EEG.trials == 1 +if EEGOUT.trials == 1 recurrence = 2; - EEG = eeg_regepochs(EEG, recurrence, [0 2]); + EEGOUT = eeg_regepochs(EEGOUT, recurrence, [0 2]); end chansel = {}; @@ -276,11 +283,17 @@ g.leadfield = EEG.dipfit.sourcemodel; end -EEG = roi_activity(EEG, 'leadfield', g.leadfield, 'headmodel', EEG.dipfit.hdmfile, ... +EEGOUT = roi_activity(EEGOUT, 'leadfield', g.leadfield, 'headmodel', EEG.dipfit.hdmfile, ... 'model', g.model, 'modelparams', g.modelparams, 'sourcemodel', sourceModelFile, ... 'sourcemodel2mni', sourceModel2MNI, 'nPCA', g.nPCA,'fooof', g.fooof, 'fooof_frange', g.fooof_frange, ... 'freqresolution', g.freqresolution, 'sourcemodelatlas', g.atlas, 'chansel', chansel, moreargs{:}); +if strcmpi(g.effectchanges, 'on') + EEG = EEGOUT; +else + EEG.roi = EEGOUT.roi; +end + if nargout > 1 for iOption = 1:2:length(options) if strcmpi(options{iOption}, 'leadfield') && isequal(options{iOption+1}, EEG.dipfit.sourcemodel) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 9bdbfd9..5675777 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -239,12 +239,22 @@ snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet nsnips = floor(EEG.trials/snip_eps); if nsnips < 1 - warning('Snippet length cannot exceed data length. Using the whole data length.') - EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); - else - diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); - if diff ~= 0 - warning(strcat(int2str(diff), ' seconds are thrown away.')); + error('Snippet length cannot exceed data length.') + end + diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); + if diff ~= 0 + warning(strcat(int2str(diff), ' seconds are thrown away.')); + end + + source_roi_data_save = EEG.roi.source_roi_data; + for isnip = 1:nsnips + roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets + EEG.roi.source_roi_data = single(roi_snip); + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection); % compute connectivity over one snippet + for fc = 1:n_conn_metrics + fc_name = g.methods{fc}; + fc_matrix = EEG.roi.(fc_name); + conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure end if strcmpi(g.firstsnippet, 'on') @@ -303,3 +313,265 @@ com = sprintf( 'EEG = pop_roi_connect(EEG, %s);', vararg2str( options )); end + + + + +% pop_roi_connect - call roi_connect to connectivity between ROIs +% +% Usage: +% EEG = pop_roi_connect(EEG, 'key', 'val', ...); +% +% Inputs: +% EEG - EEGLAB dataset containing ROI activity +% +% Optional inputs: +% 'morder' - [integer] Order of autoregressive model. Default is 20. +% 'naccu' - [integer] Number of accumulation for stats. Default is 0. +% 'methods' - [cell] Cell of strings corresponding to methods. +% 'CS' : Cross spectrum +% 'aCOH' : Coherence +% 'cCOH' : (complex-valued) coherency +% 'iCOH' : absolute value of the imaginary part of coherency +% 'GC' : Granger Causality +% 'TRGC' : Time-reversed Granger Causality +% 'wPLI' : Weighted Phase Lag Index +% 'PDC' : Partial directed coherence +% 'TRPDC' : Time-reversed partial directed coherence +% 'DTF' : Directed transfer entropy +% 'TRDTF' : Time-reversed directed transfer entropy +% 'MIM' : Multivariate Interaction Measure for each ROI +% 'MIC' : Maximized Imaginary Coherency for each ROI +% 'PAC' : Phase-amplitude coupling between ROIs +% 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. +% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast computation). Default is 'off'. +% 'errornosnippet' - ['on'|off] If data shorter than snippet length, issue error (otherwise compute). Default is 'on'. +% 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. +% 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets +% (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' +% 'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). +% If specified, the signal is zero padded accordingly. Default is 0 (means no padding). +% 'fcomb' - [struct] Frequency combination for which PAC is computed (in Hz). Must have fields 'low' and +% 'high' with fcomb.low < fcomb.high. For example, fcomb.low = 10 and fcomb.high = 50 if single +% frequencies are used. fcomb.low = [4 8] and fcomb.high = [48 50] if frequency bands are used +% (might take a long time to compute so use with caution). Default is {} (this will cause an error). +% 'bs_outopts' - [integer] Option which bispectral tensors should be stored in EEG.roi.PAC. Default is 1. +% 1 - store all tensors: b_orig, b_anti, b_orig_norm, b_anti_norm +% 2 - only store: b_orig, b_anti +% 3 - only store: b_orig_norm, b_anti_norm +% 4 - only store: b_orig, b_orig_norm +% 5 - only store: b_anti, b_anti_norm +% 'roi_selection' - [cell array of integers] Cell array of ROI indices {1, 2, 3, ...} indicating for which regions (ROIs) connectivity should be computed. +% Default is empty (in this case, connectivity will be computed for all ROIs). +% 'conn_stats' - ['on'|'off'] Run statistics on connectivity metrics. Default is 'off'. +% 'nshuf' - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 1001. +% 'freqrange' - [min max] frequency range in Hz. This is used to compute and plot p-values. Default is to plot broadband power. +% +% +% Output: +% EEG - EEGLAB dataset with field 'roi' containing connectivity info. +% +% Note: Optional inputs to roi_connectivity_process() are also accepted. +% +% Author: Arnaud Delorme, UCSD, 2019 +% +% Example +% p = fileparts(which('eeglab')); % path +% EEG = pop_roi_connect(EEG, 'headmodel', ... +% EEG.dipfit.hdmfile, 'elec2mni', EEG.dipfit.coord_transform, ... +% 'sourcemodel', fullfile(p, 'functions', 'supportfiles', ... +% 'head_modelColin27_5003_Standard-10-5-Cap339.mat'), 'sourcemodel2mni', ... +% [0 -26.6046230000 -46 0.1234625600 0 -1.5707963000 1000 1000 1000]); +% +% Use pop_roi_connectivity(EEG) to conectivity + +% Copyright (C) Arnaud Delorme, arnodelorme@gmail.com +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: +% +% 1. Redistributions of source code must retain the above copyright notice, +% this list of conditions and the following disclaimer. +% +% 2. Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +% THE POSSIBILITY OF SUCH DAMAGE. + +% TO DO - Arno +% - Centralize reading head mesh and Atlas (there might be a function in +% Fieldtrip to do that) ft_read_volume ft_read_mesh +% - Make compatible with all Fieldtrip and FSL Atlases +% - Downsampling of Atlas - check bug submitted to Fieldtrip +% - Plot inside(blue) vs outside(red) voxels for source volume + +function [EEG,com] = pop_roi_connect(EEG, varargin) + +com = ''; +if nargin < 1 + help pop_roi_connect; + return +end + +if ~isfield(EEG(1), 'roi') || ~isfield(EEG(1).roi, 'source_roi_data') + error('Cannot find ROI data - ROI data first'); +end + +if nargin < 2 + + rowg = [0.1 0.6 1 0.2]; + % uigeom = { 1 1 rowg rowg 1 rowg rowg [0.1 0.6 0.9 0.3] 1 rowg 1 [0.5 1 0.35 0.5] [0.5 1 0.35 0.5] [0.5 1 0.35 0.5] [1] [0.9 1.2 1] }; + uigeom = { [1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1] [0.2 1 0.35 0.8] [0.2 1 0.35 0.8] }; + uilist = { { 'style' 'text' 'string' 'Select connectivity measures' 'fontweight' 'bold' } ... + { 'style' 'checkbox' 'string' 'Cross-spectrum' 'tag' 'cs' 'value' 1 } {} ... + {'style' 'checkbox' 'string' '(Complex-valued) Coherency' 'tag' 'ccoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Coherence' 'tag' 'acoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Imaginary Coherency' 'tag' 'icoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Weighted Phase Lag Index' 'tag' 'wpli' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Granger Causality (GC)' 'tag' 'gc' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Time-reversed GC' 'tag' 'trgc' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Partial Directed Coherence (PDC)' 'tag' 'pdc' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Time-reversed PDC' 'tag' 'trpdc' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Directed Transfer Entropy (DTF)' 'tag' 'dtf' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Time-reversed DTF' 'tag' 'trdtf' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Multivariate Interaction Measure' 'tag' 'mim' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Maximized Imaginary Coherency' 'tag' 'mic' 'value' 0 } ... + {} ... + {} { 'style' 'text' 'string' 'Autoregressive model order' } { 'style' 'edit' 'string' '20' 'tag' 'morder' } {} ... + {} { 'style' 'text' 'string' 'Bootstrap if any (n)' } { 'style' 'edit' 'string' '' 'tag' 'naccu2' } {} }; + ... + [result,~,~,out] = inputgui('geometry', uigeom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_roi_connect'')', 'title', 'pop_roiconnect - connectivity'); + if isempty(result), return, end + + % check we have the same naccu + methods = {}; + if out.cs, methods = [ methods { 'CS' } ]; end +% if out.coh, methods = [ methods { 'COH' } ]; end + if out.ccoh, methods = [ methods { 'cCOH' } ]; end + if out.acoh, methods = [ methods { 'aCOH' } ]; end + if out.icoh, methods = [ methods { 'iCOH' } ]; end + if out.gc , methods = [ methods { 'GC' } ]; end + if out.trgc, methods = [ methods { 'TRGC' } ]; end + if out.wpli, methods = [ methods { 'wPLI' } ]; end + if out.pdc , methods = [ methods { 'PDC' } ]; end + if out.trpdc, methods = [ methods { 'TRPDC' } ]; end + if out.dtf , methods = [ methods { 'DTF' } ]; end + if out.trdtf, methods = [ methods { 'TRDTF' } ]; end + if out.mim , methods = [ methods { 'MIM' } ]; end + if out.mic, methods = [ methods { 'MIC' } ]; end + options = { ... + 'morder' str2num(out.morder) ... + 'naccu' str2num(out.naccu2) ... + 'methods' methods }; +else + options = varargin; +end + +% decode input parameters +% ----------------------- +g = finputcheck(options, ... + { 'morder' 'integer' { } 20; + 'naccu' 'integer' { } 0; + 'methods' 'cell' { } { }; + 'snippet' 'string' { 'on', 'off' } 'off'; + 'errornosnippet' 'string' { 'on', 'off' } 'on'; + 'firstsnippet' 'string' { 'on', 'off' } 'off'; + 'nepochs' 'real' {} []; + 'snip_length' 'integer' { } 60; + 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; + 'freqresolution' 'integer' { } 0; + 'fcomb' 'struct' { } struct; + 'bs_outopts' 'integer' { } 1; + 'roi_selection' 'cell' { } { }; + 'conn_stats' 'string' { } 'off'; ... + 'nshuf' 'integer' { } 1001}, 'pop_roi_connect'); +if ischar(g), error(g); end + +% process multiple datasets +% ------------------------- +if length(EEG) > 1 + eeglab_options; + if nargin < 2 + if option_storedisk + res = questdlg2( [ 'Data files on disk will be automatically overwritten.' 10 ... + 'Are you sure you want to proceed with this operation?' ], ... + 'Confirmation', 'Cancel', 'Proceed', 'Proceed'); + switch lower(res) + case 'cancel', return; + case 'proceed' + end + end + end + % find common datasets accross subjects + % and limit the number of epochs + allsubjects = { EEG.subject }; + for index = 1:length(allsubjects) + optionsSubj = options; + allinds = strmatch(allsubjects{index}, allsubjects, 'exact'); + + % datasets from the same subjects need to have the same number of epochs + TMPEEG = EEG(allinds); + if TMPEEG(1).trials > 1 + optionsSubj = [ optionsSubj { 'nepochs' min([TMPEEG.trials]) }]; + end + if nargin < 2 + [ TMPEEG, com ] = eeg_eval( 'pop_roi_connect', TMPEEG, 'warning', 'off', 'params', optionsSubj ); + else + [ TMPEEG, com ] = eeg_eval( 'pop_roi_connect', TMPEEG, 'params', optionsSubj ); + end + EEG = eeg_store(EEG, TMPEEG, allinds); + end + return +end + +% compute connectivity over snippets +n_conn_metrics = length(g.methods); % number of connectivity metrics +conn_matrices_snips = {}; +if strcmpi(g.snippet, 'on') && isempty(intersect(g.methods, {'PAC'})) && strcmpi(g.conn_stats, 'off') + + trials = size(EEG.roi.source_roi_data, 3); + snippet_length = g.snip_length; % seconds + snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet + nsnips = floor(EEG.trials/snip_eps); + if nsnips < 1 + warning('Snippet length cannot exceed data length. Using the whole data length.') + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); + else + diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); + if diff ~= 0 + warning(strcat(int2str(diff), ' seconds are thrown away.')); + end + + xxxxxxx + end + +elseif strcmpi(g.conn_stats, 'on') + EEG = roi_connstats(EEG, g.methods, g.nshuf, g.roi_selection); +else + EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, ... + 'roi_selection', g.roi_selection); +end + +if ~isempty(intersect(g.methods, {'PAC'})) + if strcmpi(g.snippet, 'on') + error('Snippet analysis for PAC has not been implemented yet.') + else + EEG = roi_pac(EEG, g.fcomb, g.bs_outopts, g.roi_selection); + end +end + +if nargout > 1 + com = sprintf( 'EEG = pop_roi_connect(EEG, %s);', vararg2str( options )); +end + diff --git a/roi_activity.m b/roi_activity.m index 01febf6..d08ff64 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -220,7 +220,11 @@ error('There must be the same number of vertices/voxels in the leadfield and source model'); end if isempty(g.chansel) - g.chansel = [1:EEG.nbchan]; + if isfield(EEG.dipfit, 'chansel') + g.chansel = EEG.dipfit.chansel; + else + g.chansel = 1:EEG.nbchan; + end else g.chansel = eeg_decodechan(EEG.chanlocs, g.chansel); end @@ -257,7 +261,7 @@ alpha = lcmv_reg*trace(C)/length(C); Cr = C + alpha*eye(nbchan); [~, P_eloreta] = lcmv(Cr, leadfield, struct('alpha', 0, 'onedim', 0)); - source_voxel_data = reshape(tmpdata(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3); + source_voxel_data = single(reshape(tmpdata(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3)); source_voxel_data = 10^3*source_voxel_data; % the units are nA*m else % transform the data to continuous so we can get an estimate for each sample @@ -328,7 +332,7 @@ % compute power using the Welch method disp('Computing ROI activity...'); - [tmpWelch,ftmp] = pwelch(tmpData(:,:), data_pnts, floor(data_pnts/2), data_pnts, EEG.srate); % ftmp should be equal frqs + [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, 0, data_pnts, EEG.srate); % ftmp should be equal frqs tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels @@ -385,6 +389,8 @@ EEG.roi.freqs = frqs; EEG.roi.nPCA = g.nPCA; EEG.roi.nROI = nROI; +EEG.roi.pnts = EEG.pnts; +EEG.roi.srate = EEG.srate; EEG.roi.atlas = cortex.Atlas; EEG.roi.srate = EEG.srate; EEG.roi.leadfield = g.leadfield; diff --git a/roi_connect.m b/roi_connect.m index 2d6bfda..b608ee7 100644 --- a/roi_connect.m +++ b/roi_connect.m @@ -114,9 +114,9 @@ tmpMethods1 = intersect(g.methods, methodset1); if ~isempty(tmpMethods1) if isempty(g.roi_selection) - conn_mult = data2sctrgcmim(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution); + conn_mult = data2sctrgcmim(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution); elseif ~isempty(g.roi_selection) - conn_mult = data2sctrgcmim(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA); + conn_mult = data2sctrgcmim(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA); EEG.roi.roi_selection = g.roi_selection; end fields = fieldnames(conn_mult); @@ -147,9 +147,9 @@ tmpMethods2 = intersect(g.methods, methodset2); if ~isempty(tmpMethods2) if isempty(g.roi_selection) - conn_mult = data2spwctrgc(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution); + conn_mult = data2spwctrgc(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution); else - conn_mult = data2spwctrgc(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA); + conn_mult = data2spwctrgc(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA); end fields = fieldnames(conn_mult); for iField = 1:length(fields) diff --git a/roi_pac.m b/roi_pac.m index 372f67a..36e849b 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -52,7 +52,7 @@ segleng = ndat; segshift = floor(ndat/2); epleng = ndat; - fs = EEG.srate; + fs = EEG.roi.srate; params = []; params.nROI = nROI; diff --git a/roi_sourceplot.m b/roi_sourceplot.m index 0a7d4ae..5a0bc67 100644 --- a/roi_sourceplot.m +++ b/roi_sourceplot.m @@ -12,7 +12,7 @@ % MNI locations. % % Required inputs: -% 'freqselect' - [real] frequency of interest or frequency range of interest. +% 'freqrange' - [real] frequency of interest or frequency range of interest. % Defaut is all frequencies. % % Example: diff --git a/test_pipes/pipeline_connectivity.m b/test_pipes/pipeline_connectivity.m index 6e9ea48..ddce455 100644 --- a/test_pipes/pipeline_connectivity.m +++ b/test_pipes/pipeline_connectivity.m @@ -25,7 +25,7 @@ EEG = pop_dipfit_settings( EEG, 'hdmfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_vol.mat'), ... 'coordformat','MNI','mrifile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_mri.mat'),... 'chanfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','elec', 'standard_1005.elc'),... - 'coord_transform',[0.83215 -15.6287 2.4114 0.081214 0.00093739 -1.5732 1.1742 1.0601 1.1485] ,'chansel',[1:32] ); + 'coord_transform',[0.83215 -15.6287 2.4114 0.081214 0.00093739 -1.5732 1.1742 1.0601 1.1485] ); EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ... 'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1); From 88de73103fc07966ab72fed0cc942fa6d3e5fae9 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Thu, 20 Jul 2023 16:28:23 -1000 Subject: [PATCH 23/25] add back the erroronsnippet and remove the changes to roi_activity --- pop_roi_connect.m | 279 +++------------------------------------------- roi_activity.m | 4 +- 2 files changed, 15 insertions(+), 268 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 5675777..92dcf31 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -32,7 +32,8 @@ % 'PAC' : Phase-amplitude coupling between ROIs % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. % 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast computation). Default is 'off'. -% 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. +% 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. +% 'erroronsnippet' - ['on'|'off'] Error if snippet too short. Default 'on'. % 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets % (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' % 'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). @@ -171,7 +172,8 @@ 'naccu' 'integer' { } 0; 'methods' 'cell' { } { }; 'snippet' 'string' { 'on', 'off' } 'off'; - 'firstsnippet' 'string' { 'on', 'off' } 'off'; + 'firstsnippet' 'string' { 'on', 'off' } 'off'; + 'erroronsnippet' 'string' { 'on', 'off' } 'off'; 'nepochs' 'real' {} []; 'snip_length' 'integer' { } 60; 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; @@ -239,7 +241,12 @@ snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet nsnips = floor(EEG.trials/snip_eps); if nsnips < 1 - error('Snippet length cannot exceed data length.') + if strcmpi(opt.erroronsnippet, 'on') + error('Snippet length cannot exceed data length.') + else + frpintf(2, 'Snippet length cannot exceed data length, using the whole data') + nsnips = 1; + end end diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); if diff ~= 0 @@ -248,7 +255,9 @@ source_roi_data_save = EEG.roi.source_roi_data; for isnip = 1:nsnips - roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets + begSnip = (isnip-1)* snip_eps + 1; + endSnip = min((isnip-1)* snip_eps + snip_eps, size(source_roi_data_save,3)); + roi_snip = source_roi_data_save(:,:, begSnip:endSnip ); % cut source data into snippets EEG.roi.source_roi_data = single(roi_snip); EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection); % compute connectivity over one snippet for fc = 1:n_conn_metrics @@ -313,265 +322,3 @@ com = sprintf( 'EEG = pop_roi_connect(EEG, %s);', vararg2str( options )); end - - - - -% pop_roi_connect - call roi_connect to connectivity between ROIs -% -% Usage: -% EEG = pop_roi_connect(EEG, 'key', 'val', ...); -% -% Inputs: -% EEG - EEGLAB dataset containing ROI activity -% -% Optional inputs: -% 'morder' - [integer] Order of autoregressive model. Default is 20. -% 'naccu' - [integer] Number of accumulation for stats. Default is 0. -% 'methods' - [cell] Cell of strings corresponding to methods. -% 'CS' : Cross spectrum -% 'aCOH' : Coherence -% 'cCOH' : (complex-valued) coherency -% 'iCOH' : absolute value of the imaginary part of coherency -% 'GC' : Granger Causality -% 'TRGC' : Time-reversed Granger Causality -% 'wPLI' : Weighted Phase Lag Index -% 'PDC' : Partial directed coherence -% 'TRPDC' : Time-reversed partial directed coherence -% 'DTF' : Directed transfer entropy -% 'TRDTF' : Time-reversed directed transfer entropy -% 'MIM' : Multivariate Interaction Measure for each ROI -% 'MIC' : Maximized Imaginary Coherency for each ROI -% 'PAC' : Phase-amplitude coupling between ROIs -% 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. -% 'firstsnippet' - ['on'|off] Only use the first snippet (useful for fast computation). Default is 'off'. -% 'errornosnippet' - ['on'|off] If data shorter than snippet length, issue error (otherwise compute). Default is 'on'. -% 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds. -% 'fcsave_format' - ['mean_snips'|'all_snips'] Option to save mean over snippets -% (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.' -% 'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). -% If specified, the signal is zero padded accordingly. Default is 0 (means no padding). -% 'fcomb' - [struct] Frequency combination for which PAC is computed (in Hz). Must have fields 'low' and -% 'high' with fcomb.low < fcomb.high. For example, fcomb.low = 10 and fcomb.high = 50 if single -% frequencies are used. fcomb.low = [4 8] and fcomb.high = [48 50] if frequency bands are used -% (might take a long time to compute so use with caution). Default is {} (this will cause an error). -% 'bs_outopts' - [integer] Option which bispectral tensors should be stored in EEG.roi.PAC. Default is 1. -% 1 - store all tensors: b_orig, b_anti, b_orig_norm, b_anti_norm -% 2 - only store: b_orig, b_anti -% 3 - only store: b_orig_norm, b_anti_norm -% 4 - only store: b_orig, b_orig_norm -% 5 - only store: b_anti, b_anti_norm -% 'roi_selection' - [cell array of integers] Cell array of ROI indices {1, 2, 3, ...} indicating for which regions (ROIs) connectivity should be computed. -% Default is empty (in this case, connectivity will be computed for all ROIs). -% 'conn_stats' - ['on'|'off'] Run statistics on connectivity metrics. Default is 'off'. -% 'nshuf' - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 1001. -% 'freqrange' - [min max] frequency range in Hz. This is used to compute and plot p-values. Default is to plot broadband power. -% -% -% Output: -% EEG - EEGLAB dataset with field 'roi' containing connectivity info. -% -% Note: Optional inputs to roi_connectivity_process() are also accepted. -% -% Author: Arnaud Delorme, UCSD, 2019 -% -% Example -% p = fileparts(which('eeglab')); % path -% EEG = pop_roi_connect(EEG, 'headmodel', ... -% EEG.dipfit.hdmfile, 'elec2mni', EEG.dipfit.coord_transform, ... -% 'sourcemodel', fullfile(p, 'functions', 'supportfiles', ... -% 'head_modelColin27_5003_Standard-10-5-Cap339.mat'), 'sourcemodel2mni', ... -% [0 -26.6046230000 -46 0.1234625600 0 -1.5707963000 1000 1000 1000]); -% -% Use pop_roi_connectivity(EEG) to conectivity - -% Copyright (C) Arnaud Delorme, arnodelorme@gmail.com -% -% Redistribution and use in source and binary forms, with or without -% modification, are permitted provided that the following conditions are met: -% -% 1. Redistributions of source code must retain the above copyright notice, -% this list of conditions and the following disclaimer. -% -% 2. Redistributions in binary form must reproduce the above copyright notice, -% this list of conditions and the following disclaimer in the documentation -% and/or other materials provided with the distribution. -% -% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF -% THE POSSIBILITY OF SUCH DAMAGE. - -% TO DO - Arno -% - Centralize reading head mesh and Atlas (there might be a function in -% Fieldtrip to do that) ft_read_volume ft_read_mesh -% - Make compatible with all Fieldtrip and FSL Atlases -% - Downsampling of Atlas - check bug submitted to Fieldtrip -% - Plot inside(blue) vs outside(red) voxels for source volume - -function [EEG,com] = pop_roi_connect(EEG, varargin) - -com = ''; -if nargin < 1 - help pop_roi_connect; - return -end - -if ~isfield(EEG(1), 'roi') || ~isfield(EEG(1).roi, 'source_roi_data') - error('Cannot find ROI data - ROI data first'); -end - -if nargin < 2 - - rowg = [0.1 0.6 1 0.2]; - % uigeom = { 1 1 rowg rowg 1 rowg rowg [0.1 0.6 0.9 0.3] 1 rowg 1 [0.5 1 0.35 0.5] [0.5 1 0.35 0.5] [0.5 1 0.35 0.5] [1] [0.9 1.2 1] }; - uigeom = { [1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1.2 1] [1] [0.2 1 0.35 0.8] [0.2 1 0.35 0.8] }; - uilist = { { 'style' 'text' 'string' 'Select connectivity measures' 'fontweight' 'bold' } ... - { 'style' 'checkbox' 'string' 'Cross-spectrum' 'tag' 'cs' 'value' 1 } {} ... - {'style' 'checkbox' 'string' '(Complex-valued) Coherency' 'tag' 'ccoh' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Coherence' 'tag' 'acoh' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Imaginary Coherency' 'tag' 'icoh' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Weighted Phase Lag Index' 'tag' 'wpli' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Granger Causality (GC)' 'tag' 'gc' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Time-reversed GC' 'tag' 'trgc' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Partial Directed Coherence (PDC)' 'tag' 'pdc' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Time-reversed PDC' 'tag' 'trpdc' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Directed Transfer Entropy (DTF)' 'tag' 'dtf' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Time-reversed DTF' 'tag' 'trdtf' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Multivariate Interaction Measure' 'tag' 'mim' 'value' 0 } ... - { 'style' 'checkbox' 'string' 'Maximized Imaginary Coherency' 'tag' 'mic' 'value' 0 } ... - {} ... - {} { 'style' 'text' 'string' 'Autoregressive model order' } { 'style' 'edit' 'string' '20' 'tag' 'morder' } {} ... - {} { 'style' 'text' 'string' 'Bootstrap if any (n)' } { 'style' 'edit' 'string' '' 'tag' 'naccu2' } {} }; - ... - [result,~,~,out] = inputgui('geometry', uigeom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_roi_connect'')', 'title', 'pop_roiconnect - connectivity'); - if isempty(result), return, end - - % check we have the same naccu - methods = {}; - if out.cs, methods = [ methods { 'CS' } ]; end -% if out.coh, methods = [ methods { 'COH' } ]; end - if out.ccoh, methods = [ methods { 'cCOH' } ]; end - if out.acoh, methods = [ methods { 'aCOH' } ]; end - if out.icoh, methods = [ methods { 'iCOH' } ]; end - if out.gc , methods = [ methods { 'GC' } ]; end - if out.trgc, methods = [ methods { 'TRGC' } ]; end - if out.wpli, methods = [ methods { 'wPLI' } ]; end - if out.pdc , methods = [ methods { 'PDC' } ]; end - if out.trpdc, methods = [ methods { 'TRPDC' } ]; end - if out.dtf , methods = [ methods { 'DTF' } ]; end - if out.trdtf, methods = [ methods { 'TRDTF' } ]; end - if out.mim , methods = [ methods { 'MIM' } ]; end - if out.mic, methods = [ methods { 'MIC' } ]; end - options = { ... - 'morder' str2num(out.morder) ... - 'naccu' str2num(out.naccu2) ... - 'methods' methods }; -else - options = varargin; -end - -% decode input parameters -% ----------------------- -g = finputcheck(options, ... - { 'morder' 'integer' { } 20; - 'naccu' 'integer' { } 0; - 'methods' 'cell' { } { }; - 'snippet' 'string' { 'on', 'off' } 'off'; - 'errornosnippet' 'string' { 'on', 'off' } 'on'; - 'firstsnippet' 'string' { 'on', 'off' } 'off'; - 'nepochs' 'real' {} []; - 'snip_length' 'integer' { } 60; - 'fcsave_format' 'string' { 'mean_snips', 'all_snips'} 'mean_snips'; - 'freqresolution' 'integer' { } 0; - 'fcomb' 'struct' { } struct; - 'bs_outopts' 'integer' { } 1; - 'roi_selection' 'cell' { } { }; - 'conn_stats' 'string' { } 'off'; ... - 'nshuf' 'integer' { } 1001}, 'pop_roi_connect'); -if ischar(g), error(g); end - -% process multiple datasets -% ------------------------- -if length(EEG) > 1 - eeglab_options; - if nargin < 2 - if option_storedisk - res = questdlg2( [ 'Data files on disk will be automatically overwritten.' 10 ... - 'Are you sure you want to proceed with this operation?' ], ... - 'Confirmation', 'Cancel', 'Proceed', 'Proceed'); - switch lower(res) - case 'cancel', return; - case 'proceed' - end - end - end - % find common datasets accross subjects - % and limit the number of epochs - allsubjects = { EEG.subject }; - for index = 1:length(allsubjects) - optionsSubj = options; - allinds = strmatch(allsubjects{index}, allsubjects, 'exact'); - - % datasets from the same subjects need to have the same number of epochs - TMPEEG = EEG(allinds); - if TMPEEG(1).trials > 1 - optionsSubj = [ optionsSubj { 'nepochs' min([TMPEEG.trials]) }]; - end - if nargin < 2 - [ TMPEEG, com ] = eeg_eval( 'pop_roi_connect', TMPEEG, 'warning', 'off', 'params', optionsSubj ); - else - [ TMPEEG, com ] = eeg_eval( 'pop_roi_connect', TMPEEG, 'params', optionsSubj ); - end - EEG = eeg_store(EEG, TMPEEG, allinds); - end - return -end - -% compute connectivity over snippets -n_conn_metrics = length(g.methods); % number of connectivity metrics -conn_matrices_snips = {}; -if strcmpi(g.snippet, 'on') && isempty(intersect(g.methods, {'PAC'})) && strcmpi(g.conn_stats, 'off') - - trials = size(EEG.roi.source_roi_data, 3); - snippet_length = g.snip_length; % seconds - snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet - nsnips = floor(EEG.trials/snip_eps); - if nsnips < 1 - warning('Snippet length cannot exceed data length. Using the whole data length.') - EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); - else - diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps); - if diff ~= 0 - warning(strcat(int2str(diff), ' seconds are thrown away.')); - end - - xxxxxxx - end - -elseif strcmpi(g.conn_stats, 'on') - EEG = roi_connstats(EEG, g.methods, g.nshuf, g.roi_selection); -else - EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, ... - 'roi_selection', g.roi_selection); -end - -if ~isempty(intersect(g.methods, {'PAC'})) - if strcmpi(g.snippet, 'on') - error('Snippet analysis for PAC has not been implemented yet.') - else - EEG = roi_pac(EEG, g.fcomb, g.bs_outopts, g.roi_selection); - end -end - -if nargout > 1 - com = sprintf( 'EEG = pop_roi_connect(EEG, %s);', vararg2str( options )); -end - diff --git a/roi_activity.m b/roi_activity.m index d08ff64..bf7f0e4 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -261,7 +261,7 @@ alpha = lcmv_reg*trace(C)/length(C); Cr = C + alpha*eye(nbchan); [~, P_eloreta] = lcmv(Cr, leadfield, struct('alpha', 0, 'onedim', 0)); - source_voxel_data = single(reshape(tmpdata(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3)); + source_voxel_data = reshape(tmpdata(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3); source_voxel_data = 10^3*source_voxel_data; % the units are nA*m else % transform the data to continuous so we can get an estimate for each sample @@ -332,7 +332,7 @@ % compute power using the Welch method disp('Computing ROI activity...'); - [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, 0, data_pnts, EEG.srate); % ftmp should be equal frqs + [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, data_pnts/2, data_pnts, EEG.srate); % ftmp should be equal frqs tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels From 5beca6f8ec4c3794c5f307f89355c8b582c5b2b2 Mon Sep 17 00:00:00 2001 From: Arnaud Delorme Date: Thu, 20 Jul 2023 16:29:40 -1000 Subject: [PATCH 24/25] add back the erroronsnippet and remove the changes to roi_activity --- roi_activity.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_activity.m b/roi_activity.m index bf7f0e4..5c277e2 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -332,7 +332,7 @@ % compute power using the Welch method disp('Computing ROI activity...'); - [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, data_pnts/2, data_pnts, EEG.srate); % ftmp should be equal frqs + [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, floor(data_pnts/2), data_pnts, EEG.srate); % ftmp should be equal frqs tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3)); tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3 tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels From 6be4b885cf5c7e0e091a90c48e10ed896f4de114 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 21 Jul 2023 10:45:34 +0200 Subject: [PATCH 25/25] add error message if p-values are all the same --- pop_roi_statsplot.m | 18 ++++++++++++------ roi_connect.m | 2 +- test_pipes/test_mim_stats.m | 23 ++++------------------- 3 files changed, 17 insertions(+), 26 deletions(-) diff --git a/pop_roi_statsplot.m b/pop_roi_statsplot.m index 649b1da..ef3864e 100644 --- a/pop_roi_statsplot.m +++ b/pop_roi_statsplot.m @@ -66,16 +66,22 @@ netFC = squeeze(mean(matrix, 2)); FC_pn = sum(netFC(:, 1) < netFC(:, 2:end), 2)./(size(matrix, 3) - 1); - % use FDR-correction for multiple comparison's correction - [p_fdr, ~] = fdr(FC_pn, 0.05); + % use FDR-correction for multiple comparison's correction + alpha = 0.05; % add this as a parameter? + [p_fdr, ~] = fdr(FC_pn, alpha); FC_pn(FC_pn > p_fdr) = 1; % plot load cm17; load cortex; - FC_pn(FC_pn==0) = 1 / (size(netMIM, 2) - 1); % 1 / nshuf + FC_pn(FC_pn==0) = 1 / (size(netFC, 2) - 1); % 1 / nshuf data = -log10(FC_pn); - allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3); - h = textsc(title, 'title'); - set(h, 'fontsize', 20); + try + allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3); + h = textsc(title, 'title'); + set(h, 'fontsize', 20); + catch + warning('There are no "significant" p-values to be plotted. These are the p-values:') + disp(FC_pn) + end end \ No newline at end of file diff --git a/roi_connect.m b/roi_connect.m index 2d6bfda..2731aad 100644 --- a/roi_connect.m +++ b/roi_connect.m @@ -137,7 +137,7 @@ elseif strcmpi(tmpMethods1{iMethods}, 'GC') || strcmpi(tmpMethods1{iMethods}, 'TRGC') TRGCnet = EEG.roi.(tmpMethods1{iMethods})(:, :, 1) - EEG.roi.(tmpMethods1{iMethods})(:, :, 2); EEG.roi.(tmpMethods1{iMethods}) = get_connect_mat( TRGCnet, nROI, -1); - else + else % wPLI measure = rm_components(EEG.roi.(tmpMethods1{iMethods}), EEG.roi.nPCA, tmpMethods1{iMethods}); % only keep the first principal component EEG.roi.(tmpMethods1{iMethods}) = measure; end diff --git a/test_pipes/test_mim_stats.m b/test_pipes/test_mim_stats.m index a4e037e..4660e30 100644 --- a/test_pipes/test_mim_stats.m +++ b/test_pipes/test_mim_stats.m @@ -8,7 +8,7 @@ EEG = pop_resample( EEG, 100); EEG = pop_epoch( EEG, { }, [-0.5 1.5], 'newname', 'EEG Data epochs epochs', 'epochinfo', 'yes'); EEG = pop_select( EEG, 'trial',1:30); -[ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG); +% [ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG); eeglab redraw; EEG = pop_dipfit_settings( EEG, 'hdmfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_vol.mat'), ... @@ -22,21 +22,6 @@ EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3); %% Create null distribution -EEG = pop_roi_connect(EEG, 'methods', {'COH', 'MIM'}, 'conn_stats', 'on', 'nshuf', 3); - -% ROI selection -MIM = squeeze(mean(EEG.roi.MIM, 1)); % broadband for now - -% generate p-values by comparing true MIM to null distribution -netMIM = squeeze(mean(MIM, 2)); -MIM_pn = sum(netMIM(:,1) < netMIM(:,2:end),2)./(size(MIM,3)-1); - -% plot -load cm17; -load cortex; -MIM_pn(MIM_pn==0) = 1 / (size(netMIM, 2)-1); % 1 / nshuf -data = -log10(MIM_pn); -% allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3); -allplots_cortex_BS(cortex_highres, data, [0 5], cm17a ,'-log(p)', 0.3); - - +% EEG = pop_roi_connect(EEG, 'methods', {'COH', 'MIM'}, 'conn_stats', 'on', 'nshuf', 3); +load('test_pipes/MIM_shuf.mat') +pop_roi_statsplot(EEG, 'measure', 'MIM', 'freqrange', [8 13]); \ No newline at end of file