From 521f8bfde72d4b63ae776b9fc4b88d3a9faa8250 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 12:11:27 +0200 Subject: [PATCH 1/8] bugfix in plotting --- pop_roi_connectplot.m | 26 ++++++++++++++++---------- test_pipes/test_roi_selection.m | 3 +++ 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 974cfd3..3499481 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -318,12 +318,14 @@ if isempty(g.measure) error('You must define a measure to plot'); end - - if ~isempty(S.roi_selection) - warning("Plotting options ('region', 'hemisphere', 'grouphemispheres') are disabled when ROIs were explicitely selected.") - g.region = 'all'; - g.hemisphere = 'all'; - g.grouphemispheres = 'all'; + + if isfield(S, 'roi_selection') + if ~isempty(S.roi_selection) + warning("Plotting options ('region', 'hemisphere', 'grouphemispheres') are disabled when ROIs were explicitely selected.") + g.region = 'all'; + g.hemisphere = 'all'; + g.grouphemispheres = 'all'; + end end % colormap @@ -541,8 +543,10 @@ labels = cellstr(labels); % remove region labels that were not selected - if ~isempty(EEG.roi.roi_selection) - labels = labels(cell2mat(EEG.roi.roi_selection)); + if isfield(EEG.roi, 'roi_selection') + if ~isempty(EEG.roi.roi_selection) + labels = labels(cell2mat(EEG.roi.roi_selection)); + end end end @@ -570,8 +574,10 @@ roi_loc = strrep(roi_loc, 'R', ''); % remove regions that were not selected - if ~isempty(EEG.roi.roi_selection) - roi_loc = roi_loc(cell2mat(EEG.roi.roi_selection)); + if isfield(EEG.roi, 'roi_selection') + if ~isempty(EEG.roi.roi_selection) + roi_loc = roi_loc(cell2mat(EEG.roi.roi_selection)); + end end try diff --git a/test_pipes/test_roi_selection.m b/test_pipes/test_roi_selection.m index e2ddc30..b46567b 100644 --- a/test_pipes/test_roi_selection.m +++ b/test_pipes/test_roi_selection.m @@ -59,5 +59,8 @@ cmp_roi_selection(roi_selection, squeeze(EEG3.roi.MIM(10, :, :)), squeeze(EEG4.roi.MIM(10, :, :))) %% Plot connectivity +pop_roi_connectplot(EEG3, 'measure', 'mim', 'plotcortex', 'off', 'plotmatrix', 'on', 'freqrange', [8 13]); +pop_roi_connectplot(EEG3, 'measure', 'mim', 'plotcortex', 'off', 'plotmatrix', 'on', 'freqrange', [8 13], 'region', 'central'); % region option disabled + pop_roi_connectplot(EEG4, 'measure', 'mim', 'plotcortex', 'off', 'plotmatrix', 'on', 'freqrange', [8 13]); pop_roi_connectplot(EEG4, 'measure', 'mim', 'plotcortex', 'off', 'plotmatrix', 'on', 'freqrange', [8 13], 'region', 'central'); % region option disabled \ No newline at end of file From fb76b0b3a7c8a58d05cf6894f9eb237883132616 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 12:22:40 +0200 Subject: [PATCH 2/8] extend documentation --- roi_pac.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_pac.m b/roi_pac.m index 7284ffe..661543c 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -2,7 +2,7 @@ % Wrapper function for pac_bispec(). % % Usage: -% EEG = roi_pac(EEG, fcomb, bs_outopts); +% EEG = roi_pac(EEG, fcomb, bs_outopts, 'key', 'val'); % % Inputs: % EEG - EEGLAB dataset with ROI activity computed. From 58df343f5c403885f87db4ef7dadbbbc641989ec Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 12:32:30 +0200 Subject: [PATCH 3/8] cosmetic changes in docstring --- roi_pac.m | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/roi_pac.m b/roi_pac.m index 661543c..fe4d667 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -5,19 +5,19 @@ % EEG = roi_pac(EEG, fcomb, bs_outopts, 'key', 'val'); % % Inputs: -% EEG - EEGLAB dataset with ROI activity computed. -% '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. +% EEG - EEGLAB dataset with ROI activity computed. +% 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 all (set to EEG.roi.nROI). +% 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). % % Output: % EEG - EEG structure with EEG.roi field updated and now containing From 194355e1102de2c941822ed0a5ed3dd5aab9824b Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 15:16:13 +0200 Subject: [PATCH 4/8] another change in docstring --- roi_pac.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/roi_pac.m b/roi_pac.m index fe4d667..372f67a 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -2,7 +2,7 @@ % Wrapper function for pac_bispec(). % % Usage: -% EEG = roi_pac(EEG, fcomb, bs_outopts, 'key', 'val'); +% EEG = roi_pac(EEG, fcomb, bs_outopts); % % Inputs: % EEG - EEGLAB dataset with ROI activity computed. From f10366bdcac5eb8dae563cd1692fc323a0195abf Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 16:30:45 +0200 Subject: [PATCH 5/8] add error message for pac on snippets --- pop_roi_connect.m | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 106a34e..471c48f 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -48,6 +48,7 @@ % 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'. % % Output: % EEG - EEGLAB dataset with field 'roi' containing connectivity info. @@ -170,7 +171,8 @@ 'freqresolution' 'integer' { } 0; 'fcomb' 'struct' { } struct; 'bs_outopts' 'integer' { } 1; - 'roi_selection' 'cell' { } { } }, 'pop_roi_connect'); + 'roi_selection' 'cell' { } { }; + 'conn_stats' 'string' { } 'off'}, 'pop_roi_connect'); if ischar(g), error(g); end % process multiple datasets @@ -223,7 +225,8 @@ % compute connectivity over snippets n_conn_metrics = length(options{2}); % number of connectivity metrics conn_matrices_snips = {}; -if strcmpi(g.snippet, 'on') +if strcmpi(g.snippet, 'on') && isempty(intersect(g.methods, {'PAC'})) + 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); @@ -273,7 +276,11 @@ end if ~isempty(intersect(g.methods, {'PAC'})) - EEG = roi_pac(EEG, g.fcomb, g.bs_outopts, g.roi_selection); + 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 From 05dd1d1386aded26c4ba86a24a26733e6f57b486 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Fri, 16 Jun 2023 17:23:18 +0200 Subject: [PATCH 6/8] started adding MIM statistics --- libs/pellegrini/fp_cpsd_mt.m | 61 +++++++++++++++++ libs/pellegrini/fp_cpsd_mt_matlab2015.m | 61 +++++++++++++++++ libs/pellegrini/fp_cpsd_mt_matlab2019.m | 58 ++++++++++++++++ libs/pellegrini/fp_cpsd_welch.m | 26 +++++++ libs/pellegrini/fp_fc_shuffletest.m | 23 +++++++ libs/pellegrini/fp_npcs2inds.m | 19 ++++++ libs/pellegrini/fp_plot_fc_shuffletest.m | 36 ++++++++++ libs/pellegrini/fp_tsdata_to_cpsd.m | 87 ++++++++++++++++++++++++ libs/pellegrini/fp_unwrap_conn.m | 50 ++++++++++++++ libs/pellegrini/shuffle_MIM.m | 69 +++++++++++++++++++ pop_roi_connect.m | 8 ++- roi_connstats.m | 59 ++++++++++++++++ 12 files changed, 555 insertions(+), 2 deletions(-) create mode 100644 libs/pellegrini/fp_cpsd_mt.m create mode 100644 libs/pellegrini/fp_cpsd_mt_matlab2015.m create mode 100644 libs/pellegrini/fp_cpsd_mt_matlab2019.m create mode 100644 libs/pellegrini/fp_cpsd_welch.m create mode 100644 libs/pellegrini/fp_fc_shuffletest.m create mode 100644 libs/pellegrini/fp_npcs2inds.m create mode 100644 libs/pellegrini/fp_plot_fc_shuffletest.m create mode 100644 libs/pellegrini/fp_tsdata_to_cpsd.m create mode 100644 libs/pellegrini/fp_unwrap_conn.m create mode 100644 libs/pellegrini/shuffle_MIM.m create mode 100644 roi_connstats.m diff --git a/libs/pellegrini/fp_cpsd_mt.m b/libs/pellegrini/fp_cpsd_mt.m new file mode 100644 index 0000000..3973147 --- /dev/null +++ b/libs/pellegrini/fp_cpsd_mt.m @@ -0,0 +1,61 @@ +function S = fp_cpsd_mt(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +ind_pow = intersect(ind_1, ind_2); +nfft = 2*(h-1); + +n1 = numel(ind_1); +n2 = numel(ind_2); +n3 = numel(ind_pow); +S = complex(zeros(h,n1,n2)); +pow = zeros(h,n3); + +winstep = window-noverlap; +ntapers = size(taparray,2); + +% compute tapered periodogram with FFT + +for k = 1:nchunks + + XSEG1 = X1((1:window) + (k-1)*winstep,:); + XSEG2 = X2((1:window) + (k-1)*winstep,:); + + % compute periodogram + P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P1 = P1(1:h,:,:); + P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P2 = P2(1:h,:,:); + + % now make cross-products of them to fill cross-spectrum matrix + + for ii = 1:n1 + o = ind_1(ii); + for jj = 1:n2 + oo = ind_2(jj); + S(:,ii,jj) = S(:,ii,jj) + mean(P1(:,:,o) .* conj(P2(:,:,oo)),2); + end + end + + if ~isempty(ind_pow) + for pp = 1:n3 + u = ind_pow(pp); + pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); + end + end + +end + +S = S/nchunks; + +if ~isempty(ind_pow) + pow = pow/nchunks; + + for pp = 1:n3 + clear a b + a = find(ind_1 == ind_pow(pp)); + b = find(ind_2 == ind_pow(pp)); + S(:,a,b) = pow(:,pp); + end +end + diff --git a/libs/pellegrini/fp_cpsd_mt_matlab2015.m b/libs/pellegrini/fp_cpsd_mt_matlab2015.m new file mode 100644 index 0000000..3973147 --- /dev/null +++ b/libs/pellegrini/fp_cpsd_mt_matlab2015.m @@ -0,0 +1,61 @@ +function S = fp_cpsd_mt(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +ind_pow = intersect(ind_1, ind_2); +nfft = 2*(h-1); + +n1 = numel(ind_1); +n2 = numel(ind_2); +n3 = numel(ind_pow); +S = complex(zeros(h,n1,n2)); +pow = zeros(h,n3); + +winstep = window-noverlap; +ntapers = size(taparray,2); + +% compute tapered periodogram with FFT + +for k = 1:nchunks + + XSEG1 = X1((1:window) + (k-1)*winstep,:); + XSEG2 = X2((1:window) + (k-1)*winstep,:); + + % compute periodogram + P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P1 = P1(1:h,:,:); + P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P2 = P2(1:h,:,:); + + % now make cross-products of them to fill cross-spectrum matrix + + for ii = 1:n1 + o = ind_1(ii); + for jj = 1:n2 + oo = ind_2(jj); + S(:,ii,jj) = S(:,ii,jj) + mean(P1(:,:,o) .* conj(P2(:,:,oo)),2); + end + end + + if ~isempty(ind_pow) + for pp = 1:n3 + u = ind_pow(pp); + pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); + end + end + +end + +S = S/nchunks; + +if ~isempty(ind_pow) + pow = pow/nchunks; + + for pp = 1:n3 + clear a b + a = find(ind_1 == ind_pow(pp)); + b = find(ind_2 == ind_pow(pp)); + S(:,a,b) = pow(:,pp); + end +end + diff --git a/libs/pellegrini/fp_cpsd_mt_matlab2019.m b/libs/pellegrini/fp_cpsd_mt_matlab2019.m new file mode 100644 index 0000000..9552a73 --- /dev/null +++ b/libs/pellegrini/fp_cpsd_mt_matlab2019.m @@ -0,0 +1,58 @@ +function S = fp_cpsd_mt_matlab2019(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +ind_pow = intersect(ind_1, ind_2); +nfft = 2*(h-1); + +n1 = numel(ind_1); +n2 = numel(ind_2); +n3 = numel(ind_pow); +S = complex(zeros(h,n1,n2)); +pow = zeros(h,n3); + +winstep = window-noverlap; +ntapers = size(taparray,2); + +% compute tapered periodogram with FFT + +for k = 1:nchunks + + XSEG1 = X1((1:window) + (k-1)*winstep,:); + XSEG2 = X2((1:window) + (k-1)*winstep,:); + + % compute periodogram + P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P1 = P1(1:h,:,:); + P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); + P2 = P2(1:h,:,:); + + % now make cross-products of them to fill cross-spectrum matrix + + for ii = 1:n1 + o = ind_1(ii); + S(:,ii,:) = S(:,ii,:) + mean(repmat(P1(:,:,o), 1, 1, length(ind_2)) .* conj(P2(:,:,ind_2)),2); + end + + if ~isempty(ind_pow) + for pp = 1:n3 + u = ind_pow(pp); + pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); + end + end + +end + +S = S/nchunks; + +if ~isempty(ind_pow) + pow = pow/nchunks; + + for pp = 1:n3 + clear a b + a = find(ind_1 == ind_pow(pp)); + b = find(ind_2 == ind_pow(pp)); + S(:,a,b) = pow(:,pp); + end +end + diff --git a/libs/pellegrini/fp_cpsd_welch.m b/libs/pellegrini/fp_cpsd_welch.m new file mode 100644 index 0000000..a3b8096 --- /dev/null +++ b/libs/pellegrini/fp_cpsd_welch.m @@ -0,0 +1,26 @@ +function S = fp_cpsd_welch(X_1, X_2,ind_1,ind_2,h,window,noverlap) +%X_1 and X_2 can be the same - then the usual CS is +%calculated. Otherwise X_2 should contain data from another trial than +%X_1. ind_1 and ind_2 contain channels of interest. The output S will be of +%size length(ind_1) x length(ind_2) x nfreq. + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +ind_pow = intersect(ind_1, ind_2); +nfft = 2*(h-1); + +n1 = numel(ind_1); +n2 = numel(ind_2); +S = complex(zeros(n1,n2,h)); + +for ii = 1:n1 + o = ind_1(ii); + S(ii,:,:) = transpose(cpsd(X_1(:,o),X_2(:,ind_2),window,noverlap,nfft)); % cross-spectra + if ismember(o,ind_pow) + clear b + b = find(ind_2 == o); + S(ii,b,:) = transpose(cpsd(X_1(:,o),X_1(:,o),window,noverlap,nfft)); + end + +end + \ No newline at end of file diff --git a/libs/pellegrini/fp_fc_shuffletest.m b/libs/pellegrini/fp_fc_shuffletest.m new file mode 100644 index 0000000..072efcb --- /dev/null +++ b/libs/pellegrini/fp_fc_shuffletest.m @@ -0,0 +1,23 @@ +function fp_fc_shuffletest(isb) +%calculates true MIM and generates MIM null distribution +%isb = subject index +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +DIRIN = './Data_MI/'; + +%subjects with high performance classification +subs = [3 4 5 8 9 11 12 14 15 16 17 18 19 21 22 23 25 27 28 29 30 31 33 34 35 37]; +nshuf = 1001; %first shuffle = true MIM, then 1000 shuffles + +%% load preprocessed EEG +sub = ['vp' num2str(subs(isb))]; +load([DIRIN sub '.mat']) + +%% shuffle +%one shuffle ~ 1 min on local machine +%first shuffle is true MIM +npcs = repmat(3,1,nroi); +MIM_s = fp_shuffle_MIM(data,npcs,fs,filt1, nshuf); + +%% save +save([DIRIN 'RDE_shuf_' num2str(isb) '.mat'],'-v7.3') diff --git a/libs/pellegrini/fp_npcs2inds.m b/libs/pellegrini/fp_npcs2inds.m new file mode 100644 index 0000000..65c96ec --- /dev/null +++ b/libs/pellegrini/fp_npcs2inds.m @@ -0,0 +1,19 @@ +function [inds, PCA_inds] = fp_npcs2inds(npcs) + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +%% +beg_inds = cumsum([1 npcs(1:end-1)]); +end_inds = cumsum([npcs]); + +for iroi = 1:numel(npcs) + PCA_inds{iroi} = beg_inds(iroi):end_inds(iroi); +end + +inds = {}; ninds = 0; +for iroi = 1:numel(npcs) + for jroi = (iroi+1):numel(npcs) + inds{ninds+1} = {PCA_inds{iroi}, PCA_inds{jroi}}; + ninds = ninds + 1; + end +end \ No newline at end of file diff --git a/libs/pellegrini/fp_plot_fc_shuffletest.m b/libs/pellegrini/fp_plot_fc_shuffletest.m new file mode 100644 index 0000000..6e1fbbb --- /dev/null +++ b/libs/pellegrini/fp_plot_fc_shuffletest.m @@ -0,0 +1,36 @@ +function fp_plot_fc_shuffletest +%Function that generates p-values from the MIM null distribution and plots +%them +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +DIRIN = '~/Dropbox/Franziska/Data_MEG_Project/RDE_shuffletest/right_MI/'; +DIRFIG = '~/Desktop/'; + +nsub = 26; + +%% generate p-values by comparing true MIM to null distribution +for isb = 1:nsub + in = load([DIRIN 'RDE_shuf_' num2str(isb) '.mat']); + + %average over one region dimension to obtain netMIM + MIM_s = squeeze(mean(in.MIM_s,2)); + MIM_pn(:,isb) = sum(MIM_s(:,1)< MIM_s(:,2:end),2)./(size(in.MIM_s,3)-1); +end + +%% Use Stouffer's method to aggregate p-values +nroi = size(MIM_pn,1); +for iroi = 1:nroi + MIM_pn_s(iroi) = fp_stouffer(squeeze(MIM_pn(iroi,:))); +end + +%% Use FDR-correction for multiple comparison's correction +[p_fdr, ~] = fdr(MIM_pn_s, 0.05); +MIM_pn_s(MIM_pn_s>p_fdr)=1; + +%% Plot +load cm17; +load('bs_results.mat'); % load cortex +MIM_pn_s(MIM_pn_s==0) = 0.001; % 1/nshuf +data = -log10(MIM_pn_s); +allplots_cortex_BS(cortex_highres,data, [35 52], cm17a ,'-log(p)', 0.3,... + [DIRFIG 'netMIM_right']) \ No newline at end of file diff --git a/libs/pellegrini/fp_tsdata_to_cpsd.m b/libs/pellegrini/fp_tsdata_to_cpsd.m new file mode 100644 index 0000000..8ef4dc4 --- /dev/null +++ b/libs/pellegrini/fp_tsdata_to_cpsd.m @@ -0,0 +1,87 @@ +function S = fp_tsdata_to_cpsd(X,fres,method,ind_1, ind_2, id_trials_1, id_trials_2, window,noverlap,nw,ntapers) +% Estimate cross-power spectral density from time series data between +% channels ind_1 and channels ind_2. +% Shuffle id_trials_2 for surrogate images. + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +%% +[n_chan,n_times,n_trials] = size(X); +X = demean(X); +X = permute(X,[2 1 3]); % transpose row, col (works for single-trial data too) +nfft = 2*fres; + +if ~exist('window','var') + window = min(n_times,nfft); % according to Chronux ... by default Matlab 'pwelch' splits data into 8 overlapping segments +end +assert(window <= n_times,'window cannot be longer than data'); + +if ~exist('overlap','var') + noverlap = round(window/2); +end +assert(noverlap < window,'overlap must be shorter than window'); + +if nargin < 3 || isempty(method) + method = 'WELCH'; +end + +if ~isequal(sort(ind_1),sort(unique(ind_1))) || ~isequal(sort(ind_2), sort(unique(ind_2))) + error('ind_1 and ind_2 must be unique') +end + +if strcmpi(method,'MT') + + nchunks = floor(((n_times-noverlap)/(window-noverlap))); % FFT chunks per channel + + if nargin < 10 || isempty(nw) + nw = 3; + end + + if nargin < 11 || isempty(ntapers) + ntapers = 2*nw -1; + end + + tapers = dpss(window,nw,ntapers,'calc'); % Slepian sequences: tapers is a matrix of size window x ntapers + taparray = tapers(:,:,ones(1,n_chan)); + + S = 0; + for r = 1:n_trials % works for single-trial too + + %trial shuffling + x_original = X(:,:,id_trials_1(r)); + x_perm = X(:,:,id_trials_2(r)); + + clear s + s = fp_cpsd_mt_matlab2019(x_original,x_perm,ind_1, ind_2,fres+1,window,noverlap,nchunks,taparray); + S = S+s; + + end + S = permute(S,[2 3 1])/n_trials; + + +elseif strcmpi(method,'WELCH') + + S = 0; + + for r = 1:n_trials + %trial shuffling + x_original = X(:,:,id_trials_1(r)); + x_perm = X(:,:,id_trials_2(r)); + + clear s + s = fp_cpsd_welch(x_original,x_perm,ind_1,ind_2,fres+1,window,noverlap); + S = S + s; + end + S = S/n_trials; + +else + error('unknown method ''%s''',method); +end + + + + + + + + diff --git a/libs/pellegrini/fp_unwrap_conn.m b/libs/pellegrini/fp_unwrap_conn.m new file mode 100644 index 0000000..007ad05 --- /dev/null +++ b/libs/pellegrini/fp_unwrap_conn.m @@ -0,0 +1,50 @@ +function [MIM_, MIC_, GC_, DIFFGC_, iCOH_, aCOH_] = fp_unwrap_conn(conn,nroi,filt1,PCA_inds) + +% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + +%% +%initialize output variables +MIM_ = []; +MIC_ = []; +GC_=[]; +DIFFGC_ = []; +iCOH_ =[]; +aCOH_ = []; + +% filt1 +% filt1.band_inds + +iinds = 0; +for iroi = 1:nroi + for jroi = (iroi+1):nroi + iinds = iinds + 1; + + if isfield(conn,'MIM') + MIM_(iroi, jroi) = mean(conn.MIM(filt1.band_inds, iinds)); + MIM_(jroi,iroi) = MIM_(iroi,jroi); + end + if isfield(conn,'MIC') + MIC_(iroi, jroi) = mean(conn.MIC(filt1.band_inds, iinds)); + MIC_(jroi,iroi) = MIC_(iroi,jroi); + end + if isfield(conn,'TRGC') + DIFFGC_(iroi,jroi) = mean(squeeze(conn.TRGC(filt1.band_inds,iinds,1) - conn.TRGC(filt1.band_inds,iinds,2))); + DIFFGC_(jroi,iroi) = -DIFFGC_(iroi,jroi); + end + if isfield(conn,'GC') + GC_(iroi,jroi) = mean(squeeze(conn.GC(filt1.band_inds,iinds,1) - conn.GC(filt1.band_inds,iinds,2))); + GC_(jroi,iroi) = -GC_(iroi,jroi); + end + end +end + +if isfield(conn,'COH') + for iroi = 1:nroi + for jroi = 1:nroi + iCOH_(iroi, jroi) = mean(mean(mean(abs(imag(conn.COH(filt1.band_inds, PCA_inds{iroi}, PCA_inds{jroi}))),1), 2), 3); + iCOH_(jroi,iroi) = iCOH_(iroi,jroi); + aCOH_(iroi, jroi) = mean(mean(mean(abs(conn.COH(filt1.band_inds, PCA_inds{iroi}, PCA_inds{jroi})),1), 2), 3); + aCOH_(jroi,iroi) = aCOH_(iroi,jroi); + end + end +end \ No newline at end of file diff --git a/libs/pellegrini/shuffle_MIM.m b/libs/pellegrini/shuffle_MIM.m new file mode 100644 index 0000000..132b6ce --- /dev/null +++ b/libs/pellegrini/shuffle_MIM.m @@ -0,0 +1,69 @@ +function MIM_s = shuffle_MIM(data,npcs,fres, nshuf) + %data is chan x l_epo x trials + % Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe + + [nchan,l_epo,nepo]= size(data); + + [inds, PCA_inds] = fp_npcs2inds(npcs); + ninds = length(inds); + + warning('One iteration takes about 90 seconds.') + fprintf('Progress of %d:', nshuf); + for ishuf = 1:nshuf %one iteration takes ~90 sec on my local laptop + if mod(ishuf, 10) == 0 + fprintf('%d', ishuf); + elseif mod(ishuf, 2) == 0 + fprintf('.'); + end + + %shuffle trials + if ishuf == 1 + shuf_inds = 1:nepo; %true MIM values + else + shuf_inds = randperm(nepo); + end + + clear MIM2 CS COH2 + CS = fp_tsdata_to_cpsd(data, fres, 'WELCH', 1:nchan, 1:nchan,1:nepo,shuf_inds); + + for ifreq = 1:size(CS,3) + clear pow + pow = real(diag(CS(:,:,ifreq))); + COH2(:,:,ifreq) = CS(:,:,ifreq)./ sqrt(pow*pow'); + end + + % loop over sender/receiver combinations to compute time-reversed GC + for iind = 1:ninds + if ~isequal(inds{iind}{1}, inds{iind}{2}) + %ind configuration + subset = [inds{iind}{1} inds{iind}{2}]; + subinds = {1:length(inds{iind}{1}), length(inds{iind}{1}) + (1:length(inds{iind}{2}))}; + + %MIC and MIM + [~ , MIM2(:, iind)] = roi_mim2(COH2(subset, subset, :), subinds{1}, subinds{2}); + end + end + + % extract measures out of the conn struct + clear conn + conn.MIM = MIM2; + conn.inds = inds; + nroi = nchan/npcs(1); + MIM_s(:, :, :, ishuf) = get_connect_mat(conn.MIM, nroi, +1); + % [MIM_s(:,:,ishuf), ~, ~, ~, ~, ~] = fp_unwrap_conn(conn,nroi,filt1,PCA_inds); + end +end + +function measure = get_connect_mat( measureOri, nROI, signVal) +% create a ROI x ROI connectivity matrix, if needed +% TRGCmat(f, ii, jj) is net TRGC from jj to ii + measure = []; + iinds = 0; + for iroi = 1:nROI + for jroi = (iroi+1):nROI + iinds = iinds + 1; + measure(:, iroi, jroi) = signVal * measureOri(:, iinds); + measure(:, jroi, iroi) = measureOri(:, iinds); + end + end +end diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 471c48f..b3c0cbf 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -49,6 +49,7 @@ % '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. % % Output: % EEG - EEGLAB dataset with field 'roi' containing connectivity info. @@ -172,7 +173,8 @@ 'fcomb' 'struct' { } struct; 'bs_outopts' 'integer' { } 1; 'roi_selection' 'cell' { } { }; - 'conn_stats' 'string' { } 'off'}, 'pop_roi_connect'); + 'conn_stats' 'string' { } 'off'; ... + 'nshuf' 'integer' { } 1001}, 'pop_roi_connect'); if ischar(g), error(g); end % process multiple datasets @@ -225,7 +227,7 @@ % compute connectivity over snippets n_conn_metrics = length(options{2}); % number of connectivity metrics conn_matrices_snips = {}; -if strcmpi(g.snippet, 'on') && isempty(intersect(g.methods, {'PAC'})) +if strcmpi(g.snippet, 'on') && isempty(intersect(g.methods, {'PAC'})) && strcmpi(g.conn_stats, 'off') snippet_length = g.snip_length; % seconds snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet @@ -270,6 +272,8 @@ EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct end 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); diff --git a/roi_connstats.m b/roi_connstats.m new file mode 100644 index 0000000..436fb15 --- /dev/null +++ b/roi_connstats.m @@ -0,0 +1,59 @@ +% roi_connstats() - Compute connectivity between ROIs and run statistics. +% +% Usage: +% EEG = roi_connstats(EEG, fcomb, nshuf); +% +% Inputs: +% EEG - EEGLAB dataset with ROI activity computed. +% 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 +% nshuf' - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 101. +% 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). + +function EEG = roi_connstats(EEG, methods, nshuf, roi_selection) + + 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'); + else + data = EEG.roi.source_roi_data; + end + + if ~isempty(intersect(methods, {'COH'})) + warning("'COH' is not supported anymore and will be replaced with aCOH (coherence). " + ... + "Please double-check with the documentation if this is what you want.") + coh_idx = strcmpi(methods, 'COH'); + g.methods{coh_idx} = 'aCOH'; + end + + if length(methods) == 1 && ~strcmpi(methods{1}, 'MIM' ) + error('Statistics are only supported for MIM, more methods will be added.') + else + warning('Statistics are only supported for MIM, more methods will be added.') + end + + methodset1 = { 'MIM' }; % the remaining methods will be added eventually + tmpMethods1 = intersect(methods, methodset1); + if ~isempty(tmpMethods1) + npcs = repmat(EEG.roi.nPCA, 1, EEG.roi.nROI); + MIM_s = shuffle_MIM(data, npcs, EEG.srate, nshuf); % (nfreq, nROI, nROI, nshuf) + EEG.roi.MIM = MIM_s; + end +end \ No newline at end of file From 77ea4135a5ae9759b01a9d11cb76a17d015315d4 Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Wed, 28 Jun 2023 15:41:29 +0200 Subject: [PATCH 7/8] started adding function to compute p-values --- test_pipes/test_mim_stats.m | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 test_pipes/test_mim_stats.m diff --git a/test_pipes/test_mim_stats.m b/test_pipes/test_mim_stats.m new file mode 100644 index 0000000..a4e037e --- /dev/null +++ b/test_pipes/test_mim_stats.m @@ -0,0 +1,42 @@ +% Test shuffling algorithm that calculates true MIM and generates MIM null distribution. +%% Run pipeline +clear +eeglab + +eeglabp = fileparts(which('eeglab.m')); +EEG = pop_loadset('filename','eeglab_data_epochs_ica.set','filepath',fullfile(eeglabp, 'sample_data/')); +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); +eeglab redraw; + +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] ); + +EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'functions','supportfiles','head_modelColin27_5003_Standard-10-5-Cap339.mat'), ... + 'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1); + +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); + + From b560ca29914a54a5ee917ef21448b806cd85451e Mon Sep 17 00:00:00 2001 From: Hiyeri <48970646+Hiyeri@users.noreply.github.com> Date: Wed, 5 Jul 2023 16:11:29 +0200 Subject: [PATCH 8/8] bugfix in GUI --- pop_roi_connect.m | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index b3c0cbf..228240c 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -50,6 +50,8 @@ % 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. @@ -115,10 +117,12 @@ 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] [0.2 1 0.35 0.8] [0.2 1 0.35 0.8] }; + 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' 'Coherence' 'tag' 'coh' 'value' 0 } ... + {'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 } ... @@ -138,7 +142,7 @@ % check we have the same naccu methods = {}; if out.cs, methods = [ methods { 'CS' } ]; end - if out.coh, methods = [ methods { 'COH' } ]; 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