Skip to content

Commit

Permalink
add preprocessing scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@FeHoff
FeHoff committed Aug 29, 2024
0 parents commit 6a3e25c
Show file tree
Hide file tree
Showing 41 changed files with 7,323 additions and 0 deletions.
208 changes: 208 additions & 0 deletions CAT125sub_all.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,208 @@
function CAT125sub_all(sub_raw, xsub, TMPdir, targetdir, force)
%%%
% CAT12 segmentation using SHOOTING & DARTEL
%
% script based on:
% $Id: cat_defaults_expert.m 895 2016-03-11 16:07:09Z gaser $
% $Id: cat_defaults.m 1183 2017-09-08 16:45:08Z dahnke $
%
% Felix Hoffstaedter
% INM-7 - Brain and Behaviour
% Data and Platforms
% March 2018 - Research Centre Juelich
%
% CATsub_all( <subjectfolder>, <subjectname>, <TEMPdirectory>, <targetdirectory>, <force processing> )
%%%

warning off all

addpath /data/BnB2/TOOLS/spm12_cat12.5

% --------------- start script only if T1w exists ---------------
if exist(sub_raw,'file')

% --- how many subjects are processed in parallel [in window = 0 ] --
jobs = 0; % a job will use min 4 cores for i.e. denoising
% --- how many cores are use for fractal dimension estimation --
cores = 0; % 1 if multiple subjects are processed
% ------- smoothing of modulated volume data im mm [5 & 8 default] --
vbmFWHM = [5 8];
% ---------- estimate fractal dimensions - SLOW ! ------------------
FD = 1;
% ---------- smoothing of cortical thickness in mm [>=15] -----------
thkFWHM = [0 15];
% ---------- smoothing of folding parameters in mm [>=25] -----------
srfFWHM = [0 25];
%--------------- load cat12 expert defaults -------------------------
def = fullfile(spm_file(which('CAT125sub_all'), 'path'), 'cat125_defaults_bnb');
try global cat; if cat.extopts.expertgui == 0; spm fmri; cat12(def); end
catch; spm fmri; cat12(def); end

sub_nii = fullfile(TMPdir, [xsub '.nii']); % Subjects T1 file for VBM

% --- delete old data if force = true ---
if force == 1
try rmdir(TMPdir,'s'); end
end

if ~exist(sub_nii,'file')
% ---------- copy & extract nifti file to CAT dir -----------------
mkdir(TMPdir);

if strfind(sub_raw,'run-02')
% --- realignment and mean image creation if multiple T1w acquired ---
% ----- write loop for more than 2 acquisitions per scan session -----
copyfile(sub_raw, fullfile(TMPdir, [xsub '-1.nii.gz']));
gunzip(fullfile(TMPdir, [xsub '-1.nii.gz']));
delete(fullfile(TMPdir, [xsub '-1.nii.gz']));
sub_raw2 = strrep(sub_raw,'run-01','run-02');
copyfile(sub_raw2, fullfile(TMPdir, [xsub '-2.nii.gz']));
gunzip(fullfile(TMPdir, [xsub '-2.nii.gz']));
delete(fullfile(TMPdir, [xsub '-2.nii.gz']));
matlabbatch = {};
matlabbatch{1}.spm.spatial.coreg.estimate.ref = {fullfile(TMPdir, [xsub '-1.nii'])};
matlabbatch{1}.spm.spatial.coreg.estimate.source = {fullfile(TMPdir, [xsub '-2.nii'])};
spm_jobman('run', matlabbatch)
matlabbatch = {};
matlabbatch{1}.spm.util.imcalc.input = {fullfile(TMPdir, [xsub '-1.nii'])
fullfile(TMPdir, [xsub '-2.nii'])};
matlabbatch{1}.spm.util.imcalc.output = [xsub '.nii'];
matlabbatch{1}.spm.util.imcalc.outdir = {TMPdir};
matlabbatch{1}.spm.util.imcalc.expression = 'mean(X)';
matlabbatch{1}.spm.util.imcalc.options.dmtx = 1;
spm_jobman('run', matlabbatch)
else
if strcmp(spm_file(sub_raw,'ext'),'gz')
copyfile(sub_raw, fullfile(TMPdir, [xsub '.nii.gz']));
gunzip(fullfile(TMPdir, [xsub '.nii.gz']));
delete(fullfile(TMPdir, [xsub '.nii.gz']));
elseif strcmp(spm_file(sub_raw,'ext'),'nii')
copyfile(sub_raw, sub_nii);
else
return
end
end

