megconnectome pipeline v2.0

Contents.m

@@ -1,5 +1,5 @@
 % megconnectome
-% Version 1.0 www.humanconnectome.org 24-Feb-2014
+% Version 2.0 www.humanconnectome.org 10-Jul-2014
 
 % Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
 %

analysis_functions/hcp_ICA_RMEG_classification.m

@@ -54,11 +54,9 @@
 cfgin.saveres  = ft_getopt(options, 'saveres');
 cfgin.saveformat     = ft_getopt(options, 'saveformat');
 
-% if isempty(cfgin.dataset)  , cfgin.dataset = '0'; end
 if isempty(cfgin.ref_channels)  , cfgin.ref_channels = {'E31', 'E32'}; end
 if isempty(cfgin.bpfreq) , cfgin.bpfreq = iteration(1,1).comp(1,1).bandpass ; end
 if isempty(cfgin.bsfreq) , cfgin.bsfreq = iteration(1,1).comp(1,1).bandstop ; end
-% if isempty(cfgin.grad) ,  sens=ft_read_sens(cfgin.dataset) ; cfgin.grad = sens ; end
 if isempty(cfgin.showbestit), cfgin.showbestit = 'no'; end
 if isempty(cfgin.plottype)  , cfgin.plottype = 'summary'; end
 if isempty(cfgin.store_bestcomp)    , cfgin.store_bestcomp = 'yes';             end
@@ -70,11 +68,6 @@
 
 if isfield(cfgin,'ref_channels')
 
-    %     if ischar(ref_data)
-    %         options=ft_setopt(options,'channels',cfgin.ref_channels)
-    %     ref_data = hcp_ICA_preprocessing(ref_data, options);
-    %     end
-
     nref= size(ref_data.trial{1},1);
     ntrl= size(ref_data.trial,2)
 
@@ -100,7 +93,7 @@
     end
     imgname=[subject '_icaclass_refch'];
     hcp_write_figure(imgname, h);
-    clear h; % FIXME should this be clear or close?
+    clear h;
 
     %   ----- Time Course of Power of Electric Channels -----------------------
 
@@ -133,7 +126,6 @@
     selec = sqrt(freq_ref_data.powspctrm);
 end
 
-%%   Performing IC CLASSIFICATION
 
 n_iter=size(iteration,2);
 for j=1:n_iter
@@ -162,14 +154,6 @@
     %----> IC POWER TIME COURSE <----
     %--------------------------------
 
-    %     cfg = [];
-    %
-    %     cfg.bpfilter      = 'yes'; %'no' or 'yes'  bandpass filter (default = 'no')
-    %     cfg.bpfreq        = cfgin.bpfreq;
-    %     cfg.hilbert='abs';
-    %
-    %     pow_data                 = ft_preprocessing(cfg,comp_iter);
-    %
     pow_data=comp_iter.pow_data;
     temp_struc(j).pow_data=pow_data;
     pIC = cell2mat(pow_data.trial).^2;
@@ -187,7 +171,6 @@
 
         IC=cell2mat(comp_iter.trial);
         if isfield(cfgin,'ref_channels')
-            %             elec=ref_data.trial{1};
             if nref > 0
 
                 %   ----- Correlation between the ICs and electric channels
@@ -218,31 +201,13 @@
 
         %   ----- Fit with 1/f spectrum, flat spectrum "quantification", IC time kurtosis
 
-
-        %!FROM GIORGOS:
-        % Replaced below the fittype function with hcp_fitspectrum, so the compiled
-        % version can run on WashU, according to instructions from the help
-        % of the function i.e.:
-        % ------------------------------------------
-        % Example:
-        %   x = (1:100)';
-        %   y = 2./x + 3 + 0.1*rand(size(x));
-        %   [fitx, goodness] = hcp_fitspectrum(x, y)
-        %
-        % This replaces the following use of the fit() function from
-        % the Mathworks curve fitting toolbox
-        %   g = fittype('abs(a)*1/x + abs(b)')
-        %   [fitx, goodness] = fit(x, y, g)
-        %--------------------------------------------
-        %g = fittype('abs(a)*1/x + b'); % COMMENTED OUT BY GIORGOS
         frq = [cfgin.bpfreq(1,1) 13];
         frq2 = [2 78];
         vec = find(F > frq(1) & F < frq(2));
         vec2 = find(F > frq2(1) & F < frq2(2));
         for ix = 1:Nc
             sspettro = mspettro(ix,:);
             sspettro_fit=sspettro./sspettro(vec(1,1));
