homotopy_denoising.m 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286
  1. % homotopy_denoising.m
  2. % Within-subjects:
  3. % 1. Linear trends
  4. % 2. Quadratic trends
  5. % 3. 6 rigid-body motion parameters
  6. % 4. 5 principal components of CSF signals
  7. % 5. Body temperature
  8. % 6. Isoflurane dose
  9. % 7. Global signal regression (optional)
  10. % 8. Bandpass filtering
  11. % Between-subjects:
  12. % 1. Age
  13. % 2. Mean Body Temperature
  14. % 3. Mean Isoflurane Dose
  15. % 4. ROI pair Euclidean distance
  16. % 5. ROI volume
  17. % Note: requires Tools for NIfTI and ANALYZE image
  18. % (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
  19. % Use smoothed or unsmoothed data for ROI extraction?
  20. useSmoothedData = false;
  21. % Use homotopic-SC as well as homotopic-nSC?
  22. useSC = true;
  23. % Regress out global signal?
  24. useGlobalSignal = false;
  25. %% LOAD DATA
  26. % Collect file paths to normalized functional images:
  27. functPath = '';
  28. nscan = 19;
  29. timeSeries = cell(1,nscan); % 5 normals, 4 scans each (except normal #5, 3 scans)
  30. for scan = 1:nscan
  31. timeSeries{scan} = cell(1,175);
  32. for timepoint = 6:180 % skip first 5 volumes (T1 equilibration)
  33. if useSmoothedData
  34. timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
  35. sprintf('s_w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
  36. else
  37. timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
  38. sprintf('w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
  39. end
  40. end
  41. end
  42. % Load ROIs:
  43. ROI_path = '';
  44. SongSystem = fullfile(ROI_path,'SongSystemROIs.nii'); % n = 20
  45. mask = load_nii(fullfile(ROI_path,'new_temp_brain_mask.nii'));
  46. nROI = 20;
  47. % Get ROI Subscripts/Coordinates:
  48. rois = load_nii(SongSystem);
  49. roi_subscripts = cell(1,nROI);
  50. for i = 1:nROI
  51. [x, y, z] = ind2sub(size(rois.img),find(rois.img(:)==i));
  52. roi_subscripts{i} = [x,y,z];
  53. end
  54. % Get CSF Subscripts/Coordinates:
  55. CSF_mask = load_nii(fullfile(ROI_path, 'CSF_mask.nii'));
  56. [CSF_x,CSF_y,CSF_z] = ind2sub(size(CSF_mask.img), find(CSF_mask.img(:)>0));
  57. n_CSF = length(CSF_x);
  58. % Extract ROI Data & CSF Signals:
  59. ROI_dat = cell(1,nscan); CSF = cell(1,nscan); globalSignal = zeros(175,nscan); % Initialize
  60. hWait = waitbar(0,'Extracting functional data...');
  61. for scan = 1:nscan
  62. waitbar(scan/nscan, hWait);
  63. ROI_dat{scan} = zeros(175, nROI); CSF_dat = zeros(175, n_CSF); % Initialize
  64. for j = 1:175 % iterate through time (j)
  65. img = load_nii(timeSeries{scan}{j});
  66. globalSignal(j, scan) = mean(img.img(mask.img(:)>0)); % Average whole brain signal for time point j
  67. for k = 1:nROI % iterate through ROIs (k)
  68. nvox = size(roi_subscripts{k}, 1); % # voxels in ROI
  69. roi_holder = zeros(1, nvox);
  70. for l = 1:nvox % iterate through voxels of ROI
  71. roi_holder(l) = img.img(roi_subscripts{k}(l,1), roi_subscripts{k}(l,2), roi_subscripts{k}(l,3));
  72. end
  73. ROI_dat{scan}(j,k) = mean(roi_holder); % average voxels of ROI
  74. end
  75. for k = 1:n_CSF % iterate through CSF voxels (k)
  76. CSF_dat(j,k) = img.img(CSF_x(k), CSF_y(k), CSF_z(k));
  77. end
  78. end
  79. [~, PCs] = pca(CSF_dat); % pca across voxels of CSF mask
  80. CSF{scan} = PCs(:,1:5); % save PC's 1-5
  81. end; delete(hWait)
  82. % Load Motion Data:
  83. motionPath = '';
  84. motion_params = cell(1,nscan);
  85. for i = 1:nscan
  86. load(fullfile(motionPath, sprintf('motion_scan_%02g.mat',i)))
  87. motion_params{i} = motion;
  88. end
  89. % Load Physiology Data (Temperature (degrees Celsius), Isoflurane (%)):
  90. physio_path = '';
  91. physio_params = cell(1,nscan);
  92. physio_times = [60:60:540,576]; % 60 seconds to 576 seconds
  93. TR = 1.0671*3; scan_times = TR:TR:TR*175;
  94. for i = 1:nscan
  95. tbl = readtable(fullfile(physio_path,'All_Normal_Condition_Physiology_Records.xlsx'),'Sheet',i);
  96. pp_temp = spline(physio_times,tbl.Temperature); % piecewise polynomial of cubic spline
  97. pred_temp = ppval(pp_temp,scan_times); % evaluate PP at x
  98. pp_iso = spline(physio_times,tbl.Isoflurane); % piecewise polynomial of cubic spline
  99. pred_iso = ppval(pp_iso,scan_times); % evaluate PP at x
  100. physio_params{i} = [pred_temp',pred_iso'];
  101. end
  102. % Labels:
  103. SongSystem_labels = {'L Area X','R Area X','L LMAN','R LMAN','L RA','R RA',...
  104. 'L HVC','R HVC','L Aud. Forebrain','R Aud. Forebrain','L Field L',...
  105. 'R Field L','L Diencephalon','R Diencephalon','L Medulla','R Medulla',...
  106. 'L Uva','R Uva','L DM','R DM'};
  107. ages = [25, 45, 65, 89, 24, 45, 65, 90, 25, 46, 67, 89, 24, 46, 66, 90, 24, 45, 67];
  108. % Get ROI-level tSNR Data:
  109. tSNR_mat = zeros(nscan, nROI);
  110. for i = 1:nscan
  111. for j = 1:nROI
  112. mean_sig = mean(ROI_dat{i}(:,j));
  113. noise_SD = std(ROI_dat{i}(:,j));
  114. tSNR_mat(i,j) = mean_sig ./ noise_SD;
  115. end
  116. end
  117. % Get Mean Physiological Parameters:
  118. temps = zeros(19,1); iso = zeros(19,1);
  119. for i = 1:19
  120. temps(i) = mean(physio_params{i}(:,1));
  121. iso(i) = mean(physio_params{i}(:,2));
  122. end
  123. % Get Mean Volumes and Euclidean Distances:
  124. vox_dim = rois.hdr.dime.pixdim(2:4);
  125. vox_volume = prod(vox_dim);
  126. m_volumes = zeros(nROI); euclideans = zeros(nROI);
  127. % Get ROI Centroids:
  128. centroids = zeros(nROI,3);
  129. for i = 1:nROI
  130. [x,y,z] = ind2sub(size(rois.img),find(rois.img==i));
  131. centroids(i,:) = [mean(x),mean(y),mean(z)].*vox_dim;
  132. end
  133. for i = 1:nROI
  134. for j = 1:nROI
  135. m_volumes(i,j) = mean([sum(rois.img(:)==i)*vox_volume,sum(rois.img(:)==j)*vox_volume]);
  136. euclideans(i,j) = sqrt(sum((centroids(i,:)-centroids(j,:)).^2));
  137. end
  138. end
  139. %% RUN DENOISING
  140. % Labels for FC:
  141. labels = {'L_Area_X','R_Area_X','L_LMAN','R_LMAN','L_RA','R_RA',...
  142. 'L_HVC','R_HVC','L_Aud_Forebrain','R_Aud_Forebrain','L_Field_L',...
  143. 'R_Field_L','L_Diencephalon','R_Diencephalon','L_Medulla','R_Medulla',...
  144. 'L_Uva','R_Uva','L_DM','R_DM'};
  145. useScans = 1:19;
  146. if useSC
  147. nROI = 20;
  148. else
  149. nROI = 12;
  150. end
  151. % Preprocess ROI Data / Within-Subjects Nuisance Regression:
  152. Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2; intercept = ones(175,1);
  153. hWait = waitbar(0,'Denoising...');
  154. for scan = useScans
  155. if ishandle(hWait); waitbar(scan/nscan, hWait); end
  156. % Get nuisance regression covariates:
  157. if all(physio_params{scan}(:,2)==mode(physio_params{scan}(:,2))) % iso is a constant, so don't need
  158. if useGlobalSignal
  159. nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1), globalSignal(:,scan)];
  160. else
  161. nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1)];
  162. end
  163. else
  164. if useGlobalSignal
  165. nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}, globalSignal(:,scan)];
  166. else
  167. nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}];
  168. end
  169. end
  170. % Perform nuisance regression, retain residuals:
  171. for j = 1:nROI
  172. [~, ~, ROI_dat{scan}(:,j)] = regress(ROI_dat{scan}(:,j), nuisance);
  173. end
  174. % Perform bandpass filtering:
  175. [ROI_dat{scan},~,~] = bst_bandpass_filtfilt(ROI_dat{scan}', Fs, .008, .1, 0, 'iir');
  176. ROI_dat{scan} = ROI_dat{scan}';
  177. end
  178. if ishandle(hWait); delete(hWait); end
  179. % Calculate Correlation Matrix:
  180. zmats = zeros(nROI, nROI, nscan);
  181. for scan = useScans
  182. zmats(:,:,scan) = fisherz(corr(ROI_dat{scan}(:,1:nROI)));
  183. end
  184. % Between-Subjects Nuisance Regression:
  185. nscan = length(useScans);
  186. numconn = nROI*(nROI-1)/2*nscan;
  187. birds = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5];
  188. scan = zeros(numconn,1); bird = zeros(numconn,1);
  189. FC = zeros(numconn,1); age = zeros(numconn,1);
  190. connection_type = zeros(numconn,1); euclidean = zeros(numconn,1);
  191. m_volume = zeros(numconn,1); mean_temp = zeros(numconn,1);
  192. mean_iso = zeros(numconn,1); conn_label = cell(numconn,1); count = 0;
  193. for i = 1:nROI-1
  194. for j = i+1:nROI
  195. FC(count+1:count+nscan) = squeeze(zmats(i,j,:));
  196. conn_label(count+1:count+nscan) = {[labels{i},'_',labels{j}]};
  197. age(count+1:count+nscan) = ages(useScans);
  198. bird(count+1:count+nscan) = birds(useScans);
  199. scan(count+1:count+nscan) = useScans;
  200. euclidean(count+1:count+nscan) = repmat(euclideans(i,j),nscan,1);
  201. m_volume(count+1:count+nscan) = repmat(m_volumes(i,j),nscan,1);
  202. if ((mod(i,2) && ~mod(j,2)) || (~mod(i,2) && mod(j,2))) && ~(mod(i,2) && j==(i+1)) % 4: heterotopic (contralateral)
  203. connection_type(count+1:count+nscan) = 4*ones(nscan,1);
  204. elseif mod(i,2) && j==(i+1)
  205. if i>12 || j>12
  206. connection_type(count+1:count+nscan) = 3*ones(nscan,1); % 3: homotopic-SC
  207. else
  208. connection_type(count+1:count+nscan) = 2*ones(nscan,1); % 2: homotopic-nSC
  209. end
  210. elseif (mod(i,2) && mod(j,2)) || (~mod(i,2) && ~mod(j,2)) % 1: ipsilateral
  211. connection_type(count+1:count+nscan) = ones(nscan,1);
  212. end
  213. mean_temp(count+1:count+nscan) = temps(useScans);
  214. mean_iso(count+1:count+nscan) = iso(useScans);
  215. count = count+nscan;
  216. end
  217. end
  218. bird = categorical(bird); scan = categorical(scan);
  219. if useSC
  220. connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic-nSC','Homotopic-SC','Heterotopic'});
  221. else
  222. connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic','Homotopic','Heterotopic'});
  223. end
  224. % Mixed Model:
  225. FC_data = table(FC, bird, scan, conn_label, connection_type, age, mean_temp, ...
  226. mean_iso, euclidean, m_volume, 'VariableNames', {'FC', 'Bird', 'Scan', ...
  227. 'Connection', 'Type', 'Age', 'Temp', 'Iso', 'Euclidean', 'Volume'});
  228. % Select Between-subjects nuisance parameters for formula:
  229. FC_formula = 'FC ~ Type*Age + Type*Iso + Temp + Euclidean + Volume + (1|Bird)';
  230. % Run LME random intercept model:
  231. lme = fitlme(FC_data, FC_formula, 'StartMethod','random', 'verbose',true, 'CheckHessian',true);
  232. lme.display
  233. % Get random effect coefficients:
  234. [B, Bnames, stats] = randomEffects(lme) %#ok
  235. % Check residuals:
  236. figure; lme.plotResiduals
  237. % Get marginal (fixed-effects only) vs. conditional (w/ random) R2
  238. lme.Rsquared
  239. resid_marginal = residuals(lme, 'Conditional',false); % fixed-effects only
  240. resid_conditional = residuals(lme, 'Conditional',true); % fixed-effects only
  241. R2_marginal = 1 - sum(resid_marginal.^2) / lme.SST
  242. R2_conditional = 1 - sum(resid_conditional.^2) / lme.SST