DMN_comp.m

clear all
close all
clc

%% adding path

path = genpath('/Users/iricchi/Documents/MIP_LabImmersion');
addpath(path)
clc

%% load and run
p = '/Users/iricchi/Documents/MIP_LabImmersion/d/';

to_load = {'WM_opt_Ls', 'S_opt_Ls',  'RS_opt_Ls', ...
           'RM_opt_Ls', 'M_opt_Ls',  'L_opt_Ls', 'E_opt_Ls'};
titles = {'Working Memory',  'Social',  'Resting State', ...
           'Relational Memory', 'Motor', 'Language', 'Emotion'};
       
to_loadmu = {'WM_opt_mus', 'S_opt_mus',  'RS_opt_mus', ...
           'RM_opt_mus', 'M_opt_mus',  'L_opt_mus', 'E_opt_mus'};

ind_dmn = [4, 1, 2, 1, 4, 1, 1];
res_spec_dmn = zeros(length(to_load));
res_cos_dmn = zeros(length(to_load));
cor_tot = zeros(length(to_load));
p_values = zeros(length(to_load));
for i= 1:length(to_load)
    for j = 1:length(to_load)
        ftl1 = [p,to_load{i}];
        Ls1 = load(ftl1).Ls;
        ftl2 = [p,to_load{j}];
        Ls2 = load(ftl2).Ls;
        k1 = ind_dmn(i);
        k2 = ind_dmn(j);
        L1 = Ls1(:,:,k1);
        L2 = Ls2(:,:,k2);
        W1 = diag(diag(L1)) - L1;
        W1(abs(W1)<1e-6) = 0;
        W2 = diag(diag(L2)) - L2;
        W2(abs(W2)<1e-6) = 0; 
        % normalize every matrix
        w1 = normalize(W1, 'range');
        w2 = normalize(W2, 'range');
        
        tr_1 = triu(w1,1);
        tr_2 = triu(w2,1);
        
        tr1 = tr_1.';
        m = tril(true(size(tr_1)),1);
        v1  = tr1(m).';
        tr2 = tr_2.';
        m = tril(true(size(tr_2)),1);
        v2  = tr2(m).';
        
        
        % cor = corrcoef(v1,v2);
        [cor, pv] = corr(v1',v2');
        cor_tot(i,j) = cor;
        
        ftl1 = [p,to_loadmu{i}];
        m1 = load(ftl1).mus;
        ftl2 = [p,to_loadmu{j}];
        m2 = load(ftl2).mus;
        
        corr_means(i,j) = corr(m1(:,k1),m2(:,k2));
%         % t-test to check for significance
%         [h,pv] = ttest(v1,v2);
%         p_values(i,j) = pv;
        
        [eigv1, d1] =  eig(w1);
        d1 = diag(d1);
        
        [eigv2, d2] =  eig(w2);
        d2 = diag(d2);
        rotatemat = rri_bootprocrust(eigv1, eigv2);
        
        Z = eigv2.* d2.*rotatemat; % Z is your transformed eigenvector 
        rot_eigvals = sqrt(sum(Z.^2))'; % rotate eigenvalues             
        %norm_cos = normalize(cos,'range');
        res_cos_dmn(i,j) = getCosineSimilarity(d1,rot_eigvals);
        res_spec_dmn(i,j) = euclidean_dist(d1,rot_eigvals);
        
        
    end
    
%     figure()
%     imagesc(w1)
%     title(titles{i})
end    


%% spec dist between each DMN in every dataset
% generate figure 5 D

to_load = {'WM_opt_Ls', 'S_opt_Ls',  'RS_opt_Ls', ...
           'RM_opt_Ls', 'M_opt_Ls',  'L_opt_Ls', 'E_opt_Ls'};
titles = {'Working Memory',  'Social',  'Resting State', ...
           'Relational Memory', 'Motor', 'Language', 'Emotion'};
labels = {'WM', 'S', 'RS','RM','M','L','E'};


% structural connectome
SC = load('GLMM_subjects_SC_AAL90.mat').SC;
SCave = mean(SC,3);
sc = SCave+SCave';
mask_diag = logical(eye(size(sc)));
sc(mask_diag)=0;
sc = normalize(sc, 'range');


[eigvs, evals] = eig(sc);
evals = diag(evals); 

