Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/Computations_ISFC_Excursions.m
b90a4357cb97master

Last Change
Raw File

Computations_ISFC_Excursions.m
View Options

%% This script runs the ISFC computations and ISFC transient extraction on
% a set of typically developing subjects, for a movie-watching paradigm
%
% Code written by Thomas Bolton, last checked on July 24th 2018
%% 1. Data loading, paths generation and parameters
% Adds the path towards all utilities
addpath(genpath(fullfile(pwd,'Utilities')));
%%%%%%%%%%%%%%%%%% To be entered by the user %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name of the atlas file
Atlas_name = 'Craddock_300.nii';
% Name of the data file
Data_name = 'JOVE_Full_Data.mat';
% Name of the file containing the output data
Save_name = 'JOVE_Results';
% TR for the considered dataset
TR=2.0;
% Subjects that were included in the atlasing step a priori (i.e., name of
% the related subject folders)
Subjects_OI = {'S1','S9','S10','S14','S16','S17','S19','S20','S21','S23','S25',...
'S26','S27','S31','S32'};
Ses_prefix = 'SES';
% Threshold to use for motion selection (scrubbing), according to Power's
% criterion; advised value is 0.5 [mm]
Mot_Threshold = 0.5;
% Percentage frame that can maximally be scrubbed before the subject is
% taken out; advised value is 10 [%], can go up to 30 [%]
TScrub = 10;
% Size of the sliding windows (in TRs); advised value is W = 10 [TR]
W = 10;
% Step between two windows (in TRs); advised value is Step = 1 [TR]
Step = 1;
% Indices of the sessions, from the ones to analyze, combining RS and
% movie-watching subparts
combined_sessions = [1,2];
% Timing (in samples from the data matrix) of the movie-watching and RS
% subparts for those sessions
MW_timing = 3:178;
RS_timing = 194:340;
% Number of bootstrapping folds to use for ISFC computations
n_boot = 250;
% Number of reference subjects to use for the bootstraping process
n_refsubj = 6;
% alpha-value to use for ISFC thresholding
Alpha = 0.01;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Path towards the atlas
PathToAtlas = fullfile(pwd,'Utilities',Atlas_name);
% Path towards the folder will all my subjects
PathToSubjFolder = fullfile(pwd,'Subjects');
% Loads the data (called 'TCS' from the previous script)
load(fullfile(pwd,'Results',Data_name));
% Number of different sessions acquired per subject
n_sessions = size(TCS,1);
% For each session, we assign a label
Subject_Label = repmat(Subjects_OI,size(TCS,1),1);
Subject_Label = Subject_Label(:);
% Please switch to true if desiring to plot some figures of the process
is_plot = 1;
% Please switch to true if desiring to plot some text about the process
is_text = 1;
%% 2. Removal of bad quality sessions
% Number of subjects
n_subj = length(Subjects_OI);
% Label for session (1 or 2)
Session_Label = repmat(1:n_sessions,1,length(Subjects_OI));
% Computation of motion for all subjects
[~,~,PScrub] = ComputeFD(PathToSubjFolder,Subjects_OI,Mot_Threshold,Ses_prefix);
% We keep subjects if they show less than 30% of time points with large
% framewise displacement
ToKeep = (PScrub <= TScrub);
if is_text
disp([num2str(sum(~ToKeep)),'/',num2str(length(ToKeep)),' sessions to discard...']);
end
if is_plot
figure;
set(gcf,'color','w');
bar(1:length(PScrub),PScrub);
set(gca,'Box','off');
xlabel('Sessions');
ylabel('Percentage scrubbed frames [%]');
xlim([0.5 length(PScrub)+0.5]);
ylim([0 50]);
hold on;
plot([1 length(PScrub)],[TScrub TScrub],'k--');
end
Session_Label = Session_Label(ToKeep);
Subject_Label = Subject_Label(ToKeep);
% I convert my 3 x 30 data into a 1 x 90 format in order to keep only the
% sessions of interest
TCS_all = reshape(TCS,1,max(Session_Label)*size(TCS,2));
TCS_all = TCS_all(ToKeep);
clear PScrub
% Removes the regions in which at least one of the subjects has NaN values
% (i.e. the region did not exist). We use all the sessions for this step
[TCS_ip,idx_removedROIs] = RemoveNaNRegions(TCS_all);
if is_text
disp([num2str(length(idx_removedROIs)),' regions removed from the atlas...']);
end
% Filtering of the time courses at cutoff 1/W [s]
TCS_rear = cell(1,length(TCS_ip));
for i = 1:length(TCS_ip)
TCS_rear{i} = y_IdealFilter(TCS_ip{i}', TR, [1/(W*TR) 0]);
TCS_rear{i} = TCS_rear{i}';
end
%% 3. Atlas-related details
% We want a vector indexing the ROIs, going from 1 to the number of regions
% kept (not NaN)
ROI_vec = 1:size(TCS_all{1},1);
ROI_vec(idx_removedROIs) = [];
ROI_vec = 1:length(ROI_vec);
% Number of regions and of connections
n_ROI = length(ROI_vec);
n_conn = n_ROI*(n_ROI-1)/2;
if is_text
disp(['Total number of regions: ',num2str(n_ROI),'...']);
disp(['Total number of connections: ',num2str(n_conn),'...']);
end
% Converts the region vector into a connection vector (an example entry is
% '1-13' between ROI 1 and ROI 13)
Con_vec = GetConVec(ROI_vec);
% I gather the labels of the regions involved in each connection (for the
% above example, C1(i) = 1 and C2(i) = 13)
[C1,C2] = ParseConVec(Con_vec);
% Creation of the atlas codebook and of 3D center of gravity information,
% with COG of size 3 x n_ROI
[CodeBook,COG] = CreateCodeBook(PathToAtlas,ROI_vec,idx_removedROIs);
% Euclidean distance between all connections considered
for i = 1:n_conn
D_euclidean(i,1) = sqrt(sum((squeeze(COG(:,C1(i)))-squeeze(COG(:,C2(i)))).^2));
end
if is_plot
figure;
set(gcf,'color','w');
hist(D_euclidean,1000);
xlim([0 200]);
set(gca,'Box','off');
xlabel('Distance between regions [mm]');
ylabel('Number of connections [-]');
end
%% 4. Extraction of resting-state and movie-watching portions
idx_RSMW = 1;
idx_RS = 1;
% We iterate through all the preprocessed data
for i = 1:length(TCS_rear)
% tmp contains the time courses (n_ROI x n_TPs) for one subject
tmp = TCS_rear{i};
% If the current recording is a combined stimulus + resting-state case,
% then we sample both subparts
switch (ismember(Session_Label(i),combined_sessions))
case 0
TC_RS2{idx_RS} = tmp;
idx_RS = idx_RS + 1;
case 1
TC_RS{idx_RSMW} = tmp(:,RS_timing);
TC_MW{idx_RSMW} = tmp(:,MW_timing);
idx_RSMW = idx_RSMW + 1;
end
clear tmp
end
if is_text
disp(['Number of retained pure resting-state sessions: ',num2str(length(TC_RS2)),'...']);
disp(['Number of retained stimulus-related sessions: ',num2str(length(TC_RS)),'...']);
end
%% 5. Computing ISFC estimates on the data
% First, we wish to solely perform ISFC compuations on the movie-watching
% subparts; we create the ISFC_MW cell array, each cell about to contain
% the outputs for one subject (n_conn x n_slidingwindow_timepoint); then,
% we do the same for the resting-state subparts from movie sessions, and
% finally for the 'pure resting-state' sessions
if is_text
disp('Computing ISFC on stimulus-related sessions (stimulus segments)...');
end
[ISFC_MW] = Perform_Bootstrapping(TC_MW,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
if is_text
disp('Computing ISFC on stimulus-related sessions (resting-state segments)...');
end
[ISFC_RS] = Perform_Bootstrapping(TC_RS,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
if is_text
disp('Computing ISFC on purely resting-state sessions...');
end
[ISFC_RS2] = Perform_Bootstrapping(TC_RS2,Session_Label(~ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);
%% 6. Significance thresholding of connectivity pairs
if is_text
disp('Thresholding ISFC time courses...');
end
% First, I want to determine when in time, and for what connections, I
% observe a significant excursion above from the mean connectivity
% SignMat has size n_conn x n_TP x n_sessions, and contains:
% -1 if there is a negative excursion
% 0 if there is no excursion
% +1 if there is a positive excursion
SignMat = FindExcursions(ISFC_MW,ISFC_RS,ISFC_RS2,Alpha);
% Saves the data
save(fullfile(pwd,'Results',Save_name),'CodeBook','TC_RS','TC_MW',...
'TC_RS2','ISFC_MW','ISFC_RS','ISFC_RS2','SignMat','Alpha','W','Step',...
'n_boot','n_refsubj','TR','Mot_Threshold','TScrub','-v7.3');
c4science · Help