-            %[fitx,goodness] = fit(F(vec)',sspettro_fit(vec)',g); % COMMENTED OUT BY GIORGOS
             [fitx,goodness] = hcp_fitspectrum(F(vec)',sspettro_fit(vec)');
             if(fitx.a>0)
                 G(ix) = goodness.rsquare;
@@ -261,8 +226,7 @@
         thres_d = 0.5;  % PSD 1/f
         thres_e = 2.02; % Spectrum Flat
         thres_f = 15;   % IC Time Kurtosis
-        %     thres_g = 0.1;
-        %     thres_h = 15;
+
 
         %----------------------------------------------------------
         %       classification: 0=artifact   1=brain signal

analysis_functions/hcp_ICA_freq.m

@@ -28,32 +28,47 @@
 %   saveres:    'yes' or 'no' save the resutls in a file (default = 'no')
 %   saveformat: image file format (default 'fig')
 %   grad:       fieldtrip gradiometer structure
+%
+% See also HCP_ICA_PLOT HCP_ICA_PLOTCLASSIFICATION HCP_ICA_RMEG_CLASSIFICATION
+
+% Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
+%
+% This file is part of megconnectome.
+%
+% megconnectome is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% megconnectome is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with megconnectome.  If not, see <http://www.gnu.org/licenses/>.
 
 switch nargin
     case 3
         % reconstruct the component time courses from the data and the unmixing
         % matrix
-        %     comp.time     = datain.time;
         cfg           = [];
         cfg.unmixing  = comp.unmixing;
         cfg.topolabel = comp.topolabel;
         cfg.order     = comp.order;
         comp          = ft_componentanalysis(cfg, datain);
         comp.order    = cfg.order;
-        %     comp.trial{1} = comp.unmixing*datain.trial{1};
     case 2
         % check whether IC time courses are in the first input
         if ~isfield(comp,'trial') , error('IC time course not defined, data matrix required'); end
     case 1
         error('too few input arguments');
-        % FIXME in principle this should work because then all options should
-        % be set to default within the function
         if ~isfield(comp,'trial') , error('IC time course not defined, data matrix required'); end
     otherwise
         error('too many input arguments');
 end
 
-cfgin.doplot   = ft_getopt(options, 'doplot', 'no'); % FIXME I think the plotting should not be done within this function
+cfgin.doplot   = ft_getopt(options, 'doplot', 'no');
 cfgin.grad     = ft_getopt(options, 'grad');
 
 if ~isempty(cfgin.grad),
@@ -76,10 +91,7 @@
 freq_comp    = ft_freqanalysis(cfg,comp_seg);
 comp.freq_comp=freq_comp;
 
-%--------------------------------
 %----> IC POWER TIME COURSE <----
-%--------------------------------
-
 cfg               = [];
 cfg.bpfilter      = 'no'; %'no' or 'yes'  bandpass filter (default = 'no')
 cfg.bpfreq        = [1 150];
@@ -94,5 +106,5 @@
 comp.pow_data = pow_data;
 
 if strcmp(cfgin.doplot,'yes')
-    hcp_ICA_plot(comp,options);
+    hcp_ICA_plot(comp,options); % plot results
 end

analysis_functions/hcp_ICA_plot.m

@@ -7,7 +7,7 @@ function hcp_ICA_plot(comp_freq, options)
 % two different layouts (see cfg.plottype)
 %
 % Use as