ind_dmn = [4, 1, 2, 1, 4, 1, 1];
res_sdiff = zeros(length(ind_dmn),1);
vars = zeros(length(ind_dmn),1);
cors = zeros(length(ind_dmn),1);
for i= 1:length(to_load)
    ftl = [p,to_load{i}];
    Ls = load(ftl).Ls;
    
    L = Ls(:,:,ind_dmn(i));
    W = diag(diag(L)) - L;
    W(abs(W)<1e-6) = 0;  % 25th percentile
    % normalize every matrix
    w = normalize(W, 'range');
    r = corrcoef((triu(w,1)),(triu(sc,1)));
    disp(['k = ', num2str(i), ', corr = ', num2str(r(2))])
    cors(i) = r(2);
    [eigvf, evalf] = eig(w);
    evalf = diag(evalf);
    rotatemat = rri_bootprocrust(eigvs, eigvf); 

    Z = eigvf .* evalf .* rotatemat; % Z transformed eigenvalues
    rot_eigvals = sqrt(sum(Z.^2))'; % rotate eigenvalues     

    sp_diff = (rot_eigvals - evals).^2;
    res_sdiff(i) = sum(sp_diff);
    % vars(i) = var(sp_diff);
    
end
%%
spdiff = normalize(res_sdiff,'range');
%spdiff = res_sdiff;
[sv, si] = sort(spdiff);
% stds = sqrt(vars);
% errlow = sv-stds;
% errhigh = sv+stds;
% er = errorbar(1:length(sv),sv,errlow,errhigh);    
figure()
b = bar(sv, 0.7,'red','LineWidth',2);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom', 'FontSize', 14)
ylabel("Norm spect diff %")
%ylim([0,112])
xticklabels(labels(si))

% hold on
% 
% er.Color = [0 0 0];                            
% er.LineStyle = 'none'; 

%% 

[sv, si] = sort(res_sdiff);
stds = sqrt(vars);
figure()
b = bar(sv, 0.7,'red','LineWidth',2);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom', 'FontSize', 14)
xticklabels(labels(si))
ylabel("Norm spect diff")

%%
[sv, ssi] = sort(cors,'descend');
figure()
b = bar(1-sv, 0.7,'red','LineWidth',1.5);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom')
ylabel("Norm spect diff")
ylim([0,0.6])
xticklabels(labels(si))

%% Means correlations

% trying to replicate Anjali's matrix (fig 4c)
% this is correlating everything with everything
p = '/Users/iricchi/Documents/MIP_LabImmersion/data_plot_nets/';

to_load = {'opt_k_9_task_W_means', 'opt_k_8_task_S_means',  'RS_opt_mus', ...
           'opt_k_8_task_R_means', 'opt_k_4_task_M_means',  'opt_k_7_task_L_means',...
           'opt_k_5_task_E_means'};
              
titles = {'Working Memory',  'Social',  'Resting State', ...
           'Relational Memory', 'Motor', 'Language', 'Emotion'};

labels = {'WM', 'S', 'RS','RM','M','L','E'};

%ind_dmn = [3,2,2,1,4,7,2];

mu_all = zeros(90,7);
p_vals = zeros(7,1);
for i= 1:length(to_load)
    ftl = [p,to_load{i}];
    mus = load(ftl).mus;
    
    mu = mus(:,ind_dmn(i));
    mu_all(:,i) = mu;
    
end

cos_sim = zeros(7,7);
for i=1:length(to_load)
    for j=1:length(to_load)
        cos_sim(i,j) = getCosineSimilarity(mu_all(:,i),mu_all(:,j)); 
    end
end

[sorted, idx] = sort(sum(cos_sim,1));
figure()
imagesc(cos_sim(idx,idx))
colorbar()
colormap("jet")

xticklabels(labels(idx))
yticklabels(labels(idx))


[corr_means, pv] = corr(mu_all);
[sorted, idx] = sort(sum(corr_means,1));

figure()
imagesc(corr_means(idx,idx))
colorbar()
colormap("jet")

xticklabels(labels(idx))
yticklabels(labels(idx))

% % apply t-test
% [h, p_vals(1)] = ttest(mu_all(:,1), mu_all(:,2));
% [h, p_vals(2)] = ttest(mu_all(:,2), mu_all(:,3));
% [h, p_vals(3)] = ttest(mu_all(:,3), mu_all(:,4));
% [h, p_vals(4)] = ttest(mu_all(:,4), mu_all(:,5));
% [h, p_vals(5)] = ttest(mu_all(:,5), mu_all(:,6));
% [h, p_vals(6)] = ttest(mu_all(:,6), mu_all(:,7));
% 
% % apply anova
% [p,tbl,stats] = anova1(mu_all);
% p = kruskalwallis(mu_all);




%% spec dist between each DMN in every dataset
% generate figure 5 D
p = '/Users/iricchi/Documents/MIP_LabImmersion/data_plot_nets/';

