fmri_pls_analysis.m

绝对经典,老外制作的功能强大的matlab实现PLS_TOOBOX

function fmri_pls_analysis(varargin)
%
%  USAGE: 
%          fmri_pls_analysis
%    or    fmri_pls_analysis(SessionProfiles,ContrastFile, ...
%                              num_perm,num_boot,grp_analysis_flg,output_file)
%
%    Apply PLS on the fMRI data based on the session information saved in
%    the "sessionFile".  Assume session information, and st_datamat have 
%    been created.
%
%
%  INPUT:
%    SessionProfiles - a cell array, one element per group.  Each element 
%               in the array is another cell array contains the names of 
%               session profiles for the group.
%    ContrastFile - the contrast file to be used to generate the design matrix.
%                   Helmert matrix will be specified using the string of
%                   'HELMERT'.  If contrast is empty, use deviation from 
%                   grand mean for the contrast.
%    num_perm - number of permutations to be performed.
%    num_boot - number of bootstrap resampling to be performed.
%    grp_analysis_flg - flag indicates group analysis.
%    output_file - (optional) the name of output file
%
%                  
%  OUTPUT FILE:
%         - file stores the information of the PLS result.
%
%  NOTE:
%    To create session information, use 'pls_input_session_info'.
%    To create st_datamat, use 'fmri_gen_datamat', 'fmri_combine_coords',
%       and 'fmri_gen_st_datamat'.
%
%  Script needed: fmri_perm_test.m, fmri_deviation_perm_test.m
%
%   -- Created July 2001 by Wilkin Chau, Rotman Research Institute
%

  singledatamat = 0;		% init singledatamat to false

  if (nargin == 0),
     [SessionProfiles,ContrastFile,num_perm,group_analysis] = options_query;
  else
     SessionProfiles = varargin{1};
     ContrastFile = varargin{2};
     num_perm = varargin{3};
     num_boot = varargin{4};
     Clim = varargin{5};
     posthoc = varargin{6};
     save_datamat = varargin{7};
     group_analysis = varargin{8};
     cond_selection = varargin{9};
     behavname = varargin{10};
     behavdata = varargin{11};
     behavdata_lst = varargin{12};
     bscan = varargin{13};

     if (nargin > 13)
        output_file = varargin{14};
        for_batch = 1;
     else
        for_batch = 0;
     end;
  end;

  session_files_timestamp = SessionProfiles;
  datamat_files_timestamp = SessionProfiles;

  change_timestamp = 0;

  for i = 1:length(SessionProfiles)
     for j = 1:length(SessionProfiles{i})
        tmp = dir(SessionProfiles{i}{j});
        session_files_timestamp{i}{j} = tmp.date;

        load(SessionProfiles{i}{j},'session_info');
        datamat_prefix = session_info.datamat_prefix;

        if findstr('BfMRIsession.mat', SessionProfiles{i}{j})
           datamat_file = [datamat_prefix,'_BfMRIdatamat.mat'];
        else
           datamat_file = [datamat_prefix,'_fMRIdatamat.mat'];
        end


        warning off;
        load(datamat_file, 'singleprecision');
        warning on;

        if exist('singleprecision','var') & singleprecision
           singledatamat = 1;
        end


        tmp = dir(datamat_file);
        datamat_files_timestamp{i}{j} = tmp.date;

        if datenum(session_files_timestamp{i}{j}) > datenum(datamat_files_timestamp{i}{j})
           change_timestamp = 1;
        end
     end
  end

  if change_timestamp
     msg1 = ['One or more datamat files are older than their session files, '];
     msg2 = 'If you believe that the session files is just touched (e.g. due to copy) but not modified, you can click "Ignore All".';
     msg3 = 'Otherwise, please click "Stop", and re-generate the datamat file.';

     quest = questdlg({msg1 '' msg2 '' msg3 ''}, 'Choose','Proceed All','Stop','Stop');

     if strcmp(quest,'Stop')
        return;
     end
  end

  %  save results
  %
  if exist('output_file','var') & ~isempty(output_file)
    resultFile = output_file;
  else

    fn = SessionProfiles{1}{1};
    load(fn,'session_info');
    datamat_prefix = session_info.datamat_prefix;

    if findstr('BfMRIsession.mat', fn)
       [result_file,result_path] = ...
		uiputfile([datamat_prefix,'_BfMRIresult.mat'],'Saving PLS Result');
    else
       [result_file,result_path] = ...
		uiputfile([datamat_prefix,'_fMRIresult.mat'],'Saving PLS Result');
    end

    if isequal(result_file,0)			% Cancel was clicked
