% Analyses - Transcranial direct current stimulation reverses stroke-induced network alterations in mice
% Version 1.1
% 20230202
% - additional QM measures
dataPath = "D:\Labor\Projekte\Stroke_tDCS_rsfmri\B_rsfmri_stroke_tDCS_mice";
cd (dataPath)
if ~exist("results",'dir')
mkdir("results\")
end
savePath = fullfile(dataPath,'results');
addpath(".\scripts\toolbox\BCT\")
addpath("D:\Labor\Analysis\Matlab")
addpath("D:\Labor\Analysis\Toolbox\fMRI\FSLNets")
% load data
cd (fullfile(dataPath,"store_mat\"))
feat = readtable("PT_rsfmri_PT_tDCS_features.xlsx","VariableNamingRule","preserve");
load("PT_tDCS_Mat.mat");
codings = readtable("PT_tDCS_ID_list.xlsx","VariableNamingRule","preserve");
load("labels_mouse.mat")
% Tmaster_old = readtable("PT_tDCS_quali_stat.xlsx","VariableNamingRule","preserve");
Tmaster = readtable("PT_tDCS_quali_stat_rate20230201.xlsx","VariableNamingRule","preserve");
% Graphmaster = readtable ('PT_tDCS_graph_parameter.xlsx',"VariableNamingRule","preserve");
load("baseline_networks.mat");
load("PT_tDCS_mNET_all_raw.mat");
Tsub = readtable("PT_tDCS_Data_Table.xlsx",'VariableNamingRule','preserve');
data = readtable("PT_tDCS_graph_mean.xlsx",'VariableNamingRule','preserve');
Graph = readtable("Graph_data.xlsx");
% Setup Style
colors(1,:) = [.7 .7 .7];
colors(2,:) = hex2rgb('#D95319');
colors(3,:) = hex2rgb('#4DBEEE');
ax_name = 'Times';
ax_size = 16;
plotSize = [300 100 700 600];
group = unique(feat.group);
group = group([3 1 2]);
time = ["baseline","Day 3","Day14","Day 28"];
% define groupings
groups = unique(codings.group);
timepoints = unique(Tmaster.ses);
timelabel = {'baseline','Day 3','Day 14','Day 28'};
%%
% add group labels
for ix = 1:size(Tmaster,1)
Tmaster(ix,'group') = table2cell(codings(strcmp(codings.ID,Mat(ix).animal),'group'));
Graphmaster(ix,'group') = table2cell(codings(strcmp(codings.ID,Mat(ix).animal),'group'));
end
% adapt data
% % Tmaster(find(Tmaster.tSNR < 20 & (strcmp(Tmaster.qualy,'y')|strcmp(Tmaster.qualy,'d'))),'qualy') = {'n'};
% % Tmaster(find(strcmp(Tmaster.qualy,'n') & Tmaster.DVARS<0.4),"qualy")= {'x'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 2 group mean networks sorted by modularity
% GROUP DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
group_mean = nan(98,98,length(timepoints),length(groups));
group_delta = nan(98,98,length(timepoints),length(groups));
% group_idx = nan(max(max(n_good)),length(groups)*length(timepoints));
for gx = 1:length(groups)
for tx = 1:length(timepoints)
sidx = (gx-1)*length(timepoints)+tx;
ptx = find(strcmp(Tsub.group,groups(gx)) & strcmp(Tsub.ses,timepoints(tx)));
group_idx(1:length(ptx),sidx) = ptx;
% mNet(:,:,tx,gx) = reshape([Matsub(ptx).mNet_scrub_GSR],98,98,[]);
group_mean(:,:,tx,gx) = squeeze(mean(mNet(:,:,ptx),3,'omitnan'));
group_delta(:,:,tx,gx) = group_mean(:,:,tx,gx)-group_mean(:,:,1,gx);
end
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODULARITY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find modules of baseline networks
W = squeeze(mean(squeeze(group_mean(:,:,1,:)),3,'omitnan'));
if ~exist("M",'var')
% % Iterative community finetuning.
% W is the input connection matrix.
n = size(W,1); % number of nodes
M = 1:n; % initial community affiliations
Q0 = -1; Q1 = 0; % initialize modularity values
while Q1-Q0>1e-5 % while modularity increases
Q0 = Q1; % perform community detection
[M, Q1] = community_louvain(W, 1, M,'negative_sym');
end
end
% [M, Q1] = community_louvain(W, 1.1, M,'negative_sym')
[X,Y,On] = grid_communities(M);
x1 = unique(X);
x1 = x1(~isnan(x1));
y1 = unique(Y);
y1 = y1(~isnan(y1));
cd(savePath)
lx = labels.mNet(On);
for mx = 1:length(M)
modules{mx,1} = M(mx);
modules{mx,2} = lx{mx};
end
[~,ix] = sort(cell2mat(modules(:,1)));
Tm = cell2table(modules(ix,:),'VariableNames',{'ID','Region'});
writetable(Tm,'PT_tDCS_modules.xlsx')
if ~exist(fullfile(savePath,'Mat'),"dir")
mkdir(fullfile(savePath,'Mat'));
end
cd (fullfile(savePath,'Mat'))
% SupplFig3: Baseline Mat of Modules (annotated)
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(W(On,On),[],[-.4 0.8],0,labels.mNet(On));
axis square;
hold on
plot(X,Y,'k--','linewidth',2)
for xl = 2:length(x1)-1
plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
end
set(gca,'FontSize', 7)
xticklabels('')
t = title('baseline subnetworks defined by maximized modularity');
t.FontSize = 14;
c = colorbar;
c.FontSize = 14;
c.Label.String= 'functional connectivity (Z-scored)';
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'right';
c.Label.Position = [5 0];
saveas(gcf,'PT_tDCS_baseline_modules.png')
saveas(gcf,'PT_tDCS_baseline_modules.fig')
%%
if ~exist(fullfile(savePath,'Mat/RawGroupMat'),"dir")
mkdir(fullfile(savePath,'Mat/RawGroupMat'));
end
cd (fullfile(savePath,'Mat/RawGroupMat'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot group matrices sorted by modularity
for gx = 1:length(groups)
for tx = 1:length(timepoints)
%
% figure('units','normalized','outerposition',[0 0 1 1])
% plotmat_values_rsfMRI(group_mean(On,On,tx,gx),[],[-.4 .8],0,[]);
% axis square;
% hold on
% plot(X,Y,'k--','linewidth',2)
% for xl = 2:length(x1)-1
% plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
% plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
% end
% set(gca,'FontSize',ax_size,'FontName',ax_name)
% xticklabels('')
% yticklabels('')
% c = colorbar;
% c.FontSize = 14;
% c.Label.String= 'functional connectivity (Z-scored)';
% c.Label.Rotation = 270;
% c.Label.HorizontalAlignment = 'right'
% c.Label.Position = [5 -0.1]
% % title(strrep(strcat('Raw:_',groups{gx},'_@_',timepoints(tx)),'_',' '));
%
% saveas(gca,char(strcat('PT_tDCS_Raw_mNET_',groups(gx),'_',timepoints(tx),'.png')))
% saveas(gca,char(strcat('PT_tDCS_Raw_mNET_',groups(gx),'_',timepoints(tx),'.fig')))
% close
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(group_delta(On,On,tx,gx),[],[-.4 .4],0,[]);
axis square;
hold on
plot(X,Y,'k--','linewidth',2)
for xl = 2:length(x1)-1
plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
end
set(gca,'FontSize',ax_size,'FontName',ax_name)
xticklabels('')
yticklabels('')
ylabel(groups(gx))
c = colorbar;
c.FontSize = 20;
c.Label.String= strrep(strcat('delta FC (',timelabel(tx),'_- baseline)'),'_',' ');
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'right';
c.Label.Position = [5 -0.25];
% title(strrep(strcat('Delta to baseline:_',groups{gx},'_@_',timepoints(tx)),'_',' '));
saveas(gca,char(strcat('PT_tDCS_Delta_mNET_',groups(gx),'_',timepoints(tx),'.png')))
saveas(gca,char(strcat('PT_tDCS_Delta_mNET_',groups(gx),'_',timepoints(tx),'.fig')))
% pause(2)
close
end
end
%%
% plot tDCS delta stroke
group_DD = group_delta(:,:,:,3)-group_delta(:,:,:,2);%-group_delta(:,:,:,1);
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(group_DD(On,On,3),[],[-.4 .4],0,labels.mNet);
axis square;
hold on
plot(X,Y,'k--','linewidth',2)
for xl = 2:length(x1)-1
plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
end
set(gca,'FontSize',ax_size,'FontName',ax_name)
xticklabels('')
yticklabels('')
ylabel('')
c = colorbar;
c.FontSize = 20;
c.Label.String= 'delta FC Day 14 (PT tDCS - PT sham)';
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'right';
c.Label.Position = [5 -0.3];
cd (savePath)
saveas(gca,char('Delta_mNET_tDCS-sham_D14.png'))
saveas(gca,char('Delta_mNET_tDCS-sham_D14.fig'))
%%
% delta tDCS - sham sorted by hemispheres
group_dx = group_delta(:,:,:,3)-group_delta(:,:,:,2);
side = {'left','right'}
for ix = 1:2
stx = (ix-1)*48+1;
inp_data = group_dx(stx:stx+49,:,3);
Mx = M(stx:stx+49);
[~,Yx,Ox] = grid_communities(Mx);
plotmat_values_rsfMRI(inp_data(Ox,On),[],[-.4 .4],0,[]);
hold on
plot(X,Yx,'k--','linewidth',2)
title(side(ix))
end
%%
% compute delta to baseline and delta to control
D2_mat = squeeze(group_delta(:,:,:,3)-group_delta(:,:,:,2));
figure('units','normalized','outerposition',[0 0 1 1])
for tx=1:size(D2_mat,3)-1
subplot(1,3,tx)
plotmat_values_rsfMRI(D2_mat(On,On,tx+1),[],[-.3 .3],0,[]);
axis square;
hold on
plot(X,Y,'k--','linewidth',2)
for xl = 2:length(x1)-1
plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
end
set(gca,'FontSize',14)
title(timelabel(tx(2:end)))
end
%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fig 2
% plot delta to baseline of sham, tDCS, control
figure('units','normalized','outerposition',[0 0 1 1])
% prepare module labels
for xm = 1:length(unique(M))
xmean (xm) = mean([x1(xm+1),x1(xm)]);
end
ymean = zeros(1,length(xmean));
modules = {'default mode like','thalamus related','motor related','sensory related'};
%plot sham D14
dataMat = squeeze(group_delta(:,:,3,2));
plotmat_values_rsfMRI(deltaC(On,On),[],[-.4 .4],0,[]);
axis square;
hold on
plot(X,Y,'k--','linewidth',2)
for xl = 2:length(x1)-1
plot(get(gca,'xlim'),[y1(xl) y1(xl)],'k--','LineWidth',1);
plot([x1(xl) x1(xl)],get(gca,'ylim'),'k--','LineWidth',1);
end
set(gca,'FontSize',14)
title(timelabel(tx(2:end)))
%%
prm = {Gsub.strenght,Gsub.modularity,Gsub.SW,Gsub.CPL,Gsub.CC};
for gx = 1:length(prm)
tbx = prm{gx};
for ix = 1:size(group_idx,2)
ix_sub = group_idx(~isnan(group_idx(:,ix)),ix);
for ixs = 1:length(ix_sub)
end
end
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2 D - E
% Graph data
% get group means
cd (savePath)
if ~exist('Graph','dir')
mkdir("Graph")
end
cd ("Graph\")
parameter = Graph.Properties.VariableNames;
parameter = parameter(3:end);
groups = unique(Graph.group);
times = unique(Graph.time);
for gx = 1:length(parameter)
data = array2table(table2array(Graph(:,parameter(gx))),"VariableNames","value");
data.group = Graph.group;
data.time = Graph.time;
% two ANOVA of time and group for rotating beam errors
% [p,tbl,~,~] = anovan(data.value,{data.group data.time},'model',2,'varnames',{'Group','Time'});
% stat_table = cell2table(tbl(2:end,:));
% stat_table.Properties.VariableNames = tbl(1,:);
%
% writetable(stat_table,char(strcat('PT_tDCS_graph_',parameter(gx),'_anova2_stat.xlsx')))
% close all
% if p(3)<0.05
for tx = 1:length(times)
subdata = data(data.time == times(tx),:);
disp(strcat('Analyzing tp',num2str(times(tx))));
[ph, tbl_post, statistics] = anova1(subdata.value,subdata.group,'off');
stat_table = cell2table(tbl_post(2:end,:));
stat_table.Properties.VariableNames = tbl_post(1,:);
writetable(stat_table,char(strcat('PT_tDCS_graph_',parameter(gx),'_',num2str(times(tx)),'_anovaTimes_stat.xlsx')))
if ph < 0.05
[results,~,~,gnames] = multcompare(statistics,'Display','off');
ctbl = array2table(results,"VariableNames", ...
["Group","Control Group","Lower Limit","Difference","Upper Limit","P-value"]);
ctbl.("Group") = gnames(ctbl.("Group"));
ctbl.("Control Group") = gnames(ctbl.("Control Group"));
writetable(ctbl,char(strcat('PT_tDCS_graph_',parameter(gx),'_',num2str(times(tx)),'_anovaTimesPost_stat.xlsx')))
end
end
end
%%
% plot boxplots of graph measures
%
% compute delta of graph measures to baseline and control
feat_name = feat.Properties.VariableNames;
graph_measures = ["SW","CPL","CC"];
for gx = 1:3
rf = ~cellfun(@isempty,strfind(feat_name,graph_measures{gx}));
subdata = feat(:,rf);
graph_feat = feat_name(rf);
for tx = 1:length(time)
C_name(tx) = strcat(graph_measures{gx},'_',time{tx});
control_mean = squeeze(mean(table2array(subdata(strcmp("control",feat.group),tx)),'omitnan'));
for cx = 1:size(subdata,1)
ref_mean = squeeze(mean(table2array(subdata(strcmp(feat{cx,"group"},feat.group),1)),'omitnan'));
value = table2array(subdata(cx,tx));
base_value = table2array(subdata(cx,1));
if tx == 1
graph_delta(cx,tx) = value - ref_mean;
graph_delta_delta(cx,tx) = graph_delta(cx,tx) - control_mean;
else
if ~isnan(value)
if ~isnan(base_value)
graph_delta(cx,tx) = value - base_value;
else
graph_delta(cx,tx) = value - ref_mean;
end
graph_delta_delta(cx,tx) = graph_delta(cx,tx) - control_mean;
else
graph_delta(cx,tx) = nan;
graph_delta_delta(cx,tx) = nan;
end
end
end
end
T_graph_delta = array2table(graph_delta,"VariableNames",C_name);
T_graph_delta.group = feat.group;
T_graph_delta_delta = array2table(graph_delta_delta,"VariableNames",C_name);
T_graph_delta_delta.group = feat.group;
cd (savePath)
writetable(T_graph_delta_delta,'PT_tDCS_graph_delta_delta_CC.xlsx')
data_long = stack(T_graph_delta_delta,C_name,"IndexVariableName","group","NewDataVariableName","value");
figure(OuterPosition=plotSize)
ax1 = nexttile;
b = boxchart(ax1,table2array(T_graph_delta_delta),'GroupByColor',feat.group);
for bx = 1:length(b)
b(bx).BoxFaceColor = color(bx,:);
b(bx).BoxEdgeColor = 'k';
end
xticklabels('')
ylabel('time to righting reflex(s)')
legend(groupnames(1:2),'Location',"northwest");
legend('boxoff')
set(gca,'FontName',ax_name,'FontSize',ax_size)
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE 3: sensorimotor subnetwork
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd 'D:\Labor\Projekte\Stroke_tDCS_rsfmri\B_rsfmri_stroke_tDCS_mice\data';
load labels.mat;
sm_nw = group_delta(labels.nr14,labels.nr14,:,:);
raw_data = nan(98,98,4,3,13);
ind_delta = nan(98,98,4,3,13);
for gx = 1:length(groups)
for tx = 1:length(timepoints)
ls_subj = group_idx(:,(gx-1)*4+tx);
ls_subj = ls_subj(~isnan(ls_subj));
for j = 1:length(ls_subj)
raw_data(:,:,tx,gx,j) = mNet(:,:,ls_subj(j));
ind_delta(:,:,tx,gx,j) = mNet(:,:,ls_subj(j)) - group_mean(:,:,1,gx);
end
end
end
nx=2; % Anzahl an tests; Bonferroni Korrektur?
zvalp=abs(norminv(0.05/nx));
zval_5 = abs(norminv(0.05/nx)); % for p=0.05
zval_1 = abs(norminv(0.01/nx)); % for p=0.01
zval_01 = abs(norminv(0.001/nx)); % for p=0.01
for i=50000:-1:1
p=i/100000;
pval(i)=abs(norminv(p));
end
% number of permutations
n_permutes = 1000;
input_data = squeeze(ind_delta(labels.nr14,labels.nr14,3,[2,3],:));
A = squeeze(input_data(:,:,1,:));
B = squeeze(input_data(:,:,2,:));
%% raw mat
% stroke mat
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(squeeze(mean(A,3,'omitnan')),[],[-.4 .4],0,labels.ROI14_short);
axis square;
hold on
axis square;
plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
plot(get(gca,'xlim'),[ min(get(gca,'ylim')) max(get(gca,'ylim'))],'k-','LineWidth',2)
title ('PT sham')
set(gca,'FontSize',ax_size,'FontName',ax_name)
c = colorbar;
c.FontSize = 20;
c.Label.String= strrep(strcat('delta FC to baseline'),'_',' ');
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'center';
c.Label.Position = [5 0];
cd(savePath)
saveas(gca,char(strcat('SM_PT_sham.png')))
saveas(gca,char(strcat('SM_PT_sham.fig')))
% tDCS mat
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(squeeze(mean(B,3,'omitnan')),[],[-.4 .4],0,labels.ROI14_short);
axis square;
hold on
axis square;
plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
plot(get(gca,'xlim'),[ min(get(gca,'ylim')) max(get(gca,'ylim'))],'k-','LineWidth',2)
title ('PT tDCS')
set(gca,'FontSize',ax_size,'FontName',ax_name)
c = colorbar;
c.FontSize = 20;
c.Label.String= strrep(strcat('delta FC to baseline'),'_',' ');
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'center';
c.Label.Position = [5 0];
cd(savePath)
saveas(gca,char(strcat('SM_PT_tDCS.png')))
saveas(gca,char(strcat('SM_PT_tDCS.fig')))
%% delta mat
eval = squeeze(mean(B,3,'omitnan'))-squeeze(mean(A,3,'omitnan'));
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(eval,[],[-.4 .4],0,labels.ROI14_short);
axis square;
hold on
axis square;
plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
plot(get(gca,'xlim'),[ min(get(gca,'ylim')) max(get(gca,'ylim'))],'k-','LineWidth',2)
title ('uncorrected')
set(gca,'FontSize',ax_size,'FontName',ax_name)
c = colorbar;
c.FontSize = 20;
c.Label.String= strrep(strcat('delta FC (PT tDCS - PT sham) and delta to baseline'),'_',' ');
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'right';
c.Label.Position = [5 -0.5];
cd(savePath)
saveas(gca,char(strcat('SM_delta_all.png')))
saveas(gca,char(strcat('SM_delta_all.fig')))
% close
% end
% title(strrep(strcat('Delta to baseline:_',groups{gx},'_@_',timepoints(tx)),'_',' '));
conc_data = squeeze(reshape(input_data,14,14,1,[]));
for permi=1:n_permutes
randorder = randperm(size(conc_data,3));
temp = conc_data(:,:,randorder);
perm_map(:,:,permi) = squeeze(mean(temp(:,:,1:13),3,'omitnan'))-squeeze(mean(temp(:,:,14:26),3,'omitnan'));
end
clear ('zstat','mean_h0','std_h0');
mean_h0 = squeeze(mean(perm_map,3));
std_h0 = std(perm_map,1,3);
zstat = abs((eval-mean_h0))./(std_h0);
idx=abs(zstat)>zval_5;
for i=1:size(zstat,1)
for k=1:size(zstat,2)
[c idx] = min(abs(pval-zstat(i,k)));
pstat(i,k)=idx/100000;
end
end
FDR=5;
fr=FDR/100;
[hidx, crit_p, adj_ci_cvrg, adj_p]=fdr_bh(pstat,fr);
% FDR corrected significance map
mat=eval.*hidx;
mat(mat==0)=nan;
figure('units','normalized','outerposition',[0 0 1 1])
plotmat_values_rsfMRI(mat,[],[-.4 .4],0,labels.ROI14_short);
hold on
axis square;
plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
plot(get(gca,'xlim'),[ min(get(gca,'ylim')) max(get(gca,'ylim'))],'k-','LineWidth',2)
set(gca,'FontSize',ax_size,'FontName',ax_name)
title('FDR corrected (q = 0.05)')
yl = get(gca,'ylim');
xl = get(gca,'xlim');
for xy = 1:14
offset = 0.5+xy;
plot(xl,[offset offset],'k-','LineWidth',0.5)
plot([offset offset],yl,'k-','LineWidth',0.5)
end
c = colorbar;
c.FontSize = 20;
c.Label.String= strrep(strcat('delta FC (PT tDCS - PT sham) and delta to baseline'),'_',' ');
c.Label.Rotation = 270;
c.Label.HorizontalAlignment = 'right';
c.Label.Position = [5 -0.5];
cd(savePath)
saveas(gca,char(strcat('SM_delta_FDR.png')))
saveas(gca,char(strcat('SM_delta_FDR.fig')))
% close
% end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE 4: correltation Delta neuroscore Day 14 to change in graph parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd (fullfile(dataPath,"store_mat\"))
parameters = readtable("Corr_Graph_Neuroscore.xlsx",'VariableNamingRule','preserve');
% Figure 4 a/b (correlation of initial change in graph parameter with outcome)
graph_val = {'Delta_CPL_Day3','Delta_CC_D3'};
xlab = {'Delta norm. charac. path length (CPL) Day 3','Delta norm. clustering coefficient (CC) Day 3'};
for gx = 1:length(graph_val)
xval = graph_val(gx);
yval = "Delta_Neuroscore_Day21";
data = parameters((isnan(table2array(parameters(:,xval)))+isnan(table2array(parameters(:,yval))))==0,:);
groups = unique(data.group);
[r,p] = corr(table2array(data(:,xval)),table2array(data(:,yval)),'type','Spearman');
tbl = [r p];
tbl_names = {'Spearman','pval'};
corr_stat = array2table(tbl,'VariableNames',tbl_names);
cd (savePath)
writetable(corr_stat,strcat('PT_tDCS_Fig4ab',num2str(gx),'.xlsx'))
figure(OuterPosition=plotSize)
gscr = gscatter(table2array(data(:,xval)),table2array(data(:,yval)),data.group,colors(2:3,:));
for lx = 1:length(gscr)
gscr(lx).MarkerSize = 50;
end
hold on
mdl = fitlm(table2array(data(:,xval)),table2array(data(:,yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = [0 0 0];
end
% ylim([-10 60])
ylabel('Recovery Day 21 (% of initial deficit)')
xlabel(xlab(gx))
set(gca,'FontName',ax_name,'FontSize',ax_size)
title('')
legend(groups,'Location','northwest')
saveas(gcf,char(strcat('PT_tDCS_Fig4ab_',num2str(gx),'.png')))
saveas(gcf,char(strcat('PT_tDCS_Fig4ab_',num2str(gx),'.fig')))
end
%%
% figure 4 c+d
% correlation of change in graph parameters (SW/CPL) at Day 14 with
% recovery at Day 14
sub = 0; % define if plot with subgroups (sub = 1), else plot with whole group
graph_val = {'Delta_CPL_Day14-Day3','Delta_SW_Day14-Day3','Delta_CC_Day14-Day3'};
xlab = {'Delta norm. charac. path length (CPL) Day 14 - Day 3','Delta small worldness (SW) Day 14 - Day 3','Delta norm. clustering coeff. (CCc) Day 14 - Day 3'};
for gx = 1:length(graph_val)
xval = graph_val(gx);
yval = "Delta_Neuroscore_Day14";
data = parameters((isnan(table2array(parameters(:,xval)))+isnan(table2array(parameters(:,yval))))==0,:);
groups = unique(data.group);
if sub == 0
[r,p] = corr(table2array(data(:,xval)),table2array(data(:,yval)),'type','Spearman');
tbl = [r p];
tbl_names = {'Spearman','pval'};
corr_stat = array2table(tbl,'VariableNames',tbl_names)
else
for sbx = 1:length(groups)
subdata = data(strcmp(data.group,groups(sbx)),:);
[r,p] = corr(table2array(subdata(:,xval)),table2array(subdata(:,yval)),'type','Pearson');
tbl(sbx,:) = [r p];
end
tbl_names = {'Spearman','pval'};
corr_stat = array2table(tbl,'VariableNames',tbl_names,'RowNames',groups)
end
cd (savePath)
writetable(corr_stat,char(strcat('PT_tDCS_Fig4_',xval,'_sub_',num2str(sub),'.xlsx')))
figure(OuterPosition=plotSize)
gscr = gscatter(table2array(data(:,xval)),table2array(data(:,yval)),data.group,colors(2:3,:));
for lx = 1:length(gscr)
gscr(lx).MarkerSize = 50;
end
hold on
if sub == 0
mdl = fitlm(table2array(data(:,xval)),table2array(data(:,yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = [0 0 0];
end
else
cbx = ['rb'];
for sbx = 1:length(groups)
hold on
subdata = data(strcmp(data.group,groups(sbx)),:);
mdl = fitlm(table2array(subdata(:,xval)),table2array(subdata(:,yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = cbx(sbx);
end
end
end
grid
ylim([0 70])
ylabel('Recovery Day 14 (% of initial deficit)')
xlabel(xlab(gx))
set(gca,'FontName',ax_name,'FontSize',ax_size)
title('')
legend(groups,'Location','northeast')
saveas(gcf,char(strcat('PT_tDCS_Fig4cd_',xval,'_sub_',num2str(sub),'.png')))
saveas(gcf,char(strcat('PT_tDCS_Fig4cd_',xval,'_sub_',num2str(sub),'.fig')))
end
%%
% figure 5 a
% correlation of deficit Day 14 with baseline CPL
% figure 5 b
% correlation of intial deficit with baseline CPL
values = parameters.Properties.VariableNames;
func_val = {'Delta_Neuroscore_Day14','initial_Deficit'};
ylab = {'Recovery Day 14 (% of initial deficit)','initial motor deficit'};
for gx = 1:length(graph_val)
yval = func_val(gx);
xval = "CPL_base";
data = parameters((isnan(table2array(parameters(:,xval)))+isnan(table2array(parameters(:,yval))))==0,:);
groups = unique(data.group);
[r,p] = corr(table2array(data(:,xval)),table2array(data(:,yval)),'type','Spearman');
tbl = [r p];
tbl_names = {'Spearman','pval'};
corr_stat = array2table(tbl,'VariableNames',tbl_names);
cd (savePath)
writetable(corr_stat,strcat('PT_tDCS_Fig5',num2str(gx),'.xlsx'))
figure(OuterPosition=plotSize)
gscr = gscatter(table2array(data(:,xval)),table2array(data(:,yval)),data.group,colors(2:3,:));
for lx = 1:length(gscr)
gscr(lx).MarkerSize = 50;
end
hold on
if gx == 2
mdl = fitlm(table2array(data(:,xval)),table2array(data(:,yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = [0 0 0];
end
legend(groups,'Location','northwest')
else
mdl = fitlm(table2array(data(strcmp(data.group,'PT sham'),xval)),table2array(data(strcmp(data.group,'PT sham'),yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = [1 0 0];
end
mdl = fitlm(table2array(data(strcmp(data.group,'PT tDCS'),xval)),table2array(data(strcmp(data.group,'PT tDCS'),yval)));
p = plot(mdl);
delete(p(1))
for px = 2:4
p(px).LineWidth = 1.25;
p(px).Color = [0 0 1];
end
legend({'PT sham','PT tDCS **'},'Location','northwest')
end
% ylim([-10 60])
ylabel(ylab(gx))
xlabel('Norm. charac. path length (CPL) baseline')
set(gca,'FontName',ax_name,'FontSize',ax_size)
title('')
saveas(gcf,char(strcat('PT_tDCS_Fig5_',num2str(gx),'.png')))
saveas(gcf,char(strcat('PT_tDCS_Fig5_',num2str(gx),'.fig')))
end
%% Quality Statistics
% Count data set quality per group
qm_labels = unique(Tmaster.qualy);
n=zeros(1,length(qm_labels));
for qi = 1:length(qm_labels)
n(1,qi) = nnz(strcmp(Tmaster.qualy,qm_labels(qi)));
end
% qm_lx = {'no data','unknown','registration error','exessive motion','distorsions','passed'};
qm_lx = {'registration error','exessive motion','distorsions','passed'};
% count available datasets
n_good = zeros(1,length(timepoints));
for gx = 1:length(groups)
for qx = 1:length(qm_labels)
n(gx+1,qx) = sum(strcmp(Tmaster.group,groups(gx)) & (strcmp(Tmaster.qualy,qm_labels(qx))));
end
end
qrow = {'All','Control','PT sham','PT tDCS'};
Q_tab = array2table(n,'VariableNames',qm_lx,'RowNames',qrow);
cd(fullfile(savePath,'QC'))
writetable(Q_tab,'PT_tDCS_N_scans_all.xlsx')
writetable(N_tab,'PT_tDCS_N_scans_groups.xlsx')
% select networks with good registration and FD<0.05
% Tsub = Tmaster((strcmp(Tmaster.qualy,'y')|strcmp(Tmaster.qualy,'d'))&(Tmaster.tSNR>20),:);
% Matsub = Mat((strcmp(Tmaster.qualy,'y')|strcmp(Tmaster.qualy,'d'))&(Tmaster.tSNR>20));
% Gsub = Graphmaster((strcmp(Tmaster.qualy,'y')|strcmp(Tmaster.qualy,'d'))&(Tmaster.tSNR>20),:);
for gx = 1:length(groups)
for tx = 1:length(timepoints)
group_list((gx-1)*length(groups)+tx) = strcat(groups(gx),'_',timepoints(tx));
n_good(gx,tx) = sum(...
strcmp(Tmaster.group,groups(gx)) ...
& strcmp(Tmaster.ses,timepoints(tx))...
& strcmp(Tmaster.qualy,'y'));
end
end
r_names = {'Control (N = 12)','PT-sham (N = 13)','PT-tDCS (N = 13)'};
N_tab = array2table(n_good,'VariableNames',timepoints,'RowNames',r_names);
% Prepare tables for printing
% Get the table #1 in string form.
T2String = evalc('disp(Q_tab)');
% Use TeX Markup for bold formatting and underscores.
T2String = strrep(T2String,'','\bf');
T2String = strrep(T2String,'','\rm');
T2String = strrep(T2String,'_','\_');
% Get the table #2 in string form.
TString = evalc('disp(N_tab)');
% Use TeX Markup for bold formatting and underscores.
TString = strrep(TString,'','\bf');
TString = strrep(TString,'','\rm');
TString = strrep(TString,'_','\_');
% Get a fixed-width font.
FixedWidth = get(0,'FixedWidthFontName');
%%
% suppl Fig 1 (available datasets)
figure('units','normalized','outerposition',[0 0 1 1])
annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
'FontName',FixedWidth,'Units','Normalized','Linestyle','none','Position',[0.1 -0.3 1 1]);
annotation(gcf,'Textbox','String',T2String,'Interpreter','Tex',...
'FontName',FixedWidth,'Units','Normalized','Linestyle','none','Position',[0.1 -0.1 1 1]);
subplot(2,2,3)
px = pie(n(1,:),[],qm_lx);
subplot(2,2,[2,4])
b = bar(categorical(qm_lx),n(2:end,:)');
for bx = 1:length(b)
b(bx).FaceColor = colors(bx,:);
end
legend(groups,'Location','northeast')
legend('boxoff')
saveas(gcf,'PT_tDCS_qm_overview.png')
saveas(gcf,'PT_tDCS_qm_overview.fig')
%%
% suppl. Fig 2 (QC meassures)
parameter = {'SNR','tSNR','DVARS','FD'};
stat_pm = {'mean','std',};
stats = groupsummary(Tsub,"group",stat_pm,parameter);
stats(4,:) = groupsummary(Tsub,"ID",stat_pm,parameter);
stats(4,"group") = {'All'};
% Get the table in string form.
TString = evalc('disp(stats)');
% Use TeX Markup for bold formatting and underscores.
TString = strrep(TString,'','\bf');
TString = strrep(TString,'','\rm');
TString = strrep(TString,'_','\_');
% Get a fixed-width font.
figure('units','normalized','outerposition',[0 0 1 1])
annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
'FontName',FixedWidth,'Units','Normalized','Linestyle','none','Position',[0.1 -0.1 1 1]);
tlx = tiledlayout(3,4);
for kx = 1:length(parameter)
nexttile(kx+4,[2,1])
b = boxchart(table2array(Tsub(:,parameter(kx))),'GroupByColor',Tsub.group);
for bx = 1:length(b)
b(bx).BoxFaceColor = colors(bx,:);
b(bx).BoxEdgeColor = 'k';
b(bx).MarkerColor = colors(bx,:);
end
xticklabels('')
ylabel(parameter(kx))
cx = get(gca,'ylim');
ylim([0 cx(2)])
set(gca,'FontName',ax_name,'FontSize',ax_size)
end
l=legend(b,'NumColumns',3);
legend('boxoff')
l.Layout.Tile = 'South';
saveas(gcf,'PT_tDCS_qm_parameters.png')
saveas(gcf,'PT_tDCS_qm_parameters.fig')