Figure1_Degree_Centrality_Caulation&Surface_Rendering.m

% Script for calculating the Degree Centrality maps and the Euclidean Distance between the children and adolescents groups.
% Author: Dong HaoMing
%% path setting
clear all;close all;

Fsdir = ['/opt/software/freesurfer'];
addpath(genpath(Fsdir))
ccs_dir = ['/opt/software/CCS/ccs_develop/'];
addpath(genpath(ccs_dir))
npypath = ['/opt/github/npy-matlab-master/'];
addpath(genpath(npypath))

outdir = ['/dir/surfFC/output/'];
pngdir = [outdir '/Figure/'];
mkdir(pngdir)
agegroup = {'child','adolescent'};

%% surface structure loading for rendering
s_lh1 = SurfStatReadSurf([Fsdir '/subjects/fsaverage5/surf/lh.inflated']);
s_rh1 = SurfStatReadSurf([Fsdir '/subjects/fsaverage5/surf/rh.inflated']);
s_lh2 = SurfStatReadSurf([Fsdir '/subjects/fsaverage5/surf/lh.mid']);
s_rh2 = SurfStatReadSurf([Fsdir '/subjects/fsaverage5/surf/rh.mid']);
s_lh.tri = s_lh1.tri;
s_lh.coord = s_lh1.coord*0.7 + s_lh2.coord*0.3;
s_rh.tri = s_rh1.tri;
s_rh.coord = s_rh1.coord*0.7 + s_rh2.coord*0.3;
fcb = ['/path/colormap.png']; % load your colormap for surface rendering
cmap = ccs_mkcolormap(fcb);

%% brainmask loading (if you have)
maskdir = ['/home/donghm/Projects/Subcortical_Gradient/'];
lhmask = load_nifti([maskdir '/brainmask_lh_fs.nii.gz']);
lhmask = lhmask.vol;
rhmask = load_nifti([fmaskdir '/brainmask_rh_fs.nii.gz']);
rhmask = rhmask.vol;
brainmask = [lhmask;rhmask];
surfnum = sum(brainmask);

%% FC matrix loading
% Here the FC matrixes were generated in Dong et al., 2021 PNAS. 
FCdir = ['/path_of_your_FCmatrix/']; 
fcFCz = [FCdir '/child_FCz.npy']; % FC matrix of children group
tmpcFCz = readNPY(fcFCz);
faFCz = [FCdir '/adolescent_FCz.npy']; % FC matrix of adolescents group
tmpaFCz = readNPY(faFCz);
% theresholding 90 percentile for each row of the FC matrix 
for vv = 1:size(tmpaFCz,1)
    if mod(vv,500)==0
        disp([num2str(vv*100/size(tmpaFCz,1)) '% completed'])
    end
    PaFCz = prctile(tmpaFCz(vv,:),90);
    PcFCz = prctile(tmpcFCz(vv,:),90);
    tmpaFCz(vv,tmpaFCz(vv,:)<PaFCz) = 0;
    tmpcFCz(vv,tmpcFCz(vv,:)<PcFCz) = 0;
end

% Binarize the FC matrix to caculate degree centrality
% The FC matrix was asymmetric after thresholding and the connenction numbers (degree) are identical for each row. So the matrix is summed for each colume to caculate the degree centrality.    
tmpaFCz(tmpaFCz<0)=0; 
tmpcFCz(tmpcFCz<0)=0;
tmpaFCz(tmpaFCz>0)=1;
tmpcFCz(tmpcFCz>0)=1;
DCa= sum(tmpaFCz,1);
DCc= sum(tmpcFCz,1);
clear tmpaFCz tmpcFCz 

%% Rendering degree centrality map 
% Surface map can be upsampled to fsaverage surface to obtain a higher resolution, here fsaverage5 is used as an example 
surfmap = zeros(size(brainmask)); 
surfmap(brainmask == 1) = DCa;
tmplh = surfmap(1:10242);
tmprh = surfmap(10243:end);

