diff --git a/BANetworkROIs_v4.txt b/BANetworkROIs_v4.txt new file mode 100644 index 0000000..e180dca --- /dev/null +++ b/BANetworkROIs_v4.txt @@ -0,0 +1,19 @@ +Cingulo_opercular_network DNM Salience_network Executive_network Working_memory_network Dorsal_attention_network +13L 23L 13L 8L 6R 17R +42L 31L 24L 9L 8R 18R +8L 39L 32L 10L 9R 19R +9L 24L 33L 46L 7R 7R +10L 25L 13R 3L 39R 8R +ThalamusL 32L 24R 5L 40R 17L +13R 33L 32R 7L 6L 18L +42R 23R 33R 39L 8L 19L +8R 31R 40L 9L 7L +9R 39R 8R 7L 8L +10R 24R 9R 39L +ThaRamusR 25R 10R 40L +39R 32R 46R +40R 33R 3R +39L 5R +40L 7R + 39R + 40R diff --git a/NGNetworkROIs_to_BA.m b/NGNetworkROIs_to_BA.m new file mode 100644 index 0000000..f2c0a72 --- /dev/null +++ b/NGNetworkROIs_to_BA.m @@ -0,0 +1,40 @@ +clear +allAreas = readtable('NGNetworkROIs_area_definition_v2.txt', 'delimiter', char(9)); +networks = readtable('NGNetworkROIs_v4.txt', 'delimiter', char(9)); +networksNew = networks; +networksNew(2:end,:) = []; + +networkNames = fieldnames(networks); +allAreasNames = fieldnames(allAreas); +allAreasNames = allAreasNames(1:end-3); +for iNet = 1:length(networkNames)-3 + + count = 1; + areas = networks.(networkNames{iNet}); + addAreaList = {}; + + % get the list of areas + for iArea = 1:length(areas) + areaNameTmp = areas{iArea}; + + if ~isempty(areaNameTmp) + if areaNameTmp(end) ~= 'L' && areaNameTmp(end) ~= 'R' + areaNameTmp1 = [ areaNameTmp '_R' ]; + areaNameTmp2 = [ areaNameTmp '_L' ]; + colPos1 = strmatch(areaNameTmp1, allAreasNames, 'exact'); + colPos2 = strmatch(areaNameTmp2, allAreasNames, 'exact'); + addAreaList = { addAreaList{:} allAreas{:,colPos1}{:} allAreas{:,colPos2}{:} }; + else + colPos = strmatch(areaNameTmp, allAreasNames, 'exact'); + addAreaList = { addAreaList{:} allAreas{:,colPos}{:} }; + end + end + end + + % write the list of areas + addAreaList(cellfun(@isempty, addAreaList)) = []; + for iArea = 1:length(addAreaList) + networksNew(iArea,iNet) = { addAreaList{iArea} }; + end +end +writetable(networksNew, 'BANetworkROIs_v4.txt', 'delimiter', char(9)); diff --git a/libs/haufe/data2cs_event.m b/libs/haufe/data2cs_event.m index ba36475..7ea4364 100644 --- a/libs/haufe/data2cs_event.m +++ b/libs/haufe/data2cs_event.m @@ -36,8 +36,6 @@ para=[]; end -maxfreqbin=min([maxfreqbin,floor(segleng/2)+1]); - segave=0; mydetrend=0; proj=[]; @@ -57,6 +55,11 @@ desired_nfreq = para.freqresolution; end +if desired_nfreq == 0 + maxfreqbin = min([maxfreqbin,floor(segleng/2)+1]); +else + maxfreqbin = desired_nfreq + 1; +end fres = segleng/2; [ndum,npat]=size(proj); @@ -168,7 +171,7 @@ end % compute wPLI -wpli = wpli_numer ./ wpli_denom; +wpli = abs(wpli_numer) ./ wpli_denom; wpli = permute(wpli, [3, 1, 2]); for f=1:maxfreqbin diff --git a/pac_bispec.m b/pac_bispec.m index 6ee0c71..63054dd 100644 --- a/pac_bispec.m +++ b/pac_bispec.m @@ -3,7 +3,7 @@ % are shown in Zandvoort and Nolte, 2021. % % Inputs: -% data - nchan x ntimepoints x ntrials ROI data, requires prior source reconstruction +% data - (nchan x ntimepoints x ntrials) or (nchan x ntimepoints * ntrials) ROI data, requires prior source reconstruction % params - Set containing parameters for bispectrum computation % % Outputs: @@ -12,7 +12,8 @@ % b_orig_norm - ROI x ROI bicoherence (normalized by threenorm) % b_anti_norm - ROI x ROI antisymmetrized bicoherence (normalized by threenorm) % -% Copyright (C) Franziska Pellegrini, franziska.pellegrini@charite.de +% Copyright (C) Franziska Pellegrini, franziska.pellegrini@charite.de, +% Tien Dung Nguyen, tien-dung.nguyen@charite.de function [b_orig, b_anti, b_orig_norm,b_anti_norm] = pac_bispec(data, params) @@ -38,7 +39,11 @@ [m, n] = ndgrid(frqs_low, frqs_high); frqs_combs = [m(:),n(:)]; n_combs = size(frqs_combs, 1); -warning('PAC is going to be estimated on %d frequency pair(s).', n_combs); +if n_combs > 50 + % according to our test simulations, the computation time scales linearly with the number of frequency pairs times 2, assuming no other ongoing CPU-heavy processes + time_est = 2 * n_combs; + warning('PAC is going to be estimated on %d frequency pair(s). Estimated time: %d seconds', n_combs, time_est); +end freqinds_low = zeros(n_combs, 2); freqinds_up = zeros(n_combs, 2); @@ -70,7 +75,7 @@ [bs_low,~] = data2bs_event(X(:,:)', segleng, segshift, epleng, freqinds_low); biv_orig_low = ([abs(bs_low(1, 2, 2, :)) abs(bs_low(2, 1, 1, :))]); xx = bs_low - permute(bs_low, [2 1 3, 4]); - biv_anti_low =([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))]); + biv_anti_low = ([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))]); % normalized by threenorm [RTP_low,~] = data2bs_threenorm(X(:,:)', segleng, segshift, epleng, freqinds_low); diff --git a/pop_roi_activity.m b/pop_roi_activity.m index e019314..10c6537 100644 --- a/pop_roi_activity.m +++ b/pop_roi_activity.m @@ -12,14 +12,14 @@ % 'sourcemodel' - [string] source model file % % Optional inputs: -% 'epochlen' - [float] epoch length (default is 2 seconds). ROIconnect -% has not been tested with other epoch lenghts. -% '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. Default is 'off'. -% 'fooof_frange' - [''] FOOOF fitting range. Default is [1 30] like in the -% example. +% 'epochlen' - [float] epoch length (default is 2 seconds). ROIconnect +% has not been tested with other epoch lenghts. +% '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 (the method 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 % specified, the signal is zero padded accordingly. % Default is 0 (means no padding). @@ -258,8 +258,10 @@ EEG = eeg_regepochs(EEG, recurrence, [0 2]); end +chansel = {}; if isempty(g.leadfield) g.leadfield = EEG.dipfit.sourcemodel; + chansel = EEG.dipfit.sourcemodel.label; end if isstruct(g.leadfield) && isfield(g.leadfield, 'file') sourceModelFile = g.leadfield.file; @@ -276,7 +278,7 @@ EEG = roi_activity(EEG, '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, moreargs{:}); + 'freqresolution', g.freqresolution, 'sourcemodelatlas', g.atlas, 'chansel', chansel, moreargs{:}); if nargout > 1 for iOption = 1:2:length(options) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 91fe18f..98018a5 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -251,7 +251,11 @@ if strcmpi(g.fcsave_format, 'all_snips') EEG.roi.(fc_name) = reshaped; else - mean_conn = squeeze(mean(reshaped, 1)); + 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 end diff --git a/roi_activity.m b/roi_activity.m index 35c8387..6b72bf2 100644 --- a/roi_activity.m +++ b/roi_activity.m @@ -32,11 +32,12 @@ % 'roiactivity' - ['on'|'off'] compute ROI activity. Default is on. If % you just need voxel activity, you can set this option to % ' off'. +% 'chansel' - [cell array of string] channel selection. Default is all. % 'exportvoxact' - ['on'|'off'] export voxel activity in EEG.roi.source_voxel_data % These array are huge, so the default is 'off'. -% 'fooof' - ['on'|'off'] enable FOOOF analysis. Default is 'off'. -% 'fooof_frange' - [ ] FOOOF fitting range. Default is [1 30] like in the -% example. +% 'fooof' - ['on'|'off'] enable FOOOF analysis (the method 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 % specified, the signal is zero padded accordingly. % Default is 0 (means no padding). @@ -94,6 +95,7 @@ 'model' 'string' { 'eLoretaFieldtrip' 'lcmvFieldtrip' 'eLoreta' 'lcmv' } 'lcmv'; 'nPCA' 'integer' { } 3; 'downsample' 'integer' { } 1; + 'chansel' 'cell' { } {}; 'roiactivity' 'string' { 'on' 'off' } 'on'; 'channelpower' 'string' { 'on' 'off' } 'off'; 'exportvoxact' 'string' { 'on' 'off' } 'off'; @@ -217,8 +219,13 @@ if ~isequal(nvox, nvox2) error('There must be the same number of vertices/voxels in the leadfield and source model'); end -if ~isequal(size(leadfield,1), EEG.nbchan) - error('There must be the same number of channels in the leadfield and in the dataset'); +if isempty(g.chansel) + g.chansel = [1:EEG.nbchan]; +else + g.chansel = eeg_decodechan(EEG.chanlocs, g.chansel); +end +if ~isequal(size(leadfield,1), length(g.chansel)) + error('There must be the same number of channels in the leadfield and in the list of selected channels'); end fres = EEG.pnts/2; @@ -227,11 +234,12 @@ frqs = sfreqs(fres, EEG.srate); % common average reference transform -H = eye(EEG.nbchan) - ones(EEG.nbchan) ./ EEG.nbchan; +nbchan = length(g.chansel); +H = eye(nbchan) - ones(nbchan) ./ nbchan; % apply to data and leadfield -EEG.data = reshape(H*EEG.data(:, :), EEG.nbchan, EEG.pnts, EEG.trials); -leadfield = reshape(H*leadfield(:, :), EEG.nbchan, nvox, 3); +tmpdata = reshape(H*EEG.data(g.chansel, :), nbchan, EEG.pnts, EEG.trials); +leadfield = reshape(H*leadfield(:, :), nbchan, nvox, 3); %% source reconstruction if strcmpi(g.model, 'eLoreta') @@ -240,16 +248,16 @@ P_eloreta = mkfilt_eloreta_v2(leadfield, g.modelparams{:}); % project to source space - source_voxel_data = reshape(EEG.data(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3); + source_voxel_data = reshape(tmpdata(:, :)'*P_eloreta(:, :), EEG.pnts*EEG.trials, nvox, 3); elseif strcmpi(g.model, 'LCMV') - C = cov(EEG.data(:, :)'); + C = cov(tmpdata(:, :)'); if length(g.modelparams) == 1 lcmv_reg = g.modelparams{1}; end alpha = lcmv_reg*trace(C)/length(C); - Cr = C + alpha*eye(EEG.nbchan); + Cr = C + alpha*eye(nbchan); [~, P_eloreta] = lcmv(Cr, leadfield, struct('alpha', 0, 'onedim', 0)); - source_voxel_data = reshape(EEG.data(:, :)'*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 @@ -392,7 +400,7 @@ % get channel power for comparison if strcmpi(g.channelpower, 'on') - tmpdata = permute(EEG.data, [2 1 3]); % pnts trials channels + tmpdata = permute(EEG.data(g.chansel, :, :), [2 1 3]); % pnts trials channels tmpdata = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3)); [tmpWelch,ftmp] = pwelch(tmpdata, data_pnts, data_pnts/2, data_pnts, data_pnts/2); % ftmp should be equal frqs tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.nbchan, EEG.trials); diff --git a/roi_definenetwork.m b/roi_definenetwork.m index f092974..f0ab7c1 100644 --- a/roi_definenetwork.m +++ b/roi_definenetwork.m @@ -177,7 +177,7 @@ % define networks networks = []; colNames = fieldnames(roiTable); -allLabels = { EEG.roi.atlas.Scouts.Label }; +allLabels = lower({ EEG.roi.atlas.Scouts.Label }); for iCol = 1:size(roiTable,2) % scan columns networks(end+1).name = colNames{iCol}; @@ -189,23 +189,18 @@ for iRow = 1:size(roiTable,1) val = roiTable{iRow, iCol}{1}; if ~isempty(val) - indTmp1 = strmatch(val, allLabels, 'exact'); - indTmp2 = strmatch([ 'Brodmann area ' val], allLabels, 'exact'); + indTmp1 = strmatch(lower(val), allLabels, 'exact'); + indTmp2 = strmatch(lower([ 'Brodmann area ' val]), allLabels, 'exact'); indTmp = [ indTmp1 indTmp2 ]; - if length(indTmp) == 1 + if isempty(indTmp) if strcmpi(g.ignoremissing, 'off') error('Area %s not found', val); else fprintf('Area %s not found\n', val); - if length(indTmp) > 1 indTmp = indTmp(1); end - end - else - if strcmpi(g.ignoremissing, 'off') - error('Area %s not found', val); - else - fprintf('Area %s duplicate, using the first one\n', val); - if length(indTmp) > 1 indTmp = indTmp(1); end end + elseif length(indTmp) > 1 + fprintf('Area %s duplicate, using the first one\n', val); + indTmp = indTmp(1); end inds = [ inds;indTmp ]; end diff --git a/roi_networkplot.m b/roi_networkplot.m index 02cd882..584e757 100644 --- a/roi_networkplot.m +++ b/roi_networkplot.m @@ -100,6 +100,7 @@ 'columns' 'integer' {} []; 'limits' 'float' {} []; 'plotmode' 'string' {'2D' '3D' 'both' 'none' } '2D'; + 'plotopt' 'cell' {} {}; 'filename' 'string' {} ''; 'threshold' 'float' {} 0.1; 'precomputed' 'struct' {} struct([]); @@ -167,6 +168,7 @@ else [EEG,networks,matrix] = roi_definenetwork(EEG, networks, 'addrois', g.addrois, 'connectmat', matrix, 'ignoremissing', 'on'); end +networks(cellfun(@(x)(length(x) < 2), { networks.ROI_inds })) = []; % Remove networks with less than 2 brain areas if strcmpi(g.subplots, 'on') if isempty(g.columns) @@ -237,7 +239,7 @@ % 2-D plot if strcmpi(g.plotmode, '2D') || strcmpi(g.plotmode, 'both') - plotconnectivity(networkMat(:,:), 'labels', labels, 'axis', gca, 'threshold', g.threshold(iNet), 'limits', g.limits); + plotconnectivity(networkMat(:,:), 'labels', labels, 'axis', gca, 'threshold', g.threshold(iNet), 'limits', g.limits, g.plotopt{:}); h = title(tmpTitle, 'interpreter', 'none'); pos = get(h, 'position'); set(h, 'position', pos + [0 0.1 0]); @@ -254,7 +256,7 @@ % 3-D plot if strcmpi(g.plotmode, '3D') || strcmpi(g.plotmode, 'both') - options = {'brainmovieopt' { 'moviename' '' } }; + options = {'brainmovieopt' { 'moviename' '' } g.plotopt{:} }; if ~strcmpi(g.subplots, 'on') && ~isempty(g.filename) tmpFileName = [ g.filename '_' networks(iNet).name '_3d' ]; options = { options{:} 'filename' tmpFileName }; diff --git a/test_pipes/test_pac.m b/test_pipes/test_pac.m index f3c109a..d48b2d2 100644 --- a/test_pipes/test_pac.m +++ b/test_pipes/test_pac.m @@ -33,7 +33,7 @@ %% Test bispectrum for frequency band inputs low = [4 8]; -high = [48 50]; +high = [9 10]; fcomb.low = low; fcomb.high = high;