diff --git a/JOVE_Preprocessing.m b/JOVE_Preprocessing.m index c17233b..1da88ff 100644 --- a/JOVE_Preprocessing.m +++ b/JOVE_Preprocessing.m @@ -1,257 +1,355 @@ %% This script contains the main preprocessing steps required to run the % ISFC transient excursion detection methodology. In more details, the % following steps are applied to the data: % % - Detrending of the data % - Regressing out covariates (constant,linear, quadratic terms, average % CSF and average WM signals) % - Atlasing into selected atlas % - Scrubbing % +% +% +% +% % Written by Thomas Bolton, last modified on July 30th 2018 - %% Parameters to be changed by the user %%%%%%%%% Move to before the "Subjects" folder - - - % Name of all the subject folders through which to operate Subjects_OI = {'S17'}; % Prefix of the folders containing the different sessions for a subject SES_Prefix = 'SES'; % Name of the folder containing all the subjects folders SFold_Name = 'Subjects'; % End of the file name of one functional file for a subject Func_Suffix = '001-01.hdr'; Func_Prefix = 'rf'; % File type (functional files), between 'nii' and 'hdr' Func_Type = 'hdr'; % Prefix of the folder containing the structural data Struct_Prefix = 'S'; Atlas_Prefix = 'wCraddock_300'; Save_Name = 'JOVE_Example'; % Framewise displacement threshold past which data is considered corrupted % (according to Power's criterion) Mot_Threshold = 0.5; %% Start of the preprocessing code +% Turn to true if you wish to observe different outputs from the script +is_plot = 1; + +% Turn to true if you wish to display different text summaries of the +% preprocessing steps +is_text = 1; + % Current folder (at the start of the folder tree) is added with all % subfolders Now = pwd; addpath(genpath(pwd)); % Gets to the subjects-containing folder cd(SFold_Name); % Iterates across subjects... for i = 1:length(Subjects_OI) - disp(['Preprocessing subject ',num2str(i),'/',num2str(length(Subjects_OI)),'...']); + if is_text + disp(['Preprocessing subject ',num2str(i),'/',num2str(length(Subjects_OI)),'...']); + end % Moves to the subject folder cd(Subjects_OI{i}); % Looks for the directories of interest containing the different % sessions to look at (all start with 'S') DOR = dir([SES_Prefix,'*']); % Iterates across sessions for j = 1:length(DOR) - disp(['Processing session ',num2str(j),'/',num2str(length(DOR)),'...']); + if is_text + disp(['Processing session ',num2str(j),'/',num2str(length(DOR)),'...']); + end % Gets to the session directory of interest cd(DOR(j).name); % Locates one functional file from the session in question tmp_refname = dir(['*',Func_Suffix]); % Gets the full path of that file REF = fullfile(pwd,tmp_refname(1).name); % Looks for the directory containing structural information and % segmentation results SDIR = dir([Struct_Prefix,'*']); % Goes to the structural directory cd(SDIR(1).name); % Gets the link towards the native space Craddock atlas Atlas_link = cellstr(spm_select('List',pwd,['^', Atlas_Prefix, '.*\.', 'nii', '$'])); % Makes this link a full path Atlas_link = fullfile(pwd,Atlas_link); + if is_text + disp('Creating WM and CSF masks...'); + end + % Gets the links towards the native space WM and CSF probabilistic % maps (given the c2 and c3 prefixes by SPM12) WM_link = cellstr(spm_select('List',pwd,['^' 'c2' '.*\.' 'nii' '$'])); CSF_link = cellstr(spm_select('List',pwd,['^' 'c3' '.*\.' 'nii' '$'])); % Converts the resolution to be of right size WM = mapVolumeToVolume(WM_link{1},REF); CSF = mapVolumeToVolume(CSF_link{1},REF); % Binarises the two volumes in very conservative manner WM(WM < 0.99) = 0; WM(WM >= 0.99) = 1; WM = logical(WM); CSF(CSF < 0.99) = 0; CSF(CSF >= 0.99) = 1; CSF = logical(CSF); + if is_text + disp(['There are ',num2str(sum(sum(sum(WM)))), ' WM voxels',... + ' and ',num2str(sum(sum(sum(CSF)))),' CSF voxels...']); + end + % Now goes to the functional data folder to perform the operations % themselves, i.e. detrending and regression (and atlasing) cd('..'); % Gets the links of all functional files FLinks = dir([Func_Prefix,'*',Func_Type]); % Number of frames n_frames = length(FLinks); + if is_text + disp(['There are ',num2str(n_frames),' frames...']); + end + + if is_text + disp('Constructing data matrix...'); + end + % Will now read the data for that session frame by frame, and % create the WM and CSF regressors for d = 1:n_frames tmp_V = spm_read_vols(spm_vol(FLinks(d).name)); % Regressors for WM and CSF are constructed (average of the % retained voxels) CSF_cov(d) = mean(tmp_V(CSF)); WM_cov(d) = mean(tmp_V(WM)); % Data matrix construction (n_tp x n_vox) Data(d,:) = tmp_V(:); clear tmp_V end % Demeaning of regressors CSF_cov = CSF_cov - mean(CSF_cov); WM_cov = WM_cov - mean(WM_cov); %% Detrending + if is_text + disp('Detrending...'); + end + % Data detrending Data = detrend(Data); % Construction of the matrix of regressors: constant, linear, % quadratic, CSF and WM (no motion because scrubbing is done % afterwards) X = [ones(n_frames,1),(1:n_frames)',((1:n_frames).^2)',CSF_cov',WM_cov']; - clear n_frames - clear CSF_cov - clear WM_cov - % Will contain the preprocessed data Data_preproc = nan(size(Data)); %% Regressing out covariates + if is_text + disp('Regressing out covariates...'); + end + % Regresses voxel after voxel (DPARSFA toolbox function) with also % adding back the mean as in the lab scripts for vox = 1:size(Data,2) [~,Data_preproc(:,vox)] = y_regress_ss(Data(:,vox),X); Data_preproc(:,vox) = Data_preproc(:,vox) + mean(Data(:,vox)); end clear Data clear X + %% Atlasing + + if is_text + disp('Atlasing...'); + end + % Converts atlas into right resolution, and makes it 1D ATLAS = mapVolumeToVolume(Atlas_link{1},REF); ATLAS = ATLAS(:); ATLAS(isnan(ATLAS)) = 0; - - %% Atlasing - + % Atlases the data - TC = Make_Atlasing(Data_preproc',ATLAS); + [TC,n_disc_TC,n_kept_TC] = Make_Atlasing(Data_preproc',ATLAS); clear ATLAS + %% Scrubbing + + if is_text + disp('Scrubbing...'); + end + % I now want to read motion parameters motfile_name = dir(['rp','*','.txt']); motfile_name = motfile_name(1).name; % My 6 motion parameters Mot = textread(motfile_name); Mot = Mot(:,1:6); % Converts the rotational components into [mm] Mot(:,4:6) = 50*Mot(:,4:6); % Computes FD FD = sum(abs([0 0 0 0 0 0; diff(Mot)]),2); % Locates the time points to scrub (set to 1) Scrub_mask1 = logical(FD > Mot_Threshold); % Locates the time point that immediately follows as well Scrub_mask2 = circshift(Scrub_mask1,[1,0]); Scrub_mask2(1) = 0; % Makes an OR operations (sets Scrub_mask element to 1 as soon as % one of Scrub_mask1/Scrub_mask2 is equal to 1) Scrub_mask = Scrub_mask1 | Scrub_mask2; + if is_text + disp(['There are ',num2str(100*sum(Scrub_mask)/n_frames),... + ' frames to srub...']); + end + % Finds the values at the time points that we know are not % corrupted TCon = TC(:,find(~Scrub_mask)); TCon(isnan(TCon)) = 0; % tinter is all the time points tinter = 1:length(Scrub_mask); % torig is the time points that we know torig = tinter(find(~Scrub_mask)); % We interpolate all time point values using the info that we % know (TCon at torig) TC2 = interp1(torig,TCon',tinter,'spline')'; % We want to keep the 'NaN's that we had because they are used to % remove bad regional time courses later on TC2(isnan(TC)) = NaN; + % Plotting of different aspects of the preprocessing + if is_plot + + % Plot of WM/CSF time courses that were regressed out, and + % motion parameters + figure; + set(gcf,'color','w'); + subplot(2,1,1); + plot(Mot(:,1),'color',[0 0 0.6],'LineWidth',2); + hold on; + title(['Subject ',num2str(i),', session ',num2str(j)]); + plot(Mot(:,2),'color',[0 0.4 1],'LineWidth',2); + plot(Mot(:,3),'color',[0.6 0.6 1],'LineWidth',2); + plot(Mot(:,4),'color',[0.8 0 0.2],'LineWidth',2); + plot(Mot(:,5),'color',[1 0.2 0.2],'LineWidth',2); + plot(Mot(:,6),'color',[1 0.6 0.2],'LineWidth',2); + legend('X','Y','Z','Alpha','Beta','Gamma'); + xlabel('Time'); + ylabel('Motion [mm]'); + set(gca,'Box','off'); + + subplot(2,1,2); + plot(CSF_cov,'color',[0.2 0.6 0],'LineWidth',2); + hold on; + plot(WM_cov,'color',[0.6 0.6 0.6],'LineWidth',2); + legend('CSF','WM'); + xlabel('Time'); + ylabel('Intensity [a.u.]'); + set(gca,'Box','off'); + + figure; + set(gcf,'color','w'); + subplot(2,1,1); + hold on; + set(gca,'Box','off'); + bar(100*n_disc_TC/(n_disc_TC+n_kept_TC)); + xlabel('Regions'); + ylabel('Number of discarded time courses'); + ylim([0 110]); + subplot(2,1,2); + hold on; + set(gca,'Box','off'); + bar(100*n_kept_TC/(n_disc_TC+n_kept_TC)); + xlabel('Regions'); + ylabel('Number of kept time courses'); + ylim([0 110]); + end + clear Mot clear torig clear tinter clear TCon clear Scrub_mask clear FD + clear CSF_cov + clear WM_cov + clear n_frames % Saves the outputs as cell arrays of time courses TCS{j,i} = TC2; cd('..'); end cd('..'); end cd('..'); % Saves the data in the 'Results' folder save(fullfile(pwd,'Results',Save_Name),'TCS'); diff --git a/Results/JOVE_Example.mat b/Results/JOVE_Example.mat index 8d49666..8dc341a 100644 Binary files a/Results/JOVE_Example.mat and b/Results/JOVE_Example.mat differ