Diff for "BatchRenderRfx" - MRC CBU Imaging Wiki
location: Diff for "BatchRenderRfx"
Differences between revisions 1 and 37 (spanning 36 versions)
Revision 1 as of 2006-10-01 20:12:08
Size: 4890
Comment:
Revision 37 as of 2008-03-12 18:40:49
Size: 9160
Editor: pc0302
Comment:
Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
Modified version Ferath's script to threshhold and print SPM5 results.
The script prints rendered activations as well as a list of coordinates.
The first six blocks of lines are the ones you need to (or are most likely to
want to) change. Originally modifi
ed by DannyMitchell and subsequently by IanNimmoSmith.
Modified version Ferath's script (attachment:BatchRenderRFX.m) to threshhold and print SPM5 results. The script prints glass brains and rendered activations as well as the first page of the results table. The lines in block (3) are the ones you need to (or are most likely to want to) change. Originally modified by DannyMitchell and subsequently by IanNimmoSmith.
Line 6: Line 3:
 1. The program needs files attachment:cls_getRes.m and attachment:cls_getSPM2.m . For CBU people the following line should find them, otherwise you need to download them and point to their path. {{{addpath /imaging/dm01/MoreTools/spm2Batch}}}. Also, you should add a path to the custom spm_render_dm, which solves some of the issues of the default spm_render command {{{addpath /imaging/dm01/MEG/aaMEG}}}. Finally, if you are using AAL, you should add the path to its toolkit, e.g. {{addpath /imaging/av02/tools/aal}}
Line 7: Line 5:
 1. Path to data, and where the results are printed
 {{{dataroot = '/imaging/dm01/TestReCon';}}}
 1. Last updated 03/08 by DannyMitchell and [http://imaging.mrc-cbu.cam.ac.uk/basewiki/AlejandroVicenteGrabovetsky AlejandroVicenteGrabovetsky]
Line 10: Line 7:
 1. Name of folder, in above directory, where analysis is
 {{{outputdir='';}}}
 1. {{{
dataroot = ''; % path to data
anadirs = {''}; % name of folder(s), in above directory, where analysis is (useful if you have used several models in analysis)
cond_dir={''}; % cell array of folder names in outputdir, containing the SPMs for each contrast
% %%(AVG) Alternatively, save the data from your 'aamod_firstlevel_contrasts.m' file onto a 'contrast_names.mat' file and do the following:
% load contrast_names
% cond_dir=cname; %cname is the array of contrast names and their directory names in your second level analysis folder
outputfilesuffix='_RFX_Rendered.ps'
%{set the following variable to 'yes' to do a 2-tailed test and show
% rendered activations and deactivations (Glass brains will
% only show positive activations). Set it to 'no' to to a 1-tailed test
% and only calculate positive activations. %}
ShowDeactivations = 'no';
Threshold=0.05; % p threshold
ExtentThreshold=0; % minimum number of contiguous voxels
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.%}
MCC='FDR'; % method for correcting for multiple comparisons : 'FWE','FDR', or 'none'
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
Style=1; % Rendering Style: NaN=old, 1=normal, <1=brighter (0.75 = light, 0.5 = more, 0.25 = lots)
% for results table
num=3; % number of maxima per cluster
dis=8; % distance among clusters (mm)
doaal = 0; % Label anatomical clusters with AAL? 1=yes 0=no
ignore_empty = 1 % (AVG) Create a file if there are no significant voxels? 1=yes 0=no}}}
Line 13: Line 32:
 1. Cell array of folder names in outputdir, containing the SPMs for each contrast
 {{{cond_dir={'EarlySmoothNoReCon'};}}}

 1. Set the following variable to 'yes' to do a 2-tailed test and show rendered activations and deactivations (assumes that positive and negative second level contrasts have already been calculated; Glass brains will
 only show positive activations). Set it to 'no' to to a 1-tailed test and only calculate positive activations.
 {{{ShowDeactivations = 'no';}}}

 1. p threshold
 {{{Threshold=0.05;}}}

 1. Minimum number of contiguous voxels
 {{{ExtentThreshold=0;}}}

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.

{{{addpath /imaging/dm01/MoreTools/spm2Batch}}}

Method for correcting for multiple comparisons : 'FWE','FDR', or 'none'
{{{MCC='FDR';}}}

The brain on which you want to render the results
{{{Brain='/imaging/local/spm/spm5/rend/render_smooth_average.mat';}}}