%%       if findstr('BfMRIsession.mat', fn)
%%          resultFile = 'BfMRIresult.mat';
%%       else
%%          resultFile = 'fMRIresult.mat';
%%       end
%
%       msg1 = ['WARNING: No file is saved.'];
%       msgbox(msg1,'Uncompleted');
%       resultFile = [];
%       disp('ERROR: Result file is not saved.');
       return;
    else
       resultFile = fullfile(result_path,result_file);
    end;
  end;


   v7 = version;
   if str2num(v7(1))<7
      singleanalysis = 0;
   else
      singleanalysis = 1;
   end

   pc = computer;
   if singleanalysis & ( strcmp(pc,'GLNXA64') | strcmp(pc,'GLNXI64') | strcmp(pc,'PCWIN64') )
      quest = questdlg({'We detected that you are running MATLAB on a 64-bit system. According to MATLAB Bug Report ID 268001, we have to convert data to double precision for Intel based system.' '' 'Is this Intel 64-bit machine?' ''}, 'Choose','No','Yes','Don''t know','Don''t know');

      if ~strcmp(quest,'No')
         singleanalysis = 0;
      end
   end;


  progress_hdl = ShowProgress('initialize');

  if isempty(ContrastFile),                         % none: use deviation
     ContrastMethod = 1;
  elseif strcmp('HELMERT',upper(ContrastFile))      % use Helmert matrix
     ContrastMethod = 2;
  elseif strcmp('BEHAV',upper(ContrastFile))        % behav
     ContrastMethod = 4;
  elseif strcmp('MULTIBLOCK',upper(ContrastFile))   % multiblock
     ContrastMethod = 5;
  else
     ContrastMethod = 3;                            % design with contrast file
  end;

  if ContrastMethod ~= 4 & ContrastMethod ~= 5
     [st_datamat,st_coords,st_dims,num_conditions,st_evt_list, ...
	st_win_size,st_voxel_size,st_origin,subj_group,behavdata,behavname, ...
		subj_name, cond_name, num_behav_subj, ...
		behavdata_lst, newdata_lst, num_subj_lst ] = ...
		concat_st_datamat(singleanalysis,SessionProfiles,progress_hdl,ContrastMethod, ...
		posthoc, cond_selection, group_analysis);
  else
     [st_datamat,st_coords,st_dims,num_conditions,st_evt_list, ...
	st_win_size,st_voxel_size,st_origin,subj_group, ...
		subj_name, cond_name, num_behav_subj, ...
		newdata_lst, num_subj_lst ] = ...
		concat_st_datamat2(singleanalysis,SessionProfiles,progress_hdl,ContrastMethod, ...
		posthoc, behavdata, cond_selection);
  end

  if isempty(st_datamat),
     return;
  end;

   if ~singleanalysis
      st_datamat = double(st_datamat);

      if ~isempty(newdata_lst)
         for g = 1:length(newdata_lst)
            newdata_lst{g} = double(newdata_lst{g});
         end
      end
   end

  if (group_analysis == 0)
     subj_group = [];		% for nongroup analysis
  end;

  %  get the contrast 
  %
  switch  (ContrastMethod) 
     case {1}
        USE_DEVIATION_PERM_TEST = 1;
	isbehav = 0;
        ContrastFile = 'NONE';
     case {2}
        USE_DEVIATION_PERM_TEST = 0;
	isbehav = 0;
        contrasts = rri_helmert_matrix(num_conditions);
     case {3}
        USE_DEVIATION_PERM_TEST = 0;
	isbehav = 0;
        if isnumeric(ContrastFile)
           contrasts = ContrastFile;
        else
           contrasts = load(ContrastFile);
        end
     case {4}
        USE_DEVIATION_PERM_TEST = 0;
	isbehav = 1;
        contrasts = behavdata;
     case {5}
        USE_DEVIATION_PERM_TEST = 0;
	isbehav = 2;
        contrasts = behavdata;
  end;

  create_ver = plsgui_vernum;

  %  start permutation PLS ...
  %
  perm_result = [];
  boot_result = [];
  if (USE_DEVIATION_PERM_TEST)		% mean

     if (num_boot > 0),
        num_rep = length(st_evt_list) / num_conditions;

