123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149 |
- % VMHC.m
- % Note: requires Tools for NIfTI and ANALYZE image
- % (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
- % Requires DPABI toolbox (http://rfmri.org/dpabi) for cluster inference
- main_path = '';
- useGlobalSignal = false;
- %% CALCULATE HOMOTOPIC VOXEL-WISE CORRELATIONS:
- % Load Mask and Get Voxel-wies Coordinates:
- % Middle 2 voxels: 138-139
- % Struct: .0703 x .0703 mm
- % Funct: .1406 x .5 mm (2 x 7.12 struct voxels equivalent)
- mask = load_nii(fullfile(main_path, 'mirrored_groupwise_template.nii'));
- ind = find(mask.img(:)>0); [x,y,z] = ind2sub(size(mask.img),ind);
- ind = x<139; numvox = sum(ind);
- x = x(ind); y = y(ind); z = z(ind);
- r_x = zeros(size(x)); for i = 1:length(x); r_x(i) = 139+(138-x(i)); end
- nscan = 19;
- % Perform denoising:
- Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2;
- intercept = ones(175,1); homotopy_img = zeros(276,256,59,nscan); % voxel-wise homotopy data
- h_wait = waitbar(0,sprintf('Loading functional series %02g...',i));
- signals = cell(1,nscan); globalSignal = zeros(175,19);
- funct_path_outer = '';
- for i = 1:nscan
- % Check which physio data is usable:
- if all(physio_params{i}(:,2)==mode(physio_params{i}(:,2))) % iso is a constant, therefore unnecessary
- physio_use = physio_params{i}(:,1);
- else
- physio_use = physio_params{i};
- end
-
- % Load Functional:
- waitbar(0, h_wait,sprintf('Loading functional series %02g...',i));
- funct_path = fullfile(funct_path_outer,sprintf('funct_%02g',i));
- funct4D = zeros(276,256,59,175);
- for j = 6:180
- waitbar(j/180,h_wait);
- funct = load_nii(fullfile(funct_path,sprintf('s_R_mirror_funct3D_bird_%02g_time_%03g.nii',[i,j])));
- globalSignal(j-5, i) = nanmean(funct.img(mask.img(:)>0));
- funct4D(:,:,:,j-5) = funct.img;
- end
-
- % Calculate correlations:
- waitbar(0,h_wait,sprintf('Calculating homotopic correlations for series %02g...',i));
- signals{i}{1} = zeros(175,numvox); signals{i}{2} = zeros(175,numvox);
- for j = 1:numvox
- if rem(j,10000)==0; waitbar(j/numvox,h_wait); end
- signals{i}{1}(:,j) = squeeze(funct4D(x(j),y(j),z(j),:));
- signals{i}{2}(:,j) = squeeze(funct4D(r_x(j),y(j),z(j),:));
- if useGlobalSignal
- [~,~,signals{i}{1}(:,j)] = regress(signals{i}{1}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use,globalSignal(:,i)]);
- [~,~,signals{i}{2}(:,j)] = regress(signals{i}{2}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use,globalSignal(:,i)]);
- else
- [~,~,signals{i}{1}(:,j)] = regress(signals{i}{1}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use]);
- [~,~,signals{i}{2}(:,j)] = regress(signals{i}{2}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use]);
- end
- end
- [signals{i}{1},~,~] = bst_bandpass_filtfilt(signals{i}{1}',Fs,.008,.1,0,'iir');
- [signals{i}{2},~,~] = bst_bandpass_filtfilt(signals{i}{2}',Fs,.008,.1,0,'iir');
- signals{i}{1} = signals{i}{1}'; signals{i}{2} = signals{i}{2}';
- r = fast_corr(signals{i}{1},signals{i}{2});
- for j = 1:numvox
- homotopy_img(x(j),y(j),z(j),i) = r(j);
- homotopy_img(r_x(j),y(j),z(j),i) = r(j);
- end
- end; close(h_wait)
- %% t-statistics across subjects
- controlCovariates = true;
- covariates = [zscore(age), zscore(temps), zscore(iso)];
- % Calculate t-stats summarizing voxel-wise homotopy across birds/scans:
- homotopy_t = zeros(size(mask.img));
- age_t = homotopy_t;
- h_wait = waitbar(0,'Calculating t-statistics...');
- for i = 1:numvox
- if rem(i,1000)==0; waitbar(i/numvox, h_wait); end
- zvec = fisherz(squeeze(homotopy_img(x(i),y(i),z(i),:)));
- if controlCovariates
- lm = fitlm(covariates, zvec);
- homotopy_t(x(i),y(i),z(i)) = lm.Coefficients.tStat(1);
- homotopy_t(r_x(i),y(i),z(i)) = lm.Coefficients.tStat(1);
- age_t(x(i),y(i),z(i)) = lm.Coefficients.tStat(2);
- age_t(r_x(i),y(i),z(i)) = lm.Coefficients.tStat(2);
- else
- [~,~,~,stats] = ttest(zvec);
- homotopy_t(x(i),y(i),z(i)) = stats.tstat;
- homotopy_t(r_x(i),y(i),z(i)) = stats.tstat;
- end
- end; close(h_wait)
- % Save t-Stats Images:
- % Middle 2 voxels: 138-139
- img = mask; img.img = homotopy_t; img.img(137:140,:,:) = 0; % Discard 2 voxels on each side of midline
- if controlCovariates
- save_nii(img,fullfile(main_path,'t_stat_homotopy_controlled.nii'))
- img2 = mask; img2.img = age_t;
- save_nii(img2, fullfile(main_path,'t_stat_age_controlled.nii'))
- else
- save_nii(img, fullfile(main_path,'t_stat_homotopy.nii'))
- end
- img.img(135:142,:,:) = 0; % Discard 4 voxels on each side of midline
- if controlCovariates
- save_nii(img,fullfile(main_path,'t_stat_homotopy_8voxgap_controlled.nii'))
- img2.img(136:141,:,:) = 0;
- save_nii(img2, fullfile(main_path,'t_stat_age_8voxgap_controlled.nii'))
- else
- save_nii(img,fullfile(main_path,'t_stat_homotopy_8voxgap.nii'))
- end
- %% Run cluster inference: Homotopy
- % Extract the significant clusters of homotopy using cluster-extent method
- % on only 1 hemisphere (hemispheres are mirror images)
- if controlCovariates
- StatsImgFile = fullfile(main_path,'t_stat_homotopy_8voxgap_controlled.nii'); % more conservative midline deletion
- OutputName = fullfile(main_path,'T_homotopy_clusters_controlled.nii');
- else
- StatsImgFile = fullfile(main_path,'t_stat_homotopy_8voxgap.nii'); % more conservative midline deletion
- OutputName = fullfile(main_path,'T_homotopy_clusters.nii');
- end
- MaskFile = fullfile(main_path, 'mirrored_groupwise_template.nii');
- IsTwoTailed = false; VoxelPThreshold = .001; ClusterPThreshold = .05;
- Flag = 'T'; Df1 = 18; Df2 = []; VoxelSize = mask.hdr.dime.pixdim(2:4);
- Header = []; dLh = [];
- [Data_Corrected, ClusterSize, Header] = y_GRF_Threshold(StatsImgFile, ...
- VoxelPThreshold,IsTwoTailed,ClusterPThreshold,OutputName,MaskFile,...
- Flag,Df1,Df2,VoxelSize,Header, dLh);
- y_ClusterReport(Data_Corrected, Header, 18) % 18 is the cluster (voxel neighbor) connectivity criterion
- % Save bilateral version of cluster-thresholded t-stats:
- if controlCovariates
- clusterT = load_nii(fullfile(main_path,'ClusterThresholded_T_homotopy_clusters_controlled.nii'));
- clusterT.img(143:276,:,:) = clusterT.img(134:-1:1,:,:); % 135:142 voxels were discarded, so start from 143:(134+142)
- save_nii(clusterT,fullfile(main_path,'T_homotopy_clusters_bilat_controlled.nii'));
- else
- clusterT = load_nii(fullfile(main_path,'ClusterThresholded_T_homotopy_clusters.nii'));
- clusterT.img(143:276,:,:) = clusterT.img(134:-1:1,:,:); % 135:142 voxels were discarded, so start from 143:(134+142)
- save_nii(clusterT,fullfile(main_path,'T_homotopy_clusters_bilat.nii'));
- end
|