% for robustness set (0 0 0) coordinate to middle of image
V = spm_vol(sub_nii);
% pre-estimated COM of MNI template
com_reference = [0 -20 -15];
Affine = eye(4);
vol = spm_read_vols(V);
avg = mean(vol(:));
avg = mean(vol(vol>avg));
% don't use background values
[x,y,z] = ind2sub(size(vol),find(vol>avg));
com = V.mat(1:3,:)*[mean(x) mean(y) mean(z) 1]';
com = com';
M = spm_get_space(V.fname);
Affine(1:3,4) = (com - com_reference)';
spm_get_space(V.fname,Affine\M);

end

% ---------- high-dimensional SHOOTING & DARTEL normalization -----------------
if ~exist(fullfile(TMPdir, 'mri', ['m0wp1' xsub '.nii']), 'file') || force == 1
matlabbatch = {};
matlabbatch{1}.spm.tools.cat.estwrite.nproc = jobs;
% data to process
matlabbatch{1}.spm.tools.cat.estwrite.data{1} = [sub_nii ',1'];
% run things
spm_jobman('run', matlabbatch)
end

% ---------- Smooth modulated gray & white matter segments --------------------
if ~exist(fullfile(TMPdir, 'mri', ['s' num2str(vbmFWHM(1)) 'm0wp1' xsub '.nii']),'file') ...
|| force == 1
matlabbatch = {};
matlabbatch{1}.spm.spatial.smooth.data = {
fullfile(TMPdir, 'mri', ['mwp1' xsub '.nii'])
fullfile(TMPdir, 'mri', ['m0wp1' xsub '.nii'])};
for i = 1:numel(vbmFWHM)
matlabbatch{1}.spm.spatial.smooth.fwhm = [vbmFWHM(i) vbmFWHM(i) vbmFWHM(i)];
matlabbatch{1}.spm.spatial.smooth.prefix = ['s' num2str(vbmFWHM(i))];
spm_jobman('run', matlabbatch)
end
end

% ---------- Estimate total intracranial volume (TIV) & -------------------
% ---------- global GM, WM, CSF, WM-hyperintensity vol --------------------
if (~exist(fullfile(TMPdir, ['TIV_' xsub '.txt']),'file') || force == 1) && ...
exist(fullfile(TMPdir, 'report', ['cat_' xsub '.xml']), 'file'); % only with existing 'report/cat_subj.mat'
matlabbatch = {};
matlabbatch{1}.spm.tools.cat.tools.calcvol.data_xml = ...
{fullfile(TMPdir, 'report', ['cat_' xsub '.xml'])};
matlabbatch{1}.spm.tools.cat.tools.calcvol.calcvol_TIV = 0;
matlabbatch{1}.spm.tools.cat.tools.calcvol.calcvol_name = [ 'TIV_' xsub '.txt'];
spm_jobman('run', matlabbatch);
movefile(fullfile(pwd, ['TIV_' xsub '.txt']), fullfile(TMPdir, ['TIV_' xsub '.txt']))
end

if exist(fullfile(TMPdir, 'surf', ['lh.central.' xsub '.gii']),'file')
% ---------- Extract gyrification, cortical complexity and sulcus depth --------------------
if ~exist(fullfile(TMPdir, 'surf', ['rh.sqrtsulc.' xsub]),'file') || force == 1
matlabbatch = {};
matlabbatch{1}.spm.tools.cat.stools.surfextract.data_surf = {
fullfile(TMPdir, 'surf', ['lh.central.' xsub '.gii'])};
matlabbatch{1}.spm.tools.cat.stools.surfextract.GI = 1;
matlabbatch{1}.spm.tools.cat.stools.surfextract.FD = FD;
matlabbatch{1}.spm.tools.cat.stools.surfextract.SD = 1;
matlabbatch{1}.spm.tools.cat.stools.surfextract.nproc = cores;
spm_jobman('run', matlabbatch)
end