-%   hcp_ICA_freq(comp_freq, options, datain)
+%   hcp_ICA_plot(comp_freq, options, datain)
 % where the input comp_freq structure should be obtained from FT_COMPONENTANALYSIS and
 % the datain structure is a fieldtrip data structure contining the channel time courses
 % used as imput of the FT_COMPONENTANALYSIS .
@@ -20,13 +20,32 @@ function hcp_ICA_plot(comp_freq, options)
 %   saveres            : 'yes' or 'no' save the resutls in a file (default = 'no')
 %   saveformat         : image file format (default 'fig')
 %   grad               : fieldtrip gradiometer structure
+%
+% See also HCP_ICA_FREQ HCP_ICA_PLOTCLASSIFICATION HCP_ICA_RMEG_CLASSIFICATION
+
+% Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
+%
+% This file is part of megconnectome.
+%
+% megconnectome is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% megconnectome is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with megconnectome.  If not, see <http://www.gnu.org/licenses/>.
 
 Nc = size(comp_freq.topo,2);
 
 cfgin.plottype   = ft_getopt(options, 'plottype',   'components'); %summary or components
-cfgin.saveres    = ft_getopt(options, 'saveres');
+cfgin.saveres    = ft_getopt(options, 'saveres'); 
 cfgin.saveformat = ft_getopt(options, 'saveformat', 'fig');
-cfgin.fileout    = ft_getopt(options, 'fileout');
+cfgin.fileout    = ft_getopt(options, 'fileout'); 
 cfgin.parameter  = ft_getopt(options, 'parameter',  'topo');
 cfgin.component  = ft_getopt(options, 'component',  1:Nc);
 cfgin.comment    = ft_getopt(options, 'comment',    'no');
@@ -125,7 +144,6 @@ function hcp_ICA_plot(comp_freq, options)
         if(~isempty(cfgin.posi))
             f = figure;
             set(f, 'visible', fvis,'on','Position', cfgin.posi);
-            %      set(f, 'Position', cfgin.posi)
         else
             f = figure;
             set(f, 'visible', fvis);

analysis_functions/hcp_ICA_plotclassification.m

@@ -1,5 +1,19 @@
 function hcp_ICA_plotclassification(comp,options_ICA,datain)
 
+% HCP_ICA_PLOTCLASSIFICATION allows the plots of Independent
+% Components automatically classified by HCP_ICA_RMEG_CLASSIFICATION function. This
+% function is able to classify Independent Components into Brain Activity
+% and ECG/EOG related artifacts provided that ECG and EOR reference
+% channels are used.
+%
+% Use as
+%   hcp_ICA_plotclassification(comp,options_ICA,datain)
+%
+%     where the input comp structure should be obtained from
+%     HCP_ICA_RMEG_CLASSIFICATION and the datain structure is a fieldtrip data
+%     structure contining the channel time courses used as imput of the
+%     FT_COMPONENTANALYSIS .
+%
 % Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
 %
 % This file is part of megconnectome.
@@ -21,7 +35,7 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
 bandpass  = ft_getopt(options_ICA, 'bandpass');
 bandstop  = ft_getopt(options_ICA, 'bandstop');
 grad      = ft_getopt(options_ICA, 'grad');
-modality      = ft_getopt(options_ICA, 'modality', 'MEG'); % GIORGOS
+modality      = ft_getopt(options_ICA, 'modality', 'MEG'); 
 
 if strcmp(modality,'MEG')
     layout='4D248.mat';
@@ -37,7 +51,7 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
 
 
 if(~isfield(comp,'pow_data'))
-    doplot='no'; % Summury Results cn be plotted also at this level without calling hcp_plotICA
+    doplot='no'; 
     options = {'doplot',doplot, 'grad', grad,'plottype','components'};
 
     [comp_freq]=hcp_ICA_freq(comp, options, datain);
@@ -46,7 +60,7 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
 end
 n_IC=size(comp.unmixing,1);
 
-frange=bandpass; %to be decided if in compo or other
+frange=bandpass;
 
 cfgin =[];
 cfgin.grad=grad;
