Browse Source

Upload files to 'code'

Elliot Layden 7 months ago
parent
commit
cd696da300
1 changed files with 149 additions and 0 deletions
  1. 149 0
      code/VMHC.m

+ 149 - 0
code/VMHC.m

@@ -0,0 +1,149 @@
+% 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