% ---------- Resample and smooth cortical thickness --------------------
if ~exist(fullfile(TMPdir, 'surf', ['s' num2str(thkFWHM(1)) '.rh.thickness.resampled.' xsub '.gii']),'file') ...
|| force == 1
matlabbatch = {};
matlabbatch{1}.spm.tools.cat.stools.surfresamp.data_surf = {
fullfile(TMPdir, 'surf', ['lh.thickness.' xsub])};
matlabbatch{1}.spm.tools.cat.stools.surfresamp.nproc = cores;
for i = 1:numel(thkFWHM)
matlabbatch{1}.spm.tools.cat.stools.surfresamp.fwhm_surf = thkFWHM(i);
matlabbatch{1}.spm.tools.cat.stools.surfresamp.merge_hemi = 0;
matlabbatch{1}.spm.tools.cat.stools.surfresamp.mesh32k = 0;
spm_jobman('run', matlabbatch)
matlabbatch{1}.spm.tools.cat.stools.surfresamp.merge_hemi = 1;
matlabbatch{1}.spm.tools.cat.stools.surfresamp.mesh32k = 1;
spm_jobman('run', matlabbatch)
end
end

% ---------- Resample and smooth gyrification and sulcus depth --------------------
if ~exist(fullfile(TMPdir, 'surf', ['s' num2str(srfFWHM(1)) '.rh.sqrtsulc.resampled.' xsub '.gii']),'file') ...
|| force == 1
matlabbatch = {};
if FD == 0
matlabbatch{1}.spm.tools.cat.stools.surfresamp.data_surf = {
fullfile(TMPdir, 'surf', ['lh.gyrification.' xsub])
fullfile(TMPdir, 'surf', ['lh.sqrtsulc.' xsub])};
else
matlabbatch{1}.spm.tools.cat.stools.surfresamp.data_surf = {
fullfile(TMPdir, 'surf', ['lh.gyrification.' xsub])
fullfile(TMPdir, 'surf', ['lh.sqrtsulc.' xsub])
fullfile(TMPdir, 'surf', ['lh.fractaldimension.' xsub])};
end
matlabbatch{1}.spm.tools.cat.stools.surfresamp.nproc = cores;
for i = 1:numel(srfFWHM)
matlabbatch{1}.spm.tools.cat.stools.surfresamp.fwhm_surf = srfFWHM(i);
matlabbatch{1}.spm.tools.cat.stools.surfresamp.merge_hemi = 0;
matlabbatch{1}.spm.tools.cat.stools.surfresamp.mesh32k = 0;
spm_jobman('run', matlabbatch)
matlabbatch{1}.spm.tools.cat.stools.surfresamp.merge_hemi = 1;
matlabbatch{1}.spm.tools.cat.stools.surfresamp.mesh32k = 1;
spm_jobman('run', matlabbatch)
end
end
end

if exist(fullfile(TMPdir, 'surf', ['s' num2str(srfFWHM(numel(srfFWHM))) ...
'.mesh.sqrtsulc.resampled_32k.' xsub '.gii']),'file')
mkdir(targetdir);
movefile(fullfile(TMPdir,'*'),targetdir);
rmdir(TMPdir);
end

end
exit
end
189 changes: 189 additions & 0 deletions RS_resmooth.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,189 @@
function RS_resmooth(sub_raw, FMdir_raw, xsub, cat_dir, targetdir, FMdir, RSdir, TR, FWHM, force)
cwd = pwd;

addpath /data/BnB2/TOOLS/spm12
addpath /data/BnB2/TOOLS/DVARS-master/

FWHM = 5;

origdir = fullfile(RSdir,'Orig'); % Unprocessed resting-state EPI folder
stdir = fullfile(RSdir,'ST'); % Slicetimed data after realigment
normdir = fullfile(RSdir,'Norm'); % Normalized resting-state EPI folder, 3D
smoothdir = fullfile(RSdir,['Smooth_' int2str(FWHM) 'mm']); % Normalized resting-state EPI folder
T1dir = cat_dir; % 3D CAT folder


if exist(fullfile(RSdir,['sw' xsub '.nii']),'file')
movefile(fullfile(RSdir,['sw' xsub '.nii']),fullfile(RSdir,['s5w' xsub '.nii']))
end

if ~exist(normdir,'dir') || force == 1
% ---------- split 4D NIFTI ---------------
mkdir(normdir);
matlabbatch = {};
matlabbatch{1}.spm.util.split.vol = {fullfile(RSdir,['w' xsub '.nii'])};
matlabbatch{1}.spm.util.split.outdir = {normdir};
spm_jobman('run',matlabbatch)
end

