Scheduled service maintenance on November 22


On Friday, November 22, 2024, between 06:00 CET and 18:00 CET, GIN services will undergo planned maintenance. Extended service interruptions should be expected. We will try to keep downtimes to a minimum, but recommend that users avoid critical tasks, large data uploads, or DOI requests during this time.

We apologize for any inconvenience.

checkMRI.m 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407
  1. % function MRI_quality(BIDS_path)
  2. % Stefan J. Blaschke, Neurology Department, University Hospital Cologne
  3. % VERSION 1.0., 31.12.2022
  4. %
  5. % matlab script to examine preprocessed fMRI dataset
  6. % display of raw imaging
  7. % overlay with registered atlas
  8. % SNR and tSNR
  9. % display of tSNR with mask outlines
  10. % DVARS and FD
  11. % - Afyouni S. & Nichols T.E., Insights and inference for DVARS, 2017
  12. % http://www.biorxiv.org/content/early/2017/04/06/125021.1
  13. % necessity of GSR (Chen et al., 2022, Magn Reson Med)
  14. % header - set up path to BIDS folder
  15. head_path = 'G:\Projects\StrokePT_tDCS_rsfmri'; %
  16. rate_data = 1;
  17. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  18. %% run script
  19. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  20. cd (head_path)
  21. addpath(".\scripts\toolbox\DVARS-master\")
  22. addpath(".\scripts\toolbox\DVARS-master\Nifti_Util\")
  23. addpath(".\scripts\toolbox\BCT\")
  24. addpath(".\scripts\toolbox\spm12")
  25. addpath(".\scripts\toolbox\fMRI-Quality-Checker-master")
  26. addpath(".\scripts\toolbox\spm12\matlabbatch")
  27. addpath(".\scripts\toolbox\nifti_utils-master\")
  28. addpath(".\scripts\toolbox\nifti_utils-master\external\NIFTI")
  29. addpath(".\scripts\toolbox\FSLNets\")
  30. addpath(".\scripts\")
  31. %%
  32. % create data table to store features
  33. varnames = {'animal','ses','ID','MRI date','file path','SNR','tSNR','DVARS','bad points','IOR','MNS','physio','qualy'};
  34. vartype = repmat({'cellstr'},length(varnames),1);
  35. Tmaster = table('Size',[100 length(varnames)],'VariableTypes',vartype,'VariableNames',varnames);
  36. parent_path = fullfile(head_path,"subjects/"); % parent_path = BIDS MRI path
  37. cd (parent_path)
  38. px = dir('sub*');
  39. %%
  40. % loop through Data
  41. %%%%%%%%%%%%%%%%%%%%
  42. %%%%%%%%%%%%%%%%%%%%
  43. idx = 1;
  44. for sub_idx = 1:length(px) % start loop subjects
  45. animal_path = fullfile(parent_path,px(sub_idx).name);
  46. cd (animal_path)
  47. tx = dir('ses*');
  48. for ses_idx = 1:length(tx) % start loop sessions
  49. %%%%%%%%%%%%%%%%%%%%
  50. %%%%%%%%%%%%%%%%%%%%
  51. session_path = fullfile(animal_path,tx(ses_idx).name);
  52. cd (session_path)
  53. [path,session,~] = fileparts(session_path);
  54. [~,animal,~] = fileparts(path);
  55. Tmaster{idx,"animal"} = {animal};
  56. Tmaster{idx,"ses"} = {session};
  57. Mat(idx).animal = animal;
  58. Mat(idx).ses = session;
  59. fprintf('Run# %u: running dataset from %s, %s...',idx,animal,session)
  60. % grab images
  61. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  62. if exist(fullfile(session_path,"anat/"))==7
  63. % a) get structural
  64. cd (fullfile(session_path,"anat/"))
  65. struct_orig = dir("*.1.nii.gz");
  66. struct_name = struct_orig(length(struct_orig)).name;
  67. struct_fullfile = fullfile(struct_orig(length(struct_orig)).folder,struct_orig(length(struct_orig)).name);
  68. struct_parts = strsplit(struct_orig(length(struct_orig)).name,'_');
  69. Tmaster(idx,'ID') = struct_parts(2);
  70. Mat(idx).ID = struct_parts(2);
  71. Tmaster(idx,'MRI date') = struct_parts(1);
  72. Mat(idx).date = struct_parts(1);
  73. % b.1) get fMRI
  74. if exist(fullfile(session_path,"func/"))==7
  75. cd (fullfile(session_path,"func/"))
  76. fmri_orig = dir("*.1.nii.gz");
  77. if length(fmri_orig)>1
  78. % for lx = 1:length(fmri_orig)
  79. % sx = niftiinfo(fmri_orig(lx).name);
  80. % volumes(lx) = sx.ImageSize(4);
  81. % end
  82. image_size = [fmri_orig.bytes];
  83. % image_size(volumes<100) = 0;
  84. fmri_orig = fmri_orig(find(image_size==max([image_size])));
  85. end
  86. ANNO_orig = dir("*AnnoSplit_rsfMRI.nii.gz");
  87. if ~isempty(fmri_orig) && ~isempty(ANNO_orig)
  88. fmri_name = fmri_orig(length(fmri_orig)).name;
  89. fMRI_fullfile = fullfile(fmri_orig(length(fmri_orig)).folder,fmri_orig(length(fmri_orig)).name);
  90. ANNO_fullfile = fullfile(ANNO_orig(length(ANNO_orig)).folder,ANNO_orig(length(ANNO_orig)).name);
  91. atlas = load_untouch_nii(ANNO_fullfile);
  92. atlas_img = imrotate(atlas.img,-90);
  93. M = atlas_img;
  94. M(M>0) = 1;
  95. % check physio regression
  96. if exist("rawMonData","dir")==7
  97. cd ("rawMonData\")
  98. if ~isempty(dir('*I32'))
  99. Tmaster{idx,'physio'} = {'y'};
  100. else
  101. Tmaster{idx,'physio'} = {'n'};
  102. end
  103. else
  104. Tmaster{idx,'physio'} = {'n'};
  105. end
  106. % b.2) get Regression File
  107. cd (fullfile(session_path,"func/"))
  108. save_path = fullfile(session_path,"func/","results/")
  109. if ~exist(save_path)
  110. mkdir(save_path)
  111. end
  112. if exist('regr','dir')==7
  113. cd regr\
  114. SFRGR_orig = dir("*_SFRGR.nii.gz");
  115. TCfile = dir("MasksTCsSplit*.txt.mat");
  116. fmri_mask = dir("*thres_mask.nii.gz");
  117. if ~isempty(TCfile) && ~isempty(SFRGR_orig)
  118. SFRGR_file = char(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
  119. TC_fullfile = fullfile(TCfile.folder,TCfile.name);
  120. Tmaster(idx,'file path') = {TC_fullfile};
  121. % load and plot SFRGR file
  122. SFRGR = load_untouch_nii(SFRGR_file);
  123. SFRGR_img = SFRGR.img;
  124. SFRGR_mean = imrotate(squeeze(mean(SFRGR_img,4)),-90);
  125. % load and add boundaries of mask
  126. mask = load_untouch_nii(fullfile(fmri_mask.folder,fmri_mask.name));
  127. mask_img = imrotate(mask.img,-90);
  128. mask_img(mask_img>0) = 1;
  129. % % % Global signal regression
  130. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  131. % RestData = gunzip(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
  132. % WB_mask = gunzip(fullfile(fmri_mask.folder,fmri_mask.name));
  133. % GNI = Rest2GNI(char(RestData),char(WB_mask),100);
  134. % if GNI >3
  135. % need_GSR = "GSR";
  136. % else
  137. % need_GSR = "no GSR";
  138. % end
  139. %
  140. % flag{1} = need_GSR;
  141. % % % SNR (spatial, temporal)
  142. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  143. info = niftiinfo(fMRI_fullfile);
  144. V1 = load_untouch_nii(fMRI_fullfile);
  145. V2 = V1.img;
  146. fmri_data = flip(flip(imrotate(V2,-90),3),2);
  147. [Ni, Nj, Nk, Nt] = size(fmri_data);
  148. func_2Dimg = reshape(fmri_data, Ni*Nj*Nk, Nt);
  149. [Nv, Nt] = size(func_2Dimg);
  150. func_mean = mean(func_2Dimg, 2,'omitnan');
  151. func_mean3D = reshape(func_mean, Ni, Nj, Nk);
  152. func_stddev = std(double(func_2Dimg), 0, 2,'omitnan');
  153. func_stddev3D = reshape(func_stddev, Ni, Nj, Nk);
  154. tSNR_2Dimg = mean(func_2Dimg, 2,'omitnan')./std(double(func_2Dimg), 0, 2,'omitnan');
  155. tSNR_3Dimg = reshape(tSNR_2Dimg, Ni, Nj, Nk);
  156. tSNR = squeeze(mean2(tSNR_3Dimg(mask_img==1)));
  157. SNR = squeeze(mean(func_mean3D(mask_img==1)))/std(func_stddev3D(mask_img==0));
  158. Tmaster{idx,'tSNR'} = {tSNR};
  159. Tmaster{idx,'SNR'} = {SNR};
  160. Mat(idx).tSNR = tSNR;
  161. Mat(idx).SNR = SNR;
  162. % cd ..\
  163. % niftiwrite(tSNR_image,strcat(fmri_name(1:end-10),'tSNR.nii'))
  164. SFRGR_img(mask_img==0)=nan;
  165. GSR = squeeze(mean(mean(mean(SFRGR_img,3,'omitnan'),2,'omitnan'),1,'omitnan'));
  166. % %
  167. % % % get DVARS
  168. fmri_BET = fmri_data;
  169. fmri_BET(mask_img==0) = nan;
  170. [X0,Y0,Z0,T0] = size(fmri_BET);
  171. I0 = prod([X0,Y0,Z0]);
  172. Y = reshape(fmri_BET,[I0,T0]);
  173. [DVARS,DVARS_Stat]=DVARSCalc(Y,'scale',1/100,'TransPower',1/3,'RDVARS','verbose',0);
  174. [V,DSE_Stat]=DSEvars(Y,'scale',1/100);
  175. Tmaster{idx,'DVARS'} = {mean(DVARS)};
  176. Dx = find(DVARS_Stat.pvals<0.05./(T0-1) & DVARS_Stat.DeltapDvar>5); %print corrupted DVARS data-points
  177. bad_points = Dx;%(Dx > 3);
  178. if isempty(bad_points)
  179. badTP = 0;
  180. else
  181. badTP = length(bad_points);
  182. end
  183. scratch = zeros(1,4*length(bad_points));
  184. for ptx = 1:length(bad_points)
  185. scratch(((ptx-1)*4+1):((ptx-1)*4+4)) = bad_points(ptx)-2:bad_points(ptx)+1;
  186. end
  187. scratch = unique(scratch(scratch>0));
  188. scratch(scratch>info.ImageSize(4)-5) = [];
  189. Tmaster{idx,"bad points"} = {badTP*100/(info.ImageSize(4)-5)};
  190. Mat(idx).badvol = {badTP*100/(info.ImageSize(4)-5)};
  191. Mat(idx).DVARS = DVARS;
  192. % getFD
  193. cd (fullfile(session_path,"func/"))
  194. if exist("txtRegrPython","dir")
  195. cd txtRegrPython\
  196. txtRegr_files = dir("*mcf*");
  197. if size(txtRegr_files,1)>1
  198. txtRegr = txtRegr_files(floor(length(txtRegr_files)/2)).name;
  199. txtRegr_table = readtable(txtRegr,'VariableNamingRule','preserve');
  200. if ~isempty(txtRegr_table)
  201. Regression_parameters = txtRegr_table(:,5:10);
  202. MP = cell2mat(table2cell(Regression_parameters));
  203. FD_measures = calculateFD(MP, .5, 0.5);
  204. FD = FD_measures.FD;
  205. Tmaster{idx,'FD'}=squeeze(mean(FD));
  206. Mat(idx).FD = FD;
  207. else
  208. FD = [];
  209. Mat(idx).FD = [];
  210. disp('Dataset has no available regression file')
  211. end
  212. else
  213. FD = [];
  214. Mat(idx).FD = [];
  215. disp('Dataset has no available regression file')
  216. end
  217. else
  218. FD = [];
  219. Mat(idx).FD = [];
  220. disp('Dataset has no available regression file')
  221. end
  222. [mNet,mNet_scrub,mNet_scrub_GSR,IOR] = grabMRI(TC_fullfile,GSR,scratch,rate_data);
  223. Tmaster{idx,'IOR'} = {IOR};
  224. Mat(idx).IOR = IOR;
  225. MNS = squeeze(mean(strengths_und(mNet)));
  226. Tmaster{idx,'MNS'} = {MNS/(size(mNet,1)-1)};
  227. Mat(idx).mNet = mNet;
  228. Mat(idx).mNet_scrub = mNet_scrub;
  229. Mat(idx).mNet_scrub_GSR = mNet_scrub_GSR;
  230. T_sub = Tmaster(idx,["animal","ses","tSNR","SNR","DVARS","bad points","IOR","MNS"]);
  231. %% plotting and rating data
  232. %%%%%%%%%%%%%%%%%%%%%%%%%%
  233. if rate_data == 1
  234. pause()
  235. saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_plots.png')))
  236. close
  237. structural_fn = char(gunzip(struct_fullfile));
  238. vc= load_untouch_nii(structural_fn);
  239. vc = imrotate(vc.img,-90);
  240. figure(2)
  241. montage(vc)
  242. clim([0 max(max(max(vc)))])
  243. pause()
  244. functional4D_fn = char(gunzip(fMRI_fullfile));
  245. vc= load_untouch_nii(functional4D_fn);
  246. vc = imrotate(squeeze(mean(vc.img,4)),-90);
  247. figure(3)
  248. montage(vc)
  249. clim([0 max(max(max(vc)))])
  250. pause()
  251. plotMRI(tSNR_3Dimg,mask_img,'hot')
  252. pause()
  253. plotMRI(SFRGR_mean,M)
  254. pause()
  255. figure('units','normalized','outerposition',[0 0 1 1])
  256. subplot(10,3,4:3:30)
  257. x = find(mNet~=0);
  258. [c,stat] = cdfplot(abs(mNet(x)));
  259. c.LineWidth = 2;
  260. c.Color = [0 0 1];
  261. hold on
  262. x = find(mNet_scrub_GSR~=0);
  263. [c,stat] = cdfplot(abs(mNet_scrub_GSR(x)));
  264. xlim([0 1])
  265. c.LineWidth = 2;
  266. c.Color = [1 0 0];
  267. legend({'raw Mat','GSR Mat'},'Location','southeast')
  268. title('')
  269. ylabel('cumulative distribution')
  270. xlabel('node strength')
  271. hold off
  272. subplot(10,3,5:3:30)
  273. mNetD = threshold_proportional(mNet,0.2);
  274. % nNetT = zeros(size(mNetD));
  275. % mNetT(mNetD~=0) = mNet_full(mNetD~=0)
  276. plotmat_values_rsfMRI(mNetD,[],[-(stat.mean*3) stat.mean*3],0,[]);
  277. hold on
  278. axis square;
  279. plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
  280. plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
  281. subplot(10,3,6:3:30)
  282. GSRx=GSR;
  283. GSRx(scratch)=nan;
  284. plot(GSRx,'r-','LineWidth',2)
  285. xlabel('timepoint')
  286. ylabel('value')
  287. title('globel signal')
  288. %Get the table in string form.
  289. TString = evalc('disp(T_sub)');
  290. % Use TeX Markup for bold formatting and underscores.
  291. TString = strrep(TString,'<strong>','\bf');
  292. TString = strrep(TString,'</strong>','\rm');
  293. TString = strrep(TString,'_','\_');
  294. % Get a fixed-width font.
  295. FixedWidth = get(0,'FixedWidthFontName');
  296. % Output the table using the annotation command.
  297. annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
  298. 'FontName',FixedWidth,'Units','Normalized','Position',[0 0 1 1]);
  299. pause()
  300. saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_parameters.png')))
  301. % dialog to rate preprocessing
  302. qm_set = 0;
  303. while qm_set == 0
  304. qm = input("was registration successful? y/n:",'s');
  305. if isempty(qm) | length(qm)>1
  306. disp('please print y or n' )
  307. elseif qm ~= 'n' && qm ~= 'y' && qm ~= 'd' && qm ~= 'g' && qm ~= 'x'
  308. disp('please print y or n' )
  309. else
  310. qm_set = 1;
  311. end
  312. end
  313. close all
  314. end
  315. else
  316. qm = 0;
  317. end % if loop no Anno?
  318. else
  319. qm = 0;
  320. end % if loop no anat?
  321. else
  322. qm = 0;
  323. end % if loop no func?
  324. else
  325. qm = 0;
  326. end % if loop no TCfile?
  327. else
  328. qm = 0;
  329. end % if loop no regression?
  330. if ~exist('qm','var')
  331. qm = '';
  332. end
  333. Tmaster{idx,'qualy'} = {qm};
  334. Mat(idx).qualy = qm;
  335. %%
  336. clearvars -except head_path idx Mat Tmaster parent_path px tx rate_data ses_idx sub_idx animal_path
  337. idx = idx +1;
  338. end % end loop sessions
  339. end % end loop subjects
  340. %%
  341. cd (head_path)
  342. writetable(Tmaster,'PT_tDCS_quali_stat_rate20230201.xlsx')
  343. save("PT_tDCS_Mat_rate_20230201","Mat")
  344. %%
  345. fprintf('idx = %u\n',idx)
  346. fprintf('subject: %u\n',sub_idx)
  347. fprintf('session: %u\n',ses_idx)
  348. % 3.1.23 11:27: idx = 206, sub 46, ses 1