Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/Matlab_Scripts/Script_Vertes_Analyses.m
257e25cb59ccmaster

Last Change
Raw File

Script_Vertes_Analyses.m
View Options

%% 6. Significance thresholding of connectivity pairs
disp('Thresholding of connectivity pairs...');
SignMat = FindExcursions(ISFC_RS,ISFC_MW,99.98,'Subject-Marginalized-Added-RS',ISFC_PRS);
SignMat_TD = FindExcursions(ISFC_RS(Group_Label==1),ISFC_MW(Group_Label==1),99.98,'Subject-Marginalized-Added-RS',ISFC_PRS(Group_Label_pureRS==1));
SignMat_AST = FindExcursions(ISFC_RS(Group_Label==2),ISFC_MW(Group_Label==2),99.98,'Subject-Marginalized-Added-RS',ISFC_PRS(Group_Label_pureRS==2));
% I want to determine how many excursions there are per time point in my
% subjects; n_excursions_time has size n_TP x n_subjects
n_excursions_time = squeeze(sum(abs(SignMat),1));
% Number of excursions past which we must be in a subject
T_exc_time = prctile(n_excursions_time,100-perc_ret_frames,1);
if is_plot > 1
adj_plot = zeros(n_ROI,n_ROI);
for i = 1:length(Con_vec2)
adj_plot(C1(i),C2(i)) = 1;
end
adj_plot = adj_plot + adj_plot';
PlotBrainGraph_Kmeans(adj_plot,sum(adj_plot,2),CodeBook,0.5,max(sum(adj_plot,2)),max(sum(adj_plot,2)),2,2,'hot','gray',1,0.05,[]);
end
%% 7. Analysis of the number of excursions
disp('Analysis of the number of excursions...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A. Number of excursions across groups
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Number of times that each connection is showing significance
% We want to include the absolute value, because we do not care about the
% direction of the change at this stage: as soon as a connection is away
% from baseline, we consider this a 'hit'
n_excursions = squeeze(sum(abs(SignMat2),2));
n_excursions_CTR = n_excursions(:,Group_Label==1);
n_excursions_ASD = n_excursions(:,Group_Label==2);
% Number of excursions in both groups throughout the brain
for i = 1:(n_CTR+n_ASD)
n_exc_total(i) = sum(n_excursions(:,i));
end
% Two-tailed t-testing of the global excursion number
[h_exc,p_exc] = ttest2(n_exc_total(Group_Label==1),n_exc_total(Group_Label==2));
% Non-parametric assessment of significance
for n_shuffle = 1:1000
% Random index vector
tmp = randperm(n_CTR+n_ASD);
% Group data
G1 = n_exc_total(tmp(1:n_CTR));
G2 = n_exc_total(tmp(n_CTR+1:end));
% p-value under the null
[~,p_exc_null(n_shuffle)] = ttest2(G1,G2);
end
% Alpha = 1% assessment
p_exc_threshold = prctile(p_exc_null,5);
if p_exc < p_exc_threshold
disp(['Significant n_excursions difference at p-value ',num2str(p_exc),'...']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B. Distribution of the number of excursions across groups
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting of the distribution of values
if is_plot > 0
figure;
set(gcf,'color','w');
set(gca,'Box','off');
histogram(mean(n_excursions_CTR,2),0:0.4:35,'FaceColor',[0 0.4 0.6],'FaceAlpha',0.5,'EdgeColor','none');
hold on;
histogram(mean(n_excursions_ASD,2),0:0.4:35,'FaceColor',[1 0.2 0.2],'FaceAlpha',0.5,'EdgeColor','none');
xlabel('Number of significant excursions');
ylabel('Counts');
set(gca,'Box','off');
end
% Kolmogorov Smirnov test to compare the two groups (average histogram of
% the populations): thos operates at the level of the cumulative
% probability distributions
[h_KS,p_KS] = kstest2(mean(n_excursions_CTR,2),mean(n_excursions_ASD,2));
if h_KS
disp(['Significant diference between CTR and ASD distributions of number of excursions, at p-value ',...
num2str(p_KS),'...']);
end
% Non-parametric assessment of the 'threshold p-value'
for n_shuffle = 1:1000
% random index vector
tmp = randperm(n_CTR+n_ASD);
% Values from both groups
G1 = n_excursions(:,tmp(1:n_CTR));
G2 = n_excursions(:,tmp(n_CTR+1:end));
% p-value under the null
[~,p_KS_null(n_shuffle)] = kstest2(mean(G1,2),mean(G2,2));
end
% One-sided test, alpha = 5%
p_KS_threshold = prctile(p_KS_null,5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C. Statistical group comparison at the connection level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjacency matrix describing the number of excursions for different
% regions (better than the vector form because we may want to visualize the
% results)
Adj_ncon = zeros(n_ROI,n_ROI,(n_CTR+n_ASD));
% We will also put our t-values in a similar format
Adj_ncon_tstat = zeros(n_ROI,n_ROI);
Adj_ncon_ANOVA = zeros(n_ROI,n_ROI,3);
% Filling the adjacency matrix for all subjects and connections
for j = 1:(n_CTR+n_ASD)
for i = 1:n_conn2
Adj_ncon(C1(i),C2(i),j) = n_excursions(i,j);
end
% I want a full matrix, so I must transpose what I have from above and
% add it
Adj_ncon(:,:,j) = Adj_ncon(:,:,j) + Adj_ncon(:,:,j)';
end
%%%%%%%%
% This is a parametric assessment with a t-test
%%%%%%%
% Now that I filled the adjacency, I can pull out my group comparison and
% fill my t-value matrix
for i = 1:n_conn2
[~,p_ncon(i),~,stats] = ttest2(Adj_ncon(C1(i),C2(i),Group_Label==1),Adj_ncon(C1(i),C2(i),Group_Label==2));
% Filling both the matrix and a vector
Adj_ncon_tstat(C1(i),C2(i)) = stats.tstat;
tvec_ncon(i) = stats.tstat;
clear stats
end
Adj_ncon_tstat = Adj_ncon_tstat + Adj_ncon_tstat';
% I can compute the significance threshold with Bonferroni correction
alpha_ncon = 0.05/n_conn2;
% I can determine what values it matches in terms of t-statistic
t_inf_ncon = tinv(alpha_ncon/2,(n_CTR+n_ASD-2));
t_sup_ncon = tinv(1-alpha_ncon/2,(n_CTR+n_ASD-2));
idx_sign_ncon_BON = find(tvec_ncon < t_inf_ncon | tvec_ncon > t_sup_ncon);
% FDR testing
[h_FDR_ncon] = fdr_bh(p_ncon,0.05,'pdep');
idx_sign_ncon_FDR = find(h_FDR_ncon == 1);
disp('Significant group differences for number of excursions (connections)');
disp('-------------------------------------------------');
disp('Bonferroni');
Con_vec2(idx_sign_ncon_BON)
disp('-------------------------------------------------');
disp('FDR');
Con_vec2(idx_sign_ncon_FDR)
disp('-------------------------------------------------');
%%%%%%%%
% This is an assessment with an ANOVA design
%%%%%%%
for i = 1:n_conn2
% Current set of values to consider
tmp_anova = squeeze(Adj_ncon(C1(i),C2(i),:));
[p_anova_ncon(:,i),table] = anovan(tmp_anova,{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
% Filling both the matrix and a vector
Adj_ncon_ANOVA(C1(i),C2(i),1) = table{2,6};
Adj_ncon_ANOVA(C1(i),C2(i),2) = table{3,6};
Adj_ncon_ANOVA(C1(i),C2(i),3) = table{4,6};
fvec_ncon(:,i) = [table{2,6};table{3,6};table{4,6}];
clear table
end
for an = 1:3
Adj_ncon_ANOVA(:,:,an) = Adj_ncon_ANOVA(:,:,an) + Adj_ncon_ANOVA(:,:,an)';
end
% I can determine what values it matches in terms of F-statistic
f_sup_ncon = finv(1-alpha_ncon,3,27);
idx_sign_ncon_ANOVA_BON{1} = find(fvec_ncon(1,:) > f_sup_ncon);
idx_sign_ncon_ANOVA_BON{2} = find(fvec_ncon(2,:) > f_sup_ncon);
idx_sign_ncon_ANOVA_BON{3} = find(fvec_ncon(3,:) > f_sup_ncon);
if is_plot > 0
figure;
set(gcf,'color','w');
imagesc(squeeze(mean(Adj_ncon(:,:,Group_Label==1),3)));
colormap(jet);
caxis([0 40]);
axis square
figure;
set(gcf,'color','w');
imagesc(squeeze(mean(Adj_ncon(:,:,Group_Label==2),3)));
colormap(jet);
caxis([0 40]);
axis square
figure;
set(gcf,'color','w');
hist(tvec_ncon,100);
hold on;
set(gca,'Box','off');
plot([t_inf_ncon,t_inf_ncon],get(gca,'YLim'),'r');
plot([t_sup_ncon,t_sup_ncon],get(gca,'YLim'),'r');
% ANOVA plots
figure;
set(gcf,'color','w');
hist(fvec_ncon(1,:),100);
hold on;
set(gca,'Box','off');
plot([f_sup_ncon,f_sup_ncon],get(gca,'YLim'),'r');
figure;
set(gcf,'color','w');
hist(fvec_ncon(2,:),100);
hold on;
set(gca,'Box','off');
plot([f_sup_ncon,f_sup_ncon],get(gca,'YLim'),'r');
figure;
set(gcf,'color','w');
hist(fvec_ncon(3,:),100);
hold on;
set(gca,'Box','off');
plot([f_sup_ncon,f_sup_ncon],get(gca,'YLim'),'r');
% I want to check the temporal evolution of those connections that
% differ across sessions
for i = 1:length(idx_sign_ncon_ANOVA_BON{1})
figure;
set(gcf,'color','w');
tmp = mean(squeeze(SignMat2(idx_sign_ncon_ANOVA_BON{1}(i),:,Group_Label==1)),2);
tmp2 = mean(squeeze(SignMat2(idx_sign_ncon_ANOVA_BON{1}(i),:,Group_Label==2)),2);
plot(19:2:351,tmp,'color',[0.8 0.6 0.2]);
set(gca,'Box','off');
hold on;
plot(19:2:351,tmp2,'color',[0.4 0.4 0.4]);
xlabel('Movie time [s]');
ylabel('Average excursions strength');
end
if is_plot > 1
%%%%%% For BON displays
% tmp_plot is the vector filled with 1 for CTR > ASD and -1 with
% ASD > CTR
tmp_plot = zeros(n_conn2,1);
tmp_plot(idx_sign_ncon_BON) = sign(tvec_ncon(idx_sign_ncon_BON))*1;
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_plot;
% With a value of 1.0, I have blue (CTR > ASD) and with a value of
% 0.97, I have red (ASD > CTR)
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
tmp_plot2(tmp_plot2 == 1) = 1.0;
tmp_plot2(tmp_plot2 == -1) = 0.97;
PlotBrainGraph_Kmeans(tmp_plot2,sum(tmp_plot2,2),CodeBook,0,max(sum(tmp_plot2,2)),max(sum(tmp_plot2,2)),2,2,'hot','lines',1,1);
%%%%%% For FDR displays
% tmp_plot is the vector filled with 1 for CTR > ASD and -1 with
% ASD > CTR
tmp_plot = zeros(n_conn2,1);
tmp_plot(idx_sign_ncon_FDR) = sign(tvec_ncon(idx_sign_ncon_FDR))*1;
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_plot;
% With a value of 1.0, I have blue (CTR > ASD) and with a value of
% 0.97, I have red (ASD > CTR)
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
tmp_plot2(tmp_plot2 == 1) = 1.0;
tmp_plot2(tmp_plot2 == -1) = 0.97;
% I use a slightly lower edge thickness here because there are MANY
% edges left !
PlotBrainGraph_Kmeans(tmp_plot2,sum(tmp_plot2,2),CodeBook,0,max(sum(tmp_plot2,2)),max(sum(tmp_plot2,2)),2,2,'hot','lines',1,0.5,[]);
%%%%%%%%% For ANOVA displays
% tmp_plot is the vector filled with 1 for CTR > ASD and -1 with
% ASD > CTR
tmp_plot = zeros(n_conn2,1);
tmp_plot(idx_sign_ncon_ANOVA_BON{1}) = sign(tvec_ncon(idx_sign_ncon_ANOVA_BON{1}))*1;
% With a value of 1.0, I have blue (CTR > ASD) and with a value of
% 0.97, I have red (ASD > CTR)
tmp_plot2 = jVecToUpperTriMat(tmp_plot,n_ROI)+jVecToUpperTriMat(tmp_plot,n_ROI)';
tmp_plot2(tmp_plot2 == 1) = 1.0;
tmp_plot2(tmp_plot2 == -1) = 0.97;
PlotBrainGraph_Kmeans(tmp_plot2,sum(tmp_plot2,2),CodeBook,0,max(sum(tmp_plot2,2)),max(sum(tmp_plot2,2)),2,2,'hot','lines',1,0.5,[]);
% tmp_plot is the vector filled with 1 for CTR > ASD and -1 with
% ASD > CTR
tmp_plot = zeros(n_conn2,1);
tmp_plot(idx_sign_ncon_ANOVA_BON{2}) = 1;
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_plot;
% With a value of 1.0, I have blue (CTR > ASD) and with a value of
% 0.97, I have red (ASD > CTR)
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
tmp_plot2(tmp_plot2 == 1) = 1.0;
tmp_plot2(tmp_plot2 == -1) = 0.97;
PlotBrainGraph_Kmeans(tmp_plot2,sum(tmp_plot2,2),CodeBook,0,max(sum(tmp_plot2,2))*2,max(sum(tmp_plot2,2)),2,2,'hot','jet',1,0.5,[]);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D. Statistical group comparison at the ROI level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I can look at the degree now
Degree_ncon = squeeze(sum(Adj_ncon,1));
% As before, I compute my t and p values
for i = 1:n_ROI
[~,p_ndeg(i),~,stats] = ttest2(Degree_ncon(i,Group_Label==1)',Degree_ncon(i,Group_Label==2)');
tvec_ndeg(i) = stats.tstat;
clear stats
end
% I can compute the significance threshold with Bonferroni correction
alpha_ndeg = 0.05/n_ROI;
% I can determine what values it matches in terms of t-statistic
t_inf_ndeg = tinv(alpha_ndeg/2,(n_CTR+n_ASD-2));
t_sup_ndeg = tinv(1-alpha_ndeg/2,(n_CTR+n_ASD-2));
idx_sign_ndeg_BON = find(tvec_ndeg < t_inf_ndeg | tvec_ndeg > t_sup_ndeg);
% FDR testing
[h_FDR_ndeg] = fdr_bh(p_ndeg,0.05,'pdep');
idx_sign_ndeg_FDR = find(h_FDR_ndeg == 1);
disp('Significant group differences for number of excursions (ROIs)');
disp('-------------------------------------------------');
disp('Bonferroni');
ROI_vec(idx_sign_ndeg_BON)
disp('-------------------------------------------------');
disp('FDR');
ROI_vec(idx_sign_ndeg_FDR)
disp('-------------------------------------------------');
% The same with an ANOVA again
% As before, I compute my t and p values
for i = 1:n_ROI
[p_anova_ndeg(:,i),table] = anovan(Degree_ncon(i,:),{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
fvec_ndeg(:,i) = [table{2,6};table{3,6};table{4,6}];
clear table
end
f_sup_ndeg = finv(1-alpha_ndeg,3,27);
idx_sign_ndeg_ANOVA_BON{1} = find(fvec_ndeg(1,:) > f_sup_ndeg);
idx_sign_ndeg_ANOVA_BON{2} = find(fvec_ndeg(2,:) > f_sup_ndeg);
idx_sign_ndeg_ANOVA_BON{3} = find(fvec_ndeg(3,:) > f_sup_ndeg);
disp('Significant group differences for number of excursions (ROIs)');
disp('-------------------------------------------------');
disp('Group difference');
ROI_vec(idx_sign_ndeg_ANOVA_BON{1})
disp('-------------------------------------------------');
disp('Session difference');
ROI_vec(idx_sign_ndeg_ANOVA_BON{2})
disp('-------------------------------------------------');
disp('Group/Session interaction');
ROI_vec(idx_sign_ndeg_ANOVA_BON{3})
if is_plot > 0
figure;
set(gcf,'color','w');
b1 = bar(mean(Degree_ncon(:,Group_Label==1),2),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(mean(Degree_ncon(:,Group_Label==2),2),'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
set(gca,'Box','off');
xlabel('ROIs');
ylabel('Degree');
figure;
set(gcf,'color','w');
hist(tvec_ndeg,20);
hold on;
set(gca,'Box','off');
plot([t_inf_ndeg,t_inf_ndeg],get(gca,'YLim'),'r');
plot([t_sup_ndeg,t_sup_ndeg],get(gca,'YLim'),'r');
if is_plot > 1
% Plots the node information (just degree of all nodes in each group)
tmp_CTR = mean(Adj_ncon(:,:,Group_Label==1),3);
tmp_ASD = mean(Adj_ncon(:,:,Group_Label==2),3);
PlotBrainGraph_Kmeans(tmp_CTR,sum(tmp_CTR,2),CodeBook,1.5,max(sum(tmp_CTR,2)),max(sum(tmp_CTR,2)),2,2,'hot','gray',1,0.5);
PlotBrainGraph_Kmeans(tmp_ASD,sum(tmp_ASD,2),CodeBook,1.5,max(sum(tmp_ASD,2)),max(sum(tmp_ASD,2)),2,2,'hot','gray',1,0.5);
% BONFERRONI NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ndeg_BON) = sign(tvec_ndeg(idx_sign_ndeg_BON));
tmp_plot(tmp_plot < 0) = 0.97;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% FDR NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ndeg_FDR) = sign(tvec_ndeg(idx_sign_ndeg_FDR));
tmp_plot(tmp_plot < 0) = 0.97;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% ANOVA NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ndeg_ANOVA_BON{1}) = sign(tvec_ndeg(idx_sign_ndeg_ANOVA_BON{1}));
tmp_plot(tmp_plot < 0) = 0.97;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E. Top connections in each group for excursion number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Indices of the connections that show most excursions in both groups, and
% of the regions as well
[~,idx_top_ncon_CTR] = sort(mean(n_excursions_CTR,2),'descend');
[~,idx_top_ndeg_CTR] = sort(mean(Degree_ncon(:,Group_Label==1),2),'descend');
[~,idx_top_ncon_ASD] = sort(mean(n_excursions_ASD,2),'descend');
[~,idx_top_ndeg_ASD] = sort(mean(Degree_ncon(:,Group_Label==2),2),'descend');
if is_plot > 1
% Plotting top connections for the CTR group
tmp_plot = zeros(n_ROI,n_ROI);
for i = 1:n_top_connections
tmp_plot(C1(idx_top_ncon_CTR(i)),C2(idx_top_ncon_CTR(i))) = 1;
end
tmp_plot = tmp_plot + tmp_plot';
PlotBrainGraph_Kmeans(tmp_plot,sum(tmp_plot,2),CodeBook,0.5,max(sum(tmp_plot,2)),max(sum(tmp_plot,2)),2,2,'hot','gray',1,0.5,[]);
% The same for the ASD group
% Plotting top connections for the CTR group
tmp_plot = zeros(n_ROI,n_ROI);
for i = 1:n_top_connections
tmp_plot(C1(idx_top_ncon_ASD(i)),C2(idx_top_ncon_ASD(i))) = 1;
end
tmp_plot = tmp_plot + tmp_plot';
PlotBrainGraph_Kmeans(tmp_plot,sum(tmp_plot,2),CodeBook,0.5,max(sum(tmp_plot,2)),max(sum(tmp_plot,2)),2,2,'hot','gray',1,0.5,[]);
end
%% 9. Graph analysis
disp('Graph analysis...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A. Computation of adjacency and degree information for both groups
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adj will contain the adjacency matrices for each time point
Adj = zeros(n_ROI,n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
% Adj_short and Adj_long contain the short-range and long-range connection
% parts
Adj_short = zeros(n_ROI,n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
Adj_long = zeros(n_ROI,n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
% Associated degree time courses for each ROI
Degree_long = nan(n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
Degree_short = nan(n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
Degree_all = nan(n_ROI,size(ISFC_MW{1},2),(n_CTR+n_ASD));
% Distance between the connections that were kept by the feature extraction
% approach
D2_euclidean = D_euclidean(idx_consistent);
% For each time point...
for t = 1:size(ISFC_MW{1},2)
% For each subject...
for s = 1:(n_CTR + n_ASD)
% I take the connectivity information
tmp = squeeze(ISFC_MW{s}(:,t));
% Thresholding of the graph
% Puts back the connections in the full vector
tmp2 = zeros(n_conn,1);
tmp2(idx_consistent) = tmp;
Adj(:,:,t,s) = jVecToUpperTriMat(tmp2,n_ROI);
Adj(:,:,t,s) = Adj(:,:,t,s) + Adj(:,:,t,s)';
Adj_short(:,:,t,s) = Adj_short(:,:,t,s) + Adj_short(:,:,t,s)';
Adj_long(:,:,t,s) = Adj_long(:,:,t,s) + Adj_long(:,:,t,s)';
% Low distance degree
Degree_short(:,t,s) = sum(squeeze(Adj_short(:,:,t,s)));
% Long distance degree
Degree_long(:,t,s) = sum(squeeze(Adj_long(:,:,t,s)));
Degree_all(:,t,s) = sum(squeeze(Adj(:,:,t,s)));
end
end
% My quantities of interest
Shortdeg = mean(Degree_short(:,:,Group_Label==1),3)-mean(Degree_short(:,:,Group_Label==2),3);
Longdeg = mean(Degree_long(:,:,Group_Label==1),3)-mean(Degree_long(:,:,Group_Label==2),3);
Alldeg = mean(Degree_all(:,:,Group_Label==1),3)-mean(Degree_all(:,:,Group_Label==2),3);
if is_plot > 0
figure;
set(gcf,'color','w');
subplot(3,1,1);
imagesc(19:2:351,1:299,mean(Degree_all(:,:,Group_Label==1),3));
xlabel('Movie time [s]');
ylabel('ROI');
colormap(jet);
caxis([0 20]);
subplot(3,1,2);
set(gcf,'color','w');
imagesc(19:2:351,1:299,mean(Degree_all(:,:,Group_Label==2),3));
xlabel('Movie time [s]');
ylabel('ROI');
colormap(jet);
caxis([0 20]);
subplot(3,1,3);
set(gcf,'color','w');
imagesc(19:2:351,1:299,Alldeg);
xlabel('Movie time [s]');
ylabel('ROI');
colormap(jet);
caxis([-15 15]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C. Significance assessment for short-range and long-range nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Degree_scon and Degree_lcon have size n_ROI x n_subjects
Degree_scon = squeeze(sum(squeeze(sum(Adj_short,1)),2));
Degree_lcon = squeeze(sum(squeeze(sum(Adj_long,1)),2));
% Computation of p-values and t-values
for i = 1:n_ROI
% Short-range
[~,p_sdeg(i),~,stats1] = ttest2(Degree_scon(i,Group_Label==1)',Degree_scon(i,Group_Label==2)');
tvec_sdeg(i) = stats1.tstat;
clear stats1
% Long-range
[~,p_ldeg(i),~,stats2] = ttest2(Degree_lcon(i,Group_Label==1)',Degree_lcon(i,Group_Label==2)');
tvec_ldeg(i) = stats2.tstat;
clear stats2
end
% Significance threshold with Bonferroni correction
alpha_sdeg = 0.05/n_ROI/2;
% I can determine what values it matches in terms of t-statistic
t_inf_sdeg = tinv(alpha_sdeg/2,(n_CTR+n_ASD-2));
t_sup_sdeg = tinv(1-alpha_sdeg/2,(n_CTR+n_ASD-2));
% Indices of significant ROIs in short and long-range cases
idx_sign_sdeg_BON = find(tvec_sdeg < t_inf_sdeg | tvec_sdeg > t_sup_sdeg);
idx_sign_ldeg_BON = find(tvec_ldeg < t_inf_sdeg | tvec_ldeg > t_sup_sdeg);
% FDR testing for both cases
[h_FDR_sdeg] = fdr_bh(p_sdeg,0.05/2,'pdep');
[h_FDR_ldeg] = fdr_bh(p_ldeg,0.05,'pdep');
% Related indices of significance
idx_sign_sdeg_FDR = find(h_FDR_sdeg == 1);
idx_sign_ldeg_FDR = find(h_FDR_ldeg == 1);
disp('Significant group differences for short-range degree');
disp('-------------------------------------------------');
disp('Bonferroni');
ROI_vec(idx_sign_sdeg_BON)
disp('-------------------------------------------------');
disp('FDR');
ROI_vec(idx_sign_sdeg_FDR)
disp('-------------------------------------------------');
disp('Significant group differences for long-range degree');
disp('-------------------------------------------------');
disp('Bonferroni');
ROI_vec(idx_sign_ldeg_BON)
disp('-------------------------------------------------');
disp('FDR');
ROI_vec(idx_sign_ldeg_FDR)
disp('-------------------------------------------------');
% SAME BUT WITH AN ANOVA DESIGN
% The same with an ANOVA again
% As before, I compute my t and p values
for i = 1:n_ROI
[p_anova_sdeg(:,i),table] = anovan(Degree_scon(i,:),{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
fvec_sdeg(:,i) = [table{2,6};table{3,6};table{4,6}];
clear table
end
f_sup_sdeg = finv(1-alpha_sdeg,3,27);
idx_sign_sdeg_ANOVA_BON{1} = find(fvec_sdeg(1,:) > f_sup_sdeg);
idx_sign_sdeg_ANOVA_BON{2} = find(fvec_sdeg(2,:) > f_sup_sdeg);
idx_sign_sdeg_ANOVA_BON{3} = find(fvec_sdeg(3,:) > f_sup_sdeg);
% As before, I compute my t and p values
for i = 1:n_ROI
[p_anova_ldeg(:,i),table] = anovan(Degree_lcon(i,:),{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
fvec_ldeg(:,i) = [table{2,6};table{3,6};table{4,6}];
clear table
end
f_sup_ldeg = finv(1-alpha_sdeg,3,27);
idx_sign_ldeg_ANOVA_BON{1} = find(fvec_ldeg(1,:) > f_sup_ldeg);
idx_sign_ldeg_ANOVA_BON{2} = find(fvec_ldeg(2,:) > f_sup_ldeg);
idx_sign_ldeg_ANOVA_BON{3} = find(fvec_ldeg(3,:) > f_sup_ldeg);
disp('Significant group differences for number of excursions (ROIs)');
disp('-------------------------------------------------');
disp('Group difference');
disp('Short-range');
ROI_vec(idx_sign_sdeg_ANOVA_BON{1})
disp('Group difference');
disp('Long-range');
ROI_vec(idx_sign_ldeg_ANOVA_BON{1})
disp('-------------------------------------------------');
disp('Session difference');
disp('Short-range');
ROI_vec(idx_sign_sdeg_ANOVA_BON{2})
disp('Session difference');
disp('Long-range');
ROI_vec(idx_sign_ldeg_ANOVA_BON{2})
disp('-------------------------------------------------');
disp('Group/session interaction');
disp('Short-range');
ROI_vec(idx_sign_sdeg_ANOVA_BON{3})
disp('Session difference');
disp('Long-range');
ROI_vec(idx_sign_ldeg_ANOVA_BON{3})
if is_plot > 0
% Short-range plots
figure;
set(gcf,'color','w');
subplot(3,1,1);
b1 = bar(mean(Degree_ncon(:,Group_Label==1),2),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(mean(Degree_ncon(:,Group_Label==2),2),'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
ylim([0 1000]);
xlabel('ROI');
ylabel('Degree');
subplot(3,1,2);
b1 = bar(mean(Degree_scon(:,Group_Label==1),2),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(mean(Degree_scon(:,Group_Label==2),2),'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
ylim([0 700]);
xlabel('ROI');
ylabel('Short-range degree');
subplot(3,1,3);
set(gcf,'color','w');
b1 = bar(mean(Degree_lcon(:,Group_Label==1),2),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(mean(Degree_lcon(:,Group_Label==2),2),'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
ylim([0 400]);
xlabel('ROI');
ylabel('Long-range degree');
% I want a look at the nodes showing significant differences
tmp_short = squeeze(sum(squeeze(sum(Adj_short,3)),1));
for i = 1:length(idx_sign_sdeg_ANOVA_BON{2})
figure;
set(gcf,'color','w');
bar([mean(tmp_short(idx_sign_sdeg_ANOVA_BON{2}(i),Session_Label == 1 & Group_Label == 1),2),...
mean(tmp_short(idx_sign_sdeg_ANOVA_BON{2}(i),Session_Label == 1 & Group_Label == 2),2)...
mean(tmp_short(idx_sign_sdeg_ANOVA_BON{2}(i),Session_Label == 2 & Group_Label == 1),2),...
mean(tmp_short(idx_sign_sdeg_ANOVA_BON{2}(i),Session_Label == 2 & Group_Label == 2),2)]);
end
% I want a look at the nodes showing significant differences
tmp_long = squeeze(sum(squeeze(sum(Adj_long,3)),1));
for i = 1:length(idx_sign_ldeg_ANOVA_BON{2})
figure;
set(gcf,'color','w');
bar([mean(tmp_long(idx_sign_ldeg_ANOVA_BON{2}(i),Session_Label == 1 & Group_Label == 1),2),...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{2}(i),Session_Label == 1 & Group_Label == 2),2)...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{2}(i),Session_Label == 2 & Group_Label == 1),2),...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{2}(i),Session_Label == 2 & Group_Label == 2),2)]);
end
for i = 1:length(idx_sign_ldeg_ANOVA_BON{3})
figure;
set(gcf,'color','w');
bar([mean(tmp_long(idx_sign_ldeg_ANOVA_BON{3}(i),Session_Label == 1 & Group_Label == 1),2),...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{3}(i),Session_Label == 1 & Group_Label == 2),2)...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{3}(i),Session_Label == 2 & Group_Label == 1),2),...
mean(tmp_long(idx_sign_ldeg_ANOVA_BON{3}(i),Session_Label == 2 & Group_Label == 2),2)]);
end
if is_plot > 1
%%%% SHORT_RANGE NODES
% This is the 'overall' short-range adjacency
tmp_short = squeeze(sum(Adj_short,3));
% Plots the node information (just degree of all nodes in each
% group) for short-range nodes
tmp_CTR = mean(tmp_short(:,:,Group_Label==1),3);
tmp_ASD = mean(tmp_short(:,:,Group_Label==2),3);
PlotBrainGraph_Kmeans(tmp_CTR,sum(tmp_CTR,2),CodeBook,1.5,max(sum(tmp_CTR,2)),max(sum(tmp_CTR,2)),2,2,'hot','gray',1,0.5,[]);
PlotBrainGraph_Kmeans(tmp_ASD,sum(tmp_ASD,2),CodeBook,1.5,max(sum(tmp_ASD,2)),max(sum(tmp_ASD,2)),2,2,'hot','gray',1,0.5,[]);
% BONFERRONI NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_sdeg_BON) = sign(tvec_sdeg(idx_sign_sdeg_BON));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% FDR NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_sdeg_FDR) = sign(tvec_sdeg(idx_sign_sdeg_FDR));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% ANOVA NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_sdeg_ANOVA_BON{1}) = sign(tvec_sdeg(idx_sign_sdeg_ANOVA_BON{1}));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_sdeg_ANOVA_BON{2}) = 1;
tmp_plot(tmp_plot < 0) = 0.6;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[0 1]);
%%% LONG-RANGE NODES
% This is the 'overall' short-range adjacency
tmp_long = squeeze(sum(Adj_long,3));
% Plots the node information (just degree of all nodes in each
% group) for short-range nodes
tmp_CTR = mean(tmp_long(:,:,Group_Label==1),3);
tmp_ASD = mean(tmp_long(:,:,Group_Label==2),3);
PlotBrainGraph_Kmeans(tmp_CTR,sum(tmp_CTR,2),CodeBook,1.5,max(sum(tmp_CTR,2)),max(sum(tmp_CTR,2)),2,2,'hot','gray',1,0.5,[]);
PlotBrainGraph_Kmeans(tmp_ASD,sum(tmp_ASD,2),CodeBook,1.5,max(sum(tmp_ASD,2)),max(sum(tmp_ASD,2)),2,2,'hot','gray',1,0.5,[]);
% BONFERRONI NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ldeg_BON) = sign(tvec_ldeg(idx_sign_ldeg_BON));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% FDR NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ldeg_FDR) = sign(tvec_ldeg(idx_sign_ldeg_FDR));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5);
% ANOVA NODES
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ldeg_ANOVA_BON{1}) = sign(tvec_ldeg(idx_sign_ldeg_ANOVA_BON{1}));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ldeg_ANOVA_BON{2}) = sign(tvec_ldeg(idx_sign_ldeg_ANOVA_BON{2}));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
% For the plots to come, I want the 'normal' nodes at 0.5, the CTR
% > ASD ones at 1, the ASD > CTR ones at -1
tmp_plot = 0.29*ones(n_ROI,1);
tmp_plot(idx_sign_ldeg_ANOVA_BON{3}) = sign(tvec_ldeg(idx_sign_ldeg_ANOVA_BON{3}));
tmp_plot(tmp_plot < 0) = 0.98;
PlotBrainGraph_Kmeans(diag(tmp_plot),tmp_plot,CodeBook,1.5,max(tmp_plot),max(tmp_plot),2,2,'lines','gray',1,0.5,[]);
end
end
%% 10. Extraction of connectivity states: k-means clustering
disp('Extraction of connectivity states through k-means clustering...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A. K-means clustering on the raw data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creation of my data matrix: I add only the frames for which there is a
% large enough number of excursions
% Will contain all the frames of interest (n_datapoints x n_dimensions)
Data_kmeans = [];
% Associated group labels
GL_kmeans = [];
% Associated subject labels
SL_kmeans = [];
TL_kmeans = [];
% For all subjects...
for s = 1:(n_CTR + n_ASD)
% For all time points...
for t = 1:size(ISFC_MW{1},2)
% Sampling of the frame of interest
tmp = ISFC_MW{s}(:,t);
% Removal of the non-conserved connections
%tmp2 = tmp(idx_consistent);
% If the particular time point survives our threshold, we add it to
% the data matrix
Data_kmeans = [Data_kmeans;tmp'];
GL_kmeans = [GL_kmeans;Group_Label(s)];
SL_kmeans = [SL_kmeans;s];
TL_kmeans = [TL_kmeans;t];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B. Analysis of the results from k-means clustering: choice of K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I want to test many values of K
K_label = [];
Kmeans_boxplot = [];
index_kmeans = 1;
% I fill in a vector of values with my reproducibility measures for all
% clusters
MEAN_JAC = [];
MED_JAC = [];
QUAR25_JAC = [];
QUAR75_JAC = [];
SEM_JAC = [];
for K = 10:5:150
disp(['Currently clustering for K=',num2str(K),'...']);
% IDX contains the indices, and centroids of the actual connection patterns
[IDX{index_kmeans},Centroids{index_kmeans}] = kmeans(Data_kmeans,K,'MaxIter',200,'Replicates',20);
% I now want to compute summarizing metrics
[Entropy{index_kmeans},J_CTR{index_kmeans},J_ASD{index_kmeans},J_all{index_kmeans}] = AssessClusterQuality(IDX{index_kmeans},GL_kmeans,SL_kmeans,TL_kmeans,n_CTR,n_ASD,size(ISFC_MW{1},2),K);
% I want to retain the best Jaccard value for each cluster
tmp_Jaccard = [J_CTR{index_kmeans};J_ASD{index_kmeans};J_all{index_kmeans}];
tmp_Jaccard_box = max(tmp_Jaccard,[],1);
MEAN_JAC = [MEAN_JAC,mean(tmp_Jaccard_box)];
SEM_JAC = [SEM_JAC,std(tmp_Jaccard_box)/sqrt(length(tmp_Jaccard_box))];
MED_JAC = [MED_JAC,median(tmp_Jaccard_box)];
QUAR25_JAC = [QUAR25_JAC,prctile(tmp_Jaccard_box,25)];
QUAR75_JAC = [QUAR75_JAC,prctile(tmp_Jaccard_box,75)];
K_label = [K_label;repmat(K,K,1)];
Kmeans_boxplot = [Kmeans_boxplot;tmp_Jaccard_box'];
index_kmeans = index_kmeans + 1;
end
% Now that I generated my clustering, I want to analyse my outcomes and I
% start with the boxplot
if is_plot > 0
figure;
set(gcf,'color','w');
boxplot(Kmeans_boxplot,K_label);
h = findobj('Tag','Box');
set(h,'color','k','LineWidth',2);
h = findobj('Tag','Median');
set(h,'color',[0.5 0.5 0.5],'LineWidth',2);
h = findobj('Tag','Outliers');
set(h,'Marker','o','MarkerFaceColor',[0.5 0.5 0.5],...
'MarkerEdgeColor',[0.5 0.5 0.5],'MarkerSize',3);
set(gca,'Box','off');
xlabel('Number of clusters');
ylabel('Jaccard index');
figure;
set(gcf,'color','w');
errorbar(10:5:150,MED_JAC,MED_JAC-QUAR25_JAC,-MED_JAC+QUAR75_JAC);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C. K-means clustering: features computed on all clusters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For each cluster, I want to generate a time course for each group, with 1
% if the cluster is in all subjects at a given time point and 0 if it is in
% no subject
idx_kmeans = 1;
for K = 10:5:150
disp(num2str(K));
%%%% What comes below is the temporal information on the clusters
% Current data to consider (indices)
tmp_index = IDX{idx_kmeans};
% I want my data to be converted into a matrix format (n_TP x n_subj)
tmp_index = reshape(tmp_index,167,(n_CTR+n_ASD));
% Reproducibility measurements and related thresholds
[TC_CTR{idx_kmeans},TC_ASD{idx_kmeans},TC_CTR_thresh{idx_kmeans},...
TC_ASD_thresh{idx_kmeans},Thresh_CTR{idx_kmeans},...
Thresh_ASD{idx_kmeans}] = AssessClusterTimeCourses(tmp_index,Group_Label);
% The same, but for sessions
[TC_SES1{idx_kmeans},TC_SES2{idx_kmeans},TC_thresh_SES1{idx_kmeans},...
TC_thresh_SES2{idx_kmeans},Thresh_SES1{idx_kmeans},Thresh_SES2{idx_kmeans}] = AssessClusterTimeCourses(tmp_index,Session_Label);
% Correlation matrix denoting the relationship between clusters
Corr_Mat{idx_kmeans} = corr(TC_CTR_thresh{idx_kmeans}',TC_ASD_thresh{idx_kmeans}');
% The same, but for covariance instead
% for j = 1:K
% for i = 1:K
% tmp = cov(TC_CTR{idx_kmeans}(j,:)',TC_ASD{idx_kmeans}(i,:)');
% tmp_ses = cov(TC_CTR_ses{idx_kmeans}(j,:)',TC_ASD_ses{idx_kmeans}(i,:)');
% Cov_Mat{idx_kmeans}(j,i) = (tmp(1,1) + tmp(2,2))/2;
% Cov_Mat_ses{idx_kmeans}(j,i) = (tmp_ses(1,1) + tmp_ses(2,2))/2;
% end
% end
% The same on thresholded data
%Corr_Mat_thresh{idx_kmeans} = corr(TC_CTR_thresh{idx_kmeans}',TC_ASD_thresh{idx_kmeans}');
%%%% What comes below is the group and spatial information on the clusters
% Of course, I also want to compute counts, so Counts is a matrix of
% size n_clusters x n_subjects
[Counts{idx_kmeans}] = ComputeClusterCounts(tmp_index);
% I want to assess possible count differences for all cases
tmp_counts = Counts{idx_kmeans};
for j = 1:K
[p_anova_kmeans_clus{idx_kmeans}(j,:),table] = anovan(tmp_counts(j,:),{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
[~,~,~,stats] = ttest2(tmp_counts(j,Group_Label==1)',tmp_counts(j,Group_Label==2)');
fvec_kmeans_clus{idx_kmeans}(j,:) = [table{2,6};table{3,6};table{4,6}];
tvec_kmeans_clus{idx_kmeans}(j) = stats.tstat;
clear table
clear stats
end
% Alpha-level for a given number of clusters (2320 is the total number
% of clusters throughout all the parcellations)
alpha_kmeans{idx_kmeans} = 0.05/2320;
% I can determine what values it matches in terms of F-statistic
f_sup_kmeans_clus{idx_kmeans} = finv(1-alpha_kmeans{idx_kmeans},3,27);
idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1} = find(squeeze(fvec_kmeans_clus{idx_kmeans}(:,1) > f_sup_kmeans_clus{idx_kmeans}));
idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,2} = find(squeeze(fvec_kmeans_clus{idx_kmeans}(:,2) > f_sup_kmeans_clus{idx_kmeans}));
idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,3} = find(squeeze(fvec_kmeans_clus{idx_kmeans}(:,3) > f_sup_kmeans_clus{idx_kmeans}));
% How many clusters are CTR or ASD ones
n_CTR_clusters(idx_kmeans) = sum(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) > 0);
n_ASD_clusters(idx_kmeans) = sum(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) < 0);
% Idiosynchracy in patterns, i.e. fraction of different patterns
% present at each time point normalized by subject
for t = 1:167
Idiosynch_CTR(t,idx_kmeans) = length(unique(tmp_index(t,Group_Label==1)))/n_CTR;
Idiosynch_ASD(t,idx_kmeans) = length(unique(tmp_index(t,Group_Label==2)))/n_ASD;
end
% Spatial similarity between all clusters
Sim_Mat_CTRASD{idx_kmeans} = ...
corr(Centroids{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) > 0),:)',...
Centroids{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) < 0),:)');
% Spatial similarity between only the CTR clusters
Sim_Mat_CTR{idx_kmeans} = ...
corr(Centroids{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) > 0),:)');
Sim_Mat_ASD{idx_kmeans} = ...
corr(Centroids{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}(tvec_kmeans_clus{idx_kmeans}(idx_sign_kmeans_clus_ANOVA_BON{idx_kmeans,1}) < 0),:)');
Sim_clusters_CTRASD_mean(idx_kmeans) = mean(jUpperTriMatToVec(Sim_Mat_CTRASD{idx_kmeans}));
Sim_clusters_CTRASD_std(idx_kmeans) = std(jUpperTriMatToVec(Sim_Mat_CTRASD{idx_kmeans}));
Sim_clusters_CTR_mean(idx_kmeans) = mean(jUpperTriMatToVec(Sim_Mat_CTR{idx_kmeans}));
Sim_clusters_CTR_std(idx_kmeans) = std(jUpperTriMatToVec(Sim_Mat_CTR{idx_kmeans}));
Sim_clusters_ASD_mean(idx_kmeans) = mean(jUpperTriMatToVec(Sim_Mat_ASD{idx_kmeans}));
Sim_clusters_ASD_std(idx_kmeans) = std(jUpperTriMatToVec(Sim_Mat_ASD{idx_kmeans}));
idx_kmeans = idx_kmeans + 1;
end
if is_plot > 0
figure;
set(gcf,'color','w');
plot(10:5:150,n_CTR_clusters,'color','b','LineWidth',2);
hold on;
plot(10:5:150,n_ASD_clusters,'color','r','LineWidth',2);
xlabel('Number of clusters for parcellation');
ylabel('Number of significant clusters in the group');
set(gca,'Box','off');
figure;
set(gcf,'color','w');
b1 = bar(10:5:150,n_CTR_clusters,'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(10:5:150,n_ASD_clusters,'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
xlabel('Number of clusters for parcellation');
ylabel('Number of significant clusters in the group');
xlim([9 151]);
set(gca,'Box','off');
figure;
set(gcf,'color','w');
errorbar(10:5:150,mean(Idiosynch_CTR,1),std(Idiosynch_CTR,[],1)/sqrt(167),std(Idiosynch_CTR,[],1)/sqrt(167),'color','b','LineWidth',2);
hold on;
errorbar(10:5:150,mean(Idiosynch_ASD,1),std(Idiosynch_ASD,[],1)/sqrt(167),std(Idiosynch_ASD,[],1)/sqrt(167),'color','r','LineWidth',2);
xlabel('Number of clusters for parcellation');
ylabel('Idiosynchracy');
set(gca,'Box','off');
figure;
set(gcf,'color','w');
errorbar(10:5:150,Sim_clusters_CTRASD_mean,Sim_clusters_CTRASD_std./sqrt((n_CTR_clusters+n_ASD_clusters).*(n_CTR_clusters+n_ASD_clusters-1)/2),Sim_clusters_CTRASD_std./sqrt((n_CTR_clusters+n_ASD_clusters).*(n_CTR_clusters+n_ASD_clusters-1)/2),'color','k','LineWidth',2);
hold on;
errorbar(10:5:150,Sim_clusters_ASD_mean,Sim_clusters_ASD_std./sqrt((n_ASD_clusters).*(n_ASD_clusters-1)/2),Sim_clusters_ASD_std./sqrt((n_ASD_clusters).*(n_ASD_clusters-1)/2),'color','r','LineWidth',2);
errorbar(10:5:150,Sim_clusters_CTR_mean,Sim_clusters_CTR_std./sqrt((n_CTR_clusters).*(n_CTR_clusters-1)/2),Sim_clusters_CTR_std./sqrt((n_CTR_clusters).*(n_CTR_clusters-1)/2),'color','b','LineWidth',2);
xlim([29 101]);
set(gca,'Box','off');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D. K-means clustering: features computed on selected cluster numbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Selected values of K
K_vector = 10:5:150;
idx_sign_K = find(K_vector == 30 | K_vector == 120);
%%% SHARED CONNECTIONS AS A FUNCTION OF COUNT DIFFERING CLUSTERS
% For the chosen clusters, plotting
for i2 = 1:length(idx_sign_K)
clear adj_kmeans
% Samples the clusters to plot
tmp_Centroids_ext = zeros(n_conn,K_vector(idx_sign_K(i2)));
tmp_Centroids_ext(idx_consistent,:) = Centroids{idx_sign_K(i2)}';
% Constructs the adjacency matrix of each cluster
for j = 1:K_vector(idx_sign_K(i2))
adj_kmeans(:,:,j) = jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI) + jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI)';
end
adj_kmeans(abs(adj_kmeans) < 0.5) = 0;
% I keep all the connections that survive thresholding
adj_kmeans(adj_kmeans ~= 0) = 1;
adj_kmeans = logical(adj_kmeans);
% Across all clusters
adj_kmeans_summary = sum(adj_kmeans,3);
% I can look for the connections that are specific to a group
idx_sign_kmeans_CTR = idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(i2),1}(tvec_kmeans_clus{idx_sign_K(i2)}(idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(i2),1})>0);
size(idx_sign_kmeans_CTR)
idx_sign_kmeans_ASD = idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(i2),1}(tvec_kmeans_clus{idx_sign_K(i2)}(idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(i2),1})<0);
size(idx_sign_kmeans_ASD)
% Only contains the CTR or the ASD-specific connections
adj_kmeans_CTR = zeros(n_ROI,n_ROI);
for j = 1:length(idx_sign_kmeans_CTR)
adj_kmeans_CTR = adj_kmeans_CTR | adj_kmeans(:,:,idx_sign_kmeans_CTR(j));
end
adj_kmeans_ASD = zeros(n_ROI,n_ROI);
for j = 1:length(idx_sign_kmeans_ASD)
adj_kmeans_ASD = adj_kmeans_ASD | adj_kmeans(:,:,idx_sign_kmeans_ASD(j));
end
% Contains what is in CTR but not in ASD
adj_kmeans_CTR_spec = adj_kmeans_CTR & ~adj_kmeans_ASD;
% Contains what is in ASD but not CTR
adj_kmeans_ASD_spec = adj_kmeans_ASD & ~adj_kmeans_CTR;
adj_kmeans_CTR = double(adj_kmeans_CTR);
adj_kmeans_ASD = double(adj_kmeans_ASD);
adj_kmeans_CTR_spec = double(adj_kmeans_CTR_spec);
adj_kmeans_ASD_spec = double(adj_kmeans_ASD_spec);
deg_kmeans_CTR_spec = sum(adj_kmeans_CTR_spec,2);
deg_kmeans_ASD_spec = sum(adj_kmeans_ASD_spec,2);
% Top nodes specific to each group
[~,idx_kmeans_topnodes_CTR_spec{i2}] = sort(deg_kmeans_CTR_spec,'descend');
[~,idx_kmeans_topnodes_ASD_spec{i2}] = sort(deg_kmeans_ASD_spec,'descend');
if is_plot > 1
% Plots the general connections across all clusters, and the
% connections specific to CTR-specific or ASD-specific clusters
%PlotBrainGraph_Kmeans(adj_kmeans_summary,sum(adj_kmeans_summary,2),CodeBook,0,max(sum(adj_kmeans_summary,2)),max(sum(adj_kmeans_summary,2)),2,2,'hot','jet',0.5,0.5,[]);
PlotBrainGraph_Kmeans(adj_kmeans_CTR_spec,sum(adj_kmeans_CTR_spec,2),CodeBook,0,max(sum(adj_kmeans_CTR_spec,2)),max(sum(adj_kmeans_CTR_spec,2)),2,2,'hot','gray',1,0.3,[]);
PlotBrainGraph_Kmeans(adj_kmeans_ASD_spec,sum(adj_kmeans_ASD_spec,2),CodeBook,0,max(sum(adj_kmeans_ASD_spec,2)),max(sum(adj_kmeans_ASD_spec,2)),2,2,'hot','gray',1,0.3,[]);
end
end
%%% SUMMARISING PLOTS OF CLUSTERS AND THEIR FEATURES
% Summarising temporal features
if is_plot > 0
for i = 1:length(idx_sign_K)
% Expression of the clusters over time for both populations
tmp_CM = colorcube(K_vector(idx_sign_K(i))+1);
figure;
set(gcf,'color','w');
title(['K=',num2str(K_vector(idx_sign_K(i)))]);
subplot(2,1,1);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(19:2:351,TC_CTR{idx_sign_K(i)}(j,:),'color',tmp_CM(j,:));
end
set(gca,'Box','off');
ylim([0 1]);
xlim([19 351]);
plot([19 351],repmat(Thresh_CTR{idx_sign_K(i)},1,2),'k--');
xlabel('Movie time [s]');
ylabel('Cluster expression (CTR)');
subplot(2,1,2);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(19:2:351,TC_ASD{idx_sign_K(i)}(j,:),'color',tmp_CM(j,:));
end
set(gca,'Box','off');
ylim([0 1]);
xlim([19 351]);
plot([19 351],repmat(Thresh_ASD{idx_sign_K(i)},1,2),'k--');
xlabel('Movie time [s]');
ylabel('Cluster expression (ASD)');
% The same for session discriminability
figure;
set(gcf,'color','w');
title(['K=',num2str(K_vector(idx_sign_K(i)))]);
subplot(2,1,1);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(19:2:351,TC_SES1{idx_sign_K(i)}(j,:),'color',tmp_CM(j,:));
end
set(gca,'Box','off');
ylim([0 1]);
xlim([19 351]);
plot([19 351],repmat(Thresh_SES1{idx_sign_K(i)},1,2),'k--');
xlabel('Movie time [s]');
ylabel('Cluster expression (SES1)');
subplot(2,1,2);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(19:2:351,TC_SES2{idx_sign_K(i)}(j,:),'color',tmp_CM(j,:));
end
set(gca,'Box','off');
ylim([0 1]);
xlim([19 351]);
plot([19 351],repmat(Thresh_SES2{idx_sign_K(i)},1,2),'k--');
xlabel('Movie time [s]');
ylabel('Cluster expression (SES2)');
% Correlation between clusters across groups
figure;
set(gcf,'color','w');
title(['K=',num2str(K_vector(idx_sign_K(i)))]);
imagesc(Corr_Mat{idx_sign_K(i)});
caxis([0 1]);
colormap(gray);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(repmat(j+0.5,1,2),[0.5 K_vector(idx_sign_K(i))+0.5],'color',[0.6 0.6 0.6],'LineWidth',0.5);
plot([0.5 K_vector(idx_sign_K(i))+0.5],repmat(j+0.5,1,2),'color',[0.6 0.6 0.6],'LineWidth',0.5);
end
xlabel('Cluster index');
ylabel('Cluster index');
figure;
set(gcf,'color','w');
title(['K=',num2str(K_vector(idx_sign_K(i)))]);
imagesc(Cov_Mat{idx_sign_K(i)});
colormap(jet);
hold on;
for j = 1:K_vector(idx_sign_K(i))
plot(repmat(j+0.5,1,2),[0.5 K_vector(idx_sign_K(i))+0.5],'color',[0.6 0.6 0.6],'LineWidth',0.5);
plot([0.5 K_vector(idx_sign_K(i))+0.5],repmat(j+0.5,1,2),'color',[0.6 0.6 0.6],'LineWidth',0.5);
end
xlabel('Cluster index');
ylabel('Cluster index');
% Temporal evolution of the clusters for subjects and groups
figure;
set(gcf,'color','w');
title(['K=',num2str(K_vector(idx_sign_K(i)))]);
tmp_plot = reshape(IDX{idx_sign_K(i)},167,(n_CTR+n_ASD));
imagesc(1:(n_CTR+n_ASD+5),19:2:351,[tmp_plot(:,Group_Label==1),K_vector(idx_sign_K(i))+1*ones(167,5),tmp_plot(:,Group_Label==2)]);
hold on;
caxis([0.5 K_vector(idx_sign_K(i))+1.5]);
colormap(colorcube(K_vector(idx_sign_K(i))+1));
xlabel('Subjects');
ylabel('Movie time [s]');
% Count differences summary for all clusters
tmp_counts = Counts{idx_sign_K(i)};
figure;
set(gcf,'color','w');
b1 = bar(mean(tmp_counts(:,Group_Label==1),2),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(mean(tmp_counts(:,Group_Label==2),2),'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
xlabel('Cluster index');
ylabel('Counts');
xlim([0 K_vector(idx_sign_K(i))+1]);
ylim([0 25]);
set(gca,'Box','off');
% Similarity matrix for the clusters (spatial similarity) for the
% CTR (top) and the ASD (bottom) cases
figure;
set(gcf,'color','w');
imagesc(Sim_Mat_CTR{idx_sign_K(i)});
caxis([0 1]);
colormap(jet);
hold on;
for j = 1:n_CTR_clusters(idx_sign_K(i))
plot(repmat(j+0.5,1,2),[0.5 n_CTR_clusters(idx_sign_K(i))+0.5],'color',[0.6 0.6 0.6],'LineWidth',0.5);
plot([0.5 n_CTR_clusters(idx_sign_K(i))+0.5],repmat(j+0.5,1,2),'color',[0.6 0.6 0.6],'LineWidth',0.5);
end
xlabel('Cluster index');
ylabel('Cluster index');
figure;
set(gcf,'color','w');
imagesc(Sim_Mat_ASD{idx_sign_K(i)});
caxis([0 1]);
colormap(jet);
hold on;
for j = 1:n_ASD_clusters(idx_sign_K(i))
plot(repmat(j+0.5,1,2),[0.5 n_ASD_clusters(idx_sign_K(i))+0.5],'color',[0.6 0.6 0.6],'LineWidth',0.5);
plot([0.5 n_ASD_clusters(idx_sign_K(i))+0.5],repmat(j+0.5,1,2),'color',[0.6 0.6 0.6],'LineWidth',0.5);
end
xlabel('Cluster index');
ylabel('Cluster index');
end
end
% PLOTS OF CLUSTERS THEMSELVES
for i = 1:length(idx_sign_K)
% Adjacency matrices of the clusters of interest (showing a
% group difference)
tmp_Centroids = Centroids{idx_sign_K(i)}';
tmp_Centroids_ext = zeros(n_conn,K_vector(idx_sign_K(i)));
tmp_Centroids_ext(idx_consistent,:) = tmp_Centroids;
% Adjacency matrices for the clusters
for j = 1:K_vector(idx_sign_K(i))
if ismember(j,idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(i),1})
if is_plot > 0
% Plotting the distribution of connections
% if tvec_kmeans_clus{idx_sign_K(i)}(j) > 0
% figure;
% set(gcf,'color','w');
% histogram(tmp_Centroids(:,j),-0.8:0.005:0.8,'FaceColor',[0 0.4 0.6],'EdgeColor','none');
% hold on;
% plot([mean(tmp_Centroids(:,j)),mean(tmp_Centroids(:,j))],get(gca,'YLim'),'k--');
% xlabel('ISFC value');
% ylabel('Counts');
% set(gca,'Box','off');
% else
% figure;
% set(gcf,'color','w');
% histogram(tmp_Centroids(:,j),-0.8:0.005:0.8,'FaceColor',[1 0.2 0.2],'EdgeColor','none');
% hold on;
% plot([mean(tmp_Centroids(:,j)),mean(tmp_Centroids(:,j))],get(gca,'YLim'),'k--');
% xlabel('ISFC value');
% ylabel('Counts');
% set(gca,'Box','off');
% end
intensity_mean(j) = mean(tmp_Centroids(:,j));
intensity_std(j) = std(tmp_Centroids(:,j));
disp(['Intensity mean and std for cluster ',num2str(j),': ',num2str(intensity_mean(j)),'+-',num2str(intensity_std(j))]);
% From 0 to 1 (most positive)
tmp_adj_kmeans = jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI) + jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI)';
% I want to extract the top nodes for each cluster
tmp_degree_kmeans = sum(abs(tmp_adj_kmeans),2);
[VAL_TOPNODES_KMEANS{i}(:,j),IDX_TOPNODES_KMEANS{i}(:,j)] = sort(tmp_degree_kmeans,'descend');
disp('Top nodes:');
IDX_TOPNODES_KMEANS{i}(VAL_TOPNODES_KMEANS{i}(:,j)>20,j)
% Makes negative values between -1 and 0
tmp_adj_kmeans_neg = tmp_adj_kmeans;
tmp_adj_kmeans_neg(tmp_adj_kmeans_neg > 0) = 0;
tmp_adj_kmeans_neg = tmp_adj_kmeans_neg/max(abs(tmp_adj_kmeans(:)));
% Makes positive values between 0 and + 1
tmp_adj_kmeans_pos = tmp_adj_kmeans;
tmp_adj_kmeans_pos(tmp_adj_kmeans_pos < 0) = 0;
tmp_adj_kmeans_pos = tmp_adj_kmeans_pos/max(abs(tmp_adj_kmeans(:)));
% Now I get something between -1 and +1
tmp_adj_plot = tmp_adj_kmeans_pos + tmp_adj_kmeans_neg;
% Now I get something thresholded between -1 and +1
tmp_adj_plot(abs(tmp_adj_plot) < T_display_kmeans) = 0;
if tvec_kmeans_clus{idx_sign_K(i)}(j) > 0
figure;
set(gcf,'color','w');
bar(1:n_ROI,sum(tmp_adj_kmeans,1),'FaceColor',[0 0.4 0.6],'EdgeColor','none');
ylim([-20 60]);
set(gca,'Box','off');
else
figure;
set(gcf,'color','w');
bar(1:n_ROI,sum(tmp_adj_kmeans,1),'FaceColor',[1 0.2 0.2],'EdgeColor','none');
ylim([-20 60]);
set(gca,'Box','off');
end
if is_plot > 1
% Plot of the cluster
PlotBrainGraph_Kmeans(tmp_adj_plot,...
sum(abs(tmp_adj_plot),2),CodeBook,0,...
max(sum(abs(tmp_adj_plot),2)),...
max(sum(abs(tmp_adj_plot),2)),2,2,...
'gray','spring',1,0.5,[-1 1]);
% Making videos for the adjacency changes in each group
if is_plot > 2
% Makes the video of the CTR and the ASD changes over time for two views:
% sagittal and horizontal
[X_CTR,X_ASD] = MakeGraphVideo(Adj,CodeBook,Group_Label,2,TC_CTR{1},TC_ASD{1},Thresh_CTR{1},Thresh_ASD{1});
[X_CTR2,X_ASD2] = MakeGraphVideo(Adj,CodeBook,Group_Label,1,TC_CTR{1},TC_ASD{1},Thresh_CTR{1},Thresh_ASD{1});
% Controls, sagittal
MyVid = VideoWriter('CTR_VIDEO_SAGITTAL','Uncompressed AVI');
MyVid.FrameRate = 0.5;
MOV_CTR = immovie(X_CTR);
open(MyVid);
writeVideo(MyVid,MOV_CTR);
close(MyVid);
% Controls, horizontal
MyVid2 = VideoWriter('CTR_VIDEO_HORIZONTAL','Uncompressed AVI');
MyVid2.FrameRate = 0.5;
MOV_CTR2 = immovie(X_CTR2);
open(MyVid2);
writeVideo(MyVid2,MOV_CTR2);
close(MyVid2);
% ASD, sagittal
MyVid3 = VideoWriter('ASD_VIDEO_SAGITTAL','Uncompressed AVI');
MyVid3.FrameRate = 0.5;
MOV_ASD = immovie(X_ASD);
open(MyVid3);
writeVideo(MyVid3,MOV_ASD);
close(MyVid3);
% ASD, horizontal
MyVid4 = VideoWriter('ASD_VIDEO_HORIZONTAL','Uncompressed AVI');
MyVid4.FrameRate = 0.5;
MOV_ASD2 = immovie(X_ASD2);
open(MyVid4);
writeVideo(MyVid4,MOV_ASD2);
close(MyVid4);
end
end
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E. K-means clustering: data-driven isolation of relevant time points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The goal in this part is to use the features computed above to select the
% most interesting time points from various perspectives
n_clusters = length(K_vector);
n_timepoints = length(19:2:351);
% Group consistancy for the CTR and the ASD group
Kmeans_consistency_CTR = zeros(n_timepoints,n_clusters);
Kmeans_consistency_ASD = zeros(n_timepoints,n_clusters);
Kmeans_consistency_SES1 = zeros(n_timepoints,n_clusters);
Kmeans_consistency_SES2 = zeros(n_timepoints,n_clusters);
% Group discriminability
Kmeans_discriminability_group = zeros(n_timepoints,n_clusters);
Kmeans_similarity_group = zeros(n_timepoints,n_clusters);
% Session discriminability
Kmeans_discriminability_session = zeros(n_timepoints,n_clusters);
Kmeans_similarity_session = zeros(n_timepoints,n_clusters);
% For each time point and each cluster configuration, I want to pick the
for i = 1:length(K_vector)
for t = 1:167
if max(TC_CTR{i}(:,t)) > Thresh_CTR{i}
Kmeans_consistency_CTR(t,i) = 1;%max(TC_CTR{i}(:,t));
tmp = find(TC_CTR{i}(:,t) == max(TC_CTR{i}(:,t)));
idx_consistency_CTR(t,i) = tmp(1);
else
idx_consistency_CTR(t,i) = 0;
end
if max(TC_ASD{i}(:,t)) > Thresh_ASD{i}
Kmeans_consistency_ASD(t,i) = 1;%max(TC_ASD{i}(:,t));
tmp = find(TC_ASD{i}(:,t) == max(TC_ASD{i}(:,t)));
idx_consistency_ASD(t,i) = tmp(1);
else
idx_consistency_ASD(t,i) = 0;
end
if max(TC_SES1{i}(:,t)) > Thresh_SES1{i}
Kmeans_consistency_SES1(t,i) = 1;%max(TC_SES1{i}(:,t));
tmp = find(TC_SES1{i}(:,t) == max(TC_SES1{i}(:,t)));
idx_consistency_SES1(t,i) = tmp(1);
else
idx_consistency_SES1(t,i) = 0;
end
if max(TC_SES2{i}(:,t)) > Thresh_SES2{i}
Kmeans_consistency_SES2(t,i) = 1;%max(TC_SES2{i}(:,t));
tmp = find(TC_SES2{i}(:,t) == max(TC_SES2{i}(:,t)));
idx_consistency_SES2(t,i) = tmp(1);
else
idx_consistency_SES2(t,i) = 0;
end
% If both the CTR and ASD groups have a cluster than stands out at
% the particular time point and parcellation degree...
if idx_consistency_CTR(t,i) > 0 && idx_consistency_ASD(t,i) > 0
% Then I fill the discriminability or similarity matrices
% appropriately
if idx_consistency_CTR(t,i) == idx_consistency_ASD(t,i)
Kmeans_similarity_group(t,i) = (Kmeans_consistency_CTR(t,i) + Kmeans_consistency_ASD(t,i))/2;
else
Kmeans_discriminability_group(t,i) = (Kmeans_consistency_CTR(t,i) + Kmeans_consistency_ASD(t,i))/2;
end
else
end
% Same for sessions
if idx_consistency_SES1(t,i) > 0 && idx_consistency_SES2(t,i) > 0
% Then I fill the discriminability or similarity matrices
% appropriately
if idx_consistency_SES1(t,i) == idx_consistency_SES2(t,i)
Kmeans_similarity_session(t,i) = (Kmeans_consistency_SES1(t,i) + Kmeans_consistency_SES2(t,i))/2;
else
Kmeans_discriminability_session(t,i) = (Kmeans_consistency_SES1(t,i) + Kmeans_consistency_SES2(t,i))/2;
end
else
end
end
end
% I can now visually locate the moments of significance
t_consistency_CTR = [2,46,58,73,98,113,148,163];
t_consistency_ASD = [1,19,33,44,100,115,141,161];
t_discriminability = [19,44,99,115,142,161];
for i = 1:length(t_discriminability)
% tmp_frame contains the frame of interest information for all the
% subjects
for s = 1:(n_CTR+n_ASD)
tmp_frame(:,s) = ISFC_MW{s}(idx_consistent,t_discriminability(i));
end
for con = 1:n_conn2
[p_anova_kmeans(i,con,:),table] = anovan(tmp_frame(con,:),{Group_Label' Session_Label'},'model',2,'varnames',{'Group_Label','Session_Label'},'display','none');
[~,~,~,stats] = ttest2(tmp_frame(con,Group_Label==1)',tmp_frame(con,Group_Label==2)');
fvec_kmeans(i,con,:) = [table{2,6};table{3,6};table{4,6}];
tvec_kmeans(i,con) = stats.tstat;
clear table
clear stats
end
end
% I can determine what values it matches in terms of F-statistic
f_sup_kmeans = finv(1-alpha_ncon/length(t_discriminability),3,27);
for i = 1:length(t_discriminability)
idx_sign_kmeans_ANOVA_BON{i,1} = find(squeeze(fvec_kmeans(i,:,1)) > f_sup_kmeans);
idx_sign_kmeans_ANOVA_BON{i,2} = find(squeeze(fvec_kmeans(i,:,2)) > f_sup_kmeans);
idx_sign_kmeans_ANOVA_BON{i,3} = find(squeeze(fvec_kmeans(i,:,3)) > f_sup_kmeans);
end
if is_plot > 0
% For each case, I plot the consistency measure in time x clusters
% fashion, and then, the 'time course of consistency' across cluster
% number, and the same across time
figure;
title('Consistency in CTR group');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_consistency_CTR);
caxis([0 1]);
colormap(jet);
figure;
title('Consistency in ASD group');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_consistency_ASD);
caxis([0 1]);
colormap(jet);
figure;
title('Similarity in groups');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_similarity_group);
caxis([0 1]);
colormap(jet);
figure;
title('Discriminability in groups');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_discriminability_group);
caxis([0 1]);
colormap(jet);
figure;
set(gcf,'color','w');
subplot(2,1,1);
plot(19:2:351,sum(logical(Kmeans_consistency_CTR),2),'color','b');
hold on;
plot(19:2:351,sum(logical(Kmeans_consistency_ASD),2),'color','r');
plot(19:2:351,sum(logical(Kmeans_similarity_group),2),'color',[0.5 0.5 0.5]);
plot(19:2:351,sum(logical(Kmeans_discriminability_group),2),'color','k');
legend('CTR consistency','ASD consistency','Group similarity','Group discriminability');
set(gca,'Box','off');
xlim([19 351]);
subplot(2,1,2);
plot(K_vector,sum(logical(Kmeans_consistency_CTR),1),'color','b');
hold on;
plot(K_vector,sum(logical(Kmeans_consistency_ASD),1),'color','r');
plot(K_vector,sum(logical(Kmeans_similarity_group),1),'color',[0.5 0.5 0.5]);
plot(K_vector,sum(logical(Kmeans_discriminability_group),1),'color','k');
legend('CTR consistency','ASD consistency','Group similarity','Group discriminability');
set(gca,'Box','off');
xlim([10 150]);
figure;
title('Consistency in SES1');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_consistency_SES1);
caxis([0 1]);
colormap(jet);
figure;
title('Consistency in SES2');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_consistency_SES2);
caxis([0 1]);
colormap(jet);
figure;
title('Similarity in sessions');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_similarity_session);
caxis([0 1]);
colormap(jet);
figure;
title('Discriminability in sessions');
set(gcf,'color','w');
imagesc(K_vector,19:2:351,Kmeans_discriminability_session);
caxis([0 1]);
colormap(jet);
figure;
set(gcf,'color','w');
subplot(2,1,1);
plot(19:2:351,sum(logical(Kmeans_consistency_SES1),2),'color','b');
hold on;
plot(19:2:351,sum(logical(Kmeans_consistency_SES2),2),'color','r');
plot(19:2:351,sum(logical(Kmeans_similarity_session),2),'color',[0.5 0.5 0.5]);
plot(19:2:351,sum(logical(Kmeans_discriminability_session),2),'color','k');
legend('SES1 consistency','SES2 consistency','SES similarity','SES discriminability');
set(gca,'Box','off');
xlim([19 351]);
subplot(2,1,2);
plot(K_vector,sum(logical(Kmeans_consistency_SES1),1),'color','b');
hold on;
plot(K_vector,sum(logical(Kmeans_consistency_SES2),1),'color','r');
plot(K_vector,sum(logical(Kmeans_similarity_session),1),'color',[0.5 0.5 0.5]);
plot(K_vector,sum(logical(Kmeans_discriminability_session),1),'color','k');
legend('SES1 consistency','SES2 consistency','SES similarity','SES discriminability');
set(gca,'Box','off');
xlim([10 150]);
if is_plot > 1
for i = 1:length(t_discriminability)
% Plotting the connections significantly different for the retained
% time points
% tmp_plot is the vector filled with 1 for CTR > ASD and -1 with
% ASD > CTR
tmp_plot = zeros(n_conn2,1);
tmp_plot(idx_sign_kmeans_ANOVA_BON{i,1}) = sign(tvec_kmeans(i,idx_sign_kmeans_ANOVA_BON{i,1}))*1;
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_plot;
% With a value of 1.0, I have blue (CTR > ASD) and with a value of
% 0.97, I have red (ASD > CTR)
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
tmp_plot2(tmp_plot2 == 1) = 1.0;
tmp_plot2(tmp_plot2 == -1) = 0.97;
PlotBrainGraph_Kmeans(tmp_plot2,sum(tmp_plot2,2),CodeBook,0,max(sum(tmp_plot2,2)),max(sum(tmp_plot2,2)),2,2,'hot','lines',1,0.5,[]);
end
for i = 1:length(t_consistency_CTR)
clear tmp_frame
for s = 1:(n_CTR+n_ASD)
tmp_frame(:,s) = ISFC_MW{s}(idx_consistent,t_consistency_CTR(i));
end
tmp_frame = tmp_frame(:,Group_Label==1);
tmp_frame = mean(tmp_frame,2);
tmp_frame(abs(tmp_frame) < 0.2) = 0;
tmp_frame = tmp_frame/max(abs(tmp_frame));
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_frame;
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
PlotBrainGraph_Kmeans(tmp_plot2,sum(abs(tmp_plot2),2),CodeBook,0,max(sum(abs(tmp_plot2),2)),max(sum(abs(tmp_plot2),2)),2,2,'hot','jet',1,0.5,[-1 1]);
end
for i = 1:length(t_consistency_ASD)
clear tmp_frame
for s = 1:(n_CTR+n_ASD)
tmp_frame(:,s) = ISFC_MW{s}(idx_consistent,t_consistency_ASD(i));
end
tmp_frame = tmp_frame(:,Group_Label==2);
tmp_frame = mean(tmp_frame,2);
tmp_frame(abs(tmp_frame) < 0.2) = 0;
tmp_frame = tmp_frame/max(abs(tmp_frame));
tmp_plot_ext = zeros(n_conn,1);
tmp_plot_ext(idx_consistent) = tmp_frame;
tmp_plot2 = jVecToUpperTriMat(tmp_plot_ext,n_ROI)+jVecToUpperTriMat(tmp_plot_ext,n_ROI)';
PlotBrainGraph_Kmeans(tmp_plot2,sum(abs(tmp_plot2),2),CodeBook,0,max(sum(abs(tmp_plot2),2)),max(sum(abs(tmp_plot2),2)),2,2,'hot','jet',1,0.5,[-1 1]);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% F. Relationship between ASD scores and clustering outcomes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creates my matrix of scores (IQ ADib ADicnv ADicv ADid ADitot1 ADitot2
% AQ)
Scores_ASD_all = [50 2 16 12 NaN 7 35 NaN;...
90 3 24 11 17 8 43 NaN;...
90 3 24 11 17 8 43 NaN;...
82 4 13 7 NaN 6 26 NaN;...
82 4 13 7 NaN 6 26 NaN;...
66 28 11 17 5 44 61 28;...
66 28 11 17 5 44 61 28;...
124 22 14 21 11 47 68 33;...
124 22 14 21 11 47 68 33;...
57 23 5 16 10 38 54 26;...
57 23 5 16 10 38 54 26;...
31 34 14 25 8 56 81 31;...
120 25 17 22 2 44 66 21;...
120 25 17 22 2 44 66 21;...
91 25 14 18 5 44 62 27;...
91 25 14 18 5 44 62 27];
Scores_ASD_all = Scores_ASD_all(:,[1:4,6:7]);
% Session label only for ASD subjects
SES_ASD = Session_Label(Group_Label == 2);
% I will do the check for all values of K
for id = 1:length(idx_sign_K)
tmp_counts = Counts{idx_sign_K(id)};
for ses = 1:2
% Counts for only ASD subjects
Counts_ASD = tmp_counts(:,Group_Label==2);
Counts_ASD = Counts_ASD(:,SES_ASD==ses);
% Gets all the significant clusters
Counts_ASD = Counts_ASD(idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(id),1},:);
% Only keeps the ones from the ASD group
Counts_ASD = Counts_ASD(tvec_kmeans_clus{idx_sign_K(id)}(idx_sign_kmeans_clus_ANOVA_BON{idx_sign_K(id),1}) < 0,:);
Scores_ASD = Scores_ASD_all(SES_ASD==ses,:)';
% pairs of subjects
Pairs_ASD = nchoosek((1:sum(SES_ASD==ses))',2);
% For each pair of subject, computation of the distance in the
% scores and in the counts
[D_Counts{ses,id},D_Scores{ses,id},D_Counts_Global{ses,id},D_Scores_Global{ses,id}] = Compute_Distances(Counts_ASD,Scores_ASD,Pairs_ASD);
end
% Assessment of the significance of the results (permutation testing)
% Number of tests carried out (2 = number of global metric computation
% types * 6 = number of scores (present in all subjects) * 2 (sessions) * number of cluster values)
n_tests = 5*1*2*1;
for ses = 1:2
for i = 1:2
for j = 1:6
[R{ses}(i,j),pval{ses}(i,j)] = corr(D_Scores{ses,id}(:,j,1),D_Counts_Global{ses,id}(:,i));
pval{ses}(i,j) = n_tests*pval{ses}(i,j);
Corr_scores{ses}(i,j,id) = R{ses}(i,j);
Pval_scores{ses}(i,j,id) = pval{ses}(i,j);
for n = 1:10000
tmp1 = D_Scores{ses,id}(:,j,1);
tmp2 = D_Counts_Global{ses,id}(:,i);
tmp1 = tmp1(randperm(length(tmp1)));
R_null{ses}(i,j,n) = corr(tmp1,tmp2);
R_thresh{ses}(i,j) = prctile(R_null{ses}(i,j,:),100-5/n_tests);
end
end
end
end
% Plotting of the significant results
if is_plot > 0
for ses = 1:2
for i = 1:2
for j = 1:6
if R{ses}(i,j) > R_thresh{ses}(i,j)
disp(['Significant relationship for session ',...
num2str(ses),', global measure ',...
num2str(i),', score ',num2str(j),...
' and parcellation order ',num2str(id)]);
figure;
set(gcf,'color','w');
scatter(D_Scores{ses,id}(:,j,1),D_Counts_Global{ses,id}(:,i));
end
end
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% G. Comparison between clusters expressed at the same time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pairs of clusters expressed at the same time in both populations
Cluster_pairs = [52,28; 58,2; 50,44;11,55];
for id = 1:length(idx_sign_K)
% Samples the clusters to plot
tmp_Centroids_ext = zeros(n_conn,K_vector(idx_sign_K(id)));
tmp_Centroids_ext(idx_consistent,:) = Centroids{idx_sign_K(id)}';
% Constructs the adjacency matrix of each cluster
for j = 1:K_vector(idx_sign_K(id))
adj_kmeans(:,:,j) = jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI) + jVecToUpperTriMat(tmp_Centroids_ext(:,j),n_ROI)';
end
% All pairs of interest
for p = 1:size(Cluster_pairs,1)
% Positive values denote the stronger connections in the controls
padj_kmeans = squeeze(adj_kmeans(:,:,Cluster_pairs(p,1)) - adj_kmeans(:,:,Cluster_pairs(p,2)));
% Makes negative values between -1 and 0
tmp_adj_kmeans_neg = padj_kmeans;
tmp_adj_kmeans_neg(tmp_adj_kmeans_neg > 0) = 0;
tmp_adj_kmeans_neg = tmp_adj_kmeans_neg/max(abs(padj_kmeans(:)));
[~,idx_clus_ASD{p}] = sort(sum(abs(tmp_adj_kmeans_neg),1),'descend');
deg_clus_ASD{p} = sum(abs(tmp_adj_kmeans_neg),1);
% Makes positive values between 0 and + 1
tmp_adj_kmeans_pos = padj_kmeans;
tmp_adj_kmeans_pos(tmp_adj_kmeans_pos < 0) = 0;
tmp_adj_kmeans_pos = tmp_adj_kmeans_pos/max(abs(padj_kmeans(:)));
[~,idx_clus_CTR{p}] = sort(sum(abs(tmp_adj_kmeans_pos),1),'descend');
deg_clus_CTR{p} = sum(abs(tmp_adj_kmeans_pos),1);
% Now I get something between -1 and +1
tmp_adj_plot = tmp_adj_kmeans_pos + tmp_adj_kmeans_neg;
% Now I get something thresholded between -1 and +1
tmp_adj_plot(abs(tmp_adj_plot) < 0.5) = 0;
figure;
set(gcf,'color','w');
b1 = bar(deg_clus_CTR{p},'FaceColor',[0 0.4 0.6],'EdgeColor','none');
set(get(b1,'Children'),'FaceAlpha',0.2);
hold on;
b2 = bar(deg_clus_ASD{p},'FaceColor','none','EdgeColor',[1 0.2 0.2]);
set(get(b2,'Children'),'FaceAlpha',0.2);
set(gca,'Box','off');
xlabel('ROIs');
ylabel('Degree');
% Plot of the cluster
PlotBrainGraph_Kmeans(tmp_adj_plot,sum(abs(tmp_adj_plot),2),CodeBook,0,max(sum(abs(tmp_adj_plot),2)),max(sum(abs(tmp_adj_plot),2)),2,2,'gray','spring',1,0.5,[-1 1]);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% H. Trajectories of the subjects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creates a new matrix with concatenated data, size n_conn x n_dp
Data_graphs = [];
for t = 1:167
for s = 1:(n_CTR+n_ASD)
Data_graphs = [Data_graphs,SignMat2(:,t,s)];
end
end
Dist_graphs = dist(Data_graphs);
Sigma_graphs = mean(jUpperTriMatToVec(Dist_graphs));
for i = 1:size(Data_graphs,2)
for j = 1:size(Data_graphs,2)
W_graphs(i,j) = exp(-(Dist_graphs(i,j)^2)./Sigma_graphs^2);
end
end
W_graphs(diag(W_graphs)) = 0;
% Computation of the normalized Laplacian of the graph, and of its
% eigenvalues (E) and eigenvectors (V)
L_graphs = compute_graph_laplacian(W_graphs);
L_graphs = L_graphs + 1e-6* speye(size(L_graphs));
[V_graphs,E_graphs] = eigs(L_graphs,4,'SM');
E_graphs = diag(E_graphs);
% Eigenvalues and eigenvectors are sorted from the smallest to the largest
% eigenvalue
[~,idx_graphs] = sort(E_graphs,'ascend');
E_graphs = E_graphs(idx_graphs);
V_graphs = V_graphs(:,idx_graphs);
%% 12. Can I find a difference in RS as well ?
disp('Assessing resting-state correlates of what was found in the ISFC analysis...');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A. Computation of the metrics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sampling the connections of interest
for i = 1:(n_CTR+n_ASD)
dFC_RS_red{i} = dFC_RS{i}(idx_consistent,:);
end
% Dictionary of building blocks that we want to use (size n_blocks x
% n_conn2)
Dico = Centroids{1};
disp('Approach 1: correlating with the building blocks');
% First approach: correlating the bulding blocks with the time courses for
% each subject
for i = 1:size(Dico,1)
for j = 1:(n_CTR+n_ASD)
% Contains the correlation time courses
RS_corr_TC(i,j,:) = corr(Dico(i,:)',dFC_RS_red{j});
% Absolute valued sum of elements
RS_corr_SAVE(i,j) = sum(abs(RS_corr_TC(i,j,:)));
end
end
disp('Approach 2: cosine distance with the building blocks');
for i = 1:size(Dico,1)
for j = 1:(n_CTR+n_ASD)
for t = 1:size(dFC_RS_red{j},2)
% Contains the cosine distance computation outcomes
RS_cos_TC(i,j,t) = pdist([Dico(i,:);dFC_RS_red{j}(:,t)'],'cosine');
end
% Absolute valued sum of elements
RS_cos_SAVE(i,j) = sum(abs(RS_cos_TC(i,j,:)));
end
end
disp('Approach 3: winner-take-all approach');
WTA_idx_corr = nan(size(dFC_RS_red{j},2),(n_CTR+n_ASD));
for j = 1:(n_CTR+n_ASD)
for t = 1:size(dFC_RS_red{j},2)
% 20 values, and I want the 'best match'
tmp = corr(Dico',dFC_RS_red{j}(:,t));
WTA_idx_corr(t,j) = find(tmp == max(tmp));
end
end
disp('Approach 4: winner-take-all approach (cosine)');
WTA_idx_cos = nan(size(dFC_RS_red{j},2),(n_CTR+n_ASD));
for j = 1:(n_CTR+n_ASD)
for t = 1:size(dFC_RS_red{j},2)
idx_best = 0;
tmp_best = realmax;
for i = 1:size(Dico,1)
% 20 values, and I want the 'best match'
tmp = pdist([Dico(i,:);dFC_RS_red{j}(:,t)'],'cosine');
if tmp < tmp_best
tmp_best = tmp;
idx_best = i;
end
WTA_idx_cos(t,j) = idx_best;
end
end
end
if is_plot > 0
figure;
set(gcf,'color','w');
imagesc(1:(n_CTR+n_ASD+5),1:size(dFC_RS_red{j},2),[WTA_idx_corr(:,Group_Label==1),21*ones(size(dFC_RS_red{j},2),5),WTA_idx_corr(:,Group_Label==2)]);
hold on;
caxis([0.5 20+1.5]);
colormap(colorcube(21));
xlabel('Subjects');
ylabel('Time');
figure;
set(gcf,'color','w');
imagesc(1:(n_CTR+n_ASD+5),1:size(dFC_RS_red{j},2),[WTA_idx_cos(:,Group_Label==1),21*ones(size(dFC_RS_red{j},2),5),WTA_idx_cos(:,Group_Label==2)]);
hold on;
caxis([0.5 20+1.5]);
colormap(colorcube(21));
xlabel('Subjects');
ylabel('Time');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% B. Assessment of statistical significance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Approach 1: t-test to compare the populations
for i = 1:size(Dico,1)
[~,pvec_RS_corr(i),~,stats] = ttest2(RS_corr_SAVE(i,Group_Label==1)',RS_corr_SAVE(i,Group_Label==2)');
tvec_RS_corr(i) = stats.tstat;
clear stats
end
% Approach 2: same thing again
for i = 1:size(Dico,1)
[~,pvec_RS_cos(i),~,stats] = ttest2(RS_cos_SAVE(i,Group_Label==1)',RS_cos_SAVE(i,Group_Label==2)');
tvec_RS_cos(i) = stats.tstat;
clear stats
end
RS_counts = ComputeClusterCounts(WTA_idx_corr);
% Approach 2: same thing again
for i = 1:size(Dico,1)
[~,pvec_RS_WTA(i),~,stats] = ttest2(RS_counts(i,Group_Label==1)',RS_counts(i,Group_Label==2)');
tvec_RS_WTA(i) = stats.tstat;
clear stats
end
RS_counts_sum1 = RS_counts(1,:) + RS_counts(3,:) + RS_counts(13,:) + RS_counts(14,:) + RS_counts(15,:) + RS_counts(16,:) + RS_counts(20,:);
RS_counts_sum2 = RS_counts(7,:) + RS_counts(12,:) + RS_counts(17,:) + RS_counts(18,:);
if is_plot > 0
figure;
set(gcf,'color','w');
b1 = bar([RS_counts_sum1(Group_Label==1),RS_counts_sum1(Group_Label==2)],'FaceColor',[0 0.4 0.6],'EdgeColor','none');
hold on;
b2 = bar([RS_counts_sum2(Group_Label==1),RS_counts_sum2(Group_Label==2)],'FaceColor','none','EdgeColor',[1 0.2 0.2]);
xlabel('Subject');
ylabel('Counts for CTR-specific or ASD-specific stuff');
xlim([0 n_CTR+n_ASD+1]);
set(gca,'Box','off');
end
c4science · Help