if exist(normdir,'dir') && ~exist(smoothdir,'dir') || force == 1
% ---------- Smooth --------------------
matlabbatch = {};
matlabbatch{1}.spm.spatial.smooth.dtype = 0;
matlabbatch{1}.spm.spatial.smooth.im = 0;
matlabbatch{1}.spm.spatial.smooth.prefix = 's';
matlabbatch{1}.spm.spatial.smooth.fwhm = [FWHM FWHM FWHM];
fils = dir(fullfile(normdir,'w*.nii')); fils = fils(~[fils.isdir]);
for file = 1:size(fils,1)
matlabbatch{1}.spm.spatial.smooth.data{file,1} = fullfile(normdir,fils(file).name);
end
spm_jobman('run',matlabbatch)
mkdir(smoothdir);
movefile(fullfile(normdir,'sw*.nii'), smoothdir)
end

if ~exist(fullfile(RSdir,['s' num2str(FWHM) 'w' xsub '.nii']),'file') && ...
exist(smoothdir,'dir')
% merge 3D to 4D Nifti
matlabbatch = {};
fils = dir(fullfile(smoothdir,['sw*.nii']));
for file = 1:size(fils,1)
matlabbatch{1}.spm.util.cat.vols{file,1} = fullfile(smoothdir,fils(file).name);
end
matlabbatch{1}.spm.util.cat.name = fullfile(RSdir,['s' num2str(FWHM) 'w' xsub '.nii']);
matlabbatch{1}.spm.util.cat.dtype = 4;
spm_jobman('run',matlabbatch);
end

% if exist(fullfile(RSdir,['s' num2str(FWHM) 'w' xsub '.nii']),'file')
% try rmdir(smoothdir,'s'); end
% end