%        boot_progress = rri_progress_ui('initialize');

        if isempty(subj_group)
           [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
              = rri_boot_check(num_rep, 1, num_boot, 1, ...
                for_batch);
%                boot_progress, for_batch);
        else
           [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
              = rri_boot_check(subj_group, num_conditions, num_boot, 1, ...
                for_batch);
%                boot_progress, for_batch);
        end;

        num_boot = new_num_boot;
     end;

     if (num_perm > 0) | (num_boot == 0),
        [brainlv,s,designlv,b_scores,d_scores,perm_result,lv_evt_list] = ...
           fmri_deviation_perm_test(st_datamat,num_conditions,st_evt_list, ...
                                   num_perm,subj_group);
     end;

     if (num_boot > 0),
        [brainlv2,s2,designlv2,b_scores2,d_scores2,boot_result,lv_evt_list2] = ...
           fmri_deviation_boot_test(st_datamat,num_conditions,st_evt_list, ...
                                   num_boot, subj_group, ...
         min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

        if num_perm == 0
           brainlv = brainlv2;
           s = s2;
           designlv = designlv2;
           b_scores = b_scores2;
           d_scores = d_scores2;
           lv_evt_list = lv_evt_list2;
           perm_result = [];
        end
     end;

     saved_info=['brainlv s designlv perm_result boot_result st_coords ', ...
  	         'st_dims lv_evt_list st_win_size st_voxel_size st_origin ', ...
		 'SessionProfiles ContrastFile b_scores d_scores ', ...
		 'subj_group num_conditions cond_name cond_selection ', ...
		 'num_subj_lst subj_name session_files_timestamp ', ...
		 'datamat_files_timestamp create_ver'];

     if save_datamat & ~isempty(brainlv)
        first = 1;
        last = 0;
%        grp_datamat = [];
        for g = 1:length(num_subj_lst)
           last = last + num_conditions*num_subj_lst(g);
           [tmp idx] = sort(st_evt_list(first:last));
           tmp = st_datamat(first:last,:);
%           grp_datamat = [grp_datamat; tmp(idx,:)];
           grp_datamat{g} = tmp(idx,:);
           first = last + 1;
        end;
        
        st_datamat = grp_datamat;
        saved_info = [saved_info, ' st_datamat'];
     end

  elseif isbehav == 2			% Multiblock Analysis

        ibehavdata_lst = behavdata_lst;

        if (num_boot > 0),
%           boot_progress = rri_progress_ui('initialize');
           [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
              = rri_boot_check(num_subj_lst, num_conditions, num_boot, 0, ...
                for_batch);
%                boot_progress, for_batch);

           num_boot = new_num_boot;
        end;

        if (num_perm > 0) | (num_boot == 0),
           [brainlv,s,designlv,behavlv,brainscores,d_scores,behavscores,lvcorrs, ...
		origpost,perm_result,behavdata,lv_evt_list,behavdata_lst,datamatcorrs_lst, ...
		b_scores,behav_row_idx] = ...
			fmri_perm_multiblock(st_datamat,contrasts,st_evt_list, ...
			ibehavdata_lst, newdata_lst, num_subj_lst, ...
			num_perm,num_conditions,num_behav_subj,posthoc,bscan);
        end;

        if (num_boot > 0),
           if num_perm == 0, origpost = []; end;
           [brainlv2,s2,designlv2,behavlv2,brainscores2,d_scores2,behavscores2,lvcorrs2, ...
		boot_result,behavdata2,lv_evt_list2,behavdata_lst2,datamatcorrs_lst2,b_scores2,behav_row_idx2] = ...
			fmri_boot_multiblock(st_datamat,contrasts,st_evt_list, ...
			ibehavdata_lst, newdata_lst, num_subj_lst, ...
			num_boot,num_conditions,num_behav_subj,Clim, ...
                        min_subj_per_group,is_boot_samples,boot_samples, ...
			new_num_boot,bscan);

           if num_perm == 0
              brainlv = brainlv2;
              s = s2;
              designlv = designlv2;
              behavlv = behavlv2;
              brainscores = brainscores2;
              d_scores = d_scores2;
              behavscores = behavscores2;
              lvcorrs = lvcorrs2;
              behavdata = behavdata2;
              lv_evt_list = lv_evt_list2;
              behavdata_lst = behavdata_lst2;
              datamatcorrs_lst = datamatcorrs_lst2;
              b_scores = b_scores2;
	      behav_row_idx = behav_row_idx2;
              perm_result = [];
           end
        end;

        ismultiblock = 1;

        saved_info=['brainlv s designlv behavlv brainscores b_scores d_scores behavscores lvcorrs ', ...
			'origpost perm_result boot_result st_coords ', ...
			'behavdata behavname datamatcorrs_lst ', ...
			'num_conditions subj_name cond_name cond_selection ', ...
			'st_dims lv_evt_list st_win_size st_voxel_size ', ...
			'subj_group behavdata_lst num_subj_lst ', ...
			'st_origin SessionProfiles ContrastFile ', ...
			'session_files_timestamp datamat_files_timestamp ', ...
			'create_ver ismultiblock bscan'];

        if save_datamat & ~isempty(brainlv)
            first = 1;
            last = 0;
%            grp_datamat = [];
            for g = 1:length(num_subj_lst)
               last = last + num_conditions*num_subj_lst(g);
               [tmp idx] = sort(st_evt_list(first:last));
               tmp = st_datamat(first:last,:);
%               grp_datamat = [grp_datamat; tmp(idx,:)];
               grp_datamat{g} = tmp(idx,:);
               first = last + 1;
            end;
        
            st_datamat = grp_datamat;
            saved_info = [saved_info, ' st_datamat'];
        end

  else					% behav & no rotate

     if (isbehav)			% Behavior Analysis

        if (num_boot > 0),
%           boot_progress = rri_progress_ui('initialize');
           [min_subj_per_group,is_boot_samples,boot_samples,new_num_boot] ...
              = rri_boot_check(num_subj_lst, num_conditions, num_boot, 0, ...
                for_batch);
%                boot_progress, for_batch);

           num_boot = new_num_boot;
        end;

        if (num_perm > 0) | (num_boot == 0),
           [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ...
		origpost,perm_result,behavdata,lv_evt_list,datamatcorrs_lst] = ...
			fmri_perm_behav(st_datamat,contrasts,st_evt_list, ...
			behavdata_lst, newdata_lst, num_subj_lst, ...
			num_perm,num_conditions,num_behav_subj,posthoc);
        end;

        if (num_boot > 0),
           if num_perm == 0, origpost = []; end;
           [brainlv2,s2,behavlv2,brainscores2,behavscores2,lvcorrs2, ...
		boot_result,behavdata,lv_evt_list2,datamatcorrs_lst2] = ...
			fmri_boot_behav(st_datamat,contrasts,st_evt_list, ...
			behavdata_lst, newdata_lst, num_subj_lst, ...
			num_boot,num_conditions,num_behav_subj,Clim, ...
                        min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);

           if num_perm == 0
              brainlv = brainlv2;
              s = s2;
              behavlv = behavlv2;
              brainscores = brainscores2;
              behavscores = behavscores2;
              lvcorrs = lvcorrs2;
              lv_evt_list = lv_evt_list2;
              datamatcorrs_lst = datamatcorrs_lst2;
              perm_result = [];
           end
        end;

        saved_info=['brainlv s behavlv brainscores behavscores lvcorrs ', ...
			'origpost perm_result boot_result st_coords ', ...
			'behavdata behavname datamatcorrs_lst ', ...
			'num_conditions subj_name cond_name cond_selection ', ...
			'st_dims lv_evt_list st_win_size st_voxel_size ', ...
			'subj_group behavdata_lst num_subj_lst ', ...
			'st_origin SessionProfiles ContrastFile ', ...
			'session_files_timestamp datamat_files_timestamp ', ...
			'create_ver'];

        if save_datamat & ~isempty(brainlv)
            first = 1;
            last = 0;
%            grp_datamat = [];
            for g = 1:length(num_subj_lst)
               last = last + num_conditions*num_subj_lst(g);
               [tmp idx] = sort(st_evt_list(first:last));
               tmp = st_datamat(first:last,:);
%               grp_datamat = [grp_datamat; tmp(idx,:)];
               grp_datamat{g} = tmp(idx,:);
               first = last + 1;
            end;