-
Notifications
You must be signed in to change notification settings - Fork 6
Running PLS in MATLAB with cross sectional data
Gabriel A. Devenyi edited this page May 1, 2023
·
9 revisions
For this example, I will be using cross-sectional DBM data (ie: the overall, as opposed to the resampled, Jacobian), specifically the absolute Jacobian mnc file (as opposed to the relative one).
To run this code, you will need tools originally from: https://www.rotman-baycrest.on.ca/index.php?section=84
If you're working at the CIC,PLS is available as a module that you can load, which will load the tools on matlab start up.
for file in ../../../derivatives/MICe_build_model_take2/ASY20190409_processed/*/stats-volumes/*_lsq6_I_lsq6_lsq12_and_nlin_inverted_displ_log_det_abs_fwhm0.2.mnc; do cp $file minc-orientation-fix/; done
for file in minc-orientation-fix/*; do minc_modify_header -dinsert "xspace:direction_cosines=1,0,0" $file; minc_modify_header -dinsert "yspace:direction_cosines=0,1,0" $file; minc_modify_header -dinsert "zspace:direction_cosines=0,0,1" $file; done
for file in minc-orientation-fix/*; do mnc2nii $file; done
%% setting up
addpath /data/chamal/projects/stephanie/cross-sect-asyn-PFF-project/analysis/PLS/plscmd/ % point to path to directory where you will run the PLS
addpath /data/chamal/projects/stephanie/cross-sect-asyn-PFF-project/analysis/PLS/plscmd/ % point to path to directory with PLS tools
% point to path to directory with PLS tools
% unnecesary now, if you module load PLS at the CIC.
% addpath /data/chamal/projects/stephanie/cross-sect-asyn-PFF-project/analysis/PLS/nii_tools/
inpath = '/data/chamal/projects/stephanie/cross-sect-asyn-PFF-project/analysis/PLS/nifti-jacobians' % point to path to directory where you have the niftis that you will run the PLS on
%% MANUALLY IMPORT DATA
% make sure to import all the columns as 'text' (select 'apply to selection to do this for all the columns)
% import as 'column vectors'
%% creating variables
allsubjects = unique(ID); % get subject ids
naidx = strcmp('NA',abs_log_full_det_blur); % find subjects that don't have jac. file
nasubj = ID(naidx); % count the number of subj that don't have jac file
subjects = setxor(nasubj,allsubjects); % make array with only the subjects that have a jac file
groups = unique(group); %make array with the group names
%% identify the variables that you want to look at with respect to the brain measures
IHC = {'sex_1',...
'GFAP1MCleft', 'GFAP1MCright', 'GFAP1SSleft1', 'GFAP1SSleft2', 'GFAP1SSright1', 'GFAP1SSright2', 'GFAPhippleft', 'GFAPhippright', 'GFAPhypleft', 'GFAPhypright', 'GFAPmedulla', 'GFAPmedulla2', 'GFAPMRNleft', 'GFAPMRNright', 'GFAPNAcleft', 'GFAPNAcright', 'GFAPPAG', 'GFAPPHN', 'GFAPponsleft', 'GFAPponsleft2', 'GFAPponsright', 'GFAPponsright2', 'GFAPSNleft', 'GFAPSNright', 'GFAPstrleft', 'GFAPstrright', 'GFAPthleft', 'GFAPthleft2', 'GFAPthright', 'GFAPthright2'};
nbehav = numel(IHC); % gets the number of variables that you are interested in (those specified above)
%% load in the DBM brain mask
mask_nii = niftiinfo(fullfile(inpath,'ASY20190409-nlin-3_statsmap-corr_mask.nii'));
mask = niftiread(mask_nii);
mask_idx = find(mask);
%% setting up variables
ngroup = numel(groups); % get the number of groups
nsubj = numel(subjects); % get the number of subjects
nvoxel = length(mask_idx); % get the number of DBM measures (ie the number of voxels)
%% creating dataframe with brain measures and variables of interest
mri_data = cell(ngroup,1); % create array to add in mri data
IHC_data = cell(ngroup,1); % create array to add in variables of interest data (in this case IHC regions)
for gg = 1:numel(groups) % runs loop for each group
tempidx= find(strcmp(groups{gg},group) & naidx ==0); % find all subjects that have a jacobians file for group assignment gg ( ie group 1, 2, etc...)
mri_data{gg} = zeros(length(tempidx),nvoxel); % for each subj in group gg, add place for all of its voxels into the dataframe
IHC_data{gg} = zeros(length(tempidx),nbehav); % for each subj in group gg, add place for all of the behaviour into the dataframe
for ii = 1:length(tempidx) % runs for the number of subjects in group gg
for jj = 1:numel(IHC) % runs for the number of IHC variables
temp = eval(IHC{jj}); % choose IHC jj
IHC_data{gg}(ii,jj) = str2double(temp(tempidx(ii))); % add to the dataframe
end
mri_path = abs_log_full_det_blur(tempidx(ii));
%nii = load_nii(mri_path);
%mri_data{gg}(ii,:) = nii.img(mask_idx);
nii = niftiinfo(mri_path);
niidata = niftiread(nii);
mri_data{gg}(ii,:) = niidata(mask_idx); % add to the dataframe
fprintf('group %i subject %i / %i done\n',gg,ii,length(tempidx))
end
end
%% permutation testing
datamat{1} = [mri_data{1}; mri_data{2}]; % create new dataframe with brain measures
option.method = 3; % 1 = mean-centering (i.e. group/condition comaprison)
% 3 = behaviour (comparing 2 sets of variables)
option.num_perm = 1000;
option.num_boot = 1000;
option.stacked_IHCdata = [IHC_data{1}; IHC_data{2}];
option.stacked_IHCdata(isnan(option.stacked_IHCdata)) = 0;
% z score manually
x = option.stacked_IHCdata;
mu = mean(x);
sigma = std(x);
z = bsxfun(@minus,x,mu);
z = bsxfun(@rdivide,z,sigma);
option.stacked_IHCdata = z;
result = pls_analysis(datamat,14,1,option); % (dataframe, #subj,#condition, option)
%% checking results
% check pvalues
result.perm_result.sprob
% percent covariance
pct_cov = (result.s .^2) / sum(result.s .^2)
%% plot covariance vs pvalue for all LVs
figure;
plot(result.perm_result.sprob,pct_cov*100, '.', 'MarkerSize',15,'LineWidth',2)
xlabel('p-value','FontSize',15)
ylabel('Percent Variance Explained','FontSize',15)
xlim([0 0.05]) % only show LVs where pvalues is less than 0.05
%% plot the variable loadings for LV1 with confidence intervals
%% LV1
upper = result.boot_result.ulcorr(:,1) - result.lvcorrs(:,1);
lower = result.lvcorrs(:,1) - result.boot_result.llcorr(:,1);
figure;
width=0.4;
barh(result.lvcorrs(:,1)); hold on
set(gca,'TickLabelInterpreter','none')
set(gca,'Ytick',1:numel(IHC),'YTickLabel',IHC,'FontSize',17);
er = errorbar(result.lvcorrs(:,1),1:length(result.lvcorrs(:,1)),lower,upper,'horizontal');
er.LineStyle = 'none';
hold off
%% LV2
upper = result.boot_result.ulcorr(:,2) - result.lvcorrs(:,2);
lower = result.lvcorrs(:,2) - result.boot_result.llcorr(:,2);
figure;
width=0.4;
barh(result.lvcorrs(:,2)); hold on
set(gca,'TickLabelInterpreter','none')
set(gca,'Ytick',1:numel(IHC),'YTickLabel',IHC,'FontSize',17);
er = errorbar(result.lvcorrs(:,2),1:length(result.lvcorrs(:,2)),lower,upper,'horizontal');
er.LineStyle = 'none';
hold off
%% LV3
upper = result.boot_result.ulcorr(:,3) - result.lvcorrs(:,3);
lower = result.lvcorrs(:,3) - result.boot_result.llcorr(:,3);
figure;
width=0.4;
barh(result.lvcorrs(:,3)); hold on
set(gca,'TickLabelInterpreter','none')
set(gca,'Ytick',1:numel(IHC),'YTickLabel',IHC,'FontSize',17);
er = errorbar(result.lvcorrs(:,3),1:length(result.lvcorrs(:,3)),lower,upper,'horizontal');
er.LineStyle = 'none';
hold off
%% plot bootstrapping histograms
figure;
histogram(result.boot_result.compare_u(:,1));
figure;
histogram(result.boot_result.compare_u(:,2));
figure;
histogram(result.boot_result.compare_u(:,3));
%% get bootstrap ratio for each LV
% For brain variables - a bootstrap ratio is calculated as the ratio of the singular vector weight of a
% brain var (ie brain score) to its standard error across bootstraps and we usually look at the 95%
% confidence interval, or like 2 standard deviations away (i.e. the tail ends)
% for LV1
bsr = result.boot_result.compare_u(:,1);
bsr_volume = mask;
bsr_volume(mask_idx) = result.boot_result.compare_u(:,1);
niftiwrite(bsr_volume, '~/Documents/PLS/LV1-brain-IHC.nii', mask_nii);
% for LV2
bsr = result.boot_result.compare_u(:,2);
bsr_volume = mask;
bsr_volume(mask_idx) = result.boot_result.compare_u(:,2);
niftiwrite(bsr_volume, '/data/chamal/projects/elisa/neonate_mia/analysis/PLS/binned_calls/LV2.nii', mask_nii);
% for LV3
bsr = result.boot_result.compare_u(:,3);
bsr_volume = mask;
bsr_volume(mask_idx) = result.boot_result.compare_u(:,3);
niftiwrite(bsr_volume, '/data/chamal/projects/elisa/neonate_mia/analysis/PLS/binned_calls/LV3.nii', mask_nii);
%% get brain and IHC scores for each LV
brainsc_tp1_LV1_demo = datamat{1} * result.u(:,1);
IHCsc_tp1_LV1_demo = z * result.v(:,1);
brainsc_tp1_LV2_demo = datamat{1} * result.u(:,2);
IHCsc_tp1_LV2_demo = z * result.v(:,2);
brainsc_tp1_LV3_demo = datamat{1} * result.u(:,3);
IHCsc_tp1_LV3_demo = z * result.v(:,3);
%% look at how the subjects load onto brain-IHC mapping with respect to their group
%for LV1
figure;
scatter(brainsc_tp1_LV1_demo(1:29),IHCsc_tp1_LV1_demo(1:29),'m','filled'); hold on
scatter(brainsc_tp1_LV1_demo(30:58),IHCsc_tp1_LV1_demo(30:58),'g','filled'); hold on
scatter(brainsc_tp1_LV1_demo(59:end),IHCsc_tp1_LV2_demo(59:end),'b','filled');
refline;You can visualize your brain outputs with the niftis you create. As for the results from a mincLM/LMER, you can use:
Display secondlevel_template0.nii.gz LV2.nii LV2.nii
The process is described in detail here (https://github.com/CoBrALab/documentation/wiki/Display-mnc-files-from-the-terminal). Set the lower threshold for your "hot"/red/positive bar to 3 and the max threshold for your "cold"/blue/negative bar to -3.