diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 328457c..32b6a7f 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -125,24 +125,28 @@ 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] }; + uigeom = { [1] [1.2 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] [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' 'Cross-spectrum (CS)' 'tag' 'cs' 'value' 1 } {} ... + {'style' 'checkbox' 'string' '(Complex-valued) Coherency (cCOH)' 'tag' 'ccoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Coherence (aCOH)' 'tag' 'acoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Imaginary Coherency (iCOH)' 'tag' 'icoh' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Weighted Phase Lag Index (wPLI)' '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' 'checkbox' 'string' 'Multivariate Interaction Measure (MIM)' 'tag' 'mim' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Maximized Imaginary Coherency (MIC)' 'tag' 'mic' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Phase-Amplitude Coupling (PAC)' 'tag' 'pac' 'value' 0 } ... + { 'style' 'checkbox' 'string' 'Time Delay Estimation (TDE)' 'tag' 'tde' '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' } {} }; + {} { 'style' 'text' 'string' 'Bootstrap if any (n)' } { 'style' 'edit' 'string' '' 'tag' 'naccu2' } {} ... + {} { 'style' 'text' 'string' 'Frequency combination in Hz (for PAC) [f1 f2]' } { 'style' 'edit' 'string' '' 'tag' 'fcomb' } {} ... + {} { 'style' 'text' 'string' 'Region selection by index (for TDE) [region1 region2]' } { 'style' 'edit' 'string' '' 'tag' 'tde_regions' } {} }; ... [result,~,~,out] = inputgui('geometry', uigeom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_roi_connect'')', 'title', 'pop_roiconnect - connectivity'); if isempty(result), return, end @@ -163,10 +167,19 @@ if out.trdtf, methods = [ methods { 'TRDTF' } ]; end if out.mim , methods = [ methods { 'MIM' } ]; end if out.mic, methods = [ methods { 'MIC' } ]; end + if out.pac, methods = [ methods { 'PAC' } ]; end + if out.tde, methods = [ methods { 'TDE' } ]; end options = { ... 'morder' str2num(out.morder) ... 'naccu' str2num(out.naccu2) ... - 'methods' methods }; + 'methods' methods ... + 'tde_regions' eval( [ '[' out.tde_regions ']' ] )}; + if ~isempty(eval( [ '[' out.fcomb ']' ] )) + out_fcomb = eval( [ '[' out.fcomb ']' ] ); + fcomb.low = out_fcomb(1); + fcomb.high = out_fcomb(2); + options = [options {'fcomb' fcomb}]; + end else options = varargin; end diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 86a2279..adef5b7 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -158,7 +158,7 @@ splot(end ).unit = 'aCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -169,7 +169,7 @@ splot(end ).unit = 'cCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -180,7 +180,7 @@ splot(end ).unit = 'iCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -191,7 +191,7 @@ splot(end ).unit = 'GC'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -202,7 +202,7 @@ splot(end ).unit = 'TRGC'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -213,7 +213,7 @@ splot(end ).unit = 'MIC'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = -1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -224,7 +224,7 @@ splot(end ).unit = 'MIM'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -235,18 +235,18 @@ splot(end ).unit = 'PAC'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end if isfield(EEG.roi, 'PAC') - splot(end+1).label = 'ROI to ROI Phase-amplitude coupling'; + splot(end+1).label = 'ROI to ROI Phase-amplitude coupling (antisymmetrized)'; splot(end ).labelshort = 'Phase-amplitude coupling'; splot(end ).acronym = 'PAC_anti'; % PAC based on antisymmetrized bicoherence splot(end ).unit = 'PAC'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = 0; + splot(end ).psd = -1; splot(end ).plot3d = plot3dFlag; end @@ -262,7 +262,7 @@ end if isfield(EEG.roi, 'TDE') - splot(end+1).label = 'ROI to ROI Time-delay estimation'; + splot(end+1).label = 'ROI to ROI Time-delay estimation (antisymmetrized)'; splot(end ).labelshort = 'Time-delay estimation'; splot(end ).acronym = 'TDE_anti'; % TDE based on antisymmetrized bispectrum splot(end ).unit = 'TDE'; % not used yet @@ -508,6 +508,7 @@ error('PAC (original bicoherence) cannot be plotted, field is missing.') end cortexPlot = mean(matrix, 2); + titleStr = sprintf('f1 = %1.1f Hz, f2 = %1.1f Hz', S.PAC.fcomb.low, S.PAC.fcomb.high); case {'pac_anti'} if isfield(S.PAC, 'b_anti_norm') @@ -518,6 +519,7 @@ error('PAC (antisymmetrized bicoherence) cannot be plotted, field is missing.') end cortexPlot = mean(matrix, 2); + titleStr = sprintf('f1 = %1.1f Hz, f2 = %1.1f Hz', S.PAC.fcomb.low, S.PAC.fcomb.high); case {'tde' 'tde_anti'} if strcmpi(g.measure, 'tde') @@ -591,7 +593,7 @@ warning('Butterfly plots (frequency x connectivity) cannot be computed for PAC because frequencies have already been specified for the computation.') else figure; plot(EEG.roi.freqs, butterflyplot, 'LineWidth', 1) - h = title([ 'ROI to ROI ' replace_underscores(g.measure) ' (' titleStr ')' ]); + h = title([ 'ROI to ROI ' replace_underscores(g.measure)]); set(h, 'fontsize', 16); xlabel('Frequency (Hz)') ylabel([replace_underscores(g.measure) ' (a.u.)']) @@ -970,8 +972,8 @@ function plot_tde(T, shift, region_X, region_Y, method) [~, peak_idx] = max(T); est_delay = shift(peak_idx); % in Hz - figure; plot(shift, T, 'LineWidth', 1) - xline(est_delay, '--r') + figure; plot(shift, T, 'black', 'LineWidth', 1) + xline(est_delay, '--r', 'LineWidth', 1) xlabel('Time (s)') ylabel('a.u.') h = title(sprintf('%s -> %s TDE | Method %d', region_X, region_Y, method)); diff --git a/pop_roi_statsplot.m b/pop_roi_statsplot.m index 865137b..ec7d806 100644 --- a/pop_roi_statsplot.m +++ b/pop_roi_statsplot.m @@ -17,11 +17,14 @@ % 'TRDTF' : Time-reversed directed transfer entropy % 'MIM' : Multivariate Interaction Measure for each ROI % 'MIC' : Maximized Imaginary Coherency for each ROI +% 'PAC' : Phase Amplitude Coupling for each ROI % 'freqrange' - [min max] frequency range or [integer] single frequency in Hz. Default is to plot broadband power. % 'alpha' - [integer] Significance level. Default is 0.05. +% 'bispec' - ['b_anti'|'b_orig'] Option to compute antisymmetric or original bispectrum. % % Author: Franziska Pellegrini, franziska.pellegrini@charite.de % Tien Dung Nguyen, tien-dung.nguyen@charite.de +% Zixuan Liu, zixuan.liu@campus.tu-berlin.de function EEG = pop_roi_statsplot(EEG, varargin) @@ -39,36 +42,51 @@ g = finputcheck(varargin, { 'measure' 'string' { } ''; 'freqrange' 'real' { } []; ... - 'alpha' 'integer' { } 0.05}, 'pop_roi_statsplot'); + 'alpha' 'integer' { } 0.05; ... + 'bispec' 'string' {'b_orig', 'b_anti'} 'b_anti'}, 'pop_roi_statsplot'); if ischar(g), error(g); end - S = EEG.roi; + % check if measure is defined. if isempty(g.measure) error('You must define a measure to plot'); end - - % extract frequency indices - if ~isempty(g.freqrange) - if length(g.freqrange) == 1 - frq_inds = find(S.freqs == g.freqrange(1)); - title = sprintf('%1.1f Hz', g.freqrange(1)); + + % adjust based on measure, PAC has one less dimension. + if strcmp(g.measure, 'PAC') % check if measure is PAC + % for PAC, check the bispectrum parameter + if isfield(EEG.roi.(g.measure), g.bispec) + matrix = EEG.roi.(g.measure).(g.bispec); % use specified bispectrum field else - 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)); + error(['The specified bispectrum field (' g.bispec ') does not exist in EEG.roi.']); 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 + % if measure is not PAC, use the EEG.roi + S = EEG.roi; + + % extract frequency indices + if ~isempty(g.freqrange) + if length(g.freqrange) == 1 + frq_inds = find(S.freqs == g.freqrange(1)); + title = sprintf('%1.1f Hz', g.freqrange(1)); + else + 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)); + end + else + frq_inds = 1:length(S.freqs); + title = 'broadband'; + end - % average over one dimension to obtain net FC, then generate p-values by comparing the true FC (first shuffle) to null distribution + % 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 + 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); diff --git a/roi_pac.m b/roi_pac.m index 43ad1a9..a1a372b 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -84,6 +84,7 @@ [b_orig, b_anti, b_orig_norm, b_anti_norm] = data2bs_pac(data, params); % options which bispectral tensors to store + EEG.roi.PAC.fcomb = fcomb; switch bs_outopts case 2 EEG.roi.PAC.b_orig = b_orig; diff --git a/test_pipes/test_pac.m b/test_pipes/test_pac.m index 9ab4874..e543f35 100644 --- a/test_pipes/test_pac.m +++ b/test_pipes/test_pac.m @@ -85,3 +85,5 @@ pop_roi_connectplot(EEG1, 'measure', 'MIM', 'plotmatrix', 'on'); pop_roi_connectplot(EEG1, 'measure', 'aCOH', 'plotmatrix', 'on'); +% Statistic test plot +pop_roi_statsplot(EEG2, 'measure', 'PAC', 'bispec', 'b_anti');