@@ -68,7 +82,7 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
 leg={'OneOverF (>0.5)  ' ; 'Specflat (<2.0)  ' ; 'Kurtosis (>15)   ' ; 'elecSig (>0.1)   '; 'elecPow (>0.25)  '; 'elecSpe (>0.95)  '};
 thrs=[0.5;2.0;15;0.1;0.25;0.95];
 class=comp.class;
-% size_s=get(0,'ScreenSize');
+
 for ix=1:n_IC
     str(1,:)={num2str(abs(comp.class.spectrum_one_over_f(1,ix)),'%5.3f')};
     str(2,:)={num2str(abs(comp.class.spectrum_flat(1,ix)),'%5.1f')};
@@ -86,36 +100,23 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
     str_ic = 'ARTIFACT';
 
     if(find(class.brain_ic==ix)), str_ic='BRAIN'; end
-%%%   In the cluster version all this settings are not required    
-%     X = 29.7;                  %# A3 paper size
-%     Y = 14.35;                  %# A3 paper size
-%     xMargin = 1;               %# left/right margins from page borders
-%     yMargin = 1;               %# bottom/top margins from page borders
-%     xSize = X - 2*xMargin;     %# figure size on paper (widht & hieght)
-%     ySize = Y - 2*yMargin;     %# figure size on paper (widht & hieght)
+
     figure('Menubar','none');
 %     
     set(gca, 'XTickLabel',[], 'YTickLabel',[], ...
         'Units','normalized', 'Position',[0 0 1 1])
-%     
-%     %# figure size on screen (50% scaled, but same aspect ratio)
-%     set(gcf, 'Units','centimeters', 'Position',[5 5 xSize ySize])
 
 %# figure size printed on paper
     set(gcf, 'visible', 'off')
     set(gcf, 'PaperUnits','inches')
-%     set(gcf, 'PaperSize',[X Y])
-%     set(gcf, 'PaperPosition',[xMargin yMargin xSize ySize])
-%     set(gcf, 'PaperOrientation','portrait')
-%%%
+
     set(gcf, 'paperposition', [1 1 20 12]);
 
     subplot(2,3,1)
     cfgin.component=ix;
     cfgin.layout=layout;
     cfgin.colorbar='yes';
     ft_topoplotIC(cfgin, comp);
-%     xlabel(['IC ' num2str(ix)],'color','k');
     title('IC sensor map','FontSize',12);
 
     subplot(2,3,4)
@@ -130,8 +131,6 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
 
     str_icnum = {['ic ' num2str(ix) ' = ' str_ic]};
 
-%     text(1.2,0.5,leg,'Units','normalized','FontSize',16)
-%     text(2,0.5,str,'Units','normalized','FontSize',16,'color',strcolor)
     text(1.2 ,0.65,leg{1},'Units','normalized','FontSize',22,'color',strcolor(1))
     text(2,0.65, str{1},'Units','normalized','FontSize',22,'color',strcolor(1))
     text(1.2,0.59,leg{2},'Units','normalized','FontSize',22,'color',strcolor(2))
@@ -145,10 +144,7 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
     text(1.2,0.335,leg{6},'Units','normalized','FontSize',22,'color',strcolor(6))
     text(2,0.335,str{6},'Units','normalized','FontSize',22,'color',strcolor(6))
 
-%     text(1.5,0.2,str_icnum,'Units','normalized','FontSize',16)
-%     text(1.5,0.1,str_ic,'Units','normalized','FontSize',16,'color','r')
     text(1.2,0.8,str_icnum,'Units','normalized','FontSize',22,'color','b')
-%     text(1.8,0.8,str_ic,'Units','normalized','FontSize',16,'color','b')
 
     subplot(2,3,[2 3])
     timeic=cell2mat(comp.time);
@@ -173,7 +169,6 @@ function hcp_ICA_plotclassification(comp,options_ICA,datain)
     disp('SAVING');
     hcp_write_figure(imgname, gcf, 'resolution', 300);
     disp('SAVED');    
-%     clear h;
     close
 end