123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414 |
- % 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,[]);
|