Diffusion ISFC clustering (master) Browse

ISFCClustering/Analysis.m
b90a4357cb97master

Analysis.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 = 512*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

ISFCClustering/Analysis.mb90a4357cb97master

Analysis.mView Options

ISFCClustering/Analysis.m
b90a4357cb97master

Analysis.m
View Options