Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/Matlab_Scripts/Script_Vertes_Computations_Mean.m
257e25cb59ccmaster

Last Change
Raw File

Script_Vertes_Computations_Mean.m
View Options

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% So, this is supposed to be the script implementing the analyses for the
% MW paper and its revisions
%
% Last updated: August 3rd 2017, Thomas Bolton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Data loading and path generation
disp('Loading data and generating paths...');
% We add the paths towards functions of interest for the analysis, namely:
% General utilities (plotting...)
addpath(genpath(fullfile(pwd,'Utils')));
% Dictionary learning
addpath(genpath(fullfile(pwd,'Utils_DL')));
% K-means clustering
addpath(genpath(fullfile(pwd,'Utils_KMeans')));
% Graph analysis
addpath(genpath(fullfile(pwd,'Utils_Graphs')));
% Graph visualization tools (toolbox advised by Prof. Vandergheynst)
addpath(genpath(fullfile(pwd,'gspbox')));
% Scripts required to visualize 3D brains with nodes and edges
addpath(genpath('/Volumes/Code/shared/static_FC_pipeline/depMat'));
addpath(genpath('/media/miplab-nas1/Code/shared/static_FC_pipeline/depMat'));
addpath(genpath('/Volumes/Code/shared/static_FC_pipeline/graph_vis'));
addpath(genpath('/media/miplab-nas1/Code/shared/static_FC_pipeline/graph_vis'));
% Atlas used in this study
addpath(genpath(fullfile(pwd,'Craddock_Atlas')));
% Path towards the atlas
PathToAtlas = fullfile(pwd,'Craddock_300.nii');
% Path towards the folder will all my subjects
PathToSubjFolder = '/Volumes/Lorena/Research projects/3. Autism_Spectrum_Disorders_Idiosyncrasy/HBM_MW_NewPrepro/Subjects';
load TCS_real
All_TCs = TCS_real;
clear TCS_real
% Is there or not a resting-state scanning session for the subjects ?
is_RS = [0 0 3 4 5 0 0 0 9 10 11 12 13 14 15 16 17 18 19 20 21 0 23 24 25 26 27 28 29 30];
%% 2. Preprocessing, reordering of nodes, exclusion of bad sessions
disp('Preprocessing...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A. Fixed analytical parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TR for the considered dataset
TR=2.0;
% Subjects that were included in the atlasing step a priori
Subjects_OI = {'S1','S2','S3','S4','S5','S6','S7','S8','S9','S10','S12',...
'S13','S14','S15','S16','S17','S18','S19','S20','S21','S22',...
'S23','S24','S25','S26','S27','S28','S29','S30','S31','S32'};
% Number of subjects
n_subj = length(Subjects_OI);
% Label for session (1 or 2)
Session_Label = repmat([1 2 3],1,length(Subjects_OI));
% Group label: 1 = CTR, 2 = ASD, NaN = no session available
Group_Label = [1 1 NaN, 2 2 NaN, 2 2 2, 2 2 2, 1 1 1, 2 2 NaN, 2 2 NaN, 2 2 NaN, 1 1 1, 1 1 1, 2 2 2, 2 2 2, 1 1 1, 2 2 2,...
1 1 1, 1 1 1, 2 2 2, 1 1 1, 1 1 1, 1 1 1, 2 2 2, 1 1 1, 2 2 NaN, 1 1 1, 1 1 1, 1 1 1, 2 2 2, 2 2 2, 2 2 2, 1 1 1, 1 1 1];
% For each subject, we assign a label (recall that there are three sessions
% per subject)
Subject_Label = {'S1','S1','S1','S2','S2','S2','S3','S3','S3','S4','S4','S4','S5','S5','S5','S6','S6','S6',...
'S7','S7','S7','S8','S8','S8','S9','S9','S9','S10','S10','S10','S12','S12','S12','S13','S13','S13','S14','S14','S14',...
'S15','S15','S15','S16','S16','S16','S17','S17','S17','S18','S18','S18','S19','S19','S19','S20','S20','S20','S21','S21','S21',...
'S22','S22','S22','S23','S23','S23','S24','S24','S24','S25','S25','S25','S26','S26','S26','S27','S27','S27',...
'S28','S28','S28','S29','S29','S29','S30','S30','S30','S31','S31','S31','S32','S32','S32'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B. Tunable analytical parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Threshold to use for motion selection (scrubbing)
MT = 0.5;
% Percentage frame that can maximally be scrubbed before the subject is
% taken out
%%%%%%%%%%%%%%% WE USE 30 PERCENT HERE, STILL A DECENT VALUE OVERALL
TScrub = 30;
% Type of connectivity measure to compute
conn='Pearson';
% Size of the sliding windows (in TRs)
W = 10;
% Step between two windows (in TRs)
Step = 1;
% Do I want to plot my information ? 0 = no, 1 = only fast plots, 2 = also
% brain graphs, 3 = also generates videos
is_plot = 0;
% Number of top connections to show in the brain graphs
n_top_connections = 100;
% Number of top ROIs to display in the brain graphs
n_top_regions = 20;
% Distance below which we say that a connection is short-range (used in the
% graph part); for now we use 80mm because it is the cutoff used in Achard
% and Bullmore, J Neuroscience
long_distance_range_threshold = 80;
short_distance_range_threshold = 80;
% Number of times we run the null process for graph extraction of time
% points with significant group differences
n_null_graph = 10;
% Percentage of frames that I want to retain per subject for k-means
perc_ret_frames = 100;
% Number of eigencomponents to look at
n_eig_patterns = 50;
% Threshold above which we consider a connection to be significant (non
% normalized, to assess shared ones across significant clusters)
T_kmeans_con = 0.4;
% Threshold for what value should be the minimal connectivity represented
% in a cluster (on normalized data from -1 to 1 separately for positive and
% negative connections)
T_display_kmeans = 0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C. Removal of corrupted sessions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computation of motion for all subjects
[~,~,PScrub] = ComputeFD(PathToSubjFolder,Subjects_OI,MT);
% After visual inspection, no case looked totally off using atlasing in the
% most recent way
ToKeepVis = logical([1 1 0, 1 1 0, 1 1 1, 1 1 1, 1 1 1, 1 1 0, 1 1 0, 1 1 0, ...
1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1,...
1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 0, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, ...
1 1 1, 1 1 1, 1 1 1]);
% We keep subjects if they show less than 30% of time points with large
% framewise displacement
ToKeepScr = (PScrub <= TScrub);
ToKeep = ToKeepVis & ToKeepScr;
Session_Label = Session_Label(ToKeep);
Subject_Label = Subject_Label(ToKeep);
Group_Label = Group_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(All_TCs,1,93);
TCS_all = TCS_all(ToKeep);
% Number of CTR and ASD subjects (computed only on the MW+RS sessions)
n_CTR = sum(Group_Label == 1 & Session_Label < 3);
n_ASD = sum(Group_Label == 2 & Session_Label < 3);
clear TCS_INFO
clear BinMat
clear PScrub
clear FD
% 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);
% Filtering of the time courses at cutoff 1/W [s]
TCS_ip_filt = cell(1,length(TCS_ip));
for i = 1:length(TCS_ip)
TCS_ip_filt{i} = y_IdealFilter(TCS_ip{i}', TR, [1/(W*TR) 0]);
TCS_ip_filt{i} = TCS_ip_filt{i}';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D. Reordering of the ROIs so that they are arranged topologically
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We want a vector indexing the ROIs, going from 1 to the number of regions
% kept (not NaN)
ROI_vec = 1:300;
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;
% 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
% I construct an adjacency matrix describing the distance between my
% regions; values close to 1 mean a very small distance, while values close
% to 0 mean a very large distance. I use the method taught by Prof. Xavier
% Bresson (former EPFL member), where sigma is the average distance across
% connections
Adj_distance = zeros(n_ROI,n_ROI);
sigma = mean(D_euclidean);
for i = 1:n_conn
Adj_distance(C1(i),C2(i)) = exp(-(D_euclidean(i)^2)./sigma^2);
end
Adj_distance = Adj_distance + Adj_distance';
% Computation of the normalized Laplacian of the graph, and of its
% eigenvalues (E) and eigenvectors (V)
L = compute_graph_laplacian(Adj_distance);
L = L + 1e-6* speye(size(L));
[V,E] = eigs(L,4,'SM');
E = diag(E);
% Eigenvalues and eigenvectors are sorted from the smallest to the largest
% eigenvalue
[~,idx] = sort(E,'ascend');
E = E(idx);
V = V(:,idx);
% The Fiedler vector gives us information about the topological arrangement
% of the nodes
Dist_manifold = V(:,2);
% IDX_ATLAS has the ROI indices sorted in a topologically relevant order
[DM2,IDX_ATLAS] = sort(Dist_manifold,'ascend');
% Of course, I then need to rearrange all that was already computed
% accordingly, to do what will follow only with the newest ROI ordering. I
% thus need to swap:
% - the codebook information
COG = COG(:,IDX_ATLAS);
for i = 1:n_ROI
CodeBook.full.center{i} = COG(:,i);
end
% - the euclidean distance information
clear D_euclidean
for i = 1:n_conn
D_euclidean(i,1) = sqrt(sum((squeeze(COG(:,C1(i)))-squeeze(COG(:,C2(i)))).^2));
end
% - the time courses
TCS_rear = cell(1,length(TCS_ip_filt));
for i = 1:length(TCS_ip_filt)
TCS_rear{i} = TCS_ip_filt{i}(IDX_ATLAS,:);
end
if is_plot > 0
G = gsp_graph(Adj_distance(IDX_ATLAS,IDX_ATLAS));
G.coords = COG';
gsp_plot_signal(G,DM2);
figure;
set(gcf,'color','w');
plot(ROI_vec,COG(2,:),'k','LineWidth',2);
xlabel('ROI index');
ylabel('MNI Y coordinate value');
set(gca,'Box','off');
figure;
set(gcf,'color','w');
imagesc(Adj_distance(IDX_ATLAS,IDX_ATLAS));
axis square
caxis([0 1]);
xlabel('ROI index');
ylabel('ROI index');
colormap(gray);
end
%% 3. Extracting resting-state and movie-watching portions
disp('Extraction of resting-state and movie-watching portions...');
% Here, we must pay particular attention to the stimulus timing, because
% there were some subject cases for which the movie did not start exactly
% at the same time
idx_RSMW = 1;
idx_RS = 1;
for i = 1:length(TCS_rear)
% tmp contains the time courses (n_ROI x n_TPs) for one subject
tmp = TCS_rear{i};
% If we are in the combined sessions...
if Session_Label(i) < 3
current_subj = Subject_Label{i};
current_ses = Session_Label(i);
current_full = [current_subj,'_SES',num2str(current_ses)];
% Depending on the subject, we extract different subportions because in
% a few cases, the movie started at a slightly different time.
%
% Because the movie started in the middle of frame 3 for most cases,
% but at the start of frames for the others, we want to interpolate
% what happens between two frames for those subcases, to have exactly
% the same time courses across subjects. We thus perform spline
% interpolation on this data. We interpolate the limited subset of
% subjects with different movie timings so that most of the data is
% left unchanged
switch current_full
% Movie from scan 4 (start) to scan 177 (end)
case {'S1_SES1','S6_SES1','S16_SES1','S25_SES1','S22_SES2'}
disp('Exception 1');
tmp2 = tmp(:,3:179);
tx = 3:179;
ty = 3.5:178.5;
tmp3 = spline(tx,tmp2,ty);
% Movie from scan 6 (start) to scan 179 (end)
case 'S22_SES1'
disp('Exception 2');
tmp2 = tmp(:,5:181);
tx = 5:181;
ty = 5.5:180.5;
tmp3 = spline(tx,tmp2,ty);
% Movie from scan 11 (start) to scan 184 (end)
case 'S24_SES1'
disp('Exception 3');
tmp2 = tmp(:,10:186);
tx = 10:186;
ty = 10.5:185.5;
tmp3 = spline(tx,tmp2,ty);
% Normal case, movie from scan 3 (half) to scan 177 (half)
otherwise
tmp3 = tmp(:,3:178);
end
TC_RS{idx_RSMW} = tmp(:,194:end);
TC_MW{idx_RSMW} = tmp(:,3:178);
idx_RSMW = idx_RSMW + 1;
clear tmp
clear tmp2
clear tmp3
else
TC_RS2{idx_RS} = tmp;
idx_RS = idx_RS + 1;
clear tmp
end
end
% For what follows, I don't care about my pure RS sessions anymore, so I
% want Subject_Label, Group_Label and Session_Label to only contain
% information about the combined sessions
Subject_Label = Subject_Label(Session_Label < 3);
Group_Label_pureRS = Group_Label(Session_Label == 3);
Group_Label = Group_Label(Session_Label < 3);
Session_Label = Session_Label(Session_Label < 3);
%% 4. Computing sliding window estimates on the data
disp('Sliding window computations...');
rmpath(genpath(fullfile(pwd,'Utils','spm12')));
% We first select only CTR or ASD sessions, because we want to apply the
% ISFC framework to each group separately
tmp_RS = TC_RS(Group_Label == 1);
tmp_MW = TC_MW(Group_Label == 1);
tmp_RS2 = TC_RS(Group_Label == 2);
tmp_MW2 = TC_MW(Group_Label == 2);
% Classical sliding window computation for each subject, for each type of
% session (for this, we don't need to separate groups, because the
% computations only depend on the particular subject presently considered)
dFC_RS = SlidingWindow(TC_RS,W,Step,'Pearson','FC');
dFC_MW = SlidingWindow(TC_MW,W,Step,'Pearson','FC');
% Now, we want to do the same computations, but comparing the subjects from
% each of our groups to a reference subgroup every time
% Number of bootstrapping folds
n_boot = 500;
% Number of reference subjects to use; ideally we want to have half covering
% session 1, and half covering session 2
n_refsubj = 6;
% I have my MW and my RS bits for the CTR and for the ASD subjects. In each
% case, what I will do is bootstrap many many times by comparing to a
% reduced reference group. And then, I will average the outcomes together.
%%%%%%
% 1. For the MW bits (TD and ASD groups)
%%%%%%
ISFC_MW1 = cell(1,sum(Group_Label == 1));
[ISFC_MW1{:}] = deal(zeros(size(dFC_MW{1})));
ISFC_MW2 = cell(1,sum(Group_Label == 2));
[ISFC_MW2{:}] = deal(zeros(size(dFC_MW{1})));
% The same, but for stationary ISFC
sISFC_MW1 = cell(1,sum(Group_Label == 1));
[sISFC_MW1{:}] = deal(zeros(270,270));
sISFC_MW2 = cell(1,sum(Group_Label == 2));
[sISFC_MW2{:}] = deal(zeros(270,270));
% To keep track of what to divide by for each subject
Denom_factor_TD = cell(1,sum(Group_Label == 1));
[Denom_factor_TD{:}] = deal(0);
Denom_factor_ASD = cell(1,sum(Group_Label == 2));
[Denom_factor_ASD{:}] = deal(0);
tmp_denom_TD = cell(1,sum(Group_Label == 1));
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD = cell(1,sum(Group_Label == 2));
[tmp_denom_ASD{:}] = deal(0);
% First we need to determine the indices of SES1 and of SES2 recordings
% for both the TD and the ASD groups
idx_TD = 1:sum(Group_Label==1);
idx_TD_SES1 = idx_TD(Session_Label(Group_Label==1)==1);
idx_TD_SES2 = idx_TD(Session_Label(Group_Label==1)==2);
idx_ASD = 1:sum(Group_Label==2);
% For each fold...
for boot = 1:n_boot
boot
% Next, we need to select n_refsubj/2 elements from SES1 and from SES2
% elements, so first we pick the indices in a random order
idx_TD_SES1_b = randperm(length(idx_TD_SES1));
idx_TD_SES2_b = randperm(length(idx_TD_SES2));
% This has the indices of all the subjects that are picked as reference
idx_ref_TD = sort([idx_TD_SES1(idx_TD_SES1_b(1:n_refsubj/2)),idx_TD_SES2(idx_TD_SES2_b(1:n_refsubj/2))],'ascend');
idx_notref_TD = 1:sum(Group_Label==1);
idx_notref_TD(idx_ref_TD) = [];
% Next, we select the subjects that make the reference group for the
% particular bootstrapping fold
tmp_ref_TD = tmp_MW(idx_ref_TD);
tmp_notref_TD = tmp_MW(idx_notref_TD);
% We can now do the computations themselves
% With ISFC, we compute the outputs for all TD and ASD subjects...
[ISFC_MW1_notref,sISFC_MW1_notref] = ISFC_ToReferenceGroup(tmp_notref_TD,tmp_ref_TD,W,Step);
[ISFC_MW2_notref,sISFC_MW2_notref] = ISFC_ToReferenceGroup(tmp_MW2,tmp_ref_TD,W,Step);
% Updates the general cell (ISFC)
tmp_cell = cell(1,sum(Group_Label==1));
[tmp_cell{:}] = deal(zeros(size(ISFC_MW1_notref{1})));
tmp_cell(idx_notref_TD) = ISFC_MW1_notref;
ISFC_MW1 = cellfun(@(x,y) x+y,tmp_cell,ISFC_MW1,'UniformOutput',0);
clear tmp_cell
% Same for sISFC
tmp_cell = cell(1,sum(Group_Label==1));
[tmp_cell{:}] = deal(zeros(270,270));
tmp_cell(idx_notref_TD) = sISFC_MW1_notref;
sISFC_MW1 = cellfun(@(x,y) x+y,tmp_cell,sISFC_MW1,'UniformOutput',0);
clear tmp_cell
% ISFC for ASD
tmp_cell = cell(1,sum(Group_Label==2));
[tmp_cell{:}] = deal(zeros(size(ISFC_MW2_notref{1})));
tmp_cell(idx_ASD) = ISFC_MW2_notref;
ISFC_MW2 = cellfun(@(x,y) x+y,tmp_cell,ISFC_MW2,'UniformOutput',0);
clear tmp_cell
% sISFC for ASD
tmp_cell = cell(1,sum(Group_Label==2));
[tmp_cell{:}] = deal(zeros(270,270));
tmp_cell(idx_ASD) = sISFC_MW2_notref;
sISFC_MW2 = cellfun(@(x,y) x+y,tmp_cell,sISFC_MW2,'UniformOutput',0);
clear tmp_cell
% Keeps track of how many times each subject has been included as a non
% reference (i.e. how many times a value has been computed for it)
tmp_denom_TD(idx_notref_TD) = {1};
Denom_factor_TD = cellfun(@(x,y) x+y,tmp_denom_TD,Denom_factor_TD,'UniformOutput',0);
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD(idx_ASD) = {1};
Denom_factor_ASD = cellfun(@(x,y) x+y,tmp_denom_ASD,Denom_factor_ASD,'UniformOutput',0);
[tmp_denom_ASD{:}] = deal(0);
end
% Averages the results; we divide by the number of times a given subject has
% been included
ISFC_MW1 = cellfun(@(x,y) x/y,ISFC_MW1,Denom_factor_TD,'UniformOutput',0);
ISFC_MW2 = cellfun(@(x,y) x/y,ISFC_MW2,Denom_factor_ASD,'UniformOutput',0);
sISFC_MW1 = cellfun(@(x,y) x/y,sISFC_MW1,Denom_factor_TD,'UniformOutput',0);
sISFC_MW2 = cellfun(@(x,y) x/y,sISFC_MW2,Denom_factor_ASD,'UniformOutput',0);
% Finally, we get back to a representation including both the TD and the
% ASD subjects
ISFC_MW = cell(1,length(dFC_MW));
ISFC_MW(Group_Label==1) = ISFC_MW1;
ISFC_MW(Group_Label==2) = ISFC_MW2;
sISFC_MW = cell(1,length(dFC_MW));
sISFC_MW(Group_Label==1) = sISFC_MW1;
sISFC_MW(Group_Label==2) = sISFC_MW2;
%%
%%%%%%
% 2. For the RS bits (TD and ASD groups)
%%%%%%
ISFC_RS1 = cell(1,sum(Group_Label == 1));
[ISFC_RS1{:}] = deal(zeros(size(dFC_RS{1})));
ISFC_RS2 = cell(1,sum(Group_Label == 2));
[ISFC_RS2{:}] = deal(zeros(size(dFC_RS{1})));
% To keep track of what to divide by for each subject
Denom_factor_TD = cell(1,sum(Group_Label == 1));
[Denom_factor_TD{:}] = deal(0);
Denom_factor_ASD = cell(1,sum(Group_Label == 2));
[Denom_factor_ASD{:}] = deal(0);
tmp_denom_TD = cell(1,sum(Group_Label == 1));
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD = cell(1,sum(Group_Label == 2));
[tmp_denom_ASD{:}] = deal(0);
% First we need to determine the indices of SES1 and of SES2 recordings
% for both the TD and the ASD groups
idx_TD = 1:sum(Group_Label==1);
idx_TD_SES1 = idx_TD(Session_Label(Group_Label==1)==1);
idx_TD_SES2 = idx_TD(Session_Label(Group_Label==1)==2);
idx_ASD = 1:sum(Group_Label==2);
% For each fold...
for boot = 1:n_boot
boot
% Next, we need to select n_refsubj/2 elements from SES1 and from SES2
% elements, so first we pick the indices in a random order
idx_TD_SES1_b = randperm(length(idx_TD_SES1));
idx_TD_SES2_b = randperm(length(idx_TD_SES2));
% This has the indices of all the subjects that are picked as reference
idx_ref_TD = sort([idx_TD_SES1(idx_TD_SES1_b(1:n_refsubj/2)),idx_TD_SES2(idx_TD_SES2_b(1:n_refsubj/2))],'ascend');
idx_notref_TD = 1:sum(Group_Label==1);
idx_notref_TD(idx_ref_TD) = [];
% Next, we select the subjects that make the reference group for the
% particular bootstrapping fold
tmp_ref_TD = tmp_RS(idx_ref_TD);
tmp_notref_TD = tmp_RS(idx_notref_TD);
% We can now do the computations themselves
% With ISFC, we compute the outputs for all TD and ASD subjects...
ISFC_RS1_notref = ISFC_ToReferenceGroup(tmp_notref_TD,tmp_ref_TD,W,Step);
ISFC_RS2_notref = ISFC_ToReferenceGroup(tmp_RS2,tmp_ref_TD,W,Step);
% Updates the general cell
tmp_cell = cell(1,sum(Group_Label==1));
[tmp_cell{:}] = deal(zeros(size(ISFC_RS1_notref{1})));
tmp_cell(idx_notref_TD) = ISFC_RS1_notref;
ISFC_RS1 = cellfun(@(x,y) x+y,tmp_cell,ISFC_RS1,'UniformOutput',0);
clear tmp_cell
tmp_cell = cell(1,sum(Group_Label==2));
[tmp_cell{:}] = deal(zeros(size(ISFC_RS2_notref{1})));
tmp_cell(idx_ASD) = ISFC_RS2_notref;
ISFC_RS2 = cellfun(@(x,y) x+y,tmp_cell,ISFC_RS2,'UniformOutput',0);
clear tmp_cell
% Keeps track of how many times each subject has been included as a non
% reference (i.e. how many times a value has been computed for it)
tmp_denom_TD(idx_notref_TD) = {1};
Denom_factor_TD = cellfun(@(x,y) x+y,tmp_denom_TD,Denom_factor_TD,'UniformOutput',0);
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD(idx_ASD) = {1};
Denom_factor_ASD = cellfun(@(x,y) x+y,tmp_denom_ASD,Denom_factor_ASD,'UniformOutput',0);
[tmp_denom_ASD{:}] = deal(0);
end
% Averages the results; we divide by the number of times a given subject has
% been included
ISFC_RS1 = cellfun(@(x,y) x/y,ISFC_RS1,Denom_factor_TD,'UniformOutput',0);
ISFC_RS2 = cellfun(@(x,y) x/y,ISFC_RS2,Denom_factor_ASD,'UniformOutput',0);
% Finally, we get back to a representation including both the TD and the
% ASD subjects
ISFC_RS = cell(1,length(dFC_RS));
ISFC_RS(Group_Label==1) = ISFC_RS1;
ISFC_RS(Group_Label==2) = ISFC_RS2;
%%
%%%%%%
% 3. For the pure RS bits (TD and ASD groups)
%%%%%%
tmp_pureRS1 = TC_RS2(Group_Label_pureRS == 1);
tmp_pureRS2 = TC_RS2(Group_Label_pureRS == 2);
ISFC_PRS1 = cell(1,sum(Group_Label_pureRS == 1));
[ISFC_PRS1{:}] = deal(zeros(size(dFC_RS{1},1),146));
ISFC_PRS2 = cell(1,sum(Group_Label_pureRS == 2));
[ISFC_PRS2{:}] = deal(zeros(size(dFC_RS{1},1),146));
% To keep track of what to divide by for each subject
Denom_factor_TD = cell(1,sum(Group_Label_pureRS == 1));
[Denom_factor_TD{:}] = deal(0);
Denom_factor_ASD = cell(1,sum(Group_Label_pureRS == 2));
[Denom_factor_ASD{:}] = deal(0);
tmp_denom_TD = cell(1,sum(Group_Label_pureRS == 1));
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD = cell(1,sum(Group_Label_pureRS == 2));
[tmp_denom_ASD{:}] = deal(0);
% First we need to determine the indices of SES1 and of SES2 recordings
% for both the TD and the ASD groups
idx_TD = 1:sum(Group_Label_pureRS==1);
idx_ASD = 1:sum(Group_Label_pureRS==2);
% For each fold...
for boot = 1:n_boot
boot
% Next, we need to select n_refsubj/2 elements from SES1 and from SES2
% elements, so first we pick the indices in a random order
idx_TD_b = randperm(length(idx_TD));
% This has the indices of all the subjects that are picked as reference
idx_ref_TD = sort(idx_TD(idx_TD_b(1:n_refsubj/2)),'ascend');
idx_notref_TD = 1:sum(Group_Label_pureRS==1);
idx_notref_TD(idx_ref_TD) = [];
% Next, we select the subjects that make the reference group for the
% particular bootstrapping fold
tmp_ref_TD = tmp_pureRS1(idx_ref_TD);
tmp_notref_TD = tmp_pureRS1(idx_notref_TD);
% We can now do the computations themselves
% With ISFC, we compute the outputs for all TD and ASD subjects...
ISFC_RS1_notref = ISFC_ToReferenceGroup(tmp_notref_TD,tmp_ref_TD,W,Step);
ISFC_RS2_notref = ISFC_ToReferenceGroup(tmp_pureRS2,tmp_ref_TD,W,Step);
% Updates the general cell
tmp_cell = cell(1,sum(Group_Label_pureRS==1));
[tmp_cell{:}] = deal(zeros(size(ISFC_RS1_notref{1})));
tmp_cell(idx_notref_TD) = ISFC_RS1_notref;
ISFC_PRS1 = cellfun(@(x,y) x+y,tmp_cell,ISFC_PRS1,'UniformOutput',0);
clear tmp_cell
tmp_cell = cell(1,sum(Group_Label_pureRS==2));
[tmp_cell{:}] = deal(zeros(size(ISFC_RS2_notref{1})));
tmp_cell(idx_ASD) = ISFC_RS2_notref;
ISFC_PRS2 = cellfun(@(x,y) x+y,tmp_cell,ISFC_PRS2,'UniformOutput',0);
clear tmp_cell
% Keeps track of how many times each subject has been included as a non
% reference (i.e. how many times a value has been computed for it)
tmp_denom_TD(idx_notref_TD) = {1};
Denom_factor_TD = cellfun(@(x,y) x+y,tmp_denom_TD,Denom_factor_TD,'UniformOutput',0);
[tmp_denom_TD{:}] = deal(0);
tmp_denom_ASD(idx_ASD) = {1};
Denom_factor_ASD = cellfun(@(x,y) x+y,tmp_denom_ASD,Denom_factor_ASD,'UniformOutput',0);
[tmp_denom_ASD{:}] = deal(0);
end
% Averages the results; we divide by the number of times a given subject has
% been included
ISFC_PRS1 = cellfun(@(x,y) x/y,ISFC_PRS1,Denom_factor_TD,'UniformOutput',0);
ISFC_PRS2 = cellfun(@(x,y) x/y,ISFC_PRS2,Denom_factor_ASD,'UniformOutput',0);
% Finally, we get back to a representation including both the TD and the
% ASD subjects
ISFC_PRS = cell(1,(sum(Group_Label_pureRS==1)+sum(Group_Label_pureRS==2)));
ISFC_PRS(Group_Label_pureRS==1) = ISFC_PRS1;
ISFC_PRS(Group_Label_pureRS==2) = ISFC_PRS2;
c4science · Help