Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/Preprocessing.m
b90a4357cb97master

Last Change
Raw File

Preprocessing.m
View Options

%% 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)
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)
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'];
% 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;
% Atlases the data
[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(n_disc_TC);
xlabel('Regions');
ylabel('Number of discarded time courses');
subplot(2,1,2);
hold on;
set(gca,'Box','off');
bar(n_kept_TC);
xlabel('Regions');
ylabel('Number of kept time courses');
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');
c4science · Help