{{{%%%%%%% Do it:
if exist(fullfile(dataroot,strcat(outputdir,'_RFX_Rendered.ps')))==2
    delete(fullfile(dataroot,strcat(outputdir,'_RFX_Rendered.ps')));
 1. {{{
%%%%%%% Do it:
for j = 1:size(anadirs,2)
    outputdir=anadirs{j};
if exist(fullfile(dataroot,strcat(outputdir,outputfilesuffix)),'file')
    delete(fullfile(dataroot,strcat(outputdir,outputfilesuffix)));
Line 42: Line 39:
disp(strcat('Printing results to:...',fullfile(dataroot,strcat(outputdir,'_RFX_Rendered.ps'))))
Line 46: Line 42:
try q=defaults.units{3}; catch; defaults.units={'mm','mm','mm'}; end;
}}}
Line 47: Line 45:
%--------------------------------------------------------------------  1. {{{
Line 49: Line 47:
Line 60: Line 57:
clear SPM}}}
Line 61: Line 59:
clear SPM

for i = 1:size(cond_dir,2)
    condir = fullfile(dataroot,outputdir,cond_dir{i});
    cd(condir);
    xSPM.spmmat = fullfile(condir,'SPM.mat'); % Get SPM path
    if ShowDeactivations(1)=='y'
        xSPM.u=0.025; % 2-tailed
        TailText='2-tailed, p<0.05';
    else
        xSPM.u=0.05; % 1-tailed
        TailText='1-tailed, p<0.05';
    end
    xSPM.k=0;
    bck_xSPM=xSPM;
    delete(gcf) % this is a bit annoying, but for some reason I had problems without it
    for i=2:-1:1
        xSPM= bck_xSPM; % Reset xSPM to default
        xSPM.Ic=i; %Set the current contrast Index
        [hReg,xSPM,SPM] = csl_getRes(xSPM);
        dat(3-i) = struct( 'XYZ', xSPM.XYZ,...
            't', xSPM.Z',...
            'mat', xSPM.M,...
            'dim', xSPM.DIM);
    end
    clear flipdat
    flipdat(1)=dat(2);
    if ShowDeactivations(1)=='y'
        flipdat(2)=dat(1);
    end
    try
        spm_render(flipdat,0.5,Brain);
        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');
    catch
 1. {{{
    for i = 1:size(cond_dir,2)
        condir = fullfile(dataroot,outputdir,cond_dir{i});
        % Checks if there is a rendered file made for this contrast and deletes the old one
        if exist(strcat(condir,outputfilesuffix),'file')
            delete(strcat(condir,outputfilesuffix));
        end;
        cd(condir);
        xSPM.spmmat = fullfile(condir,'SPM.mat'); % Get SPM path
        %%%%%%%%% modified 080308 to run 2nd level contrasts if not found
        % Set significance thresholds and run t-tests if not done already
        load(xSPM.spmmat)
        if isempty(SPM.xCon)
            SPM.xCon = spm_FcUtil('Set',cond_dir{i},'T','c',1,SPM.xX.xKXs);
        end
        if lower(ShowDeactivations(1))=='y'
            xSPM.u=Threshold/2; % 2-tailed
            TailText=strcat('2-tailed, p<',num2str(Threshold));
            if length(SPM.xCon)<2
                SPM.xCon(end+1) = spm_FcUtil('Set',strcat('Negative_',cond_dir{i}),'T','c',-1,SPM.xX.xKXs);
            end
        else
            xSPM.u=Threshold; % 1-tailed
            TailText=strcat('1-tailed, p<',num2str(Threshold));
        end
        if isempty(SPM.xCon(1).Vcon)
            spm_contrasts(SPM);
        end
        %%%%%%%%%%%%%%
        xSPM.k=ExtentThreshold;
        bck_xSPM=xSPM;
        delete(gcf) % this is a bit annoying, but for some reason I had problems without it
        clear dat
        for con=2:-1:1 % get any -ve tail first, so that glass brain ends up showing +ve tail
            xSPM=bck_xSPM; % Reset xSPM to default
            xSPM.Ic=con; %Set the current contrast Index
            if strcmp(ShowDeactivations,'no') && con==2; continue; end;
            try
                [hReg,xSPM,SPM] = csl_getRes(xSPM);
                if exist('dat','var')
                    dat(end+1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
                else
                    dat(1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
                end
            catch
                xSPM.u=Threshold; % 1-tailed
                TailText=strcat('1-tailed, p<',num2str(Threshold));
            end
        end
        dat=circshift(dat,[1 1]);
Line 96: Line 110:
            spm_render(flipdat(1),0.5,Brain);             % should work if all tails tested have significant voxels
            spm_render_dm(dat,Style,Brain);
Line 100: Line 115:
                spm_render(flipdat(2),0.5,Brain);
                ann1=annotation('textbox',[0 .94 1 .06],'Color','r','String',{strcat('RFX results (',TailText,' correction=',MCC,', RED=NEGATIVE)'),strcat('Data:...',condir),strcat('Date:...',datestr(now,0))},'edge','none');
                % just render +ve tail if no significant -ve voxels
spm_render_dm(dat(1),Style,Brain);
                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');
Line 103: Line 119:
                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');                 try
                    % just render -ve tail if no significant +ve voxels
                    dat(1)=dat(2); dat(1).t=zeros(length(dat(1).XYZ),1); % set dat(1) as empty dat(2)
                    spm_render_dm(dat,Style,Brain);
                    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');
                catch
                    % probably no significant voxels
                    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');
                    if ignore_empty == 0
                        %%(AVG) Alternatively: probably no need to draw, so do next contrast rather than print this one
                        continue;
                    end
                end
Line 106: Line 134:
        try spm_print(strcat(condir,outputfilesuffix));
        catch fprintf('\nFailed to print! (check write permissions?)'); return
        end
        TabDat=spm_list('List',xSPM,hReg,num, dis);
        tot=0;
        for r=1:size(TabDat.dat,1)
            if ~isempty(TabDat.dat{r,4}); tot=tot+TabDat.dat{r,4};end
        end
        if tot>0
            delete(ann1);
            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');
            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');
            delete(ann2);
        end
        %Following section borrowed from Mirjana Bozic's script
        if doaal == 1
            % SPM_print anatomical labels for clusters using batch-enabled AAL scripts
            gin_clusters_plabels_auto('List',xSPM);
            spm_print(strcat(condir,outputfilesuffix));
        end
        xSPM= bck_xSPM; % Reset xSPM to default
        delete(ann1);
Line 107: Line 157:
    spm_print(fullfile(dataroot,strcat(outputdir,'_RFX_Rendered.ps')));
    spm_list('List',xSPM,hReg);
    spm_print(fullfile(dataroot,strcat(outputdir,'_RFX_Rendered.ps')));
    xSPM= bck_xSPM; % Reset xSPM to default
    delete(ann1);
    fprintf('\nFinished.\nResults saved to: %s.\r',strcat(condir,outputfilesuffix));

Modified version Ferath's script (attachment:BatchRenderRFX.m) to threshhold and print SPM5 results. The script prints glass brains and rendered activations as well as the first page of the results table. The lines in block (3) are the ones you need to (or are most likely to want to) change. Originally modified by DannyMitchell and subsequently by IanNimmoSmith.

  1. The program needs files attachment:cls_getRes.m and attachment:cls_getSPM2.m . For CBU people the following line should find them, otherwise you need to download them and point to their path. addpath /imaging/dm01/MoreTools/spm2Batch. Also, you should add a path to the custom spm_render_dm, which solves some of the issues of the default spm_render command addpath /imaging/dm01/MEG/aaMEG. Finally, if you are using AAL, you should add the path to its toolkit, e.g. addpath /imaging/av02/tools/aal

  2. Last updated 03/08 by DannyMitchell and [http://imaging.mrc-cbu.cam.ac.uk/basewiki/AlejandroVicenteGrabovetsky AlejandroVicenteGrabovetsky]

  3. dataroot = ''; % path to data
    anadirs = {''}; % name of folder(s), in above directory, where analysis is (useful if you have used several models in analysis)
    cond_dir={''}; % cell array of folder names in outputdir, containing the SPMs for each contrast
    % %%(AVG) Alternatively, save the data from your 'aamod_firstlevel_contrasts.m' file onto a 'contrast_names.mat' file and do the following:
    % load contrast_names
    % cond_dir=cname; %cname is the array of contrast names and their directory names in your second level analysis folder
    outputfilesuffix='_RFX_Rendered.ps'
    %{set the following variable to 'yes' to do a 2-tailed test and show
    % rendered activations and deactivations (Glass brains will
    % only show positive activations). Set it to 'no' to to a 1-tailed test
    % and only calculate positive activations. %}
    ShowDeactivations = 'no';
    Threshold=0.05; % p threshold
    ExtentThreshold=0; % minimum number of contiguous voxels
    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.%}
    MCC='FDR'; % method for correcting for multiple comparisons : 'FWE','FDR', or 'none'
    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
    Style=1; % Rendering Style: NaN=old, 1=normal, <1=brighter (0.75 = light, 0.5 = more, 0.25 = lots)
    % for results table
    num=3;    % number of maxima per cluster
    dis=8;    % distance among clusters (mm)
    doaal = 0; % Label anatomical clusters with AAL? 1=yes 0=no
    ignore_empty = 1 % (AVG) Create a file if there are no significant voxels? 1=yes 0=no
  4. %%%%%%% Do it:
    for  j = 1:size(anadirs,2)
        outputdir=anadirs{j};
    if exist(fullfile(dataroot,strcat(outputdir,outputfilesuffix)),'file')
        delete(fullfile(dataroot,strcat(outputdir,outputfilesuffix)));
    end;
    spm_defaults;
    global defaults
    defaults.modality='FMRI';
    try q=defaults.units{3}; catch; defaults.units={'mm','mm','mm'}; end;
  5. %- xSPM is a structure for the parameters
    xSPM = struct( ...
        'swd', '', ... % full path to SPM.mat file
        'Ic', [], ... % no of contrast (or contrasts for conjunction)
        'Im', [],... % no of contrast to mask with. Empty for no masking
        'pm', [],... % masking contrast uncorrected p
        'Ex', [],... % whether masking is inclusive or exclusive
        'title', '',... % if empty results in default contrast title
        'Mcp', MCC,... % Mutiple comp method: FWE|FDR|none
        'u', Threshold,... % threshold (corrected or uncorrected, as above)
        'k', ExtentThreshold); % extent threshold
    clear SPM
  6.     for i = 1:size(cond_dir,2)
            condir = fullfile(dataroot,outputdir,cond_dir{i});
            % Checks if there is a rendered file made for this contrast and deletes the old one
            if exist(strcat(condir,outputfilesuffix),'file')
                delete(strcat(condir,outputfilesuffix));
            end;
            cd(condir);
            xSPM.spmmat = fullfile(condir,'SPM.mat'); % Get SPM path
            %%%%%%%%% modified 080308 to run 2nd level contrasts if not found
            % Set significance thresholds and run t-tests if not done already
            load(xSPM.spmmat)
            if isempty(SPM.xCon)
                SPM.xCon = spm_FcUtil('Set',cond_dir{i},'T','c',1,SPM.xX.xKXs);
            end
            if lower(ShowDeactivations(1))=='y'
                xSPM.u=Threshold/2; % 2-tailed
                TailText=strcat('2-tailed, p<',num2str(Threshold));
                if length(SPM.xCon)<2
                    SPM.xCon(end+1) = spm_FcUtil('Set',strcat('Negative_',cond_dir{i}),'T','c',-1,SPM.xX.xKXs);
                end
            else
                xSPM.u=Threshold; % 1-tailed
                TailText=strcat('1-tailed, p<',num2str(Threshold));
            end
            if isempty(SPM.xCon(1).Vcon)
                spm_contrasts(SPM);
            end
            %%%%%%%%%%%%%%
            xSPM.k=ExtentThreshold;
            bck_xSPM=xSPM;
            delete(gcf) % this is a bit annoying, but for some reason I had problems without it
            clear dat
            for con=2:-1:1 % get any -ve tail first, so that glass brain ends up showing +ve tail
                xSPM=bck_xSPM; % Reset xSPM to default
                xSPM.Ic=con; %Set the current contrast Index
                if strcmp(ShowDeactivations,'no') && con==2; continue; end;
                try
                    [hReg,xSPM,SPM] = csl_getRes(xSPM);
                    if exist('dat','var')
                        dat(end+1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
                    else
                        dat(1) = struct('XYZ', xSPM.XYZ,'t', xSPM.Z','mat', xSPM.M,'dim', xSPM.DIM);
                    end
                catch
                    xSPM.u=Threshold; % 1-tailed
                    TailText=strcat('1-tailed, p<',num2str(Threshold));
                end
            end
            dat=circshift(dat,[1 1]);
            try
                % should work if all tails tested have significant voxels
                spm_render_dm(dat,Style,Brain);
                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');
            catch
                try
                    % just render +ve tail if no significant -ve voxels
                    spm_render_dm(dat(1),Style,Brain);
                    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');
                catch
                    try
                        % just render -ve tail if no significant +ve voxels
                        dat(1)=dat(2); dat(1).t=zeros(length(dat(1).XYZ),1); % set dat(1) as empty dat(2)
                        spm_render_dm(dat,Style,Brain);
                        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');
                    catch
                        % probably no significant voxels
                        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');
                        if ignore_empty == 0
                            %%(AVG) Alternatively: probably no need to draw, so do next contrast rather than print this one
                            continue;
                        end
                    end
                end
            end
            try spm_print(strcat(condir,outputfilesuffix));
            catch fprintf('\nFailed to print! (check write permissions?)'); return
            end
            TabDat=spm_list('List',xSPM,hReg,num, dis);
            tot=0;
            for r=1:size(TabDat.dat,1)
                if ~isempty(TabDat.dat{r,4}); tot=tot+TabDat.dat{r,4};end
            end
            if tot>0
                delete(ann1);
                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');
                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');
                delete(ann2);
            end
            %Following section borrowed from Mirjana Bozic's script
            if doaal == 1
                % SPM_print anatomical labels for clusters using batch-enabled AAL scripts
                gin_clusters_plabels_auto('List',xSPM);
                spm_print(strcat(condir,outputfilesuffix));
            end
            xSPM= bck_xSPM; % Reset xSPM to default
            delete(ann1);
        end
        fprintf('\nFinished.\nResults saved to: %s.\r',strcat(condir,outputfilesuffix));
    end

CbuImaging: BatchRenderRfx (last edited 2013-03-07 21:24:07 by localhost)