% 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,'<strong>','\bf');
    T2String = strrep(T2String,'</strong>','\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,'<strong>','\bf');
    TString = strrep(TString,'</strong>','\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,'<strong>','\bf');
    TString = strrep(TString,'</strong>','\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')