Skip to content

Commit

Permalink
Merge pull request #55 from nguyen-td/feature
Browse files Browse the repository at this point in the history
Feature
  • Loading branch information
nguyen-td authored Jul 21, 2023
2 parents 522d849 + e130ce5 commit 1ccb9d4
Show file tree
Hide file tree
Showing 15 changed files with 305 additions and 126 deletions.
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,10 @@

# 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
> 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.

Expand Down
2 changes: 2 additions & 0 deletions libs/pellegrini/shuffle_MIM.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
3 changes: 2 additions & 1 deletion plotconnectivity.m
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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(:));
Expand Down
23 changes: 18 additions & 5 deletions pop_roi_activity.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -244,19 +249,21 @@
'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';
'freqresolution' 'integer' {} 0;
'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 = {};
Expand All @@ -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)
Expand Down
69 changes: 48 additions & 21 deletions pop_roi_connect.m
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,9 @@
% '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'.
% 'snip_length' - ['on'|'off'] Length of the snippets. Default is 60 seconds.
% '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.
% '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).
Expand Down Expand Up @@ -170,7 +172,9 @@
'naccu' 'integer' { } 0;
'methods' 'cell' { } { };
'snippet' 'string' { 'on', 'off' } 'off';
'nepochs' 'real' {} [ ];
'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';
'freqresolution' 'integer' { } 0;
Expand Down Expand Up @@ -237,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
Expand All @@ -246,34 +255,52 @@

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
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
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')
Expand Down
11 changes: 8 additions & 3 deletions pop_roi_connectplot.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
87 changes: 87 additions & 0 deletions pop_roi_statsplot.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
% 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
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(netFC, 2) - 1); % 1 / nshuf
data = -log10(FC_pn);
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
14 changes: 11 additions & 3 deletions roi_activity.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -328,11 +332,12 @@

% 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, 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

source_voxel_power = tmpWelch;

% fooof settings
if strcmpi(g.fooof, 'on')
f_range = g.fooof_frange; % freq range where 1/f should be fitted
Expand Down Expand Up @@ -377,12 +382,15 @@
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
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;
Expand Down
Loading

0 comments on commit 1ccb9d4

Please sign in to comment.