123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280 |
- % homotopy_denoising.m
- % Within-subjects:
- % 1. Linear trends
- % 2. Quadratic trends
- % 3. 6 rigid-body motion parameters
- % 4. 5 principal components of CSF signals
- % 5. Body temperature
- % 6. Isoflurane dose
- % 7. Global signal regression (optional)
- % 8. Bandpass filtering
- % Between-subjects:
- % 1. Age
- % 2. Mean Body Temperature
- % 3. Mean Isoflurane Dose
- % 4. ROI pair Euclidean distance
- % 5. ROI volume
- % Note: requires Tools for NIfTI and ANALYZE image
- % (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
- % Use smoothed or unsmoothed data for ROI extraction?
- useSmoothedData = false;
- % Use homotopic-SC as well as homotopic-nSC?
- useSC = true;
- % Regress out global signal?
- useGlobalSignal = false;
- %% LOAD DATA
- % Collect file paths to normalized functional images:
- functPath = '';
- nscan = 19;
- timeSeries = cell(1,nscan); % 5 normals, 4 scans each (except normal #5, 3 scans)
- for scan = 1:nscan
- timeSeries{scan} = cell(1,175);
- for timepoint = 6:180 % skip first 5 volumes (T1 equilibration)
- if useSmoothedData
- timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
- sprintf('s_w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
- else
- timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
- sprintf('w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
- end
- end
- end
- % Load ROIs:
- ROI_path = '';
- SongSystem = fullfile(ROI_path,'SongSystemROIs.nii'); % n = 20
- mask = load_nii(fullfile(ROI_path,'new_temp_brain_mask.nii'));
- nROI = 20;
- % Get ROI Subscripts/Coordinates:
- rois = load_nii(SongSystem);
- roi_subscripts = cell(1,nROI);
- for i = 1:nROI
- [x, y, z] = ind2sub(size(rois.img),find(rois.img(:)==i));
- roi_subscripts{i} = [x,y,z];
- end
- % Get CSF Subscripts/Coordinates:
- CSF_mask = load_nii(fullfile(ROI_path, 'CSF_mask.nii'));
- [CSF_x,CSF_y,CSF_z] = ind2sub(size(CSF_mask.img), find(CSF_mask.img(:)>0));
- n_CSF = length(CSF_x);
- % Extract ROI Data & CSF Signals:
- ROI_dat = cell(1,nscan); CSF = cell(1,nscan); globalSignal = zeros(175,nscan); % Initialize
- hWait = waitbar(0,'Extracting functional data...');
- for scan = 1:nscan
- waitbar(scan/nscan, hWait);
-
- ROI_dat{scan} = zeros(175, nROI); CSF_dat = zeros(175, n_CSF); % Initialize
- for j = 1:175 % iterate through time (j)
-
- img = load_nii(timeSeries{scan}{j});
- globalSignal(j, scan) = mean(img.img(mask.img(:)>0)); % Average whole brain signal for time point j
-
- for k = 1:nROI % iterate through ROIs (k)
- nvox = size(roi_subscripts{k}, 1); % # voxels in ROI
- roi_holder = zeros(1, nvox);
- for l = 1:nvox % iterate through voxels of ROI
- roi_holder(l) = img.img(roi_subscripts{k}(l,1), roi_subscripts{k}(l,2), roi_subscripts{k}(l,3));
- end
- ROI_dat{scan}(j,k) = mean(roi_holder); % average voxels of ROI
- end
-
- for k = 1:n_CSF % iterate through CSF voxels (k)
- CSF_dat(j,k) = img.img(CSF_x(k), CSF_y(k), CSF_z(k));
- end
- end
-
- [~, PCs] = pca(CSF_dat); % pca across voxels of CSF mask
- CSF{scan} = PCs(:,1:5); % save PC's 1-5
-
- end; delete(hWait)
- % Load Motion Data:
- motionPath = '';
- motion_params = cell(1,nscan);
- for i = 1:nscan
- load(fullfile(motionPath, sprintf('motion_scan_%02g.mat',i)))
- motion_params{i} = motion;
- end
- % Load Physiology Data (Temperature (degrees Celsius), Isoflurane (%)):
- physio_path = '';
- physio_params = cell(1,nscan);
- physio_times = [60:60:540,576]; % 60 seconds to 576 seconds
- TR = 1.0671*3; scan_times = TR:TR:TR*175;
- for i = 1:nscan
- tbl = readtable(fullfile(physio_path,'All_Normal_Condition_Physiology_Records.xlsx'),'Sheet',i);
- pp_temp = spline(physio_times,tbl.Temperature); % piecewise polynomial of cubic spline
- pred_temp = ppval(pp_temp,scan_times); % evaluate PP at x
- pp_iso = spline(physio_times,tbl.Isoflurane); % piecewise polynomial of cubic spline
- pred_iso = ppval(pp_iso,scan_times); % evaluate PP at x
- physio_params{i} = [pred_temp',pred_iso'];
- end
- ages = [25, 45, 65, 89, 24, 45, 65, 90, 25, 46, 67, 89, 24, 46, 66, 90, 24, 45, 67];
- % Get ROI-level tSNR Data:
- tSNR_mat = zeros(nscan, nROI);
- for i = 1:nscan
- for j = 1:nROI
- mean_sig = mean(ROI_dat{i}(:,j));
- noise_SD = std(ROI_dat{i}(:,j));
- tSNR_mat(i,j) = mean_sig ./ noise_SD;
- end
- end
- % Get Mean Physiological Parameters:
- temps = zeros(19,1); iso = zeros(19,1);
- for i = 1:19
- temps(i) = mean(physio_params{i}(:,1));
- iso(i) = mean(physio_params{i}(:,2));
- end
- % Get Mean Volumes and Euclidean Distances:
- vox_dim = rois.hdr.dime.pixdim(2:4);
- vox_volume = prod(vox_dim);
- m_volumes = zeros(nROI); euclideans = zeros(nROI);
- % Get ROI Centroids:
- centroids = zeros(nROI,3);
- for i = 1:nROI
- [x,y,z] = ind2sub(size(rois.img),find(rois.img==i));
- centroids(i,:) = [mean(x),mean(y),mean(z)].*vox_dim;
- end
- for i = 1:nROI
- for j = 1:nROI
- m_volumes(i,j) = mean([sum(rois.img(:)==i)*vox_volume,sum(rois.img(:)==j)*vox_volume]);
- euclideans(i,j) = sqrt(sum((centroids(i,:)-centroids(j,:)).^2));
- end
- end
- %% RUN DENOISING
- % Labels for FC:
- labels = {'R Area X','L Area X','R LMAN','L LMAN','R RA','L RA',...
- 'R HVC','L HVC','R Aud. Forebrain','L Aud. Forebrain','R Field L',...
- 'L Field L','R Diencephalon','L Diencephalon','R Medulla','L Medulla',...
- 'R Uva','L Uva','R DM','L DM'};
- useScans = 1:19;
- if useSC
- nROI = 20;
- else
- nROI = 12;
- end
- % Preprocess ROI Data / Within-Subjects Nuisance Regression:
- Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2; intercept = ones(175,1);
- hWait = waitbar(0,'Denoising...');
- for scan = useScans
- if ishandle(hWait); waitbar(scan/nscan, hWait); end
-
- % Get nuisance regression covariates:
- if all(physio_params{scan}(:,2)==mode(physio_params{scan}(:,2))) % iso is a constant, so don't need
- if useGlobalSignal
- nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1), globalSignal(:,scan)];
- else
- nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1)];
- end
- else
- if useGlobalSignal
- nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}, globalSignal(:,scan)];
- else
- nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}];
- end
- end
-
- % Perform nuisance regression, retain residuals:
- for j = 1:nROI
- [~, ~, ROI_dat{scan}(:,j)] = regress(ROI_dat{scan}(:,j), nuisance);
- end
-
- % Perform bandpass filtering:
- [ROI_dat{scan},~,~] = bst_bandpass_filtfilt(ROI_dat{scan}', Fs, .008, .1, 0, 'iir');
- ROI_dat{scan} = ROI_dat{scan}';
- end
- if ishandle(hWait); delete(hWait); end
- % Calculate Correlation Matrix:
- zmats = zeros(nROI, nROI, nscan);
- for scan = useScans
- zmats(:,:,scan) = fisherz(corr(ROI_dat{scan}(:,1:nROI)));
- end
- % Between-Subjects Nuisance Regression:
- nscan = length(useScans);
- numconn = nROI*(nROI-1)/2*nscan;
- birds = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5];
- scan = zeros(numconn,1); bird = zeros(numconn,1);
- FC = zeros(numconn,1); age = zeros(numconn,1);
- connection_type = zeros(numconn,1); euclidean = zeros(numconn,1);
- m_volume = zeros(numconn,1); mean_temp = zeros(numconn,1);
- mean_iso = zeros(numconn,1); conn_label = cell(numconn,1); count = 0;
- for i = 1:nROI-1
- for j = i+1:nROI
- FC(count+1:count+nscan) = squeeze(zmats(i,j,:));
- conn_label(count+1:count+nscan) = {[labels{i},'_',labels{j}]};
- age(count+1:count+nscan) = ages(useScans);
- bird(count+1:count+nscan) = birds(useScans);
- scan(count+1:count+nscan) = useScans;
- euclidean(count+1:count+nscan) = repmat(euclideans(i,j),nscan,1);
- m_volume(count+1:count+nscan) = repmat(m_volumes(i,j),nscan,1);
- if ((mod(i,2) && ~mod(j,2)) || (~mod(i,2) && mod(j,2))) && ~(mod(i,2) && j==(i+1)) % 4: heterotopic (contralateral)
- connection_type(count+1:count+nscan) = 4*ones(nscan,1);
- elseif mod(i,2) && j==(i+1)
- if i>12 || j>12
- connection_type(count+1:count+nscan) = 3*ones(nscan,1); % 3: homotopic-SC
- else
- connection_type(count+1:count+nscan) = 2*ones(nscan,1); % 2: homotopic-nSC
- end
- elseif (mod(i,2) && mod(j,2)) || (~mod(i,2) && ~mod(j,2)) % 1: ipsilateral
- connection_type(count+1:count+nscan) = ones(nscan,1);
- end
- mean_temp(count+1:count+nscan) = temps(useScans);
- mean_iso(count+1:count+nscan) = iso(useScans);
- count = count+nscan;
- end
- end
- bird = categorical(bird); scan = categorical(scan);
- if useSC
- connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic-nSC','Homotopic-SC','Heterotopic'});
- else
- connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic','Homotopic','Heterotopic'});
- end
- % Mixed Model:
- FC_data = table(FC, bird, scan, conn_label, connection_type, age, mean_temp, ...
- mean_iso, euclidean, m_volume, 'VariableNames', {'FC', 'Bird', 'Scan', ...
- 'Connection', 'Type', 'Age', 'Temp', 'Iso', 'Euclidean', 'Volume'});
- % Select Between-subjects nuisance parameters for formula:
- FC_formula = 'FC ~ Type*Age + Type*Iso + Temp + Euclidean + Volume + (1|Bird)';
- % Run LME random intercept model:
- lme = fitlme(FC_data, FC_formula, 'StartMethod','random', 'verbose',true, 'CheckHessian',true);
- lme.display
- % Get random effect coefficients:
- [B, Bnames, stats] = randomEffects(lme) %#ok
- % Check residuals:
- figure; lme.plotResiduals
- % Get marginal (fixed-effects only) vs. conditional (w/ random) R2
- lme.Rsquared
- resid_marginal = residuals(lme, 'Conditional',false); % fixed-effects only
- resid_conditional = residuals(lme, 'Conditional',true); % fixed-effects only
- R2_marginal = 1 - sum(resid_marginal.^2) / lme.SST
- R2_conditional = 1 - sum(resid_conditional.^2) / lme.SST
|