to_load = {'opt_k_9_task_W_Ls', 'opt_k_8_task_S_Ls',  'RS_opt_Ls', ...
           'opt_k_8_task_R_Ls', 'opt_k_4_task_M_Ls',  'opt_k_7_task_L_Ls', 'opt_k_5_task_E_Ls'};
       
to_load = {'W_Ls_k_4', 'S_Ls_k_3',  'RS_opt_Ls', ...
           'R_Ls_k_4', 'M_Ls_k_5',  'L_Ls_k_3', 'E_Ls_k_3'};
       
titles = {'Working Memory',  'Social',  'Resting State', ...
           'Relational Memory', 'Motor', 'Language', 'Emotion'};
labels = {'WM', 'S', 'RS','RM','M','L','E'};


% sc2 (new)
SC = load('GLMM_subjects_SC_AAL90.mat').SC;
SCave = mean(SC,3);
sc = SCave+SCave';
mask_diag = logical(eye(size(sc)));
sc(mask_diag)=0;
sc = normalize(sc, 'range');

% % sc old
% old = load('Mean_AAL116_SC_105_subjects'); % load the ALL SC
% SC = old.SC_AAL(1:90,1:90);
% %SC(abs(SC)<th) = 0;
% %sc1 = SC;
% sc = normalize(SC, 'range');
% 

[eigvs, evals] = eig(sc);
evals = diag(evals); 

ind_dmn = [2, 3, 2, 1, 5, 1, 1];
res_sdiff = zeros(length(ind_dmn),1);
vars = zeros(length(ind_dmn),1);
cors = zeros(length(ind_dmn),1);
for i= 1:length(to_load)
    ftl = [p,to_load{i}];
    Ls = load(ftl).Ls;
    
    L = Ls(:,:,ind_dmn(i));
    W = diag(diag(L)) - L;
    W(abs(W)<1e-6) = 0;  % 25th percentile
    % normalize every matrix
    w = normalize(W, 'range');
    r = corrcoef((triu(w,1)),(triu(sc,1)));
    disp(['k = ', num2str(i), ', corr = ', num2str(r(2))])
    cors(i) = r(2);
    [eigvf, evalf] = eig(w);
    evalf = diag(evalf);
    rotatemat = rri_bootprocrust(eigvs, eigvf); 

    Z = eigvf .* evalf .* rotatemat; % Z transformed eigenvalues
    rot_eigvals = sqrt(sum(Z.^2))'; % rotate eigenvalues     

    sp_diff = (rot_eigvals - evals);
    res_sdiff(i) = max(sp_diff);
    
end
%
spdiff = normalize(res_sdiff,'range');
[sv, si] = sort(spdiff);
% stds = sqrt(vars);
% errlow = sv-stds;
% errhigh = sv+stds;
% er = errorbar(1:length(sv),sv,errlow,errhigh);    
figure()
b = bar(100*sv, 0.7,'red','LineWidth',2);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom', 'FontSize', 14)
ylabel("Norm spect diff %")
ylim([0,112])
xticklabels(labels(si))

% hold on
% 
% er.Color = [0 0 0];                            
% er.LineStyle = 'none'; 

%

[sv, si] = sort(res_sdiff);
stds = sqrt(vars);
figure()
b = bar(100*sv, 0.7,'red','LineWidth',2);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom', 'FontSize', 14)
xticklabels(labels(si))
ylabel("Norm spect diff")


[sv, si] = sort(normalize(res_sdiff, 'norm'));
figure()
b = bar(sv, 0.7,'red','LineWidth',2);
xtips = b.XEndPoints;
ytips = b.YEndPoints;
label = string(b.YData);
text(xtips,ytips,label,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom', 'FontSize', 14)
xticklabels(labels(si))
ylabel("Norm spect diff")
ylim([0 0.5])


%% Neurosynth 
%%%%%%% CHANGE DIRECTORY IN NEUROSYNTH

mus_dmn = [];
for i = 1:length(ind_dmn)
    load([to_load{i}(1:end-2), 'mus'])
    mus_oi = mus(:,ind_dmn(i));
    mus_dmn = [mus_dmn, mus_oi];
end

mus_dmn = mus_dmn(:,si);
mus_dmn = zscore(mus_dmn);


for i = 1:size(mus_dmn,2)
   mu = mus_dmn(:,i);
   mu(mu<1) = 0;
   mu(mu>=1) = 1;
   name = ['dmn_mu', num2str(i), '.nii'];
   GraphToNifti(mu, name, 'AAL116_correctLR_resliced2fMRIHCP.nii')
end