Merge pull request #39 from nguyen-td/feature

Feature

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		

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));

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

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, [email protected]
+% Copyright (C) Franziska Pellegrini, [email protected],
+%               Tien Dung Nguyen, [email protected]
 
 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);

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)

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

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);

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

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 };

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;