Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/Matlab_Scripts/Script_NI_SI_SHINEY.m
257e25cb59ccmaster

Last Change
Raw File

Script_NI_SI_SHINEY.m
View Options

%% This script performs all analyses described in [ADD PAPER NAME]
%
% The content includes:
%% 0. Data loading
% Loads all the data of importance for what follows
% ISFC time courses upon movie-watching and rest for W = 10
load('Results_W10B.mat');
load('Results_W10.mat');
% Group labels (1 = TD, 2 = ASD)
load('Group_Label');
load('Group_Label_pureRS');
% Session labels
load('Session_Label');
% Number of TD and ASD sessions
load('n_CTR');
load('n_ASD');
% Framewise displacement time courses
load('FD');
% IQ scores
load('IQ');
% ASD scores
load('AQ');
% Connection indices and code book
load('C1');
load('C2');
load('CodeBook');
load('COG');
% Number of connections, time points and ROI
n_conn = 36315;
n_ROI = 270;
n_TP = 167;
% Average FD on non-scrubbed frames per subject
load('Motion_regressor');
% Euclidean distance between all connections considered
for i = 1:n_conn
D_euclidean(i,1) = sqrt(sum((squeeze(COG(:,C1(i)))-squeeze(COG(:,C2(i)))).^2));
end
%% 1. Assessment of IQ-related connections
% In this section, the goal is to determine what amount of connections
% should be retained in the following analyses
% Range of alpha-values to investigate
alpha_range = (100-[60,65,70,75,80,85,90,95,96,97,98,99,99.5,99.75,99.8,99.85,99.9])/100;
% For each, ISFC excursions are computed; then, we consider the case in
% which average FD (non scrubbed frames) is not regressed out (IOI(:,1,:)),
% or is regressed out (IOI(:,2,:)). IOI contains the indices of non
% IQ-related connections
for a = 1:length(alpha_range)
[MOI(:,a),IOI(:,1,a),IOI(:,2,a)] = ComputeIQConnections(ISFC_MW,...
ISFC_RS,ISFC_PRS,Group_Label,Group_Label_pureRS,IQ,...
Motion_regressor,alpha_range(a));
end
% Both cases
n_nonIQ_nomotreg = sum(squeeze(IOI(:,1,:)),1);
n_nonIQ_motreg = sum(squeeze(IOI(:,2,:)),1);
% Plots the results
figure;
set(gcf,'color','w');
plot(alpha_range,n_nonIQ_motreg,'color',[0.6 0.6 0.6],'LineWidth',2);
%hold on;
%plot(alpha_range,n_nonIQ_nomotreg,'color',[0.3 0.3 0.3],'LineWidth',2);
set(gca,'Box','off');
% I now want to display example slices showing the retained connections for
% three indicative alpha values
% alpha = 99.75
tmp = jVecToUpperTriMat(IOI(:,2,14),n_ROI);
tmp2 = jVecToUpperTriMat(IOI(:,2,7),n_ROI);
tmp3 = jVecToUpperTriMat(IOI(:,2,2),n_ROI);
division_factor = max([sum(tmp,2);sum(tmp2,2);sum(tmp3,2)]);
PlotBrainGraph_Kmeans_withtitle(tmp,sum(tmp,2),CodeBook,2,division_factor/4,...
division_factor,1,1,'hot','gray',1,0.1,[],'FIG1','ALPHA99DOT75_V1');
PlotBrainGraph_Kmeans_withtitle(tmp,sum(tmp,2),CodeBook,2,division_factor/4,...
division_factor,1,2,'hot','gray',1,0.1,[],'FIG1','ALPHA99DOT75_V2');
PlotBrainGraph_Kmeans_withtitle(tmp,sum(tmp,2),CodeBook,2,division_factor/4,...
division_factor,1,3,'hot','gray',1,0.1,[],'FIG1','ALPHA99DOT75_V3');
% alpha = 90
PlotBrainGraph_Kmeans_withtitle(tmp2,sum(tmp2,2),CodeBook,2,division_factor/4,...
division_factor,1,1,'hot','gray',1,0.1,[],'FIG1','ALPHA90_V1');
PlotBrainGraph_Kmeans_withtitle(tmp2,sum(tmp2,2),CodeBook,2,division_factor/4,...
division_factor,1,2,'hot','gray',1,0.1,[],'FIG1','ALPHA90_V2');
PlotBrainGraph_Kmeans_withtitle(tmp2,sum(tmp2,2),CodeBook,2,division_factor/4,...
division_factor,1,3,'hot','gray',1,0.1,[],'FIG1','ALPHA90_V3');
% alpha = 65
PlotBrainGraph_Kmeans_withtitle(tmp3,sum(tmp3,2),CodeBook,2,division_factor/4,...
division_factor,1,1,'hot','gray',1,0.1,[],'FIG1','ALPHA65_V1');
PlotBrainGraph_Kmeans_withtitle(tmp3,sum(tmp3,2),CodeBook,2,division_factor/4,...
division_factor,1,2,'hot','gray',1,0.1,[],'FIG1','ALPHA65_V2');
PlotBrainGraph_Kmeans_withtitle(tmp3,sum(tmp3,2),CodeBook,2,division_factor/4,...
division_factor,1,3,'hot','gray',1,0.1,[],'FIG1','ALPHA65_V3');
%% 2. Assessment of clustering quality
% In this section, we want to determine the optimal number of clusters to
% use for ISFC state extraction
% Will contain all the frames of interest (n_datapoints x n_dimensions)
Data_kmeans = [];
% Associated group, subject and time labels
GL_kmeans = [];
SL_kmeans = [];
TL_kmeans = [];
% Indices of the connections to consider
IDX = logical(squeeze(IOI(:,2,7)));
% 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 information of interest
tmp = ISFC_MW{s}(logical(IDX),t);
% Filling of data matrix and associated summarizing vectors
Data_kmeans = [Data_kmeans;tmp'];
GL_kmeans = [GL_kmeans;Group_Label(s)];
SL_kmeans = [SL_kmeans;s];
TL_kmeans = [TL_kmeans;t];
end
end
% As a way to evaluate the quality of a given clustering scheme, we perform
% consensus clustering
K_range = 10:5:250;
CC = CAP_ConsensusClustering({Data_kmeans'},K_range,'items',0.8,50,'correlation');
[CDF,Lorena] = ComputeClusteringQuality(CC,K_range);
figure;
set(gcf,'color','w');
bar(Lorena);
set(gca,'Box','off');
% Trying to reveal things more clearly by differentiating PAC
Daniela = diff(Lorena);
% Colormaop
CM = cbrewer('div','Spectral',20);
figure;
set(gcf,'color','w');
hold on;
for i = 1:size(Daniela,2)
plot(K_range(2:end),Daniela(:,i),'LineWidth',2,'color',CM(i,:));
end
set(gca,'Box','off');
%% 3. Computation of graph metrics
% In this section, ISGC graph theoretical metrics are computed on the
% successive sliding window ISFC measurement sets (separately for positive
% and negative ISFC elements). At the nodel level, the clustering
% coefficient, the efficiency and the betweenness centrality were computed.
% At the network level, modularity was computed. In addition, eigenvalues
% of the graph Laplacian were also computed.
for t = 1:n_TP
t
% For each subject...
for s = 1:(n_CTR + n_ASD)
% Takes the ISFC graph (full)
tmp_frame(:,s) = ISFC_MW{s}(:,t);
% Sets to zero the IQ-related connections
tmp_frame(~IDX,s) = 0;
% Converts into a matrix form
tmp_graph(:,:,s) = jVecToUpperTriMat(tmp_frame(:,s),n_ROI);
% Considers positive and negative data separately
tmpg1 = tmp_graph;
tmpg1(tmpg1 < 0) = 0;
tmpg2 = tmp_graph;
tmpg2(tmpg2 > 0) = 0;
tmpg2 = abs(tmpg2);
L1 = diag(sum(squeeze(tmpg1(:,:,s))))-squeeze(tmpg1(:,:,s));
L2 = diag(sum(squeeze(tmpg2(:,:,s))))-squeeze(tmpg2(:,:,s));
% Computes the graph metrics
EIV_plus(:,s,t) = eig(L1);
Deg_plus(:,s,t) = sum(tmpg1(:,:,s));
CC_plus(:,s,t) = clustering_coef_wu(tmpg1(:,:,s));
E_plus(:,s,t) = efficiency_wei(tmpg1(:,:,s),2);
[~,Q_plus(s,t)] = modularity_und(tmpg1(:,:,s));
Betw_plus(:,s,t) = betweenness_wei(tmpg1(:,:,s));
EIV_minus(:,s,t) = eig(L2);
Deg_minus(:,s,t) = sum(tmpg2(:,:,s));
CC_minus(:,s,t) = clustering_coef_wu(tmpg2(:,:,s));
E_minus(:,s,t) = efficiency_wei(tmpg2(:,:,s),2);
[~,Q_minus(s,t)] = modularity_und(tmpg2(:,:,s));
Betw_minus(:,s,t) = betweenness_wei(tmpg2(:,:,s));
end
clear tmp_graph
clear tmp_frame
end
%% 4. Assessment of significance for the different graph metrics
% We performed this separately for positive and negative ISFC elements,
% using the method suggested by Meskaldji and colleagues (NeuroImage,
% 2015), where p-values obtained upon Bonferroni correction are relaxed
% Eigenvalues
% The eigenvalues of the functional graphs reflect "frequencies" of the
% associated elementary building blocks. We decided to use the median
% eigenvalue as a measure of the locality/globality of processing
% We compute the median frequency
Medeig_plus = squeeze(median(EIV_plus,1));
Medeig_minus = squeeze(median(EIV_minus,1));
% For our statistical assessments, we consider 167 parallel tests. First,
% we apply Bonferroni correction, and then we relax the obtained p-values
% using the method proposed by Meskaldji et al., in order to take into
% account positive dependencies between successive sliding window
% estimates.
[P_final_eigplus] = ComputePValuesNetworkFeature(Medeig_plus,Group_Label,9999,n_TP);
[P_final_eigminus] = ComputePValuesNetworkFeature(Medeig_minus,Group_Label,9999,n_TP);
figure;
set(gcf,'color','w');
plot(time,P_final_eigplus,'LineWidth',2,'color','r');
hold on;
plot(time,P_final_eigminus,'LineWidth',2,'color','b');
% The same is performed for modularity
[P_final_Qplus] = ComputePValuesNetworkFeature(Q_plus,Group_Label,9999,n_TP);
[P_final_Qminus] = ComputePValuesNetworkFeature(Q_minus,Group_Label,9999,n_TP);
figure;
set(gcf,'color','w');
plot(time,P_final_Qplus,'LineWidth',2,'color','r');
hold on;
plot(time,P_final_Qminus,'LineWidth',2,'color','b');
%%% The following features are computed at the nodal level; this means that
%%% we need to upscale our statistical testing to incorporate both the
%%% temporal and the spatial dimensions
% Clustering coefficient
[P_final_CCplus] = ComputePValuesNodalFeatures(CC_plus,Group_Label,999,n_TP,n_ROI);
[P_final_CCminus] = ComputePValuesNodalFeatures(CC_minus,Group_Label,999,n_TP,n_ROI);
[P_final_Eplus] = ComputePValuesNodalFeatures(E_plus,Group_Label,999,n_TP,n_ROI);
[P_final_Eminus] = ComputePValuesNodalFeatures(E_minus,Group_Label,999,n_TP,n_ROI);
[P_final_Btwplus] = ComputePValuesNodalFeatures(Betw_plus,Group_Label,999,n_TP,n_ROI);
[P_final_Btwminus] = ComputePValuesNodalFeatures(Betw_minus,Group_Label,999,n_TP,n_ROI);
[P_final_Degplus] = ComputePValuesNodalFeatures(Deg_plus,Group_Label,999,n_TP,n_ROI);
[P_final_Degminus] = ComputePValuesNodalFeatures(Deg_minus,Group_Label,999,n_TP,n_ROI);
% Makes the results logical vectors
P_final_CCplus(isnan(P_final_CCplus)) = 0;
P_final_CCplus(P_final_CCplus < 0.05) = 0;
P_final_CCplus(P_final_CCplus >= 0.05) = 1;
P_final_CCplus = logical(1-P_final_CCplus);
P_final_CCminus(isnan(P_final_CCminus)) = 0;
P_final_CCminus(P_final_CCminus < 0.05) = 0;
P_final_CCminus(P_final_CCminus >= 0.05) = 1;
P_final_CCminus = logical(1-P_final_CCminus);
P_final_Eplus(isnan(P_final_Eplus)) = 0;
P_final_Eplus(P_final_Eplus < 0.05) = 0;
P_final_Eplus(P_final_Eplus >= 0.05) = 1;
P_final_Eplus = logical(1-P_final_Eplus);
P_final_Eminus(isnan(P_final_Eminus)) = 0;
P_final_Eminus(P_final_Eminus < 0.05) = 0;
P_final_Eminus(P_final_Eminus >= 0.05) = 1;
P_final_Eminus = logical(1-P_final_Eminus);
P_final_Btwplus(isnan(P_final_Btwplus)) = 0;
P_final_Btwplus(P_final_Btwplus < 0.05) = 0;
P_final_Btwplus(P_final_Btwplus >= 0.05) = 1;
P_final_Btwplus = logical(1-P_final_Btwplus);
P_final_Btwminus(isnan(P_final_Btwminus)) = 0;
P_final_Btwminus(P_final_Btwminus < 0.05) = 0;
P_final_Btwminus(P_final_Btwminus >= 0.05) = 1;
P_final_Btwminus = logical(1-P_final_Btwminus);
P_final_Degplus(isnan(P_final_Degplus)) = 0;
P_final_Degplus(P_final_Degplus < 0.05) = 0;
P_final_Degplus(P_final_Degplus >= 0.05) = 1;
P_final_Degplus = logical(1-P_final_Degplus);
P_final_Degminus(isnan(P_final_Degminus)) = 0;
P_final_Degminus(P_final_Degminus < 0.05) = 0;
P_final_Degminus(P_final_Degminus >= 0.05) = 1;
P_final_Degminus = logical(1-P_final_Degminus);
% Plots the results: summarizing and individual metrics
MakeFancyColorDisplay(P_final_CCplus',P_final_Eplus',P_final_Btwplus',...
GetsSigns(CC_plus,Group_Label),GetsSigns(E_plus,Group_Label),...
GetsSigns(Betw_plus,Group_Label),time);
MakeFancyColorDisplay(P_final_CCminus',P_final_Eminus',P_final_Btwminus',...
GetsSigns(CC_minus,Group_Label),GetsSigns(E_minus,Group_Label),...
GetsSigns(Betw_minus,Group_Label),time);
% Our Venn diagram colors
Venn = cmyk2rgb([0,0,0,70;...
0.86,100,94.81,0.2;...
100,83.06,23.79,5.1;...
74.25,8.63,100,0.34;...
3.03,97.83,7.92,0;...
75.77,21.45,0,0;...
4.57,0,88.84,0;...
0,0,0,0;...
0.86,100,94.81,20.2;...
100,83.06,23.79,25.1;...
74.25,8.63,100,20.34;...
3.03,97.83,7.92,20;...
75.77,21.45,0,20;...
4.57,0,88.84,20;...
0,0,0,20]/100);
% Clustering coefficient
figure;
set(gcf,'color','w');
plot(time,sum(P_final_CCplus'.*GetsSigns(CC_plus,Group_Label)),'LineWidth',2,'color',Venn(2,:));
set(gca,'Box','off');
hold on;
plot(time,sum(P_final_Eplus'.*GetsSigns(E_plus,Group_Label)),'LineWidth',2,'color',Venn(3,:));
plot(time,sum(P_final_Btwplus'.*GetsSigns(Betw_plus,Group_Label)),'LineWidth',2,'color',Venn(4,:));
plot(time,sum(P_final_Degplus'.*GetsSigns(Deg_plus,Group_Label)),'LineWidth',2,'color',[0.5 0.5 0.5]);
figure;
set(gcf,'color','w');
plot(time,sum(P_final_CCminus'.*GetsSigns(CC_minus,Group_Label)),'LineWidth',2,'color',Venn(2,:));
set(gca,'Box','off');
hold on;
plot(time,sum(P_final_Eminus'.*GetsSigns(E_minus,Group_Label)),'LineWidth',2,'color',Venn(3,:));
plot(time,sum(P_final_Btwminus'.*GetsSigns(Betw_minus,Group_Label)),'LineWidth',2,'color',Venn(4,:));
plot(time,sum(P_final_Degminus'.*GetsSigns(Deg_minus,Group_Label)),'LineWidth',2,'color',[0.5 0.5 0.5]);
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_CCplus.*GetsSigns(CC_plus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(2,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_CCminus.*GetsSigns(CC_minus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(2,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Eplus.*GetsSigns(E_plus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(3,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Eminus.*GetsSigns(E_minus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(3,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Btwplus.*GetsSigns(Betw_plus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(4,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Btwminus.*GetsSigns(Betw_minus,Group_Label)'),1,'LineWidth',2,'FaceColor',Venn(4,:));
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Degplus.*GetsSigns(Deg_plus,Group_Label)'),1,'LineWidth',2,'FaceColor',[0.5 0.5 0.5]);
set(gca,'Box','off');
figure;
set(gcf,'color','w');
bar(1:1.6:270*1.6,sum(P_final_Degminus.*GetsSigns(Deg_minus,Group_Label)'),1,'LineWidth',2,'FaceColor',[0.5 0.5 0.5]);
set(gca,'Box','off');
%% 5. Computation of ISFC state features of interest and group differences
% We considered three metrics: (1) the idiosyncrasy in state expression,
% (2) the homogeneity in state expression, and (3) the dissimilarity
% between populations
% Selected cluster number (from Section 2)
K = 95;
% Time vector (used for plotting)
time = (19:2:351) - 5;
% Temporal evolution of the clusters for subjects and groups
figure;
set(gcf,'color','w');
title(['K=',num2str(K)]);
tmp_plot = reshape(IDX_THECHOSENONE,n_TP,(n_CTR+n_ASD));
imagesc(1:(n_CTR+n_ASD+5),time,[tmp_plot(:,Group_Label==1),K+1*ones(167,5),tmp_plot(:,Group_Label==2)]);
hold on;
caxis([0.5 K+1.5]);
colormap(colorcube(K+1));
xlabel('Subjects');
ylabel('Movie time [s]');
%%% Homogeneity and idiosyncrasy
% First, I want to compute homogeneity and idiosyncrasy in cluster
% expression, two related variables
[HOM_TD,HOM_ASD,IDIO_TD,IDIO_ASD] = ComputeHomIdio(IDX_THECHOSENONE,...
n_TP,n_CTR,n_ASD,Group_Label);
% Extraction of significant state indices in the TD and ASD groups (the
% ones occurring at least once)
IDX_2D = reshape(IDX_THECHOSENONE,n_TP,(n_CTR+n_ASD));
for t = 1:n_TP
States_TD(t) = length(unique(IDX_2D(t,Group_Label==1)));
States_ASD(t) = length(unique(IDX_2D(t,Group_Label==2)));
end
ns_TD = median(States_TD);
ns_ASD = median(States_ASD);
% Next, I want to generate a "null surface" for the same cases
for null = 1:999
null
% For each fold, I want to generate new index time courses sampled at
% random from all possible group states
NU = nan(n_TP,n_CTR+n_ASD);
for t = 1:n_TP
NU_TD(t,:) = datasample(1:ns_TD,n_CTR);
NU_ASD(t,:) = datasample(1:ns_ASD,n_ASD);
end
NU(:,Group_Label==1) = NU_TD;
NU(:,Group_Label==2) = NU_ASD;
[HOM_TD_Null(:,null),HOM_ASD_Null(:,null),IDIO_TD_Null(:,null),IDIO_ASD_Null(:,null)] = ...
ComputeHomIdio(NU(:),n_TP,n_CTR,n_ASD,Group_Label);
clear NU
clear NU_TD
clear NU_ASD
end
CM = cbrewer('div','RdYlBu',n_TP);
% I now want to generate a summarizing figure with a background of null
% pairs, and a foreground of actual TD/ASD trajectories
figure;
set(gcf,'color','w');
f1 = scatter(HOM_TD_Null(:),IDIO_TD_Null(:),1000,[0.8 0.8 1],'s','filled');
hold on;
scatter(HOM_TD(:),IDIO_TD(:),linspace(300,60,n_TP),flipud(CM),'filled');
xlim([0,0.7]);
ylim([0.1,0.8]);
figure;
set(gcf,'color','w');
f2 = scatter(HOM_ASD_Null(:),IDIO_ASD_Null(:),1000,[1 0.8 0],'s','filled');
hold on;
scatter(HOM_ASD(:),IDIO_ASD(:),linspace(300,60,n_TP),flipud(CM),'filled');
xlim([0,0.7]);
ylim([0.1,0.8]);
% When in time do I get out of my null zone
time_sign = OutOfTheZone(HOM_TD_Null,IDIO_TD_Null,HOM_TD,IDIO_TD);
time_sign_b = OutOfTheZone(HOM_ASD_Null,IDIO_ASD_Null,HOM_ASD,IDIO_ASD);
%%% Dissimilarity in state expression across groups
[P_DISS] = ComputePValuesDissimilarity(IDX_2D,Group_Label,9999,n_TP,n_CTR,n_ASD);
% figure;
% set(gcf,'color','w');
% time = 19:2:351;
% plot(time,DG,'LineWidth',2,'color',[0 0.6 0.2]);
% hold on;
% plot(time,max(DG_null,[],2),'LineWidth',1,'color',[0.6 1 0.2]);
% plot(time,min(DG_null,[],2),'LineWidth',1,'color',[0.6 1 0.2]);
% set(gca,'Box','off');
% xlim([time(1),time(end)]);
%% Selected time points
% Time points of interest
t_OI = [54,60,68,166,204,214,218,256];
scale_overall = 1;
i = 1;
for t = 1:length(t_OI)
test_idx = find(P_final_Degplus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_DEG'],'-depsc');
test_idx = find(P_final_CCplus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_CC'],'-depsc');
test_idx = find(P_final_Eplus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_E'],'-depsc');
test_idx = find(P_final_Btwplus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_B'],'-depsc');
end
t_OI = [20,60,134,212,282];
scale_overall = 1;
i = 1;
for t = 1:length(t_OI)
test_idx = find(P_final_Degminus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_DEG_MINUS'],'-depsc');
test_idx = find(P_final_CCminus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_CC_MINUS'],'-depsc');
test_idx = find(P_final_Eminus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_E_MINUS'],'-depsc');
test_idx = find(P_final_Btwminus(time==t_OI(t),:));
toplot = zeros(270,1);
toplot(test_idx)=1;
figure;
PlotBrainGraph_Kmeans(zeros(270,270),toplot,...
CodeBook,0,scale_overall/2,scale_overall,1,i,'autumn','autumn',1,1,[]);
print(['NETWORK_TP_',num2str(t_OI(t)),'_B_MINUS'],'-depsc');
end
%% 6. PLS analysis within ASD subjects: network-level features
% Computation of behavioral features of interest: we consider three scores,
% one telling of social impairments, one of stereotypies, one of nonverbal
% communication impairment. We remove the subjects for whom all three
% scores were not available
Behavioral = AQ(Group_Label==2,[1,2,4]);
idx_toremove = isnan(sum(Behavioral,2));
Behavioral(idx_toremove,:) = [];
% Now, we create our network-level features
CCp = CreateNetworkPLSFeature(CC_plus,Group_Label);
CCm = CreateNetworkPLSFeature(CC_minus,Group_Label);
Degp = CreateNetworkPLSFeature(Deg_plus,Group_Label);
Degm = CreateNetworkPLSFeature(Deg_minus,Group_Label);
Ep = CreateNetworkPLSFeature(E_plus,Group_Label);
Em = CreateNetworkPLSFeature(E_minus,Group_Label);
Bp = CreateNetworkPLSFeature(Betw_plus,Group_Label);
Bm = CreateNetworkPLSFeature(Betw_minus,Group_Label);
%EIGp = CreateNetworkPLSFeature(EIV_plus,Group_Label);
%EIGm = CreateNetworkPLSFeature(EIV_minus,Group_Label);
% Concatenating them together
Functional_network = [CCp,CCm,Ep,Em,Bp,Bm,Degp,Degm];
%Functional_network = [CreateNetworkPLSFeature(Q_plus,Group_Label),CreateNetworkPLSFeature(Q_minus,Group_Label)];
Functional_network(idx_toremove,:) = [];
% From preliminary analyses, we know that the first two values are outliers
Behavioral([1,2],:) = [];
Functional_network([1,2],:) = [];
% Performing PLS analysis itself
[UN,VN,EVN,sprobN,BS_UN,BS_VN] = RunPLS(Functional_network,Behavioral);
% First, the behavioral saliences
figure;
set(gcf,'color','w');
bar(BS_VN(:,1));
set(gca,'Box','off');
RecoverNetworkSaliences(BS_UN(:,1),3,n_TP,n_ROI,time,CodeBook);
% Latent weights
Lf2 = Functional_network * UN(:,1);
Lb2 = Behavioral * VN(:,1);
figure;
set(gcf,'color','w');
scatter(Lf2,Lb2,30,'k','d','fill');
%% 7. PLS analysis within ASD subjects: state features
% Computation of functional features of interest (state analysis)
ISFC_OI = ISFC_MW(Group_Label==2);
STATES_OI = IDX_2D(:,Group_Label==2);
% Similarity to cluster and cluster counts for each subject
for s = 1:n_ASD
for st = 1:K
Counts(s,st) = sum(STATES_OI(:,s)==st);
end
end
% Functional_state contains our metrics obtained from state analysis: for
% each of the K ISFC states, the counts and the similarity to centroid
% values. Thus, a total of 2*K features.
Functional_state = [Counts];
Functional_state(idx_toremove,:) = [];
Functional_state([1,2],:) = [];
[US,VS,EVS,sprobS,BS_US,BS_VS] = RunPLS(Functional_state,Behavioral);
BS_US(isnan(BS_US)) = 0;
% First, the behavioral saliences
figure;
set(gcf,'color','w');
bar(BS_VS(:,1));
set(gca,'Box','off');
% Then, the state saliences (counts)
figure;
set(gcf,'color','w');
bar(BS_US(:,1));
set(gca,'Box','off');
% Latent weights
Lf = Functional_state * US(:,1);
Lb = Behavioral * VS(:,1);
figure;
set(gcf,'color','w');
scatter(Lf,Lb,30,'k','d','fill');
% Building a matrix summarizing state expression across the ASD population
for st = 1:K
% tmp has size n_TP x n_ASD
tmp = IDX_2D(:,Group_Label == 2);
tmp(:,idx_toremove) = [];
tmp(:,[1,2]) = [];
tmp(tmp ~= st) = 0;
tmp(tmp == st) = 1;
State_TC(:,st) = sum(tmp,2)/size(tmp,2);
end
% Now, I want to build a time course representative of the ISFC states
% contributing to the gradient of interest. The goal is to have a large
% value when a state that is strongly present in the saliences is expressed
% across subjects
Movie_TimeCourse = zeros(167,1);
% For each state, I weight its overall expression across subjects
% (State_TC) by its salience in the PLS analysis (BS_US)
for i = 1:K
Movie_TimeCourse = Movie_TimeCourse + BS_US(i,1) * State_TC(:,i);
end
figure;
set(gcf,'color','w');
plot(time,(Movie_TimeCourse),'LineWidth',2,'color','k');
hold on;
set(gca,'Box','off');
% I now want to generate a "brain ISFC contrast"
Contrast = zeros(n_conn,1);
for st = 1:K
Contrast(IDX) = Contrast(IDX) + BS_US(st,1) * CENTROIDS_THECHOSENONE(st,:)';
end
figure;
set(gcf,'color','w');
scatter(D_euclidean(IDX),Contrast(IDX),10,'k','o','fill');
set(gca,'Box','off');
% This is just doing the plotting: first, I take my contrast connectivities
% of interest
tmp_plot = Contrast(IDX);
% Then, I define the threshold below which a contrast value is too low to
% be displayed (i.e., the connection will not be shown)
thresh = prctile(abs(tmp_plot),90);
% tmp_plot2 has the same information, in the whole connection space
tmp_plot2 = zeros(n_conn,1);
tmp_plot2(IDX) = tmp_plot;
% I want to compute tmp_plot3 when tmp_plot2 still has all connections
% (positive and negative), because it will be used to plot the nodes in the
% display
tmp_plot3 = jVecToUpperTriMat(tmp_plot2,n_ROI);
% tmp_plot2 and tmp_plot2b contain the positive and negative-valued
% elements, respectively
tmp_plot2b = tmp_plot2;
tmp_plot2(tmp_plot2 < 0) = 0;
tmp_plot2b(tmp_plot2b > 0) = 0;
tmp_plot2b = abs(tmp_plot2b);
% Overall scaling factor used for the plotting
scale_overall = max(sum(abs(tmp_plot3)));
% What is below my threshold is not displayed
tmp_plot2(tmp_plot2 < thresh) = 0;
tmp_plot4 = jVecToUpperTriMat(tmp_plot2,n_ROI);
tmp_plot2b(tmp_plot2b < thresh) = 0;
tmp_plot4b = jVecToUpperTriMat(tmp_plot2b,n_ROI);
% Plotting itself; the only way I found to have both positive and negative
% connections plotted together is to overlay two displays
figure;
PlotBrainGraph_Kmeans(tmp_plot4,sum(abs(tmp_plot3),2),...
CodeBook,0,scale_overall/2,scale_overall,1,2,'hot','autumn',1,1,[]);
print(['PLS_CONTRAST_SHINEY_SAGPOS'],'-depsc');
figure;
PlotBrainGraph_Kmeans(tmp_plot4b,sum(abs(tmp_plot3),2),...
CodeBook,0,scale_overall/2,scale_overall,1,2,'hot','winter',1,1,[]);
print(['PLS_CONTRAST_SHINEY_SAGNEG'],'-depsc');
figure;
PlotBrainGraph_Kmeans(tmp_plot4,sum(abs(tmp_plot3),2),...
CodeBook,0,scale_overall/2,scale_overall,1,1,'hot','autumn',1,1,[]);
print(['PLS_CONTRAST_SHINEY_HORPOS'],'-depsc');
figure;
PlotBrainGraph_Kmeans(tmp_plot4b,sum(abs(tmp_plot3),2),...
CodeBook,0,scale_overall/2,scale_overall,1,1,'hot','winter',1,1,[]);
print(['PLS_CONTRAST_SHINEY_HORNEG'],'-depsc');
c4science · Help