1 % Modified version of Ferath's script to threshhold and print SPM5 second
   2 % level results.
   3 % Prints rendered activations as well as a list of coordinates.
   4 % The first 11 lines are the ones you need to (or are most likely to
   5 % want to) change.
   6 % djm 03/08
   7 
   8 dataroot = ''; % path to data
   9 anadirs = {''}; % name of folder(s), in above directory, where analysis is (useful if you have used several models in analysis)
  10 cond_dir={''}; % cell array of folder names in outputdir, containing the SPMs for each contrast
  11 % %%(AVG) Alternatively, save the data from your 'aamod_firstlevel_contrasts.m' file onto a 'contrast_names.mat' file and do the following:
  12 % load contrast_names
  13 % cond_dir=cname; %cname is the array of contrast names and their directory names in your second level analysis folder
  14 
  15 outputfilesuffix='_RFX_Rendered.ps'
  16 %{set the following variable to 'yes' to do a 2-tailed test and show
  17 % rendered activations and deactivations (assumes that positive and negative
  18 % second level contrasts have already been calculated; Glass brains will
  19 % only show positive activations). Set it to 'no' to to a 1-tailed test
  20 % and only calculate positive activations. %}
  21 ShowDeactivations = 'no';
  22 Threshold=0.05; % p threshold
  23 ExtentThreshold=0; % minimum number of contiguous voxels
  24 addpath /imaging/dm01/MoreTools/spm2Batch %{the program needs files cls_getRes.m and cls_getSPM2.m. For CBU people this line will find them, otherwise you need to download them and point to their path.%}
  25 MCC='FDR'; % method for correcting for multiple comparisons : 'FWE','FDR', or 'none'
  26 Brain='/imaging/local/spm/spm5/rend/render_smooth_average.mat'%render_single_subj.mat' %render_smooth_average.mat'; % the brain on which you want to render the results
  27 Style=1; % Rendering Style: NaN=old, 1=normal, <1=brighter (0.75 = light, 0.5 = more, 0.25 = lots)
  28 
  29 % for results table
  30 num=3;    % number of maxima per cluster
  31 dis=8;    % distance among clusters (mm)
  32 doaal = 0; % Label anatomical clusters with AAL? 1=yes 0=no
  33 ignore_empty = 1 % (AVG) Create a file if there are no significant voxels? 1=yes 0=no
  34 
  35 %%%%%%% Do it:
  36 for  j = 1:size(anadirs,2)
  37     outputdir=anadirs{j};
  38 
  39     if exist(fullfile(dataroot,strcat(outputdir,outputfilesuffix)))==2
  40         delete(fullfile(dataroot,strcat(outputdir,outputfilesuffix)));
  41     end;
  42 
  43     spm_defaults;
  44     global defaults
  45     defaults.modality='FMRI';
  46     try q=defaults.units{3}; catch; defaults.units={'mm','mm','mm'}; end; % for some reason I had to add this line 15/01/07
  47 
  48     %- xSPM is a structure for the parameters
  49 
  50 
  51     xSPM = struct( ...
  52         'swd', '', ... % full path to SPM.mat file
  53         'Ic', [], ... % no of contrast (or contrasts for conjunction)
  54         'Im', [],... % no of contrast to mask with. Empty for no masking
  55         'pm', [],... % masking contrast uncorrected p
  56         'Ex', [],... % whether masking is inclusive or exclusive
  57         'title', '',... % if empty results in default contrast title
  58         'Mcp', MCC,... % Mutiple comp method: FWE|FDR|none
  59         'u', Threshold,... % threshold (corrected or uncorrected, as above)
  60         'k', ExtentThreshold); % extent threshold
  61 
  62     clear SPM
  63 
  64     for i = 1:size(cond_dir,2)
  65         condir = fullfile(dataroot,outputdir,cond_dir{i});
  66 
  67         % Checks if there is a rendered file made for this contrast and deletes the old one
  68         if exist(strcat(condir,outputfilesuffix),'file')==2
  69             delete(strcat(condir,outputfilesuffix));
  70         end;
  71 
  72         cd(condir);
  73         xSPM.spmmat = fullfile(condir,'SPM.mat'); % Get SPM path
  74 
  75         %%%%%%%%% modified 080308 to run 2nd level contrasts if not found
  76         % Set significance thresholds and run t-tests if not done already
  77         load(xSPM.spmmat)
  78         if isempty(SPM.xCon)
  79             SPM.xCon = spm_FcUtil('Set',cond_dir{i},'T','c',1,SPM.xX.xKXs);
  80         end
  81         if lower(ShowDeactivations(1))=='y'
  82             xSPM.u=Threshold/2; % 2-tailed
  83             TailText=strcat('2-tailed, p<',num2str(Threshold));
  84             if length(SPM.xCon)<2
  85                 SPM.xCon(end+1) = spm_FcUtil('Set',strcat('Negative_',cond_dir{i}),'T','c',-1,SPM.xX.xKXs);
  86             end
  87         else
  88             xSPM.u=Threshold; % 1-tailed
  89             TailText=strcat('1-tailed, p<',num2str(Threshold));
  90         end
  91         if isempty(SPM.xCon(1).Vcon)
  92             spm_contrasts(SPM);
  93         end
  94         %%%%%%%%%%%%%%
  95 
  96         xSPM.k=ExtentThreshold;
  97         bck_xSPM=xSPM;
  98         delete(gcf) % this is a bit annoying, but for some reason I had problems without it
  99         clear dat
 100         for con=2:-1:1 % get any negative tail first, so that glass brain ends up showing positive tail
 101             xSPM=bck_xSPM; % Reset xSPM to default
 102             xSPM.Ic=con; %Set the current contrast Index
 103             if strcmp(ShowDeactivations,'no') && con==2; continue; end;
 104             try
 105                 [hReg,xSPM,SPM] = csl_getRes(xSPM);
 106 
 107                 if exist('dat','var')
 108                     dat(end+1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
 109                 else
 110                     dat(1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
 111                 end
 112             catch
 113                 xSPM.u=Threshold; % 1-tailed
 114                 TailText=strcat('1-tailed, p<',num2str(Threshold));
 115             end
 116         end
 117 
 118         dat=circshift(dat,[1 1]);
 119 
 120         try
 121             % should work if all tails tested have significant voxels
 122             spm_render_dm(dat,Style,Brain);
 123             ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,', red=positive)'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
 124         catch
 125             try
 126                 % just render +ve tail if no significant -ve voxels
 127                 spm_render_dm(dat(1),Style,Brain);
 128                 ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,', red=positive)'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
 129             catch
 130                 try
 131                     % just render -ve tail if no significant +ve voxels
 132                     dat(1)=dat(2); dat(1).t=zeros(length(dat(1).XYZ),1); % set dat(1) as empty dat(2)
 133                     spm_render_dm(dat,Style,Brain);
 134                     ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,', blue=negative)'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
 135                 catch
 136                     % probably no significant voxels
 137                     ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,') No significant effects'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
 138                     if ignore_empty == 0
 139                         %%(AVG) Alternatively: probably no need to draw, so do next contrast rather than print this one
 140                         continue;
 141                     end
 142                 end
 143             end
 144         end
 145 
 146         try spm_print(strcat(condir,outputfilesuffix));
 147         catch fprintf('\nFailed to print! (check write permissions?)'); return
 148         end
 149 
 150         TabDat=spm_list('List',xSPM,hReg,num, dis);
 151         tot=0;
 152         for r=1:size(TabDat.dat,1)
 153             if ~isempty(TabDat.dat{r,4}); tot=tot+TabDat.dat{r,4};end
 154         end
 155         if tot>0
 156             delete(ann1);
 157             ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,')'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
 158             ann2=annotation('textbox',[0 0 1 .02],'Color','r','String',strcat('First page of table only. Total number of significant voxels=',num2str(tot)),'edge','none');
 159             try spm_print(strcat(condir,outputfilesuffix));
 160             catch fprintf('\nFailed to print! (check write permissions?)'); return
 161             end
 162             delete(ann2);
 163         end
 164 
 165         %Following section borrowed from Mirjana Bozic's script
 166         if doaal == 1
 167             % SPM_print anatomical labels for clusters using batch-enabled AAL scripts
 168             gin_clusters_plabels_auto('List',xSPM);
 169             spm_print(strcat(condir,outputfilesuffix));
 170         end
 171 
 172         xSPM= bck_xSPM; % Reset xSPM to default
 173         delete(ann1);
 174     end
 175     fprintf('\nFinished.\nResults saved to: %s.\r',strcat(condir,outputfilesuffix));
 176 end