Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F97005344
Preprocessing.m
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jan 1, 13:38
Size
10 KB
Mime Type
text/x-Algol68
Expires
Fri, Jan 3, 13:38 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
23310846
Attached To
rISFCCLUS ISFC clustering
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');
Event Timeline
Log In to Comment