MC_VBM_autorun_final.m

clear all

addpath(genpath('/usr/local/spm12'));
spm_jobman('initcfg');
input         = '/DATA/FKBP51_KOxBL6/input.txt';

betdel_flag     = 1;          % DELetion of previous bet files
prep_flag       = 1;          % PREP loop
bet_flag        = 1;          % brain extraction (BET)
seg_flag        = 1;          % SEGmentation
dartemp_flag    = 1;          % DARTEL: template creation
dartwarp_flag   = 1;          % DARTEL warp: normalisation to template6
masknorm_flag   = 1;          % MASK normalisation
fgm_flag        = 1;          % MASK final_group_mask creation
fgmsum_flag     = 1;          % MASK final_group_mask_sum creation
ero_flag        = 1;          % MASK erosion (for CSF)
sm_flag         = 1;          % SMOOTHING
fgmapp_flag     = 1;          % (eroded_)final_group_mask application
vols_flag       = 1;          % Global Calculations: VOLUMES

fid = fopen(input);
tline = fgetl(fid);
i=1;
while ischar(tline)
    data_set(i,:) = tline;
    tline = fgetl(fid);
    i=i+1;
end
fclose(fid);

tmp = strfind(data_set(1,:),'/');
if tmp(1) ~= 1
    error('\ndata path should be an absolute path, starting with / \n')
end

dirname = data_set(1,1:tmp(end));
cd(dirname);

number_of_files = size(data_set,1);
aaa=[' '];
for i=1:number_of_files
        tmp = strsplit(data_set(i,:),'/');
    filename(i,:) = tmp{end};
    aaa=[aaa,'bet_',filename(i,1:end-4),'* '];
end

% delete all previous bet files related to actual list of animals we run
if betdel_flag == 1
    cmd_out = ['rm',aaa];
    system(cmd_out);
end
    
%% PREP loop
% nameset           - input names for all procedures in one array
% fgm_exp           - final_group_mask expression for its creation /step xx/

if prep_flag == 1
    nameset=cell(number_of_files,17);
    fgm_exp='';
    for i=1:number_of_files
        input_name = filename(i,:);
        nameset(i,1)=cellstr(['bet_',input_name,',1']);                          %brain extraction + segmentation?
        nameset(i,2)=cellstr(['rc1bet_',input_name,',1']);                       %Dartel
        nameset(i,3)=cellstr(['rc2bet_',input_name,',1']);                       %Dartel
        nameset(i,4)=cellstr(['u_rc1bet_',input_name(1:end-4),'_Template.nii']); %Dartel warp
        nameset(i,5)=cellstr(['mbet_',input_name]);                              %Dartel warp
        nameset(i,6)=cellstr(['c1bet_',input_name]);                             %Dartel warp
        nameset(i,7)=cellstr(['c2bet_',input_name]);                             %Dartel warp
        nameset(i,8)=cellstr(['c3bet_',input_name]);                             %Dartel warp
        nameset(i,9)=cellstr(['bet_',input_name(1:end-4),'_mask.nii']);          %Dartel warp mask norm
        nameset(i,10)=cellstr(['wbet_',input_name(1:end-4),'_mask.nii']);        %Dartel warp mask norm
        nameset(i,11)=cellstr(['mwc1bet_',input_name]);                          %Smoothing
        nameset(i,12)=cellstr(['mwc1bet_',input_name,',1']);                     %fgm app input c1
        nameset(i,13)=cellstr(['mwc2bet_',input_name,',1']);                     %fgm app input c2
        nameset(i,14)=cellstr(['mwc3bet_',input_name,',1']);                     %fgm app input c3
        nameset(i,15)=cellstr(['masked_mwc1bet_',input_name]);                   %fgm app output c1
        nameset(i,16)=cellstr(['masked_mwc2bet_',input_name]);                   %fgm app output c2
        nameset(i,17)=cellstr(['eroded_mwc3bet_',input_name]);                   %fgm app output c3
        fgm_exp=[fgm_exp,'+i',num2str(i)];
    end
        fgm_exp_sum=['(',fgm_exp(2:end),')'];
        fgm_exp=['(',fgm_exp(2:end),')==',num2str(number_of_files)];
end

%% BRAIN EXTRACTION (bet_)   
% VBM analysis pre-processing: step 1
%
% input:        sub_xxx_T1w.nii;                    /filename/
% outputs:      bet_sub_xxx_T1w.nii; bet_sub_xxx_T1w_mask.nii

if bet_flag == 1
    for i=1:number_of_files
        input_name = filename(i,:);
        output_name = ['bet_', input_name];
        cmd_out = ['bet ', input_name, ' ', output_name, ' -m'];
        system(cmd_out);
        cmd_out = ['gunzip ', output_name, '.gz'];
        system(cmd_out);
        cmd_out = ['gunzip ', output_name(1:end-4), '_mask.nii.gz'];
        system(cmd_out);        
    end
end
   
%% SEGMENTATION (c1-3) in FSL
% VBM analysis pre-procesing: step 2
% 
% write rc* files in 0.7 mm resolution
% 
% input:        bet_sub*.nii                            /nameset(;,1)/
%               /+SPM templates for GM, WM, CSF/
% outputs:      c?bet*.nii, rc?bet*.nii, mwc?bet*.nii, mbet*.nii 

if seg_flag == 1
    
    clear matlabbatch
    
    matlabbatch{1}.spm.spatial.preproc.channel.vols = nameset(:,1);
    matlabbatch{1}.spm.spatial.preproc.channel.biasreg = 1.0000e-03;
    matlabbatch{1}.spm.spatial.preproc.channel.biasfwhm = 60;
    matlabbatch{1}.spm.spatial.preproc.channel.write = [0 1];
    
    matlabbatch{1}.spm.spatial.preproc.tissue(1).tpm={'/usr/local/spm12/toolbox/spmmouse/tpm/greyr62_x10.nii,1'};
    matlabbatch{1}.spm.spatial.preproc.tissue(2).tpm={'/usr/local/spm12/toolbox/spmmouse/tpm/whiter62_x10.nii,1'};
    matlabbatch{1}.spm.spatial.preproc.tissue(3).tpm={'/usr/local/spm12/toolbox/spmmouse/tpm/csfr62_x10.nii,1'};
    matlabbatch{1}.spm.spatial.preproc.tissue(1).ngaus=1;
    matlabbatch{1}.spm.spatial.preproc.tissue(2).ngaus=1;
    matlabbatch{1}.spm.spatial.preproc.tissue(3).ngaus=2;
    matlabbatch{1}.spm.spatial.preproc.tissue(1).native=[1 1];
    matlabbatch{1}.spm.spatial.preproc.tissue(2).native=[1 1];
    matlabbatch{1}.spm.spatial.preproc.tissue(3).native=[1 1];
    matlabbatch{1}.spm.spatial.preproc.tissue(1).warped=[0 1];
    matlabbatch{1}.spm.spatial.preproc.tissue(2).warped=[0 1];
    matlabbatch{1}.spm.spatial.preproc.tissue(3).warped=[0 1];
    
    matlabbatch{1}.spm.spatial.preproc.warp.mrf = 1;
    matlabbatch{1}.spm.spatial.preproc.warp.cleanup = 1;
    matlabbatch{1}.spm.spatial.preproc.warp.reg = [0 0.0010 0.5000 0.0500 0.2000];
    matlabbatch{1}.spm.spatial.preproc.warp.affreg = '';
    matlabbatch{1}.spm.spatial.preproc.warp.fwhm = 0;
    matlabbatch{1}.spm.spatial.preproc.warp.samp = 3;
    matlabbatch{1}.spm.spatial.preproc.warp.write = [0 0];
    
    matlabbatch{1}.spm.spatial.preproc.warp.vox     =  0.7;
    
    spm_preproc_run(matlabbatch{1}.spm.spatial.preproc,'run')
end

%% DARTEL: template creation
% VBM analysis pre-processing: step 3
%
%input:         rc1bet*.nii, rc2bet*.nii                /nameset(i,2-3)/
% outputs:      u_rc1bet*.nii
%               template_?.nii /7x/

if dartemp_flag == 1
    clear matlabbatch
    
    matlabbatch{1}.spm.tools.dartel.warp.images{1}=nameset(:,2);
    matlabbatch{1}.spm.tools.dartel.warp.images{2}=nameset(:,3);
    
    matlabbatch{1}.spm.tools.dartel.warp.settings.template='Template';
    matlabbatch{1}.spm.tools.dartel.warp.settings.rform=0;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(1).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(2).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(3).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(4).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(5).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(6).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(7).its=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(1).rparam=[4.0000    2.0000    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(2).rparam=[2.0000    1.0000    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(3).rparam=[1.0000    0.5000    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(4).rparam=[0.5000    0.2500    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(5).rparam=[0.2500    0.1250    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(6).rparam=[0.2500    0.1250    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(7).rparam=[0.1000    0.0700    0.0000];
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(1).K=0;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(2).K=0;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(3).K=1;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(4).K=2;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(5).K=4;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(6).K=6;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(7).K=7;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(1).slam=16;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(2).slam=8;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(3).slam=4;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(4).slam=2;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(5).slam=1;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(6).slam=0.5000;
    matlabbatch{1}.spm.tools.dartel.warp.settings.param(7).slam=0;
    
    matlabbatch{1}.spm.tools.dartel.warp.settings.optim.lmreg=0.0100;
    matlabbatch{1}.spm.tools.dartel.warp.settings.optim.cyc=3;
    matlabbatch{1}.spm.tools.dartel.warp.settings.optim.its=3;
    
    spm_jobman('run',matlabbatch);
end

%% DARTEL warp: normalisation to template6
% VBM analysis pre-processing: step 4
%
% use u_rc1 flow field maps 
% 
% input:    u_rc1bet*_Template.nii                      /nameset(:,4)/
%           mbet*.nii                                   /nameset(:,5)/
%           c?bet*.nii                                  /nameset(:,6-8)/          
% =same procedure runs 4 times

if dartwarp_flag == 1
    
    clear matlabbatch
    
    matlabbatch{1}.spm.tools.dartel.crt_warped.flowfields=nameset(:,4);
    matlabbatch{1}.spm.tools.dartel.crt_warped.images{1}=nameset(:,5);
    matlabbatch{1}.spm.tools.dartel.crt_warped.images{2}=nameset(:,6);
    matlabbatch{1}.spm.tools.dartel.crt_warped.images{3}=nameset(:,7);
    matlabbatch{1}.spm.tools.dartel.crt_warped.images{4}=nameset(:,8);
    matlabbatch{1}.spm.tools.dartel.crt_warped.jactransf=1;
    matlabbatch{1}.spm.tools.dartel.crt_warped.K=6;
    matlabbatch{1}.spm.tools.dartel.crt_warped.interp=5;
    
    spm_jobman('run',matlabbatch);
end

%% MASK normalisation (A), final_group_mask formation (B), erosion (C)
% (A)********************************************************************
% MASK NORMALISATION
% input:        u_rc1bet*_Template.nii                      /nameset(:,4)/
%               bet*_mask.nii                               /nameset(;,9)/
% output:       wbet*_mask.nii
% = same DARTEL warp procedure as for c1-3 and mbet files, 
% but WITHOUT MODULATION /wbet instead of mwbet/

if masknorm_flag == 1
    
    clear matlabbatch
    
    matlabbatch{1}.spm.tools.dartel.crt_warped.flowfields=nameset(:,4);
    matlabbatch{1}.spm.tools.dartel.crt_warped.images{1}=nameset(:,9);
    matlabbatch{1}.spm.tools.dartel.crt_warped.jactransf=0;
    matlabbatch{1}.spm.tools.dartel.crt_warped.K=6;
    matlabbatch{1}.spm.tools.dartel.crt_warped.interp=5;
    
    spm_jobman('run',matlabbatch);
end

% (B)********************************************************************
% final_group_mask FORMATION
% input:        u_rc1bet*_Template.nii                      /nameset(:,4)/
%               wbet*_mask.nii                              /nameset(;,10)/
% output:       final_group_mask.nii

if fgm_flag == 1
    
    clear matlabbatch
    
    matlabbatch{1}.spm.util.imcalc.input=nameset(:,10);
    matlabbatch{1}.spm.util.imcalc.output='final_group_mask';
    %matlabbatch{1}.spm.util.imcalc.outdir='';
    matlabbatch{1}.spm.util.imcalc.expression=fgm_exp;
    %matlabbatch{1}.spm.util.imcalc.var='';
    matlabbatch{1}.spm.util.imcalc.options.dmtx=0;
    matlabbatch{1}.spm.util.imcalc.options.mask=0;
    matlabbatch{1}.spm.util.imcalc.options.interp=1;
    matlabbatch{1}.spm.util.imcalc.options.dtype=4;
    
    spm_jobman('run',matlabbatch);
end

% *sumed final group mask to decipher proportion of animals on regional loss of info in final_group_mask 
if fgmsum_flag == 1
    
    clear matlabbatch
    
    matlabbatch{1}.spm.util.imcalc.input=nameset(:,10);
    matlabbatch{1}.spm.util.imcalc.output='final_group_mask_sum';
    %matlabbatch{1}.spm.util.imcalc.outdir='';
    matlabbatch{1}.spm.util.imcalc.expression=fgm_exp_sum;
    %matlabbatch{1}.spm.util.imcalc.var='';
    matlabbatch{1}.spm.util.imcalc.options.dmtx=0;
    matlabbatch{1}.spm.util.imcalc.options.mask=0;
    matlabbatch{1}.spm.util.imcalc.options.interp=1;
    matlabbatch{1}.spm.util.imcalc.options.dtype=4;
    
    spm_jobman('run',matlabbatch);
end

% (C)********************************************************************
% EROSION
% input:        final_group_mask.nii
% output:       final_group_mask_ero.nii

if ero_flag == 1
    cmd_out = ['fslmaths final_group_mask.nii -ero -kernel sphere 8 final_group_mask_ero.nii'];
    system(cmd_out);
    cmd_out = ['gunzip final_group_mask_ero.nii.gz'];
    system(cmd_out);
end

%% SMOOTHING
% input:             mwc1bet_*.nii                            /nameset(:,11)/
%
% output:            smwc1bet_*.nii

if sm_flag == 1
    clear matlabbatch
    
    matlabbatch{1}.spm.spatial.smooth.data=nameset(:,11); %input images
    matlabbatch{1}.spm.spatial.smooth.fwhm=[8 8 8];       %full width at half maximum
    matlabbatch{1}.spm.spatial.smooth.dtype=0;            %data type: 0=same as input
    matlabbatch{1}.spm.spatial.smooth.im=0;               %implicit masking
    matlabbatch{1}.spm.spatial.smooth.prefix='s';         %filename prefix
    
    spm_jobman('run',matlabbatch);
end

%% fgm on mw?
% input:        mwc1bet*.nii                               /nameset(:,12)/
%               mwc2bet*.nii                               /nameset(;,13)/
%               mwc3bet*.nii                               /nameset(;,14)/
%               final_group_mask.nii; ..._mask_ero.nii     /c1 and c2; c3/
% output:       masked_mwc1*.nii                           /nameset(;,15)/
%               masked_mwc2*.nii                           /nameset(;,16)/
%               eroded_mwc3*.nii                           /nameset(;,17)/

if fgmapp_flag ==1
    clear matlabbatch
    
    %matlabbatch{1}.spm.util.imcalc.outdir=;
    matlabbatch{1}.spm.util.imcalc.expression='i1.*i2';
    %matlabbatch{1}.spm.util.imcalc.var=;
    matlabbatch{1}.spm.util.imcalc.options.dmtx=0;      %data matrix (y/n)
    matlabbatch{1}.spm.util.imcalc.options.mask=1;      %masking (non-implicit/implicit/NaN=0)
    matlabbatch{1}.spm.util.imcalc.options.interp=-5;   %interpolation (nn/3lin/xdegree sinc)
    matlabbatch{1}.spm.util.imcalc.options.dtype=16;    %data type: as input images = float32
    for j=1:3
        filetype=['mwc', num2str(j)];
        for i=1:number_of_files
            fgm_output=char(nameset(i,14+j));
            if strfind(filetype,'3')
                fgm_pair = [nameset(i,11+j);"final_group_mask_ero.nii,1"];
            else
                fgm_pair = [nameset(i,11+j);"final_group_mask.nii,1"];
            end
            matlabbatch{1}.spm.util.imcalc.input= cellstr(fgm_pair);
            matlabbatch{1}.spm.util.imcalc.output=fgm_output;
            spm_jobman('run',matlabbatch);
        end
    end
end

%% GLOBAL CALCULATION: VOLUMES
% input:        masked_mwc1bet*.nii                        /nameset(:,15)/
%               masked_mwc2bet*.nii                        /nameset(;,16)/
%               eroded_mwc3bet*.nii                        /nameset(;,17)/
% output:       Vols [n of animals, n of compartments]     

if vols_flag==1
    V = spm_vol(nameset(:,15:17));
    Vsize=size(V);
    Vols = zeros(Vsize(1),Vsize(2));
    for k=1:Vsize(2)
        for j=1:Vsize(1)
            tot = 0;
            for i=1:V{j,k}.dim(3)
                img = spm_slice_vol(V{j,k},spm_matrix([0 0 i]),V{j,k}.dim(1:2),0);
                img = img(isfinite(img)); % <-- exclude non-finite values
                tot = tot + sum(img(:));
            end
            voxvol = abs(det(V{j,k}.mat))/100^3; % volume of a voxel, in mililitres
            Vols(j,k) = tot*voxvol;
        end
    end
    save('Vols','Vols')
    dlmwrite('Vols.txt', Vols, 'delimiter','\t','newline','pc')
end
%% GLOBAL CALCULATION: FACTORIAL DESIGN