figure('visible','off')
SurfStatViewData(tmplh,s_lh1); 
colormap(cmap) % change the colormap 
fout = [ pngdir '/Adolescents_SurfDC_lh.png'];
saveas(gcf,fout,'png')
close

figure('visible','off')
SurfStatViewData(tmprh,s_rh1);
colormap(cmap)
fout = [ pngdir '/Adolescents_SurfDC.png'];
saveas(gcf,fout,'png')
close

surfmap = zeros(size(brainmask));
surfmap(brainmask == 1) = DCc;
tmplh = surfmap(1:10242);
tmprh = surfmap(10243:end);
figure('visible','off')
SurfStatViewData(tmplh,s_lh1);
colormap(cmap)
fout = [ pngdir '/Children_SurfDC.png'];
saveas(gcf,fout,'png')
close

figure('visible','off')
SurfStatViewData(tmprh,s_rh1);
colormap(cmap)
fout = [ pngdir '/Children_SurfDC.png'];
saveas(gcf,fout,'png')
close

%% caculate the Euclidean distance between the affnity matrixes of children and adolescents group
% Here, the affnity (distance) matrix was generated by caculting the cosine distance of the FC matrix, which were first generated in Dong et al., 2021 PNAS. 
postfix = ['euclidean'];
% loading the affnity matrix
fCSurfdist = [FCdir '/' agegroup{1} '_dist.npy'];
CSurfdist = readNPY(fCSurfdist);
fASurfdist = [FCdir '/' agegroup{2} '_dist.npy'];
ASurfdist = readNPY(fASurfdist);
distAC= [];
for vv = 1:size(CSurfdist,1)
    distAC(vv) = pdist2(CSurfdist(:,vv)',ASurfdist(:,vv)',postfix);
end

% caculate network-level Eculidean distance
netnum = 7;
netdir =  ['/github/Yeo_JNeurophysiol11_FreeSurfer/fsaverage5/label/']
flh7 = [netdir '/lh.Yeo2011_7Networks_N1000.annot'];
frh7 = [netdir '/rh.Yeo2011_7Networks_N1000.annot'];
[lv7, ll7, lc7] = read_annotation(flh7 );
[rv7, rl7, rc7] = read_annotation(frh7 );
bl7 = [ll7;rl7];
bl7 = bl7 .* brainmask;
bl7 = bl7(bl7>0);
NetworkED= zeros(sum(brainmask),7);
for nn = 1:7
    tmpidx = find(bl7 == lc7.table(nn+1,5));
    bary(nn) = mean(distAC(tmpidx));
    NetworkED(tmpidx,nn) = distAC(tmpidx);
    tmpidxnan = find(bl7 ~= lc7.table(nn+1,5));
    NetworkED(tmpidxnan,nn) = nan;
    erry(nn) = std(distAC(tmpidx));
end
[p,tbl,stats] = anova1(NetworkED);
c = multcompare(stats)
[sortbary,sortidx] = sort(bary);
figure('visible','off')
for nn = 1:7
    hbar(nn) = bar(nn,sortbary(nn));
    tmpxl{nn} = num2str(sortidx(nn));
    set(hbar(nn),'FaceColor',lc7.table(sortidx(nn)+1,1:3)/256);
    hold on;
end
hold on
set(gca,'XTick',[1:7],'XTickLabel',tmpxl);
errorbar(sortbary,erry(sortidx), 'k', 'Linestyle', 'None');
fout = [ pngdir '/A-C_SurfDist_' postfix 'Hist_' num2str(netnum) 'Net_Bar.png'];
saveas(gcf,fout,'png')
set(gcf, 'PaperPositionMode', 'auto', 'Color', ...
    'white', 'InvertHardCopy','off');
close

surfmap = zeros(size(brainmask));
surfmap(brainmask == 1) = distAC;
tmplh = surfmap(1:10242);
tmprh = surfmap(10243:end);

figure('visible','off')
SurfStatViewData(tmplh,s_lh1);
colormap(cmap)
fout = [ pngdir '/A-C_SurfDist_' postfix 'Surfmap_lh.png'];
saveas(gcf,fout,'png')
close

figure('visible','off')
SurfStatViewData(tmprh,s_rh1);
colormap(cmap)
fout = [ pngdir '/A-C_SurfDist_' postfix 'Surfmap_rh.png'];
saveas(gcf,fout,'png')
close

% save the areas with the largest Euclidean distance between groups
T = 90; % defining the percentage of areas kept to be examined
Tp = prctile(distcorr,T);
tmpmask = brainmask;
    
surfmap = zeros(size(brainmask));
surfmap(brainmask == 1) = distcorr;
surfmap(surfmap<Tp) = 0;
tmpmask(surfmap<Tp) = 0;
lhmaskv = tmpmask(1:10242);
rhmaskv = tmpmask(10243:end);
tmplh = surfmap(1:surfnum/2);
tmprh = surfmap(surfnum/2+1:end);
mask_outdir = ['path/to/save/mask'];
fout = [ mask_outdir '/output/SurfDist_' postfix '_T' num2str(T) '_lhP.png'];
figure('visible','off')
SurfStatViewData(tmplh,s_lh1);
saveas(gcf,fout,'png')
close
fout = [ mask_outdir '/output/SurfDist_' postfix '_T' num2str(T) '_rhP.png'];
figure('visible','off')
SurfStatViewData(tmprh,s_rh1);
saveas(gcf,fout,'png')
close
lhmaskhdr.vol = tmplh;
rhmaskhdr.vol = tmprh;
fout = [mask_outdir '/' postfix 'dist_T' num2str(T) '_lh.nii.gz']
save_nifti(lhmaskhdr,fout)
fout = [mask_outdir '/' postfix 'dist_T' num2str(T) '_rh.nii.gz']
save_nifti(rhmaskhdr,fout)

% clustering the areas above 500 mm^2 
fscommand  = ['mri_surfcluster --in ' mask_outdir '/' postfix 'dist_T' num2str(T) '_lh.nii.gz --hemi lh --o  ' mask_outdir '/euclideandist.cluster.lh.nii.gz ' ...
  '--surf white --annot aparc --sum euclideandist.lh.cluster.summary --ocn ' mask_outdir '/euclideandist.lh.cluster.ocn.' num2str(T) '_500.nii.gz ' ...
  '--subject fsaverage5 --nofixmni --minarea 500 --thmin 0.1'];
[status,results] = system(fscommand,'-echo');
fscommand  = ['mri_surfcluster --in ' mask_outdir '/' postfix 'dist_T' num2str(T) '_rh.nii.gz --hemi rh --o  ' mask_outdir '/euclideandist.cluster.rh.nii.gz ' ...
  '--surf white --annot aparc --sum euclideandist.lh.cluster.summary --ocn ' mask_outdir '/euclideandist.rh.cluster.ocn.' num2str(T) '_500.nii.gz ' ...
  '--subject fsaverage5 --nofixmni --minarea 500 --thmin 0.1'];
[status,results] = system(fscommand,'-echo');
% binarize the clusters
fscommand  = ['mri_binarize --i '  mask_outdir '/euclideandist.lh.cluster.ocn.' num2str(T) '_500.nii.gz --min 0.01 --o ' ourdir '/euclideandist.lh.cluster.bin.' num2str(T) '_500.nii.gz' ];
[status,results] = system(fscommand,'-echo');
fscommand  = ['mri_binarize --i '  mask_outdir '/euclideandist.rh.cluster.ocn.' num2str(T) '_500.nii.gz --min 0.01 --o ' ourdir '/euclideandist.rh.cluster.bin.' num2str(T) '_500.nii.gz' ];
[status,results] = system(fscommand,'-echo');