Hettwer2022_Figure3_Within_Disorder_Covariance.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Coordinated Cortical Thickness alterations Across Mental Disorders: A Transdiagnostic ENIGMA Project
% Script 3: Computes disorder-specific covariance profiles and embeds
% individual disorders within a transdiagnostic space.
%This script runs analyses presented in Figure 3 of the manuscript.
%Meike Hettwer 2022 - CNG Lab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Add paths and Toolboxes
%ENIGMA-1.1.3huch
%cbrewer
%data to load from GitHub:
%Epicenters.mat
%CT_psych_gradients.mat
%corr_dis.mat
%% Load ENIGMA summary statistics
%Schizophrenia
sum_stats = load_summary_stats('schizophrenia');
CT = sum_stats.CortThick_case_vs_controls;
CT_d_sz = CT.d_icv;
%Bipolar
sum_stats = load_summary_stats('bipolar');
CT = sum_stats.CortThick_case_vs_controls_adult;
CT_d_bipo = CT.d_icv;
%ASD
sum_stats = load_summary_stats('asd');
CT = sum_stats.CortThick_case_vs_controls_meta_analysis;
CT_d_asd = CT.d_icv;
%ADHD
sum_stats = load_summary_stats('adhd');
CT = sum_stats.CortThick_case_vs_controls_adult;
CT_d_adhd = CT.d_icv;
%MDD
sum_stats = load_summary_stats('depression');
CT = sum_stats.CortThick_case_vs_controls_adult;
CT_d_mdd = CT.d_icv;
%OCD
sum_stats = load_summary_stats('ocd');
CT = sum_stats.CortThick_case_vs_controls_adult;
CT_d_ocd = CT.d_icv;
%combine
disorders = [CT_d_bipo, CT_d_adhd, CT_d_asd, CT_d_mdd, CT_d_ocd, CT_d_sz];
disorder_names = {'BD'; 'ADHD'; 'ASD'; 'MDD'; 'OCD'; 'SCZ'};
%% Disorder-specific covariance (*Figure 3B*)
% via differences in Cohen's d values between regions
for i = 1:6
dis = disorders(:,i);
dis_diff=nan(length(dis),length(dis));
for j = 1:length(dis)%columns
for k = 1:length(dis)%rows
dis_diff(k,j) = dis(j)-dis(k); %caution: resulting matrix is asymmetrical as signs flip between rows and columns
end
end
dis_diff_all(:,:,i) = abs(dis_diff(:,:));%take abs so it's symmetric
clear dis_diff dis dis_diff_threshbin
end
%within-disorder covariance plots
color = cbrewer('seq','BuPu',50);
color(color<0) =0;
clims = [-0.2 0];
for i = 1:6
this_dis = dis_diff_all(:,:,i)*-1;
fig=figure; imagesc(this_dis, clims); title(disorder_names{i}); colormap(color);colorbar
yticks([]), xticks([]); axis square
%close all
end
%correlation between disorder-specific and transdiagnostic covariance plotted on the brain
for i = 1:6
this_dis = dis_diff_all(:,:,i);
for column = 1:68
r_temp(column,i) = corr(this_dis(:,column)*-1,corr_dis(:,column));
end
fig=figure; plot_cortical(parcel_to_surface(r_temp(:,i)),'label_text',disorder_names{i},'color_range', [-0.6 0.6])
end
%% Cross-disorder similarity matrix (*Figure 3C*)
cross_disorder_matrix = corr(disorders);
figure; imagesc(cross_disorder_matrix); set(gca, 'XTick',[1:6],'XTickLabel',disorder_names','YTick',[1:6],'YTickLabel',disorder_names,'clim',[-1 1]);colormap(flipud(cbrewer('div','RdBu',25)));colorbar
%% cluster disorders (*Figure 3D*)
tree = linkage(corr(disorders)' ,'average');
figure()
dendrogram(tree,'Labels',disorder_names);
%% Embedding of individual disorders in 2D space framed by transdiagnostic epicenters and hubs
%A) epicenters
%load HCP connectivity data
[fc_ctx, fc_ctx_labels, ~, ~] = load_fc();
[sc_ctx, sc_ctx_labels, ~, ~] = load_sc();
[~, ~, fc_sctx, fc_sctx_labels] = load_fc();
[~, ~, sc_sctx, sc_sctx_labels] = load_sc();
%load Epicenters.mat from Github !!
%cortical
epi_fc_trans = fc_ctx_trans_p_sig;
epi_fc_trans_bin = find(fc_ctx_trans_p_sig>0);
epi_sc_trans = sc_ctx_trans_p_sig;
epi_sc_trans_bin = find(sc_ctx_trans_p_sig>0);
%subcortical
epi_fc_sctx_trans = fc_sctx_epi_p_sig;
epi_fc_sctx_trans_bin = find(fc_sctx_epi_p_sig>0);
epi_sc_sctx_trans = sc_sctx_epi_p_sig;
epi_sc_sctx_trans_bin = find(sc_sctx_epi_p_sig>0);
epicenter_overlap_fc_sctx = nan(1,6);
epicenter_overlap_sc_sctx = nan(1,6);
epicenter_overlap_fc_percent = nan(1,6);
epicenter_overlap_sc_percent = nan(1,6);
epicenter_overlap_fc_sctx = nan(1,6);
epicenter_overlap_sc_sctx = nan(1,6);
for i = 1:size(disorders,2)
% Thickness maps for each disorder
fc_ctx_epi = zeros(size(fc_ctx, 1), 1);
fc_ctx_epi_p = zeros(size(fc_ctx, 1), 1);
for seed = 1:size(fc_ctx, 1)
seed_conn = fc_ctx(:, seed);
r_tmp = corrcoef(seed_conn, abs(disorders(:,i)));
fc_ctx_epi(seed) = r_tmp(1, 2);
fc_ctx_epi_p(seed) = spin_test(seed_conn, abs(disorders(:,i)), 'surface_name', 'fsa5', 'parcellation_name', ...
'aparc', 'n_rot', 1000, 'type', 'pearson');
end
% Identify cortical epicenter values (from structural connectivity)
sc_ctx_epi = zeros(size(sc_ctx, 1), 1);
sc_ctx_epi_p = zeros(size(sc_ctx, 1), 1);
for seed = 1:size(sc_ctx, 1)
seed_conn = sc_ctx(:, seed);
r_tmp = corrcoef(seed_conn, abs(disorders(:,i)));
sc_ctx_epi(seed) = r_tmp(1, 2);
sc_ctx_epi_p(seed) = spin_test(seed_conn, abs(disorders(:,i)), 'surface_name', 'fsa5', 'parcellation_name', ...
'aparc', 'n_rot', 1000, 'type', 'pearson');
end
% Select only regions with p < 0.05 (functional epicenters)
fc_ctx_epi_p_sig = zeros(length(fc_ctx_epi_p), 1);
fc_ctx_epi_p_sig(fc_ctx_epi_p < 0.05) = fc_ctx_epi(fc_ctx_epi_p<0.05);
% Selecting only regions with p < 0.05 (structural epicenters)
sc_ctx_epi_p_sig = zeros(length(sc_ctx_epi_p), 1);
sc_ctx_epi_p_sig(sc_ctx_epi_p < 0.05) = sc_ctx_epi(sc_ctx_epi_p<0.05);
%$$$$$$$$$$
%subcortical
% Epicenter mapping - subcortical seeds but cortical map!
% Identify subcortical epicenters (from functional connectivity)
fc_sctx_epi = zeros(size(fc_sctx, 1), 1);
fc_sctx_epi_p = zeros(size(fc_sctx, 1), 1);
for seed = 1:size(fc_sctx, 1)
seed_conn = fc_sctx(seed, :);
r_tmp = corrcoef(seed_conn, abs(disorders(:,i)));
fc_sctx_epi(seed) = r_tmp(1, 2);
fc_sctx_epi_p(seed) = spin_test(seed_conn, abs(disorders(:,i)), 'surface_name', 'fsa5', 'parcellation_name', ...
'aparc', 'n_rot', 1000, 'type', 'pearson');
end
% Identify subcortical epicenters (from structural connectivity)
sc_sctx_epi = zeros(size(sc_sctx, 1), 1);
sc_sctx_epi_p = zeros(size(sc_sctx, 1), 1);
for seed = 1:size(sc_sctx, 1)
seed_conn = sc_sctx(seed, :);
r_tmp = corrcoef(seed_conn, abs(disorders(:,i)));
sc_sctx_epi(seed) = r_tmp(1, 2);
sc_sctx_epi_p(seed) = spin_test(seed_conn, abs(disorders(:,i)), 'surface_name', 'fsa5', 'parcellation_name', ...
'aparc', 'n_rot', 1000, 'type', 'pearson');
end
% Threshold and Project the results on the surface brain
% functional epicenters
fc_sctx_epi_p_sig = zeros(length(fc_sctx_epi_p), 1);
fc_sctx_epi_p_sig(fc_sctx_epi_p < 0.05) = fc_sctx_epi(fc_sctx_epi_p<0.05);
% structural epicenters
sc_sctx_epi_p_sig = zeros(length(sc_sctx_epi_p), 1);
sc_sctx_epi_p_sig(sc_sctx_epi_p < 0.05) = sc_sctx_epi(sc_sctx_epi_p<0.05);
%$$$$$$$$$$
% check overlap with transdiagnostic epicenters
%cortical
fc_ctx_epi_p_sig(fc_ctx_epi_p_sig<0)=0;
sc_ctx_epi_p_sig(sc_ctx_epi_p_sig<0)=0;
functional_epicenters(:,i)= fc_ctx_epi_p_sig;
structural_epicenters(:,i)= sc_ctx_epi_p_sig;
%sub-cortical
fc_sctx_epi_p_sig(fc_sctx_epi_p_sig<0)=0;
sc_sctx_epi_p_sig(sc_sctx_epi_p_sig<0)=0;
sub_functional_epicenters(:,i)= fc_sctx_epi_p_sig;
sub_structural_epicenters(:,i)= sc_sctx_epi_p_sig;
%check disorder_specific overlap with transdiagnostic epicenters
%cortical
position_fc = find(fc_ctx_epi_p_sig>0);
position_sc = find(sc_ctx_epi_p_sig>0);
[val_fc,pos_fc]=intersect(epi_fc_trans_bin,position_fc);
[val_sc,pos_sc]=intersect(epi_sc_trans_bin,position_sc);
%sub-cortical
position_fc_sctx = find(fc_sctx_epi_p_sig>0);
position_sc_sctx = find(sc_sctx_epi_p_sig>0);
[val_fc,pos_fc]=intersect(epi_fc_sctx_trans_bin,position_fc_sctx);
[val_sc,pos_sc]=intersect(epi_sc_sctx_trans_bin,position_sc_sctx);
epicenter_overlap_fc(i) = size(intersect(epi_fc_trans_bin,position_fc),1);
epicenter_overlap_sc(i) = size(intersect(epi_sc_trans_bin,position_sc),1);
epicenter_overlap_fc_sctx(i) = size(intersect(epi_fc_sctx_trans_bin,position_fc_sctx),1);
epicenter_overlap_sc_sctx(i) = size(intersect(epi_sc_sctx_trans_bin,position_sc_sctx),1);
epicenter_overlap_fc_percent(i) = ((epicenter_overlap_fc(i)+epicenter_overlap_fc_sctx(i))/(size(epi_fc_trans_bin,1)+size(epi_fc_sctx_trans_bin,1)))*100;
epicenter_overlap_sc_percent(i) = ((epicenter_overlap_sc(i)+epicenter_overlap_sc_sctx(i))/(size(epi_sc_trans_bin,1)+size(epi_sc_sctx_trans_bin,1)))*100;
end
T_epi_overlap = table(epicenter_overlap_fc_percent', epicenter_overlap_sc_percent','RowNames',disorder_names,'VariableNames',{'Functional','Structural'});
% plot disorder_specific epicenters
%functional
for i = 1:6
fig = figure; plot_subcortical(sub_functional_epicenters(:,i),'color_range',[0 0.5],'label_text',disorder_names{i}, 'ventricles', 'False');
fig = figure; plot_cortical(parcel_to_surface(functional_epicenters(:,i)),'color_range',[0 0.5],'label_text', disorder_names{i})
end
%structural
for i = 1:6
figure; plot_subcortical(sub_structural_epicenters(:,i),'color_range',[0 0.5],'label_text',disorder_names{i}, 'ventricles', 'False');
figure; plot_cortical(parcel_to_surface(structural_epicenters(:,i)),'color_range',[0 0.5],'label_text', disorder_names{i})
end
%% Hubs vs disorder-specific maps
% Set up thresholded disorder correlation matrix and dc map
corr_dis_bin_thresh = bsxfun(@gt, corr_dis, prctile(corr_dis,80))';
% plot degree centrality map / cross-disorder covariance hub map
dis_dc = sum(corr_dis_bin_thresh);
r=nan(1,6);
p=nan(1,6);
for i = 1:6
r(i) = corr(dis_dc', abs(disorders(:,i))); %CHANGED TO ABSOLUTE IMPACT!
p(i) = spin_test(abs(disorders(:,i)), dis_dc', 'surface_name', 'fsa5', 'parcellation_name', ...
'aparc', 'n_rot', 1000, 'type', 'pearson');
end
T_transhubs = table(r',p','RowNames',disorder_names,'VariableNames',{'r','p_spin'});
%% Correlate disorder maps with Gradients, 2D plot
Grad1_dis_cor=nan(1,6);
Grad2_dis_cor=nan(1,6);
for i = 1:size(disorders,2)
Grad1_dis_cor(i) = corr(disorders(:,i),CT_psych_gradient1);
Grad2_dis_cor(i) = corr(disorders(:,i),CT_psych_gradient2);
end
%% Plot how disorder effects are driven by transdiagnostic hubs vs epicenters (*Figure 3E*)
%col = cbrewer('seq','BuPu','50')
txty = {'BD', 'ADHD', 'ASD','MDD', 'OCD', 'SCZ'};
dx = 0.07;
dy = 0.001;
figure;
subplot(1,2,1)
scatter(T_transhubs.r,T_epi_overlap.Functional,100, [0.5059 0.0588 0.4863],'filled');
set(gca,'xlim', [-1 1],'ylim', [0 100], 'XTick', [-1, -0.5 0 0.5 1],'YTick',[0 50 100],'TickDir','none','color','none','FontSize',10, 'FontName', 'Gill Sans MT');xlabel('Transdiagnostic hubs ({\itr})');ylabel('Epicenter overlap (%)'); axis square
hold on%plotting absolute correlations as direction of gradients is meaningless
% displacement so the text does not overlay the data points
%text(T_transhubs.r+dx, T_epi_overlap.Functional+dy, txty,'FontSize',10, 'FontName', 'Gill Sans MT');
hold on
scatter(T_transhubs.r,T_epi_overlap.Structural,100, [0.5451 0.4902 0.7294],'filled');
set(gca,'xlim', [-1 1],'ylim', [0 100], 'XTick', [-1, -0.5 0 0.5 1],'YTick',[0 50 100],'TickDir','none','color','none','FontSize',10, 'FontName', 'Gill Sans MT');xlabel('Transdiagnostic hubs ({\itr})');ylabel('Epicenter overlap (%)'); axis square
text(T_transhubs.r+dx, T_epi_overlap.Functional+dy, txty,'FontSize',10, 'FontName', 'Gill Sans MT');
text(T_transhubs.r+dx, T_epi_overlap.Structural+dy, txty,'FontSize',10, 'FontName', 'Gill Sans MT');
subplot(1,2,2)% gradients vs disorder maps
col = [162,181,205]/255;
scatter(Grad1_dis_cor,Grad2_dis_cor,100, col, 'filled');
hold on%
text(Grad1_dis_cor+dx, Grad2_dis_cor+dy, txty,'FontSize',10, 'FontName', 'Gill Sans MT');
set(gca,'xlim', [-1 1],'ylim', [-1 1], 'XTick', [-1, -0.5 0 0.5 1],'YTick',[-1, -0.5 0 0.5 1],'TickDir','none','color','none','FontSize',10, 'FontName', 'Gill Sans MT');xlabel('G1 ({\itr})');ylabel('G2 ({\itr})'); axis square
hold off