%%can specify an excel sheet with all your subjects listed and then read in %%the list [~,subs]=xlsread('/Volumes/data/Data/IMP/GroupAnalyses/Adults_2runs_combined/IMP_adults_successful.xlsx','Sheet1'); [~,contrasts]=xlsread('/Volumes/data/Data/IMP/GroupAnalyses/Contrasts_IMP_adult_combined.xlsx','Sheet1'); %specify pathname pathname='/Volumes/data/Data/IMP/GroupAnalyses/Adults_2runs_combined/' %read in a contrasts file (e.g. excel) so you end up with a contrasts %variable (same as other analysis) HB = {'R', 'O'}; brains = {'left', 'right', 'frontal'}; for a = 1:length(HB) hb = HB{1,a}; for b = 1:length(brains) spec_hemi = brains{1,b} ; for i=1:length(subs) fname_ginterp_betas{i,1}=strcat(pathname, subs{i,1}, '/Analysis_stepremoved_IMP/interp_beta_',spec_hemi,'Hb', hb, '.mat'); end minimumnumbers=[2]; for jjj=1:length(minimumnumbers) min_subj=minimumnumbers(1,jjj); nsubj=length(subs); dir_spm = strcat(pathname, 'IMP_group_min', num2str(min_subj),'/Hb', hb, spec_hemi, '/'); mkdir(dir_spm) [gavg_beta, group_beta, xX, index_group, nsubj_mask, brain_view, fname_interp_cov] = nirs_spm_group(fname_ginterp_betas, min_subj); SPM_nirs.nirs.level = 'group'; SPM_nirs.nirs.nsubj = nsubj; SPM_nirs.nirs.min_subj = min_subj; SPM_nirs.nirs.Hb = strcat('Hb',hb); SPM_nirs.xX = xX; SPM_nirs.nirs.brain_view = brain_view; %name,index,size SPM_nirs.nirs.index_group = index_group; SPM_nirs.nirs.nsubj_mask = nsubj_mask; SPM_nirs.nirs.fname_ginterp_cov = fname_interp_cov; SPM_nirs.nirs.fname_ginterp_beta = [dir_spm filesep 'ginterp_beta_' SPM_nirs.nirs.brain_view.name SPM_nirs.nirs.Hb '.mat']; SPM_nirs.nirs.fname_ginterp_avg_beta = [dir_spm filesep 'ginterp_avgbeta_' SPM_nirs.nirs.brain_view.name SPM_nirs.nirs.Hb '.mat']; fname_others = cellstr(spm_select('FPList', dir_spm, ['^ginterp_T.*\' SPM_nirs.nirs.Hb '.mat$'])); if strcmp(fname_others{1}, filesep) ~= 1 delete(fname_others{:}); end fname_others = cellstr(spm_select('FPList', dir_spm, ['^ginterp_F.*\' SPM_nirs.nirs.Hb '.mat$'])); if strcmp(fname_others{1}, filesep) ~= 1 delete(fname_others{:}); end field1='name'; value1=contrasts(:,1); field2='STAT'; value2='T'; field3='c'; vector = double(str2num(contrasts{1,2}))'; value3=vector; xCon=struct(field1,value1,field2,value2,field3,value3); for d=1:length(contrasts) xCon(d).c=double(str2num(contrasts{d,2}))'; end SPM_nirs.xCon=xCon; % save a SPM_nirs_group_brainview_HbX.mat in specified dir. save([dir_spm filesep 'SPM_group_' SPM_nirs.nirs.brain_view.name SPM_nirs.nirs.Hb '.mat'], 'SPM_nirs'); % save a group-average beta save(SPM_nirs.nirs.fname_ginterp_avg_beta, 'gavg_beta'); % save matrix containing individual interpolated betas save(SPM_nirs.nirs.fname_ginterp_beta, 'group_beta'); for Ic=1:length(contrasts) %------------------------------------------------------------------ % calculation of group-level T- or F-statistics %------------------------------------------------------------------ stat = []; % T- or F-statistic matrix % loading a file (interpolated T- or F-statistics) filename = [dir_spm 'ginterp_' xCon(Ic).STAT '_' num2str(Ic) '_' SPM_nirs.nirs.brain_view.name SPM_nirs.nirs.Hb '.mat']; fid = fopen(filename); if fid ~= -1 % if exists fclose(fid); try load(filename); end end %if isempty(stat) == 1 || isempty(index_maskkj) == 1 index_maskkj = []; index_groupkj = SPM_nirs.nirs.index_group; nvox = length(index_groupkj); nvox_brain = SPM_nirs.nirs.brain_view.size(1) * SPM_nirs.nirs.brain_view.size(2); nsubj = SPM_nirs.nirs.nsubj; L = SPM_nirs.nirs.nsubj_mask(index_groupkj); % # of subject on each voxel % contrast * group-level interp beta load(SPM_nirs.nirs.fname_ginterp_avg_beta); cgavg_beta = xCon(Ic).c' * gavg_beta; load(SPM_nirs.nirs.fname_ginterp_beta); % group-level residuals Xg = ones(nsubj, 1); Yg = xCon(Ic).c' * [group_beta{:}]; Yg = reshape(Yg', [nvox, nsubj])'; res = Yg - Xg * cgavg_beta; % LKC calculation disp('Calculating Lipschitz-Killing curvatures ...'); [L2] = calc_LKC(index_groupkj, SPM_nirs.nirs.brain_view.size, res, SPM_nirs.nirs.level); disp('Completed.'); r = sqrt(L2./pi); L1 = pi * r; L0 = 1; SPM_nirs.nirs.LKC{SPM_nirs.nirs.brain_view.index} = [L0 L1 L2]; SPM_nirs.xCon=xCon; df = []; % degree of freedom switch xCon(Ic).STAT case 'T' % t-statistic calculation var_s = zeros(1, nvox); for kk = 1:nsubj var_s = var_s + (xCon(Ic).c' * group_beta{kk} - cgavg_beta).^2; end var_s = var_s ./ (L-1); gsum_var = zeros(1, nvox_brain); term_df = zeros(1, nvox_brain); % calculation group-level variance for kk = 1:nsubj load(SPM_nirs.nirs.fname_ginterp_cov{kk}); index_maskkj=index_mask; indiv_cov = interp_var .* (xCon(Ic).c' * xCor * xCon(Ic).c); gsum_var(index_maskkj) = gsum_var(index_maskkj) + indiv_cov; term_df(index_maskkj) = term_df(index_maskkj) + (indiv_cov.^2)./df(2); end gsum_var = gsum_var(index_group); term_df = term_df(index_group); term = var_s .* L + gsum_var; stat = (xCon(Ic).c' * gavg_beta)./(sqrt(term)./L); term2 = ((L.^2)./(L-1)).* (var_s.^2); df(2) = sum((term.^2) ./ (term2 + term_df))./nvox; % degree of freedom case 'F' % F-statistic calculation X1o = pinv(Xg'); [trMV, trMVMV] = spm_SpUtil('trMV', X1o, 1); df(1) = trMV^2/trMVMV; R = eye(nsubj) - Xg*pinv(Xg); RVR = sum(res.^2)./trace(R); MVM = (inv(sum(X1o.^2)).* (cgavg_beta.^2))./trMV; stat = MVM./RVR; df(2) = trace(R)^2./trace(R*R); end stat = stat(:); index_maskkj = index_groupkj; save(filename, 'stat', 'df', 'index_maskkj'); %end side_hemi = SPM_nirs.nirs.brain_view.index; load([spm('dir') filesep 'rend' filesep 'render_single_subj.mat']); %%% modified (2008. 10. 13) %%% brain = rend{side_hemi}.ren; if issparse(brain), d = size(brain); B1 = spm_dctmtx(d(1),d(1)); B2 = spm_dctmtx(d(2),d(2)); brain = B1*brain*B2'; end msk = find(brain>1);brain(msk)=1; msk = find(brain<0);brain(msk)=0; brain = brain(end:-1:1,:); brain = brain * 0.5; min_stat = min(stat); max_stat = max(stat); smin_stat = max_stat - ((max_stat - min_stat)./63) * 127; sbar = linspace(smin_stat, max_stat, 128); stat_brain = ((-sbar(1) + sbar(64))/(0.5)).*brain + sbar(1); stat_brain(index_maskkj) = stat; imagesc(stat_brain); load Split colormap(split) axis off axis image hc = colorbar; set(hc, 'YLim', [sbar(65) sbar(128)]); y_tick = linspace(sbar(65), sbar(128), 5)'; set(hc, 'YTick', y_tick); set(hc, 'FontSize', 8); p_value=0.05; threshold = spm_u(p_value, df, xCon(Ic).STAT); index_over = find(stat > threshold); if isempty(index_over) == 1 disp('There is no significant voxel.'); else index_maskkj = index_maskkj(index_over); stat = stat(index_over); % min and max values of statistics for colormap control min_stat = min(stat); max_stat = max(stat); smin_stat = max_stat - ((max_stat - min_stat)./63) * 127; sbar = linspace(smin_stat, max_stat, 128); stat_brain = ((-sbar(1) + sbar(64))/(0.5)).*brain + sbar(1); stat_brain(index_maskkj) = stat; figure imagesc(stat_brain); load Split; colormap(split) axis off axis image hc = colorbar; if sbar(65)==sbar(128) disp('no sig after thresholding -- sbar issue') else set(hc, 'YLim', [sbar(65) sbar(128)]); y_tick = linspace(sbar(65), sbar(128), 5)'; set(hc, 'YTick', y_tick); set(hc, 'FontSize', 8); saveas(gcf,strcat(dir_spm,contrasts{Ic,1},'.png')); close all; end end end end end end