123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407 |
- % function MRI_quality(BIDS_path)
- % Stefan J. Blaschke, Neurology Department, University Hospital Cologne
- % VERSION 1.0., 31.12.2022
- %
- % matlab script to examine preprocessed fMRI dataset
- % display of raw imaging
- % overlay with registered atlas
- % SNR and tSNR
- % display of tSNR with mask outlines
- % DVARS and FD
- % - Afyouni S. & Nichols T.E., Insights and inference for DVARS, 2017
- % http://www.biorxiv.org/content/early/2017/04/06/125021.1
- % necessity of GSR (Chen et al., 2022, Magn Reson Med)
- % header - set up path to BIDS folder
- head_path = 'G:\Projects\StrokePT_tDCS_rsfmri'; %
- rate_data = 1;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% run script
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- cd (head_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\nifti_utils-master\external\NIFTI")
- addpath(".\scripts\toolbox\FSLNets\")
- addpath(".\scripts\")
- %%
- % create data table to store features
- varnames = {'animal','ses','ID','MRI date','file path','SNR','tSNR','DVARS','bad points','IOR','MNS','physio','qualy'};
- vartype = repmat({'cellstr'},length(varnames),1);
- Tmaster = table('Size',[100 length(varnames)],'VariableTypes',vartype,'VariableNames',varnames);
- parent_path = fullfile(head_path,"subjects/"); % parent_path = BIDS MRI path
- cd (parent_path)
- px = dir('sub*');
- %%
- % loop through Data
- %%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%
- idx = 1;
- for sub_idx = 1:length(px) % start loop subjects
- animal_path = fullfile(parent_path,px(sub_idx).name);
- cd (animal_path)
- tx = dir('ses*');
-
- for ses_idx = 1:length(tx) % start loop sessions
- %%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%
- session_path = fullfile(animal_path,tx(ses_idx).name);
- cd (session_path)
- [path,session,~] = fileparts(session_path);
- [~,animal,~] = fileparts(path);
- Tmaster{idx,"animal"} = {animal};
- Tmaster{idx,"ses"} = {session};
- Mat(idx).animal = animal;
- Mat(idx).ses = session;
- fprintf('Run# %u: running dataset from %s, %s...',idx,animal,session)
-
- % grab images
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if exist(fullfile(session_path,"anat/"))==7
- % a) get structural
- cd (fullfile(session_path,"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);
-
- struct_parts = strsplit(struct_orig(length(struct_orig)).name,'_');
- Tmaster(idx,'ID') = struct_parts(2);
- Mat(idx).ID = struct_parts(2);
- Tmaster(idx,'MRI date') = struct_parts(1);
- Mat(idx).date = struct_parts(1);
-
- % b.1) get fMRI
- if exist(fullfile(session_path,"func/"))==7
- cd (fullfile(session_path,"func/"))
- fmri_orig = dir("*.1.nii.gz");
- if length(fmri_orig)>1
- % for lx = 1:length(fmri_orig)
- % sx = niftiinfo(fmri_orig(lx).name);
- % volumes(lx) = sx.ImageSize(4);
- % end
- image_size = [fmri_orig.bytes];
- % image_size(volumes<100) = 0;
- fmri_orig = fmri_orig(find(image_size==max([image_size])));
- end
-
- ANNO_orig = dir("*AnnoSplit_rsfMRI.nii.gz");
- if ~isempty(fmri_orig) && ~isempty(ANNO_orig)
- fmri_name = fmri_orig(length(fmri_orig)).name;
- fMRI_fullfile = fullfile(fmri_orig(length(fmri_orig)).folder,fmri_orig(length(fmri_orig)).name);
- ANNO_fullfile = fullfile(ANNO_orig(length(ANNO_orig)).folder,ANNO_orig(length(ANNO_orig)).name);
- atlas = load_untouch_nii(ANNO_fullfile);
- atlas_img = imrotate(atlas.img,-90);
- M = atlas_img;
- M(M>0) = 1;
- % check physio regression
- if exist("rawMonData","dir")==7
- cd ("rawMonData\")
- if ~isempty(dir('*I32'))
- Tmaster{idx,'physio'} = {'y'};
- else
- Tmaster{idx,'physio'} = {'n'};
- end
- else
- Tmaster{idx,'physio'} = {'n'};
- end
- % b.2) get Regression File
- cd (fullfile(session_path,"func/"))
- save_path = fullfile(session_path,"func/","results/")
- if ~exist(save_path)
- mkdir(save_path)
- end
- if exist('regr','dir')==7
- cd regr\
- SFRGR_orig = dir("*_SFRGR.nii.gz");
- TCfile = dir("MasksTCsSplit*.txt.mat");
- fmri_mask = dir("*thres_mask.nii.gz");
-
- if ~isempty(TCfile) && ~isempty(SFRGR_orig)
- SFRGR_file = char(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
- TC_fullfile = fullfile(TCfile.folder,TCfile.name);
- Tmaster(idx,'file path') = {TC_fullfile};
- % load and plot SFRGR file
- SFRGR = load_untouch_nii(SFRGR_file);
- SFRGR_img = SFRGR.img;
- SFRGR_mean = imrotate(squeeze(mean(SFRGR_img,4)),-90);
-
-
- % load and add boundaries of mask
- mask = load_untouch_nii(fullfile(fmri_mask.folder,fmri_mask.name));
- mask_img = imrotate(mask.img,-90);
- 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);
- % if GNI >3
- % need_GSR = "GSR";
- % else
- % need_GSR = "no GSR";
- % end
- %
- % flag{1} = need_GSR;
- % % % SNR (spatial, temporal)
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- info = niftiinfo(fMRI_fullfile);
- V1 = load_untouch_nii(fMRI_fullfile);
- V2 = V1.img;
- fmri_data = flip(flip(imrotate(V2,-90),3),2);
- [Ni, Nj, Nk, Nt] = size(fmri_data);
- func_2Dimg = reshape(fmri_data, 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);
-
- tSNR = squeeze(mean2(tSNR_3Dimg(mask_img==1)));
- SNR = squeeze(mean(func_mean3D(mask_img==1)))/std(func_stddev3D(mask_img==0));
- Tmaster{idx,'tSNR'} = {tSNR};
- Tmaster{idx,'SNR'} = {SNR};
- Mat(idx).tSNR = tSNR;
- Mat(idx).SNR = SNR;
- % cd ..\
- % niftiwrite(tSNR_image,strcat(fmri_name(1:end-10),'tSNR.nii'))
- SFRGR_img(mask_img==0)=nan;
- GSR = squeeze(mean(mean(mean(SFRGR_img,3,'omitnan'),2,'omitnan'),1,'omitnan'));
- % %
- % % % get DVARS
- fmri_BET = fmri_data;
- fmri_BET(mask_img==0) = nan;
- [X0,Y0,Z0,T0] = size(fmri_BET);
- I0 = prod([X0,Y0,Z0]);
- Y = reshape(fmri_BET,[I0,T0]);
- [DVARS,DVARS_Stat]=DVARSCalc(Y,'scale',1/100,'TransPower',1/3,'RDVARS','verbose',0);
- [V,DSE_Stat]=DSEvars(Y,'scale',1/100);
-
- Tmaster{idx,'DVARS'} = {mean(DVARS)};
- Dx = find(DVARS_Stat.pvals<0.05./(T0-1) & DVARS_Stat.DeltapDvar>5); %print corrupted DVARS data-points
-
- bad_points = Dx;%(Dx > 3);
- if isempty(bad_points)
- badTP = 0;
- else
- badTP = length(bad_points);
- 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) = [];
- Tmaster{idx,"bad points"} = {badTP*100/(info.ImageSize(4)-5)};
- Mat(idx).badvol = {badTP*100/(info.ImageSize(4)-5)};
- Mat(idx).DVARS = DVARS;
-
- % getFD
- cd (fullfile(session_path,"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,'VariableNamingRule','preserve');
- 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;
- Tmaster{idx,'FD'}=squeeze(mean(FD));
- Mat(idx).FD = FD;
-
- else
- FD = [];
- Mat(idx).FD = [];
- disp('Dataset has no available regression file')
- end
- else
- FD = [];
- Mat(idx).FD = [];
- disp('Dataset has no available regression file')
- end
- else
- FD = [];
- Mat(idx).FD = [];
- disp('Dataset has no available regression file')
- end
-
- [mNet,mNet_scrub,mNet_scrub_GSR,IOR] = grabMRI(TC_fullfile,GSR,scratch,rate_data);
- Tmaster{idx,'IOR'} = {IOR};
- Mat(idx).IOR = IOR;
- MNS = squeeze(mean(strengths_und(mNet)));
- Tmaster{idx,'MNS'} = {MNS/(size(mNet,1)-1)};
- Mat(idx).mNet = mNet;
- Mat(idx).mNet_scrub = mNet_scrub;
- Mat(idx).mNet_scrub_GSR = mNet_scrub_GSR;
- T_sub = Tmaster(idx,["animal","ses","tSNR","SNR","DVARS","bad points","IOR","MNS"]);
- %% plotting and rating data
- %%%%%%%%%%%%%%%%%%%%%%%%%%
- if rate_data == 1
- pause()
- saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_plots.png')))
- close
- structural_fn = char(gunzip(struct_fullfile));
- vc= load_untouch_nii(structural_fn);
- vc = imrotate(vc.img,-90);
- figure(2)
- montage(vc)
- clim([0 max(max(max(vc)))])
- pause()
- functional4D_fn = char(gunzip(fMRI_fullfile));
- vc= load_untouch_nii(functional4D_fn);
- vc = imrotate(squeeze(mean(vc.img,4)),-90);
- figure(3)
- montage(vc)
- clim([0 max(max(max(vc)))])
- pause()
- plotMRI(tSNR_3Dimg,mask_img,'hot')
- pause()
- plotMRI(SFRGR_mean,M)
- pause()
- figure('units','normalized','outerposition',[0 0 1 1])
- subplot(10,3,4:3:30)
- x = find(mNet~=0);
- [c,stat] = cdfplot(abs(mNet(x)));
- c.LineWidth = 2;
- c.Color = [0 0 1];
- hold on
- x = find(mNet_scrub_GSR~=0);
- [c,stat] = cdfplot(abs(mNet_scrub_GSR(x)));
- xlim([0 1])
- c.LineWidth = 2;
- c.Color = [1 0 0];
- legend({'raw Mat','GSR Mat'},'Location','southeast')
- title('')
- ylabel('cumulative distribution')
- xlabel('node strength')
- hold off
-
- subplot(10,3,5:3:30)
- mNetD = threshold_proportional(mNet,0.2);
- % nNetT = zeros(size(mNetD));
- % mNetT(mNetD~=0) = mNet_full(mNetD~=0)
- plotmat_values_rsfMRI(mNetD,[],[-(stat.mean*3) stat.mean*3],0,[]);
- 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);
-
- subplot(10,3,6:3:30)
- GSRx=GSR;
- GSRx(scratch)=nan;
- plot(GSRx,'r-','LineWidth',2)
- xlabel('timepoint')
- ylabel('value')
- title('globel signal')
- %Get the table in string form.
- TString = evalc('disp(T_sub)');
- % 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');
- % Output the table using the annotation command.
- annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
- 'FontName',FixedWidth,'Units','Normalized','Position',[0 0 1 1]);
- pause()
- saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_parameters.png')))
- % dialog to rate preprocessing
- qm_set = 0;
- while qm_set == 0
- qm = input("was registration successful? y/n:",'s');
- if isempty(qm) | length(qm)>1
- disp('please print y or n' )
- elseif qm ~= 'n' && qm ~= 'y' && qm ~= 'd' && qm ~= 'g' && qm ~= 'x'
- disp('please print y or n' )
- else
- qm_set = 1;
- end
- end
- close all
- end
-
- else
- qm = 0;
- end % if loop no Anno?
- else
- qm = 0;
- end % if loop no anat?
- else
- qm = 0;
- end % if loop no func?
- else
- qm = 0;
- end % if loop no TCfile?
- else
- qm = 0;
- end % if loop no regression?
- if ~exist('qm','var')
- qm = '';
- end
- Tmaster{idx,'qualy'} = {qm};
- Mat(idx).qualy = qm;
- %%
- clearvars -except head_path idx Mat Tmaster parent_path px tx rate_data ses_idx sub_idx animal_path
- idx = idx +1;
- end % end loop sessions
- end % end loop subjects
- %%
- cd (head_path)
- writetable(Tmaster,'PT_tDCS_quali_stat_rate20230201.xlsx')
- save("PT_tDCS_Mat_rate_20230201","Mat")
- %%
- fprintf('idx = %u\n',idx)
- fprintf('subject: %u\n',sub_idx)
- fprintf('session: %u\n',ses_idx)
- % 3.1.23 11:27: idx = 206, sub 46, ses 1
|