% if exist(fullfile(T1dir,'mri',['wp0' xsub '.nii']),'file') && ...
% ~exist(fullfile(RSdir,['wp0' xsub '.nii']),'file')
% % resample partial volume label image to 2x2x2mm
% matlabbatch = {};
% matlabbatch{1}.spm.spatial.coreg.write.ref = ...
% {fullfile(RSdir,['w' xsub '_meanEPI.nii'])};
% matlabbatch{1}.spm.spatial.coreg.write.source = ...
% {fullfile(T1dir,'mri',['wp0' xsub '.nii'])};
% matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 0;
% matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0];
% matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0;
% matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r';
% spm_jobman('run',matlabbatch)
% movefile(fullfile(T1dir,'mri',['rwp0' xsub '.nii']),fullfile(RSdir,['wp0' xsub '.nii']))
% end
% %
% if exist(fullfile(T1dir,'mri',['mwp1' xsub '.nii']),'file') && ...
% ~exist(fullfile(RSdir,['mwp1' xsub '.nii']),'file')
% % resample modulated GM segment image to 2x2x2mm
% matlabbatch = {};
% matlabbatch{1}.spm.spatial.coreg.write.ref = ...
% {fullfile(RSdir,['w' xsub '_meanEPI.nii'])};
% matlabbatch{1}.spm.spatial.coreg.write.source = ...
% {fullfile(T1dir,'mri',['mwp1' xsub '.nii'])};
% matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 7; % 7th degree B-Spline
% matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0];
% matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0;
% matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r';
% spm_jobman('run',matlabbatch)
% movefile(fullfile(T1dir,'mri',['rmwp1' xsub '.nii']),fullfile(RSdir,['mwp1' xsub '.nii']))
% end
%
% if ~exist(fullfile(RSdir,['Counfounds_' xsub '.mat']),'file') || force == 1 || force == 2
% clear pXYZ tX tY tZ gx gx2 rp rpp A dat msk gx2 reg y
% % getting & computing movement parameter (FD & sqrt) from realignment parameter
% txt = dir(fullfile(RSdir,'rp_*.txt'));
% rp = load(fullfile(RSdir,txt.name));
% [FDts,FD_Stat]=FDCalc(rp);
% rp = rp-repmat(mean(rp),size(rp,1),1);
% drp = [zeros(1,6); diff(rp)];
% move = sqrt(sum(drp(:,1:3).^2,2));
% % compute root mean squared movement & frame wise displacement
% euler = zeros(size(rp,1),1);
% for ii=1:size(rp,1)
% euler(ii) = 50*acos((cos(drp(ii,4))*cos(drp(ii,5)) + cos(drp(ii,4))*cos(drp(ii,6)) + ...
% cos(drp(ii,5))*cos(drp(ii,6)) + sin(drp(ii,4))*sin(drp(ii,6))*sin(drp(ii,5)) - 1)/2);
% end
% RMS = sqrt(nanmean([move+euler].^2));
% tmp = sum(abs([drp(:,1:3) drp(:,4:6)*50]),2);
% FD = sqrt(nanmean(tmp.^2));
% % data file
% % Vi = spm_vol(fullfile(fsmoothdir,['sw' xsub '.nii']));
% Vi = spm_vol(fullfile(RSdir,['w' xsub '.nii']));
% % computing global, GM, WM, CSF, mean signals
% % getting partial volume mask
% msk = spm_read_vols(spm_vol(fullfile(RSdir,['wp0' xsub '.nii'])));
% % msk: min 0.0157; max 2.9961; 1=CSF, 2=GM, 3=WM
% idxGM = find(msk>1.99 & msk<2.01);
% % GMmsk = msk;
% % GMmsk(msk<1.99 | msk>2.01) = 0;
% % GMmsk = spm_erode(GMmsk); % erode once
% % % GMmsk = spm_erode(spm_erode(GMmsk)); % erode twice
% % idxGM = find(GMmsk>0);
% % idxWM = find(msk>2.99);
% idxCSF = find(msk>0.99 & msk<1.01);
% gx = nan(4,numel(Vi));
% % i.e. aCompCor does >.99 WM 2x erode & >.99 CSF NN clustering criteria
% WMmsk = msk;
% WMmsk(msk<2.99) = 0;
% WMmsk = spm_erode(WMmsk); % erode once
% % WMmsk = spm_erode(spm_erode(WMmsk)); % erode twice
% idxWM = find(WMmsk>0);
% datWM = zeros(length(idxWM),numel(Vi));
% datCSF = zeros(length(idxCSF),numel(Vi));
% datALL = zeros(length(find(msk>0)),numel(Vi));
% for i=1:numel(Vi)
% dat = spm_read_vols(Vi(i));
% gx(1,i) = mean(dat(idxGM));
% datWM(:,i) = dat(idxWM); % WM voxels only
% gx(2,i) = mean(datWM(:,i));
% datCSF(:,i) = dat(idxCSF); % CSF voxels only
% gx(3,i) = mean(datCSF(:,i));
% datALL(:,i) = dat(msk>0); % All brain voxels
% gx(4,i) = mean(dat(msk>0));
% end
%
% PCA_WM = pca(datWM - repmat(mean(datWM),size(datWM,1),1));
% PCA_CSF = pca(datCSF - repmat(mean(datCSF),size(datCSF,1),1));
%
% gx2(1,:) = gx(1,:)-mean(gx(1,:));
% gx2(2,:) = gx(2,:)-mean(gx(2,:));
% gx2(3,:) = gx(3,:)-mean(gx(3,:));
% gx2(4,:) = gx(4,:)-mean(gx(4,:));
%
% reg = [gx2' gx2'.^2 rp rp.^2 drp drp.^2 PCA_WM(:,1:7) PCA_CSF(:,1:7)];
% reg = reg-repmat(mean(reg),size(reg,1),1);
% reg = reg./repmat(std(reg),size(reg,1),1);
%
% [DVARS,DVARS_Stat]=DVARSCalc(datALL,'scale',1/100,'RDVARS');
% badTP=find(DVARS_Stat.pvals<0.05./(numel(DVARS)));
% [V,DSE_Stat]=DSEvars(datALL,'scale',1/100);
%
% % FIGURE with BOLD intensity carpet-plot
% % fMRIDiag_plot(V,DVARS_Stat,'BOLD',datALL,'ColRng',[-3 3],'FD',FDts,'AbsMov',[FD_Stat.AbsRot FD_Stat.AbsTrans]);
% % Diagnosis plot based on stored variables
% % f_hdl=figure('position',[50,500,600,600]);
% % fMRIDiag_plot(V,DVARS_Stat,'FD',FDts,'AbsMov',[FD_Stat.AbsRot FD_Stat.AbsTrans],'TickScaler',0.5,'figure',f_hdl);
%
% save(fullfile(RSdir,['Counfounds_' xsub '.mat']),'reg','rp','drp','gx',...
% 'gx2','TR','RMS','FD','DVARS','DVARS_Stat','V','DSE_Stat','FDts','FD_Stat','badTP')
%
% end

% if exist(fullfile(RSdir,['w' xsub '.nii']),'file')
% try rmdir(origdir,'s'); end
% try rmdir(stdir,'s'); end
% try rmdir(normdir,'s'); end
% try rmdir(smoothdir,'s'); end
% end

% if exist(fullfile(RSdir,['Counfounds_' xsub '.mat']),'file')
% mkdir(targetdir);
% movefile(fullfile(RSdir,'s*'),targetdir);
% rmdir(RSdir);
% end

% exit
end
Loading

0 comments on commit 6a3e25c

Please sign in to comment.