2023_Blaschke_Stroke_tDCS_rsfMRI/scripts/addendum/runCBI.m

% Define parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parent_path = "G:\Projects\CBIrep_Imaging\subjects";
AUC = 10:20;  % area under the curve for 
cut = 1;        % remove motion contaminated volumes
want2plot = 1;    % visualize raw data? (plot=1)
GSR = 1;  % perform global signal regression (yes=1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup
cd (parent_path)
addpath("..\scripts\toolbox\DVARS-master\")
addpath("..\scripts\toolbox\DVARS-master\Nifti_Util\")
addpath("..\scripts\toolbox\BCT\")
addpath("..\scripts\toolbox\spm12")
addpath("..\scripts\toolbox\fMRI-Quality-Checker-master")
addpath("..\scripts\toolbox\spm12\matlabbatch")
addpath("..\scripts\toolbox\nifti_utils-master\")
addpath("..\scripts\toolbox\FSLNets\")
addpath("..\scripts\")
addpath("D:\Labor\Analysis\fmri\nifti_utils-master\nifti_utils-master\external\NIFTI")

px = dir('sub*');

cd ..\
subjects = readtable('subject_list.xlsx');
participants = readtable('participants.csv');
part = cell2mat(table2cell(participants(:,"subject")));

master = cell(size(subjects,1),6);

for k = 1:size(subjects,1)
    ix = subjects(k,"ID");
    ix = table2array(ix);
    ix = char(ix);
    ix = ix(5:end);
    master{k,1}= char(table2array(subjects(k,"ID")));
    master{k,2} = char(table2array(subjects(k,"Time")));
    master{k,3}= char(table2array(participants(part==str2double(ix),"condition")));
end  

groups = unique(participants.condition);

mNet_raw  = zeros(size(subjects,1),98,98);
mNet_tr  = zeros(size(subjects,1),98,98);

quali_stat = cell(size(subjects,1),9);
quali_tag = ["ID";"Time";"tSNR";"SNR";"mean DVARS";...
    "bad volumes (%)";"mean FD";"max FD";"GNI"];

%% Run through data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% level of subjects
for i = 41:size(subjects,1)
 %%
    cd (parent_path)
    subjx = table2array(subjects(i,"ID"));
    tmx = table2array(subjects(i,"Time"));


quali_stat{i,1} = char(subjx);
quali_stat{i,2} = char(tmx);

    cd (char(subjx))
    cd (char(tmx))
        

% grab images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % a) get structural 
            cd (fullfile(parent_path,subjx,tmx,"anat/"))
            struct_orig = dir("*.1.nii.gz");
            struct_name = struct_orig(length(struct_orig)).name;  
            struct_fullfile = fullfile(struct_orig(length(struct_orig)).folder,struct_orig(length(struct_orig)).name); 
         
        % b) get fMRI
            cd (fullfile(parent_path,subjx,tmx,"func/"))
            fmri_orig = dir("*.1.nii.gz");
            fmri_name = fmri_orig(length(fmri_orig)).name;
            fMRI_fullfile = fullfile(fmri_orig(length(fmri_orig)).folder,fmri_orig(length(fmri_orig)).name); 
%             func_file = gunzip(fMRI_fullfile);
%             func_spm = spm_vol(char(func_file));
%             func_4Dimg = niftiread(func_file);

            %get func mask
%             fmri_mask = dir("*AnnoSplit_rsfMRI*");
%             
%             spm_image('Display', char(func_file))

        % c) get Regression File
            cd (fullfile(parent_path,subjx,tmx,"func/"))
            cd regr\
            SFRGR_orig = dir("*_SFRGR.nii.gz");
            TCfile = dir("MasksTCsSplit*.txt.mat");
            fmri_mask = dir("*thres_mask.nii.gz");
            SFRGR_file = char(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
%           fmri_mask = dir("mask.nii.gz");

            mask = load_untouch_nii(fullfile(fmri_mask.folder,fmri_mask.name));
            mask_img = mask.img;
%             mask_img = flip(mask_img,1);
%             mask_img = flip(mask_img,3);
            mask_img(mask_img>0) = 1;

% % % Global signal regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RestData = gunzip(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
WB_mask = gunzip(fullfile(fmri_mask.folder,fmri_mask.name));
GNI = Rest2GNI(char(RestData),char(WB_mask),100);
quali_stat{i,9} = num2str(GNI);

SFRGR = load_untouch_nii(SFRGR_file);
SFRGR_img = SFRGR.img;


% % % SNR (spatial, temporal)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   info = niftiinfo(fMRI_fullfile);
            V1 = load_untouch_nii(fMRI_fullfile);
            V2 = V1.img;
            [Ni, Nj, Nk, Nt] = size(V2);
            func_2Dimg = reshape(V2, Ni*Nj*Nk, Nt);
            [Nv, Nt] = size(func_2Dimg);

            func_mean = mean(func_2Dimg, 2,'omitnan');
            func_mean3D = reshape(func_mean, Ni, Nj, Nk);
            func_stddev = std(double(func_2Dimg), 0, 2,'omitnan');
            func_stddev3D = reshape(func_stddev, Ni, Nj, Nk);
            
            tSNR_2Dimg = mean(func_2Dimg, 2,'omitnan')./std(double(func_2Dimg), 0, 2,'omitnan');
            tSNR_3Dimg = reshape(tSNR_2Dimg, Ni, Nj, Nk);
            
            mask_img_flip = flip(mask_img,1);
            mask_img_flip = flip(mask_img_flip,3);
            
            tSNR = squeeze(mean2(tSNR_3Dimg(mask_img_flip==1)));
            SNR = squeeze(mean(func_mean3D(mask_img_flip==1)))/std(func_stddev3D(mask_img_flip==0));
           
            cd ..\
            niftiwrite(tSNR_image,strcat(fmri_name(1:end-10),'tSNR.nii'))

quali_stat{i,3} = num2str(tSNR);
quali_stat{i,4} = num2str(SNR);       



SFRGR_img(mask_img==0)=nan;
GSR_signal = squeeze(mean(mean(mean(SFRGR_img,3,'omitnan'),2,'omitnan'),1,'omitnan'));


     
% get DVARS
       
        X0 = size(V2,1); Y0 = size(V2,2); Z0 = size(V2,3); T0 = size(V2,4);
        I0 = prod([X0,Y0,Z0]);
        Y  = reshape(V2,[I0,T0]); clear V2 V1;
        [DVARS,DVARS_Stat]=DVARSCalc(Y,'scale',1/100,'TransPower',1/3,'RDVARS','verbose',1);
        [V,DSE_Stat]=DSEvars(Y,'scale',1/100);
        dvars_points = DVARS_Stat.pvals;
        bad_points = find(dvars_points(5:end)<0.001);

quali_stat{i,5} = squeeze(mean(DVARS));

        if isempty(bad_points)
            badTP = 0;
quali_stat{i,6} = 0;
        else
            badTP = bad_points-4;
quali_stat{i,6} = badTP/(info.ImageSize(4)-5);
        end

        scratch = zeros(1,4*length(bad_points));
        for ptx = 1:length(bad_points)
            scratch(((ptx-1)*4+1):((ptx-1)*4+4)) = bad_points(ptx)-2:bad_points(ptx)+1;
        end
        scratch = unique(scratch(scratch>0));
        scratch(scratch>info.ImageSize(4)-5) = [];

        
      
      
 % getFD
        cd (fullfile(parent_path,subjx,tmx,"func/"))
        if exist("txtRegrPython","dir")
        cd txtRegrPython\
        txtRegr_files = dir("*mcf*");
            if size(txtRegr_files,1)>1
            txtRegr = txtRegr_files(floor(length(txtRegr_files)/2)).name;
            txtRegr_table = readtable(txtRegr);
                if ~isempty(txtRegr_table)
                    Regression_parameters = txtRegr_table(:,5:10);
                    MP = cell2mat(table2cell(Regression_parameters));
                    FD_measures = calculateFD(MP, .5, 0.5);
                    FD = FD_measures.FD;
                   
                else
                    FD = [];
                    disp('Dataset has no available regression file')
                end
            else 
                FD = [];
                disp('Dataset has no available regression file')
            end
        else
                FD = [];
                disp('Dataset has no available regression file')
        end

        if ~isempty(FD)
quali_stat{i,7} = num2str(squeeze(mean(FD)));
quali_stat{i,8} = num2str(max(FD));
        end
%%
if ~isempty(TCfile)
        [mNet,mNet_Trim,SW,CC,CPL] = rsfmri_metrics(TCfile,GSR_signal,AUC,cut,scratch,want2plot,GSR);

saveas(gcf,strcat(fmri_name(1:end-10),'mNET.png'))
pause(2)


mNet_raw(i,:,:) = mNet;
mNet_tr (i,:,:) = mNet_Trim;
master{i,4} = num2str(SW);
master{i,5} = num2str(CC);
master{i,6} = num2str(CPL);
end
%%
% plot parameters
       if want2plot==1

            figure('units','normalized','outerposition',[0 0 1 1])
           
            subplot(231)

                histogram(tSNR_3Dimg(mask_img_flip==1),'Normalization','probability');
                xlabel('tSNR')
                title(strrep(strcat('mean tSNR =_',num2str(tSNR)),'_',' '))
      
            subplot(232)
                p = plot(DVARS(5:end));
                xlim([0 length(DVARS)-4])
                ylim([0 1])
                hold on
                    p1 = plot(get(gca,'xlim'),[squeeze(mean(DVARS)) squeeze(mean(DVARS))],'r-');
                hold off
                ylabel('DVARS')
                xlabel('Timepoints')
                title(strrep(strcat('mean DVARS =_',num2str(squeeze(mean(DVARS)))),'_',' '))
        
 
            subplot(233)
                h = histogram(FD,'BinEdges',0:0.01:0.5);
                ylabel('number of timepoints')
                xlabel('framewise displacemen (mm)')
                title(strrep(strcat('max. FD =_',num2str(max(FD)),'_mm'),'_',' '))
 
            for ki=1:3
            subplot(2,3,3+ki)
                ix = [5 8 11];
                imagesc(imrotate(tSNR_3Dimg(:,:,ix(ki)),-90))

                M = mask_img_flip(:,:,ix(ki)); 
%                 M = flip(flip(mask_img(:,:,20-ix(ki)),3),1); 
%                 M = mask_img(:,:,ix(ki));
                hold on
                visboundaries(imrotate(M,-90),'Color','b');
                colormap('hot')
            end
           cd (fullfile(parent_path,subjx,tmx,"func/"))
           saveas(gcf,strcat(fmri_name(1:end-10),'quality.png'))
           
       end

       %%
       pause(2)
       
%%
disp(strrep(strcat('Analysis of ID_',char(subjx),'_timepoint_',char(tmx),'_complete'),'_',' '))
disp(strrep(strcat('- tSNR:_',num2str(tSNR)),'_',' '))
disp(strrep(strcat('- SNR:_',num2str(SNR)),'_',' '))
disp(strrep(strcat('- mean DVARS:_',num2str(quali_stat{i,5})),'_',' '))
disp(strrep(strcat('- mean FD:_',num2str(quali_stat{i,7})),'_',' '))
% 
x = input('quality?');
if ~isempty(x)
    x=1;
else
    x=0;
end
qm_list{i,1}=char(subjx);
qm_list{i,2}=char(tmx);
qm_list{i,3}=num2str(x)

close all

end


head_master = ["ID","Time","group","SW","CC","CPL"];
T = cell2table(master);
T.Properties.VariableNames = head_master;

Q = cell2table(quali_stat);
Q.Properties.VariableNames = quali_tag;

cd (parent_path)
cd ..
cd analysis\

writetable(T,char(strrep(strcat('CBI_graph_',date,'.xlsx'),'-','_')));
% save(char(strrep(strcat('CBI_rawNet_',date,'.m'),'-','_')),mNet_raw);
% save(char(strrep(strcat('CBI_trimNet_',date,'.m'),'-','_')),mNet_tr);

%%
data_set = master(cellfun(@(x) ~isempty(x),master(:,3)),1:3);
group = unique(data_set(:,3));
timepoints = unique(data_set(:,2));

[X,Y,On] = grid_communities(M);
%%
C = cell2mat(quali_stat(:,5));
QM = cell2mat(quali_stat(:,3));


for gx = 1:length(group)
   figure('units','normalized','outerposition',[0 0 1 1])
    for tx = 1:length(timepoints)
    subNet = mNet_raw(...
        strcmp(master(:,3), char(group(gx))) ...
        & strcmp(master(:,2),char(timepoints(tx))) ...
        & C < 0.5 ...
        & QM > 0.25 ...
        ,:,:);
    subNet(subNet==0) = nan;

    group_subNet(:,:,tx,gx) = squeeze(mean(subNet,1,'omitnan'));
    delta_subNet(:,:,tx,gx) = squeeze(group_subNet(:,:,tx,gx)) - squeeze(group_subNet(:,:,1,gx));

    subplot(2,length(timepoints),tx)
      plotmat_values_rsfMRI(group_subNet(On,On,tx,gx),[],[-1 1],0,[]);
        hold on
        axis square;
        plot(X,Y,'k','linewidth',2)
        title(strrep(strcat('mNet_',group(gx),'_',timepoints(tx)),'_',' '));
   
    subplot(2,length(timepoints),tx+length(timepoints))
      plotmat_values_rsfMRI(delta_subNet(On,On,tx,gx),[],[-.5 .5],0,[]);
        hold on
        axis square;
        plot(X,Y,'k','linewidth',2)
        title(strrep(strcat('mNet_',group(gx),'_',timepoints(tx)),'_',' '));
    
    end
    
end


for gx = 1:length(group)
     figure('units','normalized','outerposition',[0 0 1 1])
  
    d2_subNet(:,:,:,gx) = delta_subNet(:,:,:,gx) - delta_subNet(:,:,:,4);
    
    for tx = 1:3
     subplot(1,length(timepoints),tx)
      plotmat_values_rsfMRI(d2_subNet(On,On,tx,gx),[],[-.5 .5],0,[]);
        hold on
        axis square;
        plot(X,Y,'k','linewidth',2)
        title(strrep(strcat('mNet_',group(gx),'_',timepoints(tx)),'_',' '));
    end
end

%%

for gx = 1:length(group)
%    figure('units','normalized','outerposition',[0 0 1 1])
    for tx = 1:length(timepoints)
   mnx = master(...
        strcmp(master(:,3), char(group(gx))) ...
        & strcmp(master(:,2),char(timepoints(tx))) ...
        & C<0.5...
        ,5);
   mnx(find(cellfun(@isempty,mnx)))=[];
   group_parameter(1:length(mnx),tx,gx) = cellfun(@(mnx)str2double(mnx), mnx);
   
   
    end
end

group_parameter (group_parameter==0)=nan;

mean_group = squeeze(mean(group_parameter,1,'omitnan'))

for dx = 1:length(group)
    for tx = 1:length(timepoints)
        delta_group (tx,dx) = (mean_group(tx,dx)-mean_group(1,dx))...
            - (mean_group(tx,4)-mean_group(1,4));
    end
end
delta_group


%%

      plotmat_values_rsfMRI(mNet_demean(:,:),[],[-1 1],0,[]);