diff --git a/.DS_Store b/.DS_Store index 3726ac3..228502b 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6591355 --- /dev/null +++ b/.gitignore @@ -0,0 +1,13 @@ +# Folder view configuration files +.DS_Store +Desktop.ini + +# Matlab files for this project that are too big +*.mat + +# Plot files +*.jpg +*.tif +*.ai +*.svg +*.zip \ No newline at end of file diff --git a/adj_matrices/all_adj_corr_alpha_May25.mat b/adj_matrices/all_adj_corr_alpha_May25.mat index 8cc2f5c..b496517 100644 Binary files a/adj_matrices/all_adj_corr_alpha_May25.mat and b/adj_matrices/all_adj_corr_alpha_May25.mat differ diff --git a/adj_matrices/all_adj_corr_beta_May25.mat b/adj_matrices/all_adj_corr_beta_May25.mat index f663126..a1bcd67 100644 Binary files a/adj_matrices/all_adj_corr_beta_May25.mat and b/adj_matrices/all_adj_corr_beta_May25.mat differ diff --git a/adj_matrices/all_adj_corr_delta_May25.mat b/adj_matrices/all_adj_corr_delta_May25.mat index 329730e..f6faf3a 100644 Binary files a/adj_matrices/all_adj_corr_delta_May25.mat and b/adj_matrices/all_adj_corr_delta_May25.mat differ diff --git a/adj_matrices/all_adj_corr_gamma_May25.mat b/adj_matrices/all_adj_corr_gamma_May25.mat index 1766fc6..0b880be 100644 Binary files a/adj_matrices/all_adj_corr_gamma_May25.mat and b/adj_matrices/all_adj_corr_gamma_May25.mat differ diff --git a/adj_matrices/all_adj_corr_theta_May25.mat b/adj_matrices/all_adj_corr_theta_May25.mat index 66b55db..9e7cfb1 100644 Binary files a/adj_matrices/all_adj_corr_theta_May25.mat and b/adj_matrices/all_adj_corr_theta_May25.mat differ diff --git a/atlas_render_color.txt b/atlas_render_color.txt index 32d6ec6..069a2cf 100644 --- a/atlas_render_color.txt +++ b/atlas_render_color.txt @@ -1,73 +1,73 @@ -0.81472 0.43874 0.35166 -0.90579 0.38156 0.83083 -0.12699 0.76552 0.58526 -0.91338 0.7952 0.54972 -0.12699 0.76552 0.58526 -0.91338 0.7952 0.54972 -0.63236 0.18687 0.91719 -0.09754 0.48976 0.28584 -0.63236 0.18687 0.91719 -0.09754 0.48976 0.28584 -0.2785 0.44559 0.7572 -0.54688 0.64631 0.75373 -0.2785 0.44559 0.7572 -0.54688 0.64631 0.75373 -0.2785 0.44559 0.7572 -0.54688 0.64631 0.75373 -0.74313 0.54722 0.26297 -0.39223 0.13862 0.65408 -0.95751 0.70936 0.38045 -0.96489 0.75469 0.56782 +0.93573 0.79466 0.17686 +0.45789 0.57739 0.95738 +0.24048 0.44004 0.26532 +0.7639 0.25761 0.92458 +0.24048 0.44004 0.26532 +0.7639 0.25761 0.92458 +0.75933 0.75195 0.22377 +0.74065 0.22867 0.37356 +0.75933 0.75195 0.22377 +0.74065 0.22867 0.37356 +0.74369 0.064187 0.0875 +0.10592 0.76733 0.64012 +0.74369 0.064187 0.0875 +0.10592 0.76733 0.64012 +0.74369 0.064187 0.0875 +0.10592 0.76733 0.64012 +0.12019 0.997 0.45164 +0.52505 0.22417 0.22771 +0.68156 0.6712 0.18062 +0.46326 0.71521 0.045051 1 1 1 1 1 1 -0.95751 0.70936 0.38045 -0.96489 0.75469 0.56782 -0.15761 0.27603 0.075854 -0.97059 0.6797 0.05395 -0.15761 0.27603 0.075854 -0.97059 0.6797 0.05395 -0.95717 0.6551 0.5308 -0.48538 0.16261 0.77917 -0.80028 0.119 0.93401 -0.14189 0.49836 0.12991 -0.80028 0.119 0.93401 -0.14189 0.49836 0.12991 -0.80028 0.119 0.93401 -0.14189 0.49836 0.12991 -0.42176 0.95974 0.56882 -0.91574 0.34039 0.46939 -0.79221 0.58527 0.011902 -0.95949 0.22381 0.33712 -0.42176 0.95974 0.56882 -0.91574 0.34039 0.46939 -0.65574 0.75127 0.16218 -0.035712 0.2551 0.79428 -0.65574 0.75127 0.16218 -0.035712 0.2551 0.79428 -0.65574 0.75127 0.16218 -0.035712 0.2551 0.79428 -0.84913 0.50596 0.31122 -0.93399 0.69908 0.52853 -0.84913 0.50596 0.31122 -0.93399 0.69908 0.52853 -0.84913 0.50596 0.31122 -0.93399 0.69908 0.52853 -0.67874 0.8909 0.16565 -0.75774 0.95929 0.60198 -0.74313 0.54722 0.26297 -0.39223 0.13862 0.65408 -0.65548 0.14929 0.68921 -0.17119 0.25751 0.74815 -0.65548 0.14929 0.68921 -0.17119 0.25751 0.74815 -0.70605 0.84072 0.45054 -0.031833 0.25428 0.083821 -0.70605 0.84072 0.45054 -0.031833 0.25428 0.083821 -0.27692 0.81428 0.22898 -0.046171 0.24352 0.91334 -0.27692 0.81428 0.22898 -0.046171 0.24352 0.91334 +0.68156 0.6712 0.18062 +0.46326 0.71521 0.045051 +0.21216 0.64206 0.72317 +0.098519 0.41905 0.34744 +0.21216 0.64206 0.72317 +0.098519 0.41905 0.34744 +0.82357 0.39076 0.66062 +0.17501 0.81614 0.38387 +0.16357 0.31743 0.62735 +0.66599 0.81454 0.02165 +0.16357 0.31743 0.62735 +0.66599 0.81454 0.02165 +0.16357 0.31743 0.62735 +0.66599 0.81454 0.02165 +0.89439 0.78907 0.91057 +0.51656 0.85226 0.80056 +0.7027 0.50564 0.74585 +0.15359 0.63566 0.81311 +0.89439 0.78907 0.91057 +0.51656 0.85226 0.80056 +0.95346 0.95089 0.38331 +0.54088 0.44396 0.61728 +0.95346 0.95089 0.38331 +0.54088 0.44396 0.61728 +0.95346 0.95089 0.38331 +0.54088 0.44396 0.61728 +0.67973 0.060019 0.57549 +0.036563 0.86675 0.53005 +0.67973 0.060019 0.57549 +0.036563 0.86675 0.53005 +0.67973 0.060019 0.57549 +0.036563 0.86675 0.53005 +0.8092 0.63119 0.27507 +0.74862 0.35507 0.24863 +0.12019 0.997 0.45164 +0.52505 0.22417 0.22771 +0.32583 0.65245 0.80445 +0.54645 0.60499 0.9861 +0.32583 0.65245 0.80445 +0.54645 0.60499 0.9861 +0.39888 0.38725 0.029992 +0.41509 0.14219 0.53566 +0.39888 0.38725 0.029992 +0.41509 0.14219 0.53566 +0.18074 0.025135 0.087077 +0.25539 0.42111 0.80209 +0.18074 0.025135 0.087077 +0.25539 0.42111 0.80209 1 1 1 1 1 1 1 1 1 @@ -76,15 +76,15 @@ 1 1 1 1 1 1 1 1 1 -1 1 1 -1 1 1 -0.097132 0.92926 0.15238 -0.82346 0.34998 0.82582 -0.097132 0.92926 0.15238 -0.82346 0.34998 0.82582 -0.69483 0.1966 0.53834 -0.3171 0.25108 0.99613 -0.69483 0.1966 0.53834 -0.3171 0.25108 0.99613 -0.95022 0.61604 0.078176 -0.034446 0.47329 0.44268 +0.020536 0.1841 0.98914 +0.92368 0.72578 0.066946 +0.020536 0.1841 0.98914 +0.92368 0.72578 0.066946 +0.020536 0.1841 0.98914 +0.92368 0.72578 0.066946 +0.6537 0.37036 0.9394 +0.93261 0.84156 0.018178 +0.6537 0.37036 0.9394 +0.93261 0.84156 0.018178 +0.16351 0.73423 0.68384 +0.9211 0.57103 0.78374 diff --git a/atlas_synthesis.m~ b/atlas_synthesis.m~ deleted file mode 100644 index dc789d5..0000000 --- a/atlas_synthesis.m~ +++ /dev/null @@ -1,1285 +0,0 @@ -%% Set up workspace -clear all - -% set up path - change iEEG_atlas_path to where you download the repository -iEEG_atlas_path = '/Users/jbernabei/Documents/PhD_Research/atlas_project/iEEG_atlas_dev'; - -% other paths may stay the same -metadata = readtable('data/atlas_metadata_final.xlsx'); -custom_atlas = readtable('data/custom_atlas.xlsx'); - -% anatomical atlas information -region_list = zeros(1,90); % for the 90 AAL regions we will be using -region_names = cell(1,90); -fi = fopen("localization/AAL116_WM.txt"); -for j = 1:90 - label = split(fgetl(fi)); - region_list(j) = str2double(label{3}); - region_names{j} = label{2}; -end - -% extract metadata -age_onset = metadata{:,11}; -age_surgery = metadata{:,12}; -HUP_outcome = floor(metadata{:,3}); - -% extract outcome -for s = 1:length(age_surgery) - engel_scores = metadata{s,3:5}; - engel_scores = rmmissing(engel_scores); - early_outcome(s) = floor(engel_scores(1)); % early outcome is 6 months - late_outcome(s) = floor(engel_scores(end)); % late outcome is latest time point up to 24 months -end - -therapy_field = metadata{:,6}; -implant_field = metadata{:,7}; -target_field = metadata{:,8}; -laterality_field = metadata{:,9}; -lesion_field = metadata{:,10}; -gender_field = metadata{:,13}; - -lesion_field = strcmp(lesion_field,'Lesional'); % 1 for lesional -therapy_field = strcmp(therapy_field,'Ablation'); % 1 ffor ablation - -load('data/MNI_atlas_new.mat') -load('data/HUP_atlas_May25.mat') - -all_pts = [Patient; (patient_no+106)]; % combine patient numbers for MNI & HUP - -HUP_outcome_all = zeros(length(soz_ch),1); -for i = 1:60 - HUP_outcome_all(find(patient_no==i)) = late_outcome(i); % we use outcome as late outcome -end - -clear NodesLeft -clear NodesLeftInflated -clear NodesRightInflated -clear NodesRight -clear NodesRegionLeft -clear NodesRegionRight -clear Data_N2 -clear Data_N3 -clear Data_R - -% some colors for plotting, etc -color1 = [0, 0.4470, 0.7410]; -color2 = [0.6350, 0.0780, 0.1840]; -color3 = [255, 248, 209]./255; -color4 = [103 55 155]/255; -color6 = [78 172 91]/255; - -% assign wake clips from MNI and HUP into a final data matrix -all_wake_data = [Data_W,wake_clip]; - -%% perform localization based on AAL116 atlas (which has an internal WM mask also) -all_coords = [ChannelPosition;mni_coords]; - -try [coords, all_roi, NN_ind] = nifti_values(all_coords,'localization/AAL116_WM.nii'); -catch ME - for i = 1:size(all_coords,1) - if mod(i,100) - fprintf('%d\n',i) - end - try [mni_coords1, mni_roi, NN_ind] = nifti_values(all_coords(i,:),'localization/AAL116_WM.nii'); - catch ME - mni_roi = NaN; - end - all_roi(i,1) = mni_roi; - end -end - -clear NN_ind -clear coords - -% assign AAL regions into new, aggregate regions of custom atlas. (116 -> 40 regions) -region_matrix = custom_atlas{:,2:4}; -for c = 1:length(all_roi) - old_ind = find(region_list==all_roi(c)); - try [row,col,V] = find(region_matrix==old_ind); - try new_roi(c,1) = row; - catch anyerror - new_roi(c,1) = 0; - end - catch anyerror - new_roi(c,1) = 0; - end - -end -%% define normal MNI/HUP and abnormal channels (SOZ/high-spike/RZ) -% set threshold on spike counts -spike_thresh = 24; % this is empirical, 1 spike/hour - -% find indices for which spikes are greater than threshold -spike_ind = [spike_24h>spike_thresh]; - -% define all abnormal channels -abnormal_ch = find([resected_ch+spike_ind+soz_ch]>0)+1772; - -% define all seizure onset indices -soz_ch_inds = find(soz_ch)+1772; - -% define normal HUP channels -normal_HUP_ch = find([resected_ch+spike_ind+soz_ch]==0)+1772; - -% define normal MNI channels -normal_MNI_ch = [1:1772]'; - -% define all normal channels -all_normal_ch = [normal_MNI_ch;normal_HUP_ch]; - -% define exclusive irritative zone -EIZ_ch = find(spike_ind - spike_ind.*soz_ch)+1772; - -% split EIZ into good vs poor -EIZ_good_ch = find((spike_ind - spike_ind.*soz_ch).*(HUP_outcome_all==1))+1772; -EIZ_poor_ch = find((spike_ind - spike_ind.*soz_ch).*(HUP_outcome_all>1))+1772; - -% split SOZ into good vs poor -soz_good_ch = find(soz_ch.*(HUP_outcome_all==1))+1772; -soz_poor_ch = find(soz_ch.*(HUP_outcome_all>1))+1772; - -% define non-EIZ non-SOZ poor outcome channels -poor_out_ch = find([resected_ch+spike_ind+soz_ch+(HUP_outcome_all==1)]==0)+1772; - -%% count HUP and MNI in original atlas and then new atlas -mni_roi_old = all_roi(1:1772); -hup_roi_old = all_roi(1773:end); - -for i = 1:45 - roi_count(i,1) = sum(mni_roi_old==region_list(2*i-1))+sum(mni_roi_old==region_list(2*i)); - roi_count(i,2) = sum(all_roi(normal_HUP_ch)==region_list(2*i-1))+sum(all_roi(normal_HUP_ch)==region_list(2*i)); - this_roi_name = split(string(region_names{2*i}),'_R'); % extract roi name without laterality - split_roi_name_old{i,1} = this_roi_name(1); -end - -for i = 1:20 - roi_count_new_normal(i,1) = length(intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); - roi_count_new_normal(i,2) = length(intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); - this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality - split_roi_name{i,1} = this_roi_name(1); -end - -new_counts = [roi_count_new_normal(:,1),(roi_count_new_normal(:,2))]; - -% Supplemental figure -figure(1);clf; -subplot(1,2,1) -barh(roi_count) -set(gca,'TickLabelInterpreter', 'none'); -xlabel('Number of nodes per region') -title('AAL regions') -yticks([1:45]) -yticklabels(split_roi_name_old) -subplot(1,2,2) -barh(new_counts) -set(gca,'TickLabelInterpreter', 'none'); -title('Composite regions') -yticks([1:21]) -yticklabels(split_roi_name) -xlabel('Number of nodes per region') - -%% do renderings of original and new atlas -random_colors = rand([40,3]); - -color_matrix = ones(117,3); - -for i = 1:40 - this_region = custom_atlas{i,2:4}; - for r = 1:length(rmmissing(this_region)) - color_matrix(this_region(r),:) = random_colors(i,:); - end -end - -% writes a file that can be used in BrainNet viewer to visualize -dlmwrite('atlas_render_color.txt',color_matrix(1:90,:),'delimiter',' ','precision',5) - -figure(1);clf; -imagesc(reshape(random_colors,40,1,3)) -set(gca,'TickLabelInterpreter', 'none'); -yticks([1:40]) -yticklabels(custom_atlas{:,1}) - -%% calculate power spectrum -all_wake_data = [Data_W,wake_clip]; - -[pxx,f] = pwelch(all_wake_data,200,100,[0.5:0.5:80],200); -pxx_norm = pxx./sum(pxx); - -%% Analyze power spectral density between HUP/MNI, Normal/irritative zone/SOZ - -freq_inds = [1, 8; % 0.5 - 4 Hz - 9, 16; % 4 - 8 Hz - 17, 26; % 8 - 13 Hz - 27, 60; % 13 - 30 Hz - 61, 160]; % 30 - 80 Hz - - -figure(1);clf % supplemental -figure(2);clf % supplemental -for i = 1:20 - % extract indices that correspond to a given ROI for HUP and MNI data - % (combining data from L and R hemispheres) - roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_EIZ = intersect([EIZ_good_ch;EIZ_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_soz = intersect([soz_good_ch;soz_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - - % median power spectral density across channels - median_mni = median(pxx_norm(:,roi_mni),2); - median_hup = median(pxx_norm(:,roi_hup),2); - median_normal = median(pxx_norm(:,[roi_mni;roi_hup]),2); - median_EIZ = median(pxx_norm(:,roi_EIZ),2); - median_soz = median(pxx_norm(:,roi_soz),2); - - for j = 1:5 - HUP_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_hup))'; - MNI_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_mni))'; - normal_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],[roi_mni;roi_hup]))'; - EIZ_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_EIZ))'; - soz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_soz))'; - - [p,h] = ranksum(HUP_feat,MNI_feat); % HUP vs MNI - pval_matrix_norm(i,j) = p; - - [p,h] = ranksum(normal_feat,EIZ_feat); % normal vs EIZ - pval_matrix_EIZ(i,j) = p; - - [p,h] = ranksum(normal_feat,soz_feat); % normal vs SOZ - pval_matrix_soz(i,j) = p; - - [p,h] = ranksum(EIZ_feat,soz_feat); % EIZ vs SOZ - pval_matrix_EIZ_vs_soz(i,j) = p; - - end - - % find upper and lower bounds (here 25th and 75th percentiles) - prct_25_75_mni = prctile(pxx_norm(:,roi_mni)',[25,75]); - prct_25_75_hup = prctile(pxx_norm(:,roi_hup)',[25,75]); - prct_25_75_normal = prctile(pxx_norm(:,[roi_mni;roi_hup])',[25,75]); - prct_25_75_EIZ = prctile(pxx_norm(:,roi_EIZ)',[25,75]); - prct_25_75_soz = prctile(pxx_norm(:,roi_soz)',[25,75]); - - x_axis = [0.5:0.5:80]; - - figure(1) - % 3 x 7 plot for 21 regions (aggregating across hemispheres) - subplot(4,5,i) - hold on - plot([4,4],[0,0.4],'k-') - plot([8,8],[0,0.4],'k-') - plot([13,13],[0,0.4],'k-') - plot([30,30],[0,0.4],'k-') - plot(x_axis,median_mni) - patch([x_axis fliplr(x_axis)], [prct_25_75_mni(1,:) fliplr(prct_25_75_mni(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') - plot(x_axis,median_hup) - patch([x_axis fliplr(x_axis)], [prct_25_75_hup(1,:) fliplr(prct_25_75_hup(2,:))], color6, 'FaceAlpha',0.5, 'EdgeColor','none') - - % add greek letters for frequency bands - txt_d = '\delta'; - txt_t = '\theta'; - txt_a = '\alpha'; - txt_b = '\beta'; - txt_g = '\gamma'; - txt_sig = '*'; - text(2,0.15,txt_d); - text(5.5,0.15,txt_t); - text(9.5,0.15,txt_a); - text(18,0.15,txt_b); - text(40,0.15,txt_g); - if pval_matrix_norm(i,1)<(0.05/100) - text(2,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_norm(i,2)<(0.05/100) - text(5.5,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_norm(i,3)<(0.05/100) - text(9.5,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_norm(i,4)<(0.05/100) - text(18,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_norm(i,5)<(0.05/100) - text(40,0.12,txt_sig,'FontSize', 20); - end - - % set axes, labels, subgraph titles - set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) - ylabel('Spectral Density') - xlabel('Frequency (Hz)') - this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality - split_roi_name{i,1} = this_roi_name(1); - title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title - xlim([0.5,80]) % xlimit 0.5-80 Hz - ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) - hold off - - figure(2) - subplot(4,5,i) - hold on - plot([4,4],[0,0.4],'k-') - plot([8,8],[0,0.4],'k-') - plot([13,13],[0,0.4],'k-') - plot([30,30],[0,0.4],'k-') - plot(x_axis,median_normal) - patch([x_axis fliplr(x_axis)], [prct_25_75_normal(1,:) fliplr(prct_25_75_normal(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') - plot(x_axis,median_EIZ) - patch([x_axis fliplr(x_axis)], [prct_25_75_EIZ(1,:) fliplr(prct_25_75_EIZ(2,:))], color6, 'FaceAlpha',0.5, 'EdgeColor','none') - try plot(x_axis,median_soz) - patch([x_axis fliplr(x_axis)], [prct_25_75_soz(1,:) fliplr(prct_25_75_soz(2,:))], color2, 'FaceAlpha',0.5, 'EdgeColor','none') - catch anyerror - end - - % add greek letters for frequency bands - text(2,0.15,txt_d); - text(5.5,0.15,txt_t); - text(9.5,0.15,txt_a); - text(18,0.15,txt_b); - text(40,0.15,txt_g); - if pval_matrix_EIZ(i,1)<(0.05/100) - text(2,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_EIZ(i,2)<(0.05/100) - text(5.5,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_EIZ(i,3)<(0.05/100) - text(9.5,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_EIZ(i,4)<(0.05/100) - text(18,0.12,txt_sig,'FontSize', 20); - end - if pval_matrix_EIZ(i,5)<(0.05/100) - text(40,0.12,txt_sig,'FontSize', 20); - end - - if pval_matrix_soz(i,1)<(0.05/100) - text(2,0.17,txt_sig,'FontSize', 20); - end - if pval_matrix_soz(i,2)<(0.05/100) - text(5.5,0.17,txt_sig,'FontSize', 20); - end - if pval_matrix_soz(i,3)<(0.05/100) - text(9.5,0.17,txt_sig,'FontSize', 20); - end - if pval_matrix_soz(i,4)<(0.05/100) - text(18,0.17,txt_sig,'FontSize', 20); - end - if pval_matrix_soz(i,5)<(0.05/100) - text(40,0.17,txt_sig,'FontSize', 20); - end - - % set axes, labels, subgraph titles - set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) - ylabel('Spectral Density') - xlabel('Frequency (Hz)') - this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality - split_roi_name{i,1} = this_roi_name(1); - title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title - xlim([0.5,80]) % xlimit 0.5-80 Hz - ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) - hold off - -end - -%% calculate univariate z-scores of neural activity -% power spectrum and entropy abnormality - -all_f = {'delta','theta','alpha','beta','gamma'}; % frequency bands -all_m = {'bandpower','entropy'}; % metrics - -c = 0; % c for column in z score matrix -for m = 1:2 % m for different metrics - for f = 1:5 % f for frequency bands - c = c+1; - - % first create the univariate atlas, taking only normal channels - % the univariate atlas is mean/median + stdev of a given feature - % across ROI - [mean_vec,median_vec,std_vec,feat_vals] = create_univariate_atlas(all_wake_data(:,find(all_normal_ch)),[1:40],new_roi(all_normal_ch),sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); - - % we call this function again on all of the data, but are only - % using it to extract all of the feature values on each channel - [~,~,~,feat_vals] = create_univariate_atlas(all_wake_data,[1:40],new_roi,sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); - - % now we test the univariate atlas to generate z scores - [zscores, pt_zscores, pt_ind, roi_ind] = test_univariate_atlas(median_vec, std_vec, feat_vals, all_pts, [1:40], new_roi,'patient'); - - % store the z scores - univariate_zscores(:,c) = zscores; - all_pt_zscores(:,c) = pt_zscores; - - end -end -%save('univariate_zscores.mat','univariate_zscores') - -%% mass univariate testing -% 21 ROI x 5 = 105 tests to correct for -freq_inds = [1, 8; % 0.5 - 4 Hz - 9, 16; % 4 - 8 Hz - 17, 26; % 8 - 13 Hz - 27, 60; % 13 - 30 Hz - 61, 160]; % 30 - 80 Hz - -for i = 1:20 - for j = 1:5 - roi_soz = intersect((find(soz_ch)+1772),[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_eiz = intersect((EIZ_ch),[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - soz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_soz))'; - eiz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_eiz))'; - HUP_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_hup))'; - MNI_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_mni))'; - [p,h] = ranksum(HUP_feat,MNI_feat); % HUP vs MNI - pval_matrix_norm(i,j) = p; - [p1,h] = ranksum([HUP_feat;MNI_feat],soz_feat); - pval_matrix_soz(i,j) = p1; - [p2,h] = ranksum([HUP_feat;MNI_feat],eiz_feat); - pval_matrix_eiz(i,j) = p2; - end -end - -figure(1);clf; -subplot(1,3,1) -imagesc((pval_matrix_norm.*100)<0.05) -title('HUP vs. MNI normal channels') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) -subplot(1,3,2) -imagesc((pval_matrix_eiz.*100)<0.05) -title('Irritative zone vs. composite atlas') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) -subplot(1,3,3) -imagesc((pval_matrix_soz.*100)<0.05) -title('Seizure onset zone vs. composite atlas') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) - -%% wavelet entropy validation across all regions -for i = 1:60 - for j = 1:size(all_wake_data,2) - start_inds = (i-1)*200+1; - end_inds = 200*i; - all_wentropy(i,j) = wentropy(all_wake_data((start_inds:end_inds),j),'shannon'); - end -end - -all_mean_wentropy = log(-1*median(all_wentropy)); - -r1 = all_mean_wentropy([normal_HUP_ch;normal_MNI_ch]); -r2 = all_mean_wentropy(poor_out_ch); -r3 = all_mean_wentropy(EIZ_good_ch); -r4 = all_mean_wentropy(EIZ_poor_ch); -r5 = all_mean_wentropy(soz_good_ch); -r6 = all_mean_wentropy(soz_poor_ch); - -ranksum(all_mean_wentropy(normal_HUP_ch),all_mean_wentropy(normal_MNI_ch)) - -p1 = ranksum(r1,r3) -p2 = ranksum(r3,r5) -p3 = ranksum(r1,r5) -p4 = ranksum(r2,r4) -p5 = ranksum(r4,r6) -p6 = ranksum(r2,r6) -p7 = ranksum(r1,r2) -p8 = ranksum(r3,r4) -p9 = ranksum(r5,r6) - -figure(1);clf; -hold on -scatter(ones(length(r1),1),r1,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([0.75 1.25],[median(r1) median(r1)],'k-','LineWidth',2) -scatter(2*ones(length(r2),1),r2,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([1.75 2.25],[median(r2) median(r2)],'k-','LineWidth',2) -scatter(3*ones(length(r3),1),r3,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([2.75 3.25],[median(r3) median(r3)],'k-','LineWidth',2) -scatter(4*ones(length(r4),1),r4,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([3.75 4.25],[median(r4) median(r4)],'k-','LineWidth',2) -scatter(5*ones(length(r5),1),r5,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([4.75 5.25],[median(r5) median(r5)],'k-','LineWidth',2) -scatter(6*ones(length(r6),1),r6,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([5.75 6.25],[median(r6) median(r6)],'k-','LineWidth',2) -xlim([0.5, 6.5]) -ylim([2,22]) -ylabel('Wavelet entropy') -xticks([1:6]) -xticklabels({'Engel 1, non-SOZ, non-EIZ','Engel 2+, non-SOZ, non-EIZ', 'Engel 1 EIZ','Engel 2+ EIZ','Engel 1 SOZ','Engel 2+ SOZ'}) -hold off - -%% -figure(1);clf; -for i = 1:20 - % extract indices that correspond to a given ROI for HUP and MNI data - % (combining data from L and R hemispheres) - roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_EIZ = intersect([EIZ_good_ch;EIZ_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_soz = intersect([soz_good_ch;soz_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - - entropy_normal = all_mean_wentropy([roi_mni;roi_hup]); - entropy_EIZ = all_mean_wentropy(roi_EIZ); - entropy_soz = all_mean_wentropy(roi_soz); - - pval_EIZ(i,1) = ranksum(entropy_normal,entropy_EIZ); - pval_soz(i,1) = ranksum(entropy_normal,entropy_soz); - pval_EIZ_vs_soz(i,1) = ranksum(entropy_EIZ,entropy_soz); - - % median power spectral density across channels - %entropy_mni = median(pxx_norm(:,roi_mni),2); - %median_hup = median(pxx_norm(:,roi_hup),2); - plot_cell{1,1} = entropy_normal; - plot_cell{1,2} = entropy_EIZ; - plot_cell{1,3} = entropy_soz; - - subplot(4,5,i) - hold on - violin(plot_cell,'xlabel',{'','',''},'facecolor',[color1;color6;color2],'mc',[],'medc','k');%,'edgecolor','k'); - legend('off') - txt_sig1 = '+'; - txt_sig2 = '*'; - txt_sig3 = '#'; - ylim([7 23]) - - if pval_EIZ(i,1) < (0.05./20) - plot([1,2],[20,20],'k-','LineWidth',2) - text(1.45,20.25,txt_sig2,'FontSize', 20); - end - if pval_soz(i,1) < (0.05./20) - plot([1,3],[21.75,21.75],'k-','LineWidth',2) - text(2,22,txt_sig2,'FontSize', 20); - end - if pval_EIZ_vs_soz(i,1) < (0.05./20) - plot([2,3],[20,20],'k-','LineWidth',2) - text(2.5,20.25,txt_sig2,'FontSize', 20); - end - - ylabel('Probability') - xlabel('- Log Entropy') - this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality - split_roi_name{i,1} = this_roi_name(1); - title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title - - hold off - -end - - - -%% create atlas data - > edge-wise and patient-wise should be options -% for HUP data, aggregate electrodes to ignore -% these should be poor-outcome/RZ/SOZ/>1 spike per hour - -conn_band = {'delta','theta','alpha','beta','gamma'}; -conn_type = {'corr','coh'}; - -% loop through types of connectivity and frequency bands -for c = 1:2 - for f = 1:5 - - % loop through subjects - for s = 1:max(all_pts) - s - - % extract subject data from all the data (columns are channels) - pt_data = sleep_clip;%all_wake_data(:,[all_pts==s]); - - % call functional connectivity code (200 is Hz, 1 is window size in seconds) - [pt_adj] = compute_FC(pt_data, 200, 1, conn_type{c},conn_band{f}); - - % assign median adjacency matrix into cell structure - all_pt_adj{s} = pt_adj; - - end - - % save adjacency matrices - save(sprintf('adj_matrices/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f}),'all_pt_adj') - - end -end - -%% assign atlas locations per patient into cell structure -for s = 1:max(all_pts) - pt_loc{s} = new_roi(all_pts==s)'; -end - -% indices of all abnormal channels -all_abnormal_ch = zeros(5203,1); -all_abnormal_ch(abnormal_ch) = 1; - -for s = 1:166 - - % check abnormal channels for HUP patients - if s>106 - - % abnormal channel indices - all_abn{s} = find(all_abnormal_ch(all_pts==s)); - if isempty(all_abn{s}) - all_abn{s} = []; - end - else - % if MNI patient, no abnormal channels by definition - all_abn{s} = []; - end -end - -%% functional cconnectivity atlas construction & z-score code - -conn_band = {'delta','theta','alpha','beta','gamma'}; -conn_type = {'corr','coh'}; - -% which feature -feat = 0; - -% loop through connectivity and frequency -for c = 1:2 - for f = 1:5 - - feat = feat+1; - a = 0 - - load(sprintf('adj_matrices/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f})); - - if c==2 - for s = 1:166 - all_pt_adj{s} = log(all_pt_adj{s}); - end - end - - [conn_edge, std_edge, samples_edge, sem_edge, raw_atlas_edge] = create_atlas_by_edge(all_pt_adj, pt_loc, all_abn, [1:40]', 2); - - % call for MNI patients - [conn_edge_mni, std_edge_mni, samples_edge_mni, ~, raw_atlas_edge_mni] = create_atlas_by_edge(all_pt_adj(1:106), pt_loc(1:106), all_abn(1:106), [1:40]', 1); - - % call for HUP patients - [conn_edge_hup, std_edge_hup, samples_edge_hup, ~, raw_atlas_edge_hup] = create_atlas_by_edge(all_pt_adj(107:166), pt_loc(107:166), all_abn(107:166), [1:40]', 1); - - % loop through all patients - for s = 1:166 % change this back - a = a+1; - - % do nodal space - try [native_adj_scores, corr_val] = test_native_adj(all_pt_adj{s}, pt_loc{s}, conn_edge, std_edge, [1:40]'); - - bivariate_native(feat).subj(a).data = native_adj_scores; - - catch anyerror - bivariate_native(feat).subj(a).data = NaN(length(pt_loc{s})); - end - - end - end - -end - -%% correlate connectivity between MNI and HUP -conn_edge_all = conn_edge_mni(:)+conn_edge_hup(:); -conn_edge_mni(isnan(conn_edge_all)) = []; -conn_edge_hup(isnan(conn_edge_all)) = []; - -std_edge_all = std_edge_mni(:)+std_edge_hup(:); -std_edge_mni(isnan(std_edge_all)) = []; -std_edge_hup(isnan(std_edge_all)) = []; - -figure(1);clf; -plot(exp(conn_edge_mni),exp(conn_edge_hup),'ko') -[r,p] = corr(exp(conn_edge_mni)',exp(conn_edge_hup)') -xlabel('MNI delta coherence') -ylabel('HUP delta coherence') -title('Correlation between atlas edges, r = 0.43, p = 4.6e-19') -xlim([0 0.1]) -ylim([0 0.1]) -%% process bivariate edge features into nodal features - -for f = 1:10 - this_feat = []; - abs_feat = []; - for s = 1:166 - % take 80th percentile across all edges of each node - abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),80)']; - end - - % assign into feature matrix - abs_bivariate_feats(:,f) = abs_feat; -end - -% for the single bivariate feature we want the maximum absolute Z score -% across all 10 bivariate features -single_bivariate_feat = max(abs_bivariate_feats')'; - -% assign previously calculated univariate z scores into bivariate scores -all_feat_zscores = [univariate_zscores, abs_bivariate_feats]; - -% incomplete channels are ones which are unlocalized in new atlas -incomplete_channels = find(isnan(sum(all_feat_zscores'))); -single_bivariate_feat(incomplete_channels) = []; - -%% -single_bivariate_feat_2 = single_bivariate_feat; -new_roi_2 = new_roi; -new_roi_2(incomplete_channels) = []; - -new_roi_2(new_roi_2==17) = 18; - -for ch = 1:4702 - if multi_class(ch)==1 - grouping{ch,1} = 'uninvolved'; - class_ind(ch,1) = 1; - elseif multi_class(ch)==2 - grouping{ch,1} = 'EIZ'; - class_ind(ch,1) = 2; - else - grouping{ch,1} = 'SOZ'; - class_ind(ch,1) = 3; - end -end - -univariate_feats = all_feat_zscores(:,1:10); -univariate_feats(incomplete_channels,:) = []; - -single_univariate_feat = prctile(abs(univariate_feats'),100)'; - -single_univariate_feat(isinf(single_univariate_feat)) = nanmedian(single_univariate_feat); -single_bivariate_feat_2(isinf(single_bivariate_feat_2)) = nanmedian(single_bivariate_feat_2); - -single_univariate_feat = single_univariate_feat(new_roi_2==18); -single_bivariate_feat_2 = single_bivariate_feat_2(new_roi_2==18); -grouping = grouping(new_roi_2==18); -class_ind = class_ind(new_roi_2==18); - - -figure(1);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat_2,'Group',grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,grouping,'orientation','horizontal',... - 'label',{'','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat_2,grouping,'orientation','horizontal',... - 'label', {'','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) - -d1 = computeCohen_d(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) -p1 = ranksum(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) -d2 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) -p2 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) -d3 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) -p3 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) - -d4 = computeCohen_d(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) -p4 = ranksum(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) -d5 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) -p5 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) -d6 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) -p6 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) - -[d1 d2 d3 d4 d5 d6] -[p1 p2 p3 p4 p5 p6] -%% do node abormality (max univariate+bivariate in plot) for following: -% 1. within/outside soz for lesional / non-lesional - -all_lesion = zeros(5203,1); -for i = 1:166 - if i>106 - all_lesion(all_pts==i) = lesion_field(i-106); - end -end - -within_soz_lesional = find(all_lesion.*soz_ch_all); -within_soz_nonlesional = find([all_lesion==0].*soz_ch_all); -outside_soz_lesional = find(all_lesion.*[soz_ch_all==0]); -outside_soz_nonlesional = find([all_lesion==0].*[soz_ch_all==0]); - -for ch = 1:5203 - if ismember(ch,within_soz_lesional) - lesional_grouping{ch,1} = 'within_soz_lesional'; - elseif ismember(ch,within_soz_nonlesional) - lesional_grouping{ch,1} = 'within_soz_nonlesional'; - elseif ismember(ch,outside_soz_lesional) - lesional_grouping{ch,1} = 'outside_soz_lesional'; - elseif ismember(ch,outside_soz_nonlesional) - lesional_grouping{ch,1} = 'outside_soz_nonlesional'; - end -end - -univariate_feats = all_feat_zscores(:,1:10); -single_univariate_feat = prctile(abs(univariate_feats'),100)'; -single_bivariate_feat = prctile(abs_bivariate_feats',100)'; - -uni_within_soz_lesional = single_univariate_feat(within_soz_lesional); -uni_within_soz_nonlesional = single_univariate_feat(within_soz_nonlesional); -uni_outside_soz_lesional = single_univariate_feat(outside_soz_lesional); -uni_outside_soz_nonlesional = single_univariate_feat(outside_soz_nonlesional); -bi_within_soz_lesional = single_bivariate_feat(within_soz_lesional); -bi_within_soz_nonlesional = single_bivariate_feat(within_soz_nonlesional); -bi_outside_soz_lesional = single_bivariate_feat(outside_soz_lesional); -bi_outside_soz_nonlesional = single_bivariate_feat(outside_soz_nonlesional); - -p1 = ranksum(uni_within_soz_lesional,uni_within_soz_nonlesional); -p2 = ranksum(uni_within_soz_lesional,uni_outside_soz_lesional); -p3 = ranksum(uni_within_soz_nonlesional,uni_outside_soz_nonlesional); -p4 = ranksum(uni_outside_soz_lesional,uni_outside_soz_nonlesional); -p5 = ranksum(bi_within_soz_lesional,bi_within_soz_nonlesional); -p6 = ranksum(bi_within_soz_lesional,bi_outside_soz_lesional); -p7 = ranksum(bi_within_soz_nonlesional,bi_outside_soz_nonlesional); -p8 = ranksum(bi_outside_soz_lesional,bi_outside_soz_nonlesional); - -single_univariate_feat(incomplete_channels) = []; -single_bivariate_feat(incomplete_channels) = []; -lesional_grouping(incomplete_channels) = []; - -figure(2);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',lesional_grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,lesional_grouping,'orientation','horizontal',... - 'label',{'','','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat,lesional_grouping,'orientation','horizontal',... - 'label', {'','','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) - -[p1 p2 p3 p4 p5 p6 p7 p8] -%% 2. within/outside ablation zone for good / poor outcome -all_outcome = ones(5203,1); -all_abl = NaN(5203,1); -for i = 1:166 - if i>106 - all_outcome(all_pts==i) = late_outcome(i-106); - if therapy_field(i-106) - all_abl(all_pts==i) = resect_ch_all(all_pts==i); - end - end -end - -within_abl_good = find([all_outcome==1].*[all_abl==1]); -outside_abl_good = find([all_outcome==1].*[all_abl==0]); -within_abl_poor = find([all_outcome>1].*[all_abl==1]); -outside_abl_poor = find([all_outcome>1].*[all_abl==0]); - -incomplete_channels2 = incomplete_channels; - -abl_grouping = cell(5203,1); - -for ch = 1:5203 - if ismember(ch,within_abl_good) - abl_grouping{ch,1} = 'within_abl_good'; - elseif ismember(ch,outside_abl_good) - abl_grouping{ch,1} = 'outside_abl_good'; - elseif ismember(ch,within_abl_poor) - abl_grouping{ch,1} = 'within_abl_poor'; - elseif ismember(ch,outside_abl_poor) - abl_grouping{ch,1} = 'outside_abl_poor'; - else - incomplete_channels2 = [incomplete_channels2, ch]; - end -end - -incomplete_channels2 = unique(incomplete_channels2); - -univariate_feats = all_feat_zscores(:,1:10); -single_univariate_feat = prctile(abs(univariate_feats'),100)'; -single_bivariate_feat = prctile(abs_bivariate_feats',100)'; - -uni_within_abl_good = single_univariate_feat(within_abl_good); -uni_outside_abl_good = single_univariate_feat(outside_abl_good); -uni_within_abl_poor = single_univariate_feat(within_abl_poor); -uni_outside_abl_poor = single_univariate_feat(outside_abl_poor); -bi_within_abl_good = single_bivariate_feat(within_abl_good); -bi_outside_abl_good = single_bivariate_feat(outside_abl_good); -bi_within_abl_poor = single_bivariate_feat(within_abl_poor); -bi_outside_abl_poor = single_bivariate_feat(outside_abl_poor); - -p1 = ranksum(uni_within_abl_good,uni_outside_abl_good); -p2 = ranksum(uni_within_abl_poor,uni_outside_abl_poor); -p3 = ranksum(uni_within_abl_good,uni_within_abl_poor); -p4 = ranksum(uni_outside_abl_good,uni_outside_abl_poor); -p5 = ranksum(bi_within_abl_good,bi_outside_abl_good); -p6 = ranksum(bi_within_abl_poor,bi_outside_abl_poor); -p7 = ranksum(bi_within_abl_good,bi_within_abl_poor); -p8 = ranksum(bi_outside_abl_good,bi_outside_abl_poor); - -single_univariate_feat(incomplete_channels2) = []; -single_bivariate_feat(incomplete_channels2) = []; -abl_grouping(incomplete_channels2) = []; - -figure(3);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',abl_grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,abl_grouping,'orientation','horizontal',... - 'label',{'','','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat,abl_grouping,'orientation','horizontal',... - 'label', {'','','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) - -[p1 p2 p3 p4 p5 p6 p7 p8] - -%% combine univariate and bivariate -%all_feat_zscores = [univariate_zscores, all_bivariate_feats]; -all_feat_zscores(isinf(all_feat_zscores)) = 0; -incomplete_channels = find(isnan(sum(all_feat_zscores'))); -all_feat_zscores2 = all_feat_zscores; -all_feat_zscores2(incomplete_channels,:) = []; -EIZ = [zeros(1772,1);[spike_ind - spike_ind.*soz_ch]]; -resect_ch_all = [zeros(1772,1);resected_ch]; - -% 1: patient, 2: abnormal, 3: EIZ, 4: SOZ, 5: resect, 6:outcome -all_possible_labels = [all_pts,[EIZ+[zeros(1772,1);soz_ch]],EIZ,[zeros(1772,1);soz_ch],resect_ch_all,[ones(1772,1);HUP_outcome_all]]; -multi_class_all = [~all_possible_labels(:,2)+2.*all_possible_labels(:,3)+3.*all_possible_labels(:,4)]; -all_possible_labels(incomplete_channels,:) = []; - -multi_class = [~all_possible_labels(:,2)+2.*all_possible_labels(:,3)+3.*all_possible_labels(:,4)]; - -% get only non EIZ normal versus SOZ (in good outcome patients) -non_normal_non_soz = [zeros(1772,1);[soz_ch==0]]; -non_normal_non_soz(incomplete_channels) = []; -eliminate_channels = find(all_possible_labels(:,3)==1);%find(non_normal_non_soz==1); -all_feat_zscores3 = all_feat_zscores2; -all_feat_zscores3(eliminate_channels,:) = []; -new_labels = all_possible_labels; -new_labels(eliminate_channels,:) = []; - -all_feats = all_feat_zscores2; -all_labels = new_labels(:,4); -all_pred_1 = zeros(size(all_labels)); -all_pred_2 = zeros(size(all_labels)); -all_pred_3 = zeros(size(all_labels)); - -[part] = make_xval_partition(length(all_labels), 10); -for p = 1:10 - p - X_train = all_feats(find(part~=p),:); - Y_train = all_labels(find(part~=p)); - X_test = all_feats(find(part==p),:); - Y_test = all_labels(find(part==p)); - - Mdl1 = TreeBagger(100,X_train(:,1:10),Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); - - Mdl2 = TreeBagger(100,X_train(:,11:20),Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); - - Mdl3 = TreeBagger(100,X_train,Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); - - [Yfit,pihat1] = predict(Mdl1,X_test(:,1:10)); - [Yfit,pihat2] = predict(Mdl2,X_test(:,11:20)); - [Yfit,pihat3] = predict(Mdl3,X_test); - - [~, Y_pred] = max(pihat1,[],2); - all_pred_1(find(part==p)) = pihat1(:,2); - all_pred_2(find(part==p)) = pihat2(:,2); - all_pred_3(find(part==p)) = pihat3(:,2); - - xfold_acc(p) = sum(Y_pred==Y_test)./length(Y_test); - - Y_pred = round(pihat1(:,2)); - - -end - - [X1,Y1,T,AUC1] = perfcurve(all_labels,all_pred_1,1) - - [X2,Y2,T,AUC2] = perfcurve(all_labels,all_pred_2,1) - - [X3,Y3,T,AUC3] = perfcurve(all_labels,all_pred_3,1) - - -% random forest feature importance -imp = Mdl3.OOBPermutedPredictorDeltaError; - -figure(1);clf; -bar(imp); -title('Curvature Test'); -ylabel('Predictor importance estimates'); -xticks([1:20]) - -figure(2);clf; -hold on -plot(X1,Y1,'LineWidth',2)%4 -plot(X2,Y2,'LineWidth',2)%7 -plot(X3,Y3,'LineWidth',2)%10 -legend(sprintf('Univariate - AUC: %f2',AUC1),sprintf('Bivariate - AUC: %f2',AUC2),sprintf('All - AUC: %f2',AUC3),'Location','SouthEast')%86/85/91 -xlabel('False positive rate') -ylabel('True positive rate') -title('All abnormal vs. normal channel classification') -%% render abnormalities -soz_ch_all = [zeros(1772,1);soz_ch]; -this_pt = 111; %114 res - 45/50 EIZ - 52/64, non-inv - 20,30 % 120 - res - 2., EIZ - 23./64, non-inv - 33/50. -this_pt_abn = all_pred_3(find(all_pts==this_pt)); -this_pt_soz = multi_class_all(find(all_pts==this_pt)); -this_pt_res = resect_ch_all(find(all_pts==this_pt)); -this_pt_coords = all_coords(find(all_pts==this_pt),:); -this_pt_eeg = all_wake_data(:,find(all_pts==this_pt)); -%this_pt_delta_coh_abn(isnan(this_pt_delta_coh_abn)) = 0; - -this_pt_zscores = all_feat_zscores(find(all_pts==this_pt),:); - -figure(1);clf; -hold on -for i = 1:10 - plot([1:12000],[this_pt_eeg(:,4*i)+i*175],'k-') -end -hold off - -%this_pt_soz(this_pt_soz==3) = 2; -%this_pt_soz(this_pt_res==1) = 3; - -% figure(1);clf; -% imagesc(all_feat_zscores(find(all_pts==this_pt),:)) -% colorbar -% caxis([0,5]) - -% figure(2);clf; -% subplot(1,3,1) -% imagesc(this_pt_zscores(2,:)') -% caxis([0,4]) -% colorbar -% subplot(1,3,2) -% imagesc(this_pt_zscores(23,:)') -% caxis([0,4]) -% colorbar -% subplot(1,3,3) -% imagesc(this_pt_zscores(50,:)') -% caxis([0,4]) -% colorbar -% -% figure(3); -% subplot(3,1,1) -% plot(this_pt_eeg(:,2)) -% subplot(3,1,2) -% plot(this_pt_eeg(:,23)) -% subplot(3,1,3) -% plot(this_pt_eeg(:,50)) -% - -figure(4);clf; -final_elec_matrix = [this_pt_coords,this_pt_abn,ones(size(this_pt_coords,1),1)]; -%final_elec_matrix = final_elec_matrix([2,23,50],:) -dlmwrite('render_elecs.node',final_elec_matrix,'delimiter',' ','precision',5) -BrainNet_MapCfg('BrainMesh_ICBM152_smoothed.nv','render_elecs.node','abnormality_render.mat') -delete render_elecs.node - - - -%% Examine correlations between z scores -[r,p] = corr(all_feat_zscores2) -figure(1);clf; -imagesc(r) -colorbar - -%% lesional vs non-lesional -EIZ = [zeros(1772,1);[spike_ind - spike_ind.*soz_ch]]; -all_labels = [all_pts,[EIZ+[zeros(1772,1);soz_ch]],EIZ,[zeros(1772,1);soz_ch],resect_ch_all,[ones(1772,1);HUP_outcome_all]]; - -%% get patient-level SOZ, RZ, spike indicators -all_normal_good = []; -all_normal_poor = []; -all_les_in_in = []; -all_nonles_in_in = []; -all_SOZ_in_in = []; -all_SOZ_in_out = []; -all_spike_in_in = []; -all_spike_in_out = []; - -b = 0; -ptile = 75; -for s = 107:166 - b = b+1 - a = s; - - pt_roi = pt_loc{s}; - res = resected_ch([patient_no==b]); - spike = spike_24h([patient_no==b]); - soz = soz_ch([patient_no==b]); - - roi_res = NaN*zeros(42,1); - -% for r = 1:42 -% this_roi = (r); -% roi_res(r) = mean(res(pt_roi==this_roi)); -% roi_soz(r) = mean(soz(pt_roi==this_roi)); -% roi_spike(r) = mean(spike(pt_roi==this_roi)); -% end -% -% HUP_scores(a).res = ceil(roi_res); -% HUP_scores(a).spike = (roi_spike>24); -% HUP_scores(a).soz = ceil(roi_soz); -% -% res_roi = [HUP_scores(a).res==1]; -% soz_roi = ([HUP_scores(a).soz==1]); -% spk_roi = ([HUP_scores(a).spike==1]); -% -% non_res_roi = ([HUP_scores(a).res==0]); -% non_soz_roi = ([HUP_scores(a).soz==0]); -% non_spk_roi = ([HUP_scores(a).spike==0]); -% - - % spared edges - all_normal = bivariate_native(feat).subj(a).data(find(~res),find(~res)); - HUP_normal(b,1) = nanmean(all_normal(:)); - - all_normal_good = [all_normal_good;prctile((all_normal(:)),ptile)]; - - % resected edges - %in_out_resect = bivariate_native(feat).subj(a).data(find(res),:); - in_in_resect = bivariate_native(feat).subj(a).data(find(res),find(~res)); - %HUP_resect(b,1) = nanmedian(in_out_resect(:)); % (in-out) - %HUP_resect(b,2) = nanmedian(in_in_resect(:)); % in-in - %HUP_resect(b,3) = nanmedian(test_HUP_scores(:)); % all - if lesion_field(b) - all_les_in_in = [all_les_in_in;prctile(in_in_resect(:),ptile)]; - else - all_nonles_in_in = [all_nonles_in_in;prctile(in_in_resect(:),ptile)]; - end - %all_RZ_in_out = [all_RZ_in_out;prctile((in_out_resect(:)),ptile)]; - - % soz edges (in-out) - in_out_soz = bivariate_native(feat).subj(a).data(find(soz),find(~soz)); - in_in_soz = bivariate_native(feat).subj(a).data(find(soz),find(soz)); - %HUP_soz(b,1) = nanmedian(in_out_soz(:)); % (in-out) - %HUP_soz(b,2) = nanmedian(in_in_soz(:)); % in-in - - all_SOZ_in_in = [all_SOZ_in_in;prctile(in_in_soz(:),ptile)]; - all_SOZ_in_out = [all_SOZ_in_out;prctile((in_out_soz(:)),ptile)]; - - % spike edges - try in_out_spike = bivariate_native(feat).subj(a).data(find(spike),find(~spike)); - all_spike_in_out = [all_spike_in_out;prctile((in_out_spike(:)),ptile)]; - catch anyerror - end - try in_in_spike = bivariate_native(feat).subj(a).data(find(spike),find(spike)); - all_spike_in_in = [all_spike_in_in;prctile(in_in_spike(:),ptile)]; - catch anyerror - end -% HUP_spike(a,1) = nanmedian(in_out_spike(:)); % (in-out) -% HUP_spike(a,2) = nanmedian(in_in_spike(:)); % in-in -% - - -% -end - -% -figure(1);clf; -hold on - -scatter(ones(length(all_normal_good),1),all_normal_good,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([0.75 1.25],[nanmedian(all_normal_good) nanmedian(all_normal_good)],'k-','LineWidth',2) -scatter(3*ones(length(all_SOZ_in_in),1),all_SOZ_in_in,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([2.75 3.25],[nanmedian(all_SOZ_in_in) nanmedian(all_SOZ_in_in)],'k-','LineWidth',2) -scatter(4*ones(length(all_SOZ_in_out),1),all_SOZ_in_out,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([3.75 4.25],[nanmedian(all_SOZ_in_out) nanmedian(all_SOZ_in_out)],'k-','LineWidth',2) -scatter(5*ones(length(all_spike_in_in),1),all_spike_in_in,'MarkerEdgeColor',color4,'MarkerFaceColor',color4,'jitter','on') -plot([4.75 5.25],[nanmedian(all_spike_in_in) nanmedian(all_spike_in_in)],'k-','LineWidth',2) -scatter(6*ones(length(all_spike_in_out),1),all_spike_in_out,'MarkerEdgeColor',color4,'MarkerFaceColor',color4,'jitter','on') -plot([5.75 6.25],[nanmedian(all_spike_in_out) nanmedian(all_spike_in_out)],'k-','LineWidth',2) -xlim([0.5, 6.5]) -ylim([0 6]) -hold off - -ranksum(all_normal_good,all_SOZ_in_in) -ranksum(all_normal_good,all_SOZ_in_out) -ranksum(all_normal_good,all_spike_in_in) -ranksum(all_normal_good,all_spike_in_out) - -ranksum(all_SOZ_in_in,all_SOZ_in_out) -ranksum(all_SOZ_in_in,all_spike_in_in) -ranksum(all_spike_in_in,all_spike_in_out) -ranksum(all_SOZ_in_out,all_spike_in_out) -%% boxplots baby - -num_good = sum([HUP_outcome==1]); -num_poor = sum([HUP_outcome>1]); - -good_out=HUP_normal(HUP_outcome==1); -poor_out=HUP_normal(HUP_outcome>1); - -good_in_out_soz=HUP_soz([HUP_outcome==1],1); -good_in_in_soz=HUP_soz([HUP_outcome==1],2); -poor_in_out_soz=HUP_soz([HUP_outcome>1],1); -poor_in_in_soz=HUP_soz([HUP_outcome>1],2); - -good_in_out_spike=HUP_spike([HUP_outcome==1],1); -good_in_in_spike=HUP_spike([HUP_outcome==1],2); -poor_in_out_spike=HUP_spike([HUP_outcome>1],1); -poor_in_in_spike=HUP_spike([HUP_outcome>1],2); - -all_plot_data = NaN*zeros(57,6); -all_plot_data(1:num_good,1) = good_in_in_soz; -all_plot_data(1:num_poor,2) = poor_in_in_soz; -all_plot_data(1:num_good,3) = good_in_out_soz; -all_plot_data(1:num_poor,4) = poor_in_out_soz; -all_plot_data(1:num_good,5) = good_out; -all_plot_data(1:num_poor,6) = poor_out; - -ranksum(good_in_in_soz,poor_in_in_soz) -ranksum(good_in_out_soz,poor_in_out_soz) -ranksum(good_out,poor_out) - -signrank(good_in_in_soz,good_in_out_soz) -signrank(good_in_out_soz,good_out) -signrank(good_in_in_soz,good_out) - -signrank(poor_in_in_soz,poor_in_out_soz) -signrank(poor_in_out_soz,poor_out) -signrank(poor_in_in_soz,poor_out) - -ranksum(good_in_in_spike,poor_in_in_spike) -ranksum(good_in_out_spike,poor_in_out_spike) -ranksum(good_out,poor_out) - -signrank(good_in_in_spike,good_in_out_spike) -signrank(good_in_out_spike,good_out) -signrank(good_in_in_spike,good_out) - -signrank(poor_in_in_spike,poor_in_out_spike) -signrank(poor_in_out_spike,poor_out) -signrank(poor_in_in_spike,poor_out) - -figure(1);clf; -boxplot(all_plot_data) -ylim([-2 2]) - -%% lesional vs non-lesional -r1 = HUP_soz((lesion_field==0),2); -r2 = HUP_soz(lesion_field,2); - -r3 = HUP_normal((lesion_field==0)); -r4 = HUP_normal(lesion_field); - -ranksum(r1,r2) - -figure(1);clf; -hold on -plot(ones(length(r1),1),r1,'ko') -plot(2*ones(length(r2),1),r2,'ko') -xlim([0.5 2.5]) -hold off - - diff --git a/data/atlas_metadata_final.xlsx b/data/atlas_metadata_final.xlsx index 944cc39..301f025 100644 Binary files a/data/atlas_metadata_final.xlsx and b/data/atlas_metadata_final.xlsx differ diff --git a/data/atlas_metadata_simplified.xlsx b/data/atlas_metadata_simplified.xlsx index fd9c307..f2fc460 100644 Binary files a/data/atlas_metadata_simplified.xlsx and b/data/atlas_metadata_simplified.xlsx differ diff --git a/data/custom_atlas.xlsx b/data/custom_atlas.xlsx index b791962..c16ecfa 100644 Binary files a/data/custom_atlas.xlsx and b/data/custom_atlas.xlsx differ diff --git a/data/~$atlas_metadata_final.xlsx b/data/~$atlas_metadata_final.xlsx new file mode 100644 index 0000000..19195df Binary files /dev/null and b/data/~$atlas_metadata_final.xlsx differ diff --git a/extra_code.m b/extra_code.m index de796cd..91a7f9e 100644 --- a/extra_code.m +++ b/extra_code.m @@ -380,4 +380,329 @@ plot(ones(length(r1),1),r1,'ko') plot(2*ones(length(r2),1),r2,'ko') xlim([0.5 2.5]) -hold off \ No newline at end of file +hold off + + +r1 = all_mean_wentropy([normal_HUP_ch;normal_MNI_ch]); +r2 = all_mean_wentropy(poor_out_ch); +r3 = all_mean_wentropy(EIZ_good_ch); +r4 = all_mean_wentropy(EIZ_poor_ch); +r5 = all_mean_wentropy(soz_good_ch); +r6 = all_mean_wentropy(soz_poor_ch); + +ranksum(all_mean_wentropy(normal_HUP_ch),all_mean_wentropy(normal_MNI_ch)) + +p1 = ranksum(r1,r3) +p2 = ranksum(r3,r5) +p3 = ranksum(r1,r5) +p4 = ranksum(r2,r4) +p5 = ranksum(r4,r6) +p6 = ranksum(r2,r6) +p7 = ranksum(r1,r2) +p8 = ranksum(r3,r4) +p9 = ranksum(r5,r6) + +figure(1);clf; +hold on +scatter(ones(length(r1),1),r1,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') +plot([0.75 1.25],[median(r1) median(r1)],'k-','LineWidth',2) +scatter(2*ones(length(r2),1),r2,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') +plot([1.75 2.25],[median(r2) median(r2)],'k-','LineWidth',2) +scatter(3*ones(length(r3),1),r3,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') +plot([2.75 3.25],[median(r3) median(r3)],'k-','LineWidth',2) +scatter(4*ones(length(r4),1),r4,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') +plot([3.75 4.25],[median(r4) median(r4)],'k-','LineWidth',2) +scatter(5*ones(length(r5),1),r5,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') +plot([4.75 5.25],[median(r5) median(r5)],'k-','LineWidth',2) +scatter(6*ones(length(r6),1),r6,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') +plot([5.75 6.25],[median(r6) median(r6)],'k-','LineWidth',2) +xlim([0.5, 6.5]) +ylim([2,22]) +ylabel('Wavelet entropy') +xticks([1:6]) +xticklabels({'Engel 1, non-SOZ, non-EIZ','Engel 2+, non-SOZ, non-EIZ', 'Engel 1 EIZ','Engel 2+ EIZ','Engel 1 SOZ','Engel 2+ SOZ'}) +hold off + +%% wavelet entropy validation across all regions +for i = 1:60 + for j = 1:size(all_wake_data,2) + start_inds = (i-1)*200+1; + end_inds = 200*i; + all_wentropy(i,j) = (wentropy(all_wake_data((start_inds:end_inds),j),'shannon')); + end +end + +all_mean_wentropy = log(-1*median(all_wentropy)); + + +figure(1);clf; +for i = 1:20 + % extract indices that correspond to a given ROI for HUP and MNI data + % (combining data from L and R hemispheres) + roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_EIZ = intersect([EIZ_good_ch;EIZ_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_soz = intersect([soz_good_ch;soz_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + + entropy_normal = all_mean_wentropy([roi_mni;roi_hup]); + entropy_EIZ = all_mean_wentropy(roi_EIZ); + entropy_soz = all_mean_wentropy(roi_soz); + + pval_EIZ(i,1) = ranksum(entropy_normal,entropy_EIZ); + pval_soz(i,1) = ranksum(entropy_normal,entropy_soz); + pval_EIZ_vs_soz(i,1) = ranksum(entropy_EIZ,entropy_soz); + + % median power spectral density across channels + %entropy_mni = median(pxx_norm(:,roi_mni),2); + %median_hup = median(pxx_norm(:,roi_hup),2); + plot_cell{1,1} = entropy_normal; + plot_cell{1,2} = entropy_EIZ; + plot_cell{1,3} = entropy_soz; + + subplot(4,5,i) + hold on + violin(plot_cell,'xlabel',{'','',''},'facecolor',[color1;color7;color2],'mc',[],'medc','k');%,'edgecolor','k'); + legend('off') + txt_sig1 = '+'; + txt_sig2 = '*'; + txt_sig3 = '#'; + ylim([7 23]) + + if pval_EIZ(i,1) < (0.05./20) + plot([1,2],[20,20],'k-','LineWidth',2) + text(1.45,20.25,txt_sig2,'FontSize', 20); + end + if pval_soz(i,1) < (0.05./20) + plot([1,3],[21.75,21.75],'k-','LineWidth',2) + text(2,22,txt_sig2,'FontSize', 20); + end + if pval_EIZ_vs_soz(i,1) < (0.05./20) + plot([2,3],[20,20],'k-','LineWidth',2) + text(2.5,20.25,txt_sig2,'FontSize', 20); + end + + ylabel('Log -Entropy') + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); + title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title + + hold off + +end + +%% +single_bivariate_feat_2 = single_bivariate_feat; +new_roi_2 = new_roi; +new_roi_2(incomplete_channels) = []; + +new_roi_2(new_roi_2==17) = 18; + +for ch = 1:4717 + if multi_class(ch)==1 + grouping{ch,1} = 'uninvolved'; + class_ind(ch,1) = 1; + elseif multi_class(ch)==2 + grouping{ch,1} = 'EIZ'; + class_ind(ch,1) = 2; + else + grouping{ch,1} = 'SOZ'; + class_ind(ch,1) = 3; + end +end + +univariate_feats = all_feat_zscores(:,1:5); +univariate_feats(incomplete_channels,:) = []; + +single_univariate_feat = max(abs(univariate_feats'))'; + +single_univariate_feat(isinf(single_univariate_feat)) = nanmedian(single_univariate_feat); +single_bivariate_feat_2(isinf(single_bivariate_feat_2)) = nanmedian(single_bivariate_feat_2); + +single_univariate_feat = single_univariate_feat(new_roi_2==18); +single_bivariate_feat_2 = single_bivariate_feat_2(new_roi_2==18); +grouping = grouping(new_roi_2==18); +class_ind = class_ind(new_roi_2==18); + + +figure(1);clf; +h = scatterhist(single_univariate_feat,single_bivariate_feat_2,'Group',grouping,'Marker','.','MarkerSize',12) + +clr = get(h(1),'colororder'); +boxplot(h(2),single_univariate_feat,grouping,'orientation','horizontal',... + 'label',{'','',''},'color',clr); + xlim(h(2),[0,30]) +boxplot(h(3),single_bivariate_feat_2,grouping,'orientation','horizontal',... + 'label', {'','',''},'color',clr); + xlim(h(3),[0,30]) + set(h(2:3),'XTickLabel',''); +view(h(3),[270,90]); % Rotate the Y plot +axis(h(1),'auto'); +xlim([0 10]) +ylim([0 10]) + +d1 = computeCohen_d(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) +p1 = ranksum(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) +d2 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) +p2 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) +d3 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) +p3 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) + +d4 = computeCohen_d(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) +p4 = ranksum(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) +d5 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) +p5 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) +d6 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) +p6 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) + +[d1 d2 d3 d4 d5 d6] +[p1 p2 p3 p4 p5 p6] +%% do node abormality (max univariate+bivariate in plot) for following: +% 1. within/outside soz for lesional / non-lesional + +all_lesion = zeros(5203,1); +for i = 1:166 + if i>106 + all_lesion(all_pts==i) = lesion_field(i-106); + end +end + +within_soz_lesional = find(all_lesion.*soz_ch_all); +within_soz_nonlesional = find([all_lesion==0].*soz_ch_all); +outside_soz_lesional = find(all_lesion.*[soz_ch_all==0]); +outside_soz_nonlesional = find([all_lesion==0].*[soz_ch_all==0]); + +for ch = 1:5203 + if ismember(ch,within_soz_lesional) + lesional_grouping{ch,1} = 'within_soz_lesional'; + elseif ismember(ch,within_soz_nonlesional) + lesional_grouping{ch,1} = 'within_soz_nonlesional'; + elseif ismember(ch,outside_soz_lesional) + lesional_grouping{ch,1} = 'outside_soz_lesional'; + elseif ismember(ch,outside_soz_nonlesional) + lesional_grouping{ch,1} = 'outside_soz_nonlesional'; + end +end + +univariate_feats = all_feat_zscores(:,1:10); +single_univariate_feat = prctile(abs(univariate_feats'),100)'; +single_bivariate_feat = prctile(abs_bivariate_feats',100)'; + +uni_within_soz_lesional = single_univariate_feat(within_soz_lesional); +uni_within_soz_nonlesional = single_univariate_feat(within_soz_nonlesional); +uni_outside_soz_lesional = single_univariate_feat(outside_soz_lesional); +uni_outside_soz_nonlesional = single_univariate_feat(outside_soz_nonlesional); +bi_within_soz_lesional = single_bivariate_feat(within_soz_lesional); +bi_within_soz_nonlesional = single_bivariate_feat(within_soz_nonlesional); +bi_outside_soz_lesional = single_bivariate_feat(outside_soz_lesional); +bi_outside_soz_nonlesional = single_bivariate_feat(outside_soz_nonlesional); + +p1 = ranksum(uni_within_soz_lesional,uni_within_soz_nonlesional); +p2 = ranksum(uni_within_soz_lesional,uni_outside_soz_lesional); +p3 = ranksum(uni_within_soz_nonlesional,uni_outside_soz_nonlesional); +p4 = ranksum(uni_outside_soz_lesional,uni_outside_soz_nonlesional); +p5 = ranksum(bi_within_soz_lesional,bi_within_soz_nonlesional); +p6 = ranksum(bi_within_soz_lesional,bi_outside_soz_lesional); +p7 = ranksum(bi_within_soz_nonlesional,bi_outside_soz_nonlesional); +p8 = ranksum(bi_outside_soz_lesional,bi_outside_soz_nonlesional); + +single_univariate_feat(incomplete_channels) = []; +single_bivariate_feat(incomplete_channels) = []; +lesional_grouping(incomplete_channels) = []; + +figure(2);clf; +h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',lesional_grouping,'Marker','.','MarkerSize',12) + +clr = get(h(1),'colororder'); +boxplot(h(2),single_univariate_feat,lesional_grouping,'orientation','horizontal',... + 'label',{'','','',''},'color',clr); + xlim(h(2),[0,30]) +boxplot(h(3),single_bivariate_feat,lesional_grouping,'orientation','horizontal',... + 'label', {'','','',''},'color',clr); + xlim(h(3),[0,30]) + set(h(2:3),'XTickLabel',''); +view(h(3),[270,90]); % Rotate the Y plot +axis(h(1),'auto'); +xlim([0 10]) +ylim([0 10]) + +[p1 p2 p3 p4 p5 p6 p7 p8] +%% 2. within/outside ablation zone for good / poor outcome +all_outcome = ones(5203,1); +all_abl = NaN(5203,1); +for i = 1:166 + if i>106 + all_outcome(all_pts==i) = late_outcome(i-106); + if therapy_field(i-106) + all_abl(all_pts==i) = resect_ch_all(all_pts==i); + end + end +end + +within_abl_good = find([all_outcome==1].*[all_abl==1]); +outside_abl_good = find([all_outcome==1].*[all_abl==0]); +within_abl_poor = find([all_outcome>1].*[all_abl==1]); +outside_abl_poor = find([all_outcome>1].*[all_abl==0]); + +incomplete_channels2 = incomplete_channels; + +abl_grouping = cell(5203,1); + +for ch = 1:5203 + if ismember(ch,within_abl_good) + abl_grouping{ch,1} = 'within_abl_good'; + elseif ismember(ch,outside_abl_good) + abl_grouping{ch,1} = 'outside_abl_good'; + elseif ismember(ch,within_abl_poor) + abl_grouping{ch,1} = 'within_abl_poor'; + elseif ismember(ch,outside_abl_poor) + abl_grouping{ch,1} = 'outside_abl_poor'; + else + incomplete_channels2 = [incomplete_channels2, ch]; + end +end + +incomplete_channels2 = unique(incomplete_channels2); + +univariate_feats = all_feat_zscores(:,1:10); +single_univariate_feat = prctile(abs(univariate_feats'),100)'; +single_bivariate_feat = prctile(abs_bivariate_feats',100)'; + +uni_within_abl_good = single_univariate_feat(within_abl_good); +uni_outside_abl_good = single_univariate_feat(outside_abl_good); +uni_within_abl_poor = single_univariate_feat(within_abl_poor); +uni_outside_abl_poor = single_univariate_feat(outside_abl_poor); +bi_within_abl_good = single_bivariate_feat(within_abl_good); +bi_outside_abl_good = single_bivariate_feat(outside_abl_good); +bi_within_abl_poor = single_bivariate_feat(within_abl_poor); +bi_outside_abl_poor = single_bivariate_feat(outside_abl_poor); + +p1 = ranksum(uni_within_abl_good,uni_outside_abl_good); +p2 = ranksum(uni_within_abl_poor,uni_outside_abl_poor); +p3 = ranksum(uni_within_abl_good,uni_within_abl_poor); +p4 = ranksum(uni_outside_abl_good,uni_outside_abl_poor); +p5 = ranksum(bi_within_abl_good,bi_outside_abl_good); +p6 = ranksum(bi_within_abl_poor,bi_outside_abl_poor); +p7 = ranksum(bi_within_abl_good,bi_within_abl_poor); +p8 = ranksum(bi_outside_abl_good,bi_outside_abl_poor); + +single_univariate_feat(incomplete_channels2) = []; +single_bivariate_feat(incomplete_channels2) = []; +abl_grouping(incomplete_channels2) = []; + +figure(3);clf; +h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',abl_grouping,'Marker','.','MarkerSize',12) + +clr = get(h(1),'colororder'); +boxplot(h(2),single_univariate_feat,abl_grouping,'orientation','horizontal',... + 'label',{'','','',''},'color',clr); + xlim(h(2),[0,30]) +boxplot(h(3),single_bivariate_feat,abl_grouping,'orientation','horizontal',... + 'label', {'','','',''},'color',clr); + xlim(h(3),[0,30]) + set(h(2:3),'XTickLabel',''); +view(h(3),[270,90]); % Rotate the Y plot +axis(h(1),'auto'); +xlim([0 10]) +ylim([0 10]) + +[p1 p2 p3 p4 p5 p6 p7 p8] diff --git a/final_figs/.DS_Store b/final_figs/.DS_Store new file mode 100644 index 0000000..1bb8cb7 Binary files /dev/null and b/final_figs/.DS_Store differ diff --git a/final_figs/AUPRC_comparison.svg b/final_figs/AUPRC_comparison.svg new file mode 100644 index 0000000..7201800 --- /dev/null +++ b/final_figs/AUPRC_comparison.svg @@ -0,0 +1,501 @@ + + +Engel 1Engel 2+00.20.40.60.81Area under precision-recall curveComparison between AUPRC, rank-sum p = 0.033 diff --git a/final_figs/classification_abs.svg b/final_figs/classification_abs.svg new file mode 100644 index 0000000..9f171aa --- /dev/null +++ b/final_figs/classification_abs.svg @@ -0,0 +1,214 @@ + + +00.20.40.60.81False positive rate00.10.20.30.40.50.60.70.80.91True positive rateAll abnormal vs. normal channel classificationUnivariate - AUC: 0.64Bivariate - AUC: 0.72All - AUC: 0.78 diff --git a/final_figs/curvature_abs.svg b/final_figs/curvature_abs.svg new file mode 100644 index 0000000..ad45cf8 --- /dev/null +++ b/final_figs/curvature_abs.svg @@ -0,0 +1,206 @@ + + +Curvature Test1234567891000.20.40.60.811.21.4Predictor importance estimates diff --git a/final_figs/new_figureset.zip b/final_figs/new_figureset.zip new file mode 100644 index 0000000..85c7e99 Binary files /dev/null and b/final_figs/new_figureset.zip differ diff --git a/final_figs/new_figureset/.DS_Store b/final_figs/new_figureset/.DS_Store new file mode 100644 index 0000000..d827125 Binary files /dev/null and b/final_figs/new_figureset/.DS_Store differ diff --git a/final_figs/new_figureset/Figure1_new.tif b/final_figs/new_figureset/Figure1_new.tif new file mode 100644 index 0000000..d15e77b Binary files /dev/null and b/final_figs/new_figureset/Figure1_new.tif differ diff --git a/final_figs/new_figureset/Figure2_new.tif b/final_figs/new_figureset/Figure2_new.tif new file mode 100644 index 0000000..16928fb Binary files /dev/null and b/final_figs/new_figureset/Figure2_new.tif differ diff --git a/final_figs/new_figureset/Figure3_new.tif b/final_figs/new_figureset/Figure3_new.tif new file mode 100644 index 0000000..6a246db Binary files /dev/null and b/final_figs/new_figureset/Figure3_new.tif differ diff --git a/final_figs/new_figureset/Figure4_new.tif b/final_figs/new_figureset/Figure4_new.tif new file mode 100644 index 0000000..6f9cce1 Binary files /dev/null and b/final_figs/new_figureset/Figure4_new.tif differ diff --git a/final_figs/new_figureset/Figure5_new.tif b/final_figs/new_figureset/Figure5_new.tif new file mode 100644 index 0000000..30d7457 Binary files /dev/null and b/final_figs/new_figureset/Figure5_new.tif differ diff --git a/final_figs/new_figureset/Figure6_new.tif b/final_figs/new_figureset/Figure6_new.tif new file mode 100644 index 0000000..3b4ae4e Binary files /dev/null and b/final_figs/new_figureset/Figure6_new.tif differ diff --git a/final_figs/new_figureset/Figure7_new.tif b/final_figs/new_figureset/Figure7_new.tif new file mode 100644 index 0000000..eba59e1 Binary files /dev/null and b/final_figs/new_figureset/Figure7_new.tif differ diff --git a/final_figs/new_figureset/Zscore_vs_regression_heatmaps.png b/final_figs/new_figureset/Zscore_vs_regression_heatmaps.png new file mode 100644 index 0000000..ef7198b Binary files /dev/null and b/final_figs/new_figureset/Zscore_vs_regression_heatmaps.png differ diff --git a/iEEG_analysis_final.m b/iEEG_analysis_final.m new file mode 100644 index 0000000..48b0040 --- /dev/null +++ b/iEEG_analysis_final.m @@ -0,0 +1,1139 @@ +%% Set up workspace +clear all + +% set up path - change iEEG_atlas_path to where you download the repository +iEEG_atlas_path = '/Users/jbernabei/Documents/PhD_Research/atlas_project/iEEG_atlas_dev'; + +% other paths may stay the same +metadata = readtable('data/atlas_metadata_final.xlsx'); +custom_atlas = readtable('data/custom_atlas.xlsx'); + +% anatomical atlas information +region_list = zeros(1,90); % for the 90 AAL regions we will be using +region_names = cell(1,90); +fi = fopen("localization/AAL116_WM.txt"); +for j = 1:90 + label = split(fgetl(fi)); + region_list(j) = str2double(label{3}); + region_names{j} = label{2}; +end + +% extract metadata +age_onset = metadata{:,11}; +age_surgery = metadata{:,12}; +HUP_outcome = floor(metadata{:,3}); + +% extract outcome +for s = 1:length(age_surgery) + engel_scores = metadata{s,3:5}; + engel_scores = rmmissing(engel_scores); + early_outcome(s) = floor(engel_scores(1)); % early outcome is 6 months + late_outcome(s) = floor(engel_scores(end)); % late outcome is latest time point up to 24 months. These are used in the project. +end + +% extract additional fields of metadata +therapy_field = metadata{:,6}; +implant_field = metadata{:,7}; +target_field = metadata{:,8}; +laterality_field = metadata{:,9}; +lesion_field = metadata{:,10}; +gender_field = metadata{:,13}; +lesion_field = strcmp(lesion_field,'Lesional'); % 1 for lesional +therapy_field = strcmp(therapy_field,'Ablation'); % 1 for ablation + +load('data/MNI_atlas_new.mat') +load('data/HUP_atlas_final.mat') + +all_pts = [Patient; (patient_no+106)]; % combine patient numbers for MNI & HUP + +HUP_outcome_all = zeros(length(soz_ch),1); +for i = 1:60 + HUP_outcome_all(find(patient_no==i)) = late_outcome(i); % we use outcome as latest measured outcome +end + +% clear unused data from MNI atlas +clear NodesLeft NodesLeftInflated NodesRightInflated NodesRight NodesRegionLeft NodesRegionRight Data_N2 Data_N3 Data_R + +% some colors for plotting, etc (r/g/b) +color1 = [0, 0.4470, 0.7410]; color2 = [0.6350, 0.0780, 0.1840]; color3 = [255, 248, 209]./255; color4 = [103 55 155]/255; +color6 = [78 172 91]/255; color7 = [0.9290, 0.6940, 0.1250]; color8 = [0.8500, 0.3250, 0.0980]; + +% assign wake clips from MNI and HUP into a final data matrix +all_wake_data = [Data_W,wake_clip]; + +% perform localization based on AAL116 atlas (which has an internal WM mask also) +all_coords = [ChannelPosition;mni_coords]; + +try [coords, all_roi, NN_ind] = nifti_values(all_coords,'localization/AAL116_WM.nii'); +catch ME + for i = 1:size(all_coords,1) + if mod(i,100) + fprintf('%d\n',i) + end + try [mni_coords1, mni_roi, NN_ind] = nifti_values(all_coords(i,:),'localization/AAL116_WM.nii'); + catch ME + mni_roi = NaN; + end + all_roi(i,1) = mni_roi; + end +end + +clear NN_ind +clear coords + +% assign AAL regions into new, aggregate regions of custom atlas. (116 -> 40 regions) +region_matrix = custom_atlas{:,2:4}; +for c = 1:length(all_roi) + old_ind = find(region_list==all_roi(c)); + try [row,col,V] = find(region_matrix==old_ind); + try new_roi(c,1) = row; + catch anyerror + new_roi(c,1) = 0; + end + catch anyerror + new_roi(c,1) = 0; + end + +end + +% define normal MNI/HUP and abnormal channels (SOZ/high-spike/RZ) +% set threshold on spike counts +spike_thresh = 24; % this is empirical, 1 spike/hour + +% find indices for which spikes are greater than threshold +spike_ind = [spike_24h>spike_thresh]; + +% define all abnormal channels +abnormal_ch = find([spike_ind+soz_ch+(HUP_outcome_all>1)]>0)+1772; + +% define all seizure onset indices +soz_ch_inds = find(soz_ch)+1772; + +% define normal HUP channels +normal_HUP_ch = find([spike_ind+soz_ch+(HUP_outcome_all>1)]==0)+1772; + +% define normal MNI channels +normal_MNI_ch = [1:1772]'; + +% define all normal channels +all_normal_ch = [normal_MNI_ch;normal_HUP_ch]; + +% define exclusive irritative zone +EIZ_ch = find(spike_ind - spike_ind.*soz_ch)+1772; + +% binary variables for whether channel is in exclusive irritative zone, resected, or +% seizure onset zone +EIZ_binary = [zeros(1772,1);[spike_ind - spike_ind.*soz_ch]]; +resect_binary = [zeros(1772,1);resected_ch]; +soz_binary = [zeros(1772,1);soz_ch]; + +% assign atlas locations per patient into cell structure +for s = 1:max(all_pts) + pt_loc{s} = new_roi(all_pts==s)'; +end + +% indices of all abnormal channels +all_abnormal_ch = zeros(5203,1); +all_abnormal_ch(abnormal_ch) = 1; + +% make structure of patient-level clinically abnormal channels, useful for +% atlas construction +for s = 1:166 + + % check abnormal channels for HUP patients + if s>106 + + % abnormal channel indices + all_abn{s} = find(all_abnormal_ch(all_pts==s)); + if isempty(all_abn{s}) + all_abn{s} = []; + end + else + % if MNI patient, no abnormal channels by definition + all_abn{s} = []; + end +end + +%% *** Supplemental figure 1A *** +% count HUP and MNI in original atlas and then new atlas + +mni_roi_old = all_roi(1:1772); +hup_roi_old = all_roi(1773:end); + +for i = 1:45 + roi_count(i,1) = sum(mni_roi_old==region_list(2*i-1))+sum(mni_roi_old==region_list(2*i)); + roi_count(i,2) = sum(all_roi(normal_HUP_ch)==region_list(2*i-1))+sum(all_roi(normal_HUP_ch)==region_list(2*i)); + this_roi_name = split(string(region_names{2*i}),'_R'); % extract roi name without laterality + split_roi_name_old{i,1} = this_roi_name(1); +end + +for i = 1:20 + roi_count_new_normal(i,1) = length(intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); + roi_count_new_normal(i,2) = length(intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); +end + +new_counts = [roi_count_new_normal(:,1),(roi_count_new_normal(:,2))]; + +figure(1);clf; +subplot(1,2,1) +barh(roi_count) +set(gca,'TickLabelInterpreter', 'none'); +xlabel('Number of nodes per region') +title('AAL regions') +yticks([1:45]) +yticklabels(split_roi_name_old) +subplot(1,2,2) +barh(new_counts) +set(gca,'TickLabelInterpreter', 'none'); +title('Composite regions') +yticks([1:21]) +yticklabels(split_roi_name) +xlabel('Number of nodes per region') + +%% % *** Supplemental figure 1B *** +% -> render the atlas_render_color.txt file in BrainNet viewer do renderings of original and new atlas +random_colors = rand([40,3]); + +color_matrix = ones(117,3); + +for i = 1:40 + this_region = custom_atlas{i,2:4}; + for r = 1:length(rmmissing(this_region)) + color_matrix(this_region(r),:) = random_colors(i,:); + end +end + +% writes a file that can be used in BrainNet viewer to visualize +dlmwrite('atlas_render_color.txt',color_matrix(1:90,:),'delimiter',' ','precision',5) + +% color scale for rendering +figure(1);clf; +imagesc(reshape(random_colors,40,1,3)) +set(gca,'TickLabelInterpreter', 'none'); +yticks([1:40]) +yticklabels(custom_atlas{:,1}) + +%% *** Figure 4 and Supplemental Figure 3 *** +% calculate power spectrum +all_wake_data = [Data_W,wake_clip]; + +[pxx,f] = pwelch(all_wake_data,400,200,[0.5:0.5:80],200); % original was 200/100 +pxx_norm = pxx./sum(pxx); + +% Analyze power spectral density between HUP/MNI, Normal/irritative zone/SOZ +% change this to be on a per-patient level rather than a per-node level + +freq_inds = [1, 8; % 0.5 - 4 Hz + 9, 16; % 4 - 8 Hz + 17, 26; % 8 - 13 Hz + 27, 60; % 13 - 30 Hz + 61, 160]; % 30 - 80 Hz + + +figure(1);clf % supplemental +figure(2);clf % supplemental +for i = 1:20 + % extract indices that correspond to a given ROI for HUP and MNI data + % (combining data from L and R hemispheres) + roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_EIZ = intersect([EIZ_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_soz = intersect([soz_ch_inds],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + + % median power spectral density across channels + median_mni = median(pxx_norm(:,roi_mni),2); + median_hup = median(pxx_norm(:,roi_hup),2); + median_normal = median(pxx_norm(:,[roi_mni;roi_hup]),2); + median_EIZ = median(pxx_norm(:,roi_EIZ),2); + median_soz = median(pxx_norm(:,roi_soz),2); + + for j = 1:5 + HUP_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_hup))'; + MNI_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_mni))'; + normal_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],[roi_mni;roi_hup]))'; + EIZ_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_EIZ))'; + soz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_soz))'; + + [p,h] = ranksum(HUP_feat,MNI_feat); % HUP vs MNI + pval_matrix_norm(i,j) = p; + + [p,h] = ranksum(normal_feat,EIZ_feat); % normal vs EIZ + pval_matrix_EIZ(i,j) = p; + + [p,h] = ranksum(normal_feat,soz_feat); % normal vs SOZ + pval_matrix_soz(i,j) = p; + + [p,h] = ranksum(EIZ_feat,soz_feat); % EIZ vs SOZ + pval_matrix_EIZ_vs_soz(i,j) = p; + + end + + % find upper and lower bounds (here 25th and 75th percentiles) + prct_25_75_mni = prctile(pxx_norm(:,roi_mni)',[25,75]); + prct_25_75_hup = prctile(pxx_norm(:,roi_hup)',[25,75]); + prct_25_75_normal = prctile(pxx_norm(:,[roi_mni;roi_hup])',[25,75]); + prct_25_75_EIZ = prctile(pxx_norm(:,roi_EIZ)',[25,75]); + prct_25_75_soz = prctile(pxx_norm(:,roi_soz)',[25,75]); + + x_axis = [0.5:0.5:80]; + + figure(1) + % 3 x 7 plot for 21 regions (aggregating across hemispheres) + subplot(4,5,i) + hold on + plot([4,4],[0,0.4],'k-') + plot([8,8],[0,0.4],'k-') + plot([13,13],[0,0.4],'k-') + plot([30,30],[0,0.4],'k-') + plot(x_axis,median_mni) + patch([x_axis fliplr(x_axis)], [prct_25_75_mni(1,:) fliplr(prct_25_75_mni(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') + plot(x_axis,median_hup) + patch([x_axis fliplr(x_axis)], [prct_25_75_hup(1,:) fliplr(prct_25_75_hup(2,:))], color6, 'FaceAlpha',0.5, 'EdgeColor','none') + + % add greek letters for frequency bands + txt_d = '\delta'; + txt_t = '\theta'; + txt_a = '\alpha'; + txt_b = '\beta'; + txt_g = '\gamma'; + txt_sig = '*'; + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); + if pval_matrix_norm(i,1)<(0.05/100) + text(2,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,2)<(0.05/100) + text(5.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,3)<(0.05/100) + text(9.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,4)<(0.05/100) + text(18,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,5)<(0.05/100) + text(40,0.12,txt_sig,'FontSize', 20); + end + + % set axes, labels, subgraph titles + set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) + ylabel('Spectral Density') + xlabel('Frequency (Hz)') + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); + title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title + xlim([0.5,80]) % xlimit 0.5-80 Hz + ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) + hold off + + figure(2) + subplot(4,5,i) + hold on + plot([4,4],[0,0.4],'k-') + plot([8,8],[0,0.4],'k-') + plot([13,13],[0,0.4],'k-') + plot([30,30],[0,0.4],'k-') + plot(x_axis,median_normal) + patch([x_axis fliplr(x_axis)], [prct_25_75_normal(1,:) fliplr(prct_25_75_normal(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') + plot(x_axis,median_EIZ,'Color',color7) + patch([x_axis fliplr(x_axis)], [prct_25_75_EIZ(1,:) fliplr(prct_25_75_EIZ(2,:))], color7, 'FaceAlpha',0.5, 'EdgeColor','none') + try plot(x_axis,median_soz,'Color',color2) + patch([x_axis fliplr(x_axis)], [prct_25_75_soz(1,:) fliplr(prct_25_75_soz(2,:))], color2, 'FaceAlpha',0.5, 'EdgeColor','none') + catch anyerror + end + + % add greek letters for frequency bands + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); + if pval_matrix_EIZ(i,1)<(0.05/100) + text(2,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,2)<(0.05/100) + text(5.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,3)<(0.05/100) + text(9.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,4)<(0.05/100) + text(18,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,5)<(0.05/100) + text(40,0.12,txt_sig,'FontSize', 20); + end + + if pval_matrix_soz(i,1)<(0.05/100) + text(2,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,2)<(0.05/100) + text(5.5,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,3)<(0.05/100) + text(9.5,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,4)<(0.05/100) + text(18,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,5)<(0.05/100) + text(40,0.17,txt_sig,'FontSize', 20); + end + + % set axes, labels, subgraph titles + set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) + ylabel('Spectral Density') + xlabel('Frequency (Hz)') + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); + title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title + xlim([0.5,80]) % xlimit 0.5-80 Hz + ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) + hold off + +end + +%% calculate univariate z-scores of neural activity +% power spectrum and entropy abnormality + +all_f = {'delta','theta','alpha','beta','gamma'}; % frequency bands +all_m = {'bandpower'};%,'entropy'}; % metrics + +c = 0; % c for column in z score matrix +for m = 1 % m for different metrics + for f = 1:5 % f for frequency bands + c = c+1; + + % first create the univariate atlas, taking only normal channels + % the univariate atlas is mean/median + stdev of a given feature + % across ROI + [mean_vec,median_vec,std_vec,feat_vals] = create_univariate_atlas(all_wake_data(:,find(all_normal_ch)),[1:40],new_roi(all_normal_ch),sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); + + % we call this function again on all of the data, but are only + % using it to extract all of the feature values on each channel + [~,~,~,feat_vals] = create_univariate_atlas(all_wake_data,[1:40],new_roi,sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); + + % now we test the univariate atlas to generate z scores + [zscores, pt_zscores, pt_ind, roi_ind] = test_univariate_atlas(median_vec, std_vec, feat_vals, all_pts, [1:40], new_roi,'patient'); + + % store the z scores + univariate_zscores(:,c) = zscores; + all_pt_zscores(:,c) = pt_zscores; + + end +end + +%save('univariate_zscores.mat','univariate_zscores') + + +%% create atlas data - > edge-wise and patient-wise should be options +% for HUP data, aggregate electrodes to ignore +% these should be poor-outcome/RZ/SOZ/>1 spike per hour + +conn_band = {'delta','theta','alpha','beta','gamma'}; +conn_type = {'coh'}; + +% loop through types of connectivity and frequency bands +for c = 1 + for f = 2:5 + + % loop through subjects + for s = 107:max(all_pts) + s + + % extract subject data from all the data (columns are channels) + pt_data = all_wake_data(:,[all_pts==s]); + + % call functional connectivity code (200 is Hz, 1 is window size in seconds) + [pt_adj] = compute_FC(pt_data, 200, 1, conn_type{c},conn_band{f}); + + % assign median adjacency matrix into cell structure + all_pt_adj{s} = pt_adj; + + end + + % save adjacency matrices + %save(sprintf('adj_matrices_rev2/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f}),'all_pt_adj') + + end +end + +%% functional connectivity atlas construction & z-score code + +conn_band = {'delta','theta','alpha','beta','gamma'}; +conn_type = {'coh'}; + +% which feature +feat = 0; + +% loop through connectivity and frequency +for c = 1 + for f = 1:5 + + feat = feat+1; + a = 0 + + load(sprintf('adj_matrices/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f})); + + [conn_edge, std_edge, samples_edge, sem_edge, raw_atlas_edge] = create_atlas_by_edge(all_pt_adj, pt_loc, all_abn, [1:40]', 2); + + % call for MNI patients + [conn_edge_mni, std_edge_mni, samples_edge_mni, ~, raw_atlas_edge_mni] = create_atlas_by_edge(all_pt_adj(1:106), pt_loc(1:106), all_abn(1:106), [1:40]', 1); + + % call for HUP patients + [conn_edge_hup, std_edge_hup, samples_edge_hup, ~, raw_atlas_edge_hup] = create_atlas_by_edge(all_pt_adj(107:166), pt_loc(107:166), all_abn(107:166), [1:40]', 1); + + this_feat_conn(feat).mni = conn_edge_mni; + this_feat_conn(feat).hup = conn_edge_hup; + + % loop through all patients + for s = 1:166 % change this back + a = a+1; + + % do nodal space + try [native_adj_scores, corr_val] = test_native_adj(all_pt_adj{s}, pt_loc{s}, conn_edge, std_edge, [1:40]'); + + bivariate_native(feat).subj(a).data = native_adj_scores; + + catch anyerror + bivariate_native(feat).subj(a).data = NaN(length(pt_loc{s})); + end + + end + end + +end + +%% process bivariate edge features into nodal features + +for f = 1:5 + this_feat = []; + abs_feat = []; + for s = 1:166 + % take 75th percentile across all edges of each node + %abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),75)']; + abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),75)']; % test with 50th percentile + end + + % assign into feature matrix + abs_bivariate_feats(:,f) = abs_feat; +end + +% for the single bivariate feature we want the maximum absolute Z score +% across all 10 bivariate features +[single_bivariate_feat] = median(abs_bivariate_feats')'; + +% assign previously calculated univariate z scores into bivariate scores +all_feat_zscores = [abs(univariate_zscores), abs_bivariate_feats]; +all_feat_zscores(isinf(all_feat_zscores)) = 0; % zero out any z-scores that were infinite due to zero variance in edge estimates + +% incomplete channels are ones which are unlocalized in new atlas +incomplete_channels = find(isnan(sum(all_feat_zscores'))); + +% remove incomplete channels from univariate features +univariate_feats = all_feat_zscores(:,1:5); +univariate_feats(incomplete_channels,:) = []; + + +%% set up machine learning features and labels +% define channels in which sparse atlas coverage yields NaN's +incomplete_channels = find(isnan(sum(all_feat_zscores'))); +complete_channels = find(~isnan(sum(all_feat_zscores'))); + +% start with channel indices for ML features as all channels +ch_indices = [1:5203]'; + +% remove channels for which features are incomplete from further analysis +zscores_complete = all_feat_zscores; +zscores_complete(incomplete_channels,:) = []; % remove channels that are incomplete +ch_indices(incomplete_channels) = []; % remove channel indices that are incomplete + +% define binary variable for which channels are abnormal +abn_ch_binary = zeros(5203,1); +abn_ch_binary(abnormal_ch) = 1; + +% define a matrix containing different labels. The binary variables +% correspond to the following labels +% column 1: patient number +% column 2: whether the channel is normal or abnormal +% column 3: the exclusive irritative zone +% column 4: the seizure onset zone +% column 5: whether or not the channel was resected +% column 6: surgical outcome +% column 7: whether the channel was both resected and seizure onset zone +all_class_labels = [all_pts,abn_ch_binary,EIZ_binary,soz_binary,resect_binary,[ones(1772,1);HUP_outcome_all],[[zeros(1772,1);soz_ch].*resect_binary]]; +all_class_labels(incomplete_channels,:) = []; + +% we want only resected, soz channels vs normal channels +normal_channels = find(all_class_labels(:,2)==0); +resected_soz_channels = find(all_class_labels(:,7)==1); +abn_ch_list = find(all_class_labels(:,2)==1); + +% extract features for these channels from all features and define labels +ML_features = zscores_complete([normal_channels;resected_soz_channels],:); +ML_chs = ch_indices([normal_channels;resected_soz_channels]); % original channel indices for data we are using +ML_labels = all_class_labels([normal_channels;resected_soz_channels],7); % 0 for normal 1 for resected-SOZ +ML_pt_inds = all_class_labels([normal_channels;resected_soz_channels],1); % patient indices for data we are using + +%% *** Figure 6A and 6B *** +% random forest and feature importance + +clear part all_pred_1 all_pred_2 all_pred_3 + +num_parts = 10; + +[pt_part] = make_xval_partition(166, num_parts); + +% cross-validate based on patients +for p = 1:num_parts + indices = find(pt_part==p) + for i = 1:length(indices) + part(find(indices(i)==ML_pt_inds(:,1))) = p; + + part2(find(indices(i)==all_class_labels(:,1))) = p; + end +end + +% Define classes, 0 = non seizure, 1 = seizure +ClassNames = [0,1]; + +% Define costs -> penalty here is 5 +cost.ClassNames = ClassNames; +cost.ClassificationCosts = [0 1; 1 0]; + +for p = 1:num_parts + p + X_train = ML_features(find(part~=p),:); + Y_train = ML_labels(find(part~=p)); + X_test = ML_features(find(part==p),:); + Y_test = ML_labels(find(part==p)); + + X_analysis = zscores_complete(find(part2==p),:); + + Mdl1 = TreeBagger(100,X_train(:,1:5),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + Mdl2 = TreeBagger(100,X_train(:,6:10),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + Mdl3 = TreeBagger(100,X_train,Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + [Yfit,pihat1] = predict(Mdl1,X_test(:,1:5)); + [Yfit,pihat2] = predict(Mdl2,X_test(:,6:10)); + [Yfit,pihat3] = predict(Mdl3,X_test); + [Yfit, pihat_analysis] = predict(Mdl3,X_analysis); + [Yfit, pihat_analysis2] = predict(Mdl1,X_analysis(:,1:5)); + [Yfit, pihat_analysis3] = predict(Mdl2,X_analysis(:,6:10)); + + [~, Y_pred] = max(pihat1,[],2); + all_pred_1(find(part==p)) = pihat1(:,2); + all_pred_2(find(part==p)) = pihat2(:,2); + all_pred_3(find(part==p)) = pihat3(:,2); + all_pred_analysis(find(part2==p),1) = pihat_analysis(:,2); + all_pred_analysis2(find(part2==p),1) = pihat_analysis2(:,2); + all_pred_analysis3(find(part2==p),1) = pihat_analysis3(:,2); + + xfold_acc(p) = sum(Y_pred==Y_test)./length(Y_test); + + Y_pred = round(pihat1(:,2)); + + +end + + [X1,Y1,T,AUC1] = perfcurve(ML_labels,all_pred_1,1); + [X2,Y2,T,AUC2] = perfcurve(ML_labels,all_pred_2,1); + [X3,Y3,T,AUC3] = perfcurve(ML_labels,all_pred_3,1); + + AUC1 + AUC2 + AUC3 + + Mdl3 = TreeBagger(100,ML_features,ML_labels,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + + +% random forest feature importance +imp = Mdl3.OOBPermutedPredictorDeltaError; + +figure(1);clf; +bar(imp); +title('Curvature Test'); +ylabel('Predictor importance estimates'); +xticks([1:10]) + +figure(2);clf; +hold on +plot(X1,Y1,'LineWidth',2,'Color',color4)%4 +plot(X2,Y2,'LineWidth',2,'Color',color7)%7 +plot(X3,Y3,'LineWidth',2,'Color',color1)%10 +legend(sprintf('Univariate - AUC: %.2f',AUC1),sprintf('Bivariate - AUC: %.2f',AUC2),sprintf('All - AUC: %.2f',AUC3),'Location','SouthEast')%86/85/91 +xlabel('False positive rate') +ylabel('True positive rate') +title('All abnormal vs. normal channel classification') + +%% *** Figure 5A *** +% mean Z-scores for Mesial Temporal Lobe (regions 17-20) in SOZ vs IZ vs uninvolved at nodal level +new_roi_2 = new_roi; +new_roi_2(incomplete_channels) = []; + +single_bivariate_feat_2 = single_bivariate_feat; +single_bivariate_feat_2(incomplete_channels) = []; + +MTL_uninvolved = []; +MTL_EIZ = []; +MTL_SOZ = []; + +for pt = 1:166 + uninv_pt = find([all_class_labels(:,1)==pt].*[all_class_labels(:,2)==0].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + EIZ_pt = find([all_class_labels(:,1)==pt].*[all_class_labels(:,3)==1].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + SOZ_pt = find([all_class_labels(:,1)==pt].*[all_class_labels(:,4)==1].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + + MTL_uninvolved = [MTL_uninvolved; mean(single_bivariate_feat_2(uninv_pt))]; + MTL_EIZ = [MTL_EIZ; mean(single_bivariate_feat_2(EIZ_pt))]; + MTL_SOZ = [MTL_SOZ; mean(single_bivariate_feat_2(SOZ_pt))]; +end + +mean_uninv = MTL_uninvolved; +mean_EIZ = MTL_EIZ; +mean_SOZ = MTL_SOZ; + +plot_matrix = padcat(mean_uninv,mean_EIZ,mean_SOZ); +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','EIZ','SOZ'}) + +p1 = ranksum(mean_uninv,mean_EIZ,'tail','left') +cd1 = computeCohen_d(mean_uninv,mean_EIZ) +p2 = ranksum(mean_uninv,mean_SOZ,'tail','left') +cd2 = computeCohen_d(mean_uninv,mean_SOZ) +p3 = ranksum(mean_EIZ,mean_SOZ,'tail','left') +cd3 = computeCohen_d(mean_SOZ,mean_EIZ) + +%% need to map back onto all channels and render +all_pred = zeros(5203,1); +all_pred(incomplete_channels) = 0; +all_pred(complete_channels) = all_pred_analysis; + +all_pred2 = zeros(5203,1); +all_pred2(incomplete_channels) = 0; +all_pred2(complete_channels) = all_pred_analysis2; + +all_pred3 = zeros(5203,1); +all_pred3(incomplete_channels) = 0; +all_pred3(complete_channels) = all_pred_analysis3; + +%% *** Figure 7 *** +% render abnormalities +soz_ch_all = [zeros(1772,1);soz_ch]; +this_pt = 116;%116/122/154 / 154 %114 res - 45/50 EIZ - 52/64, non-inv - 20,30 % 120 - res - 2., EIZ - 23./64, non-inv - 33/50. +this_pt_abn = all_pred(find(all_pts==this_pt)); +this_pt_soz = soz_ch_all(find(all_pts==this_pt)); +this_pt_eiz = EIZ_binary(find(all_pts==this_pt)); +this_pt_res = resect_binary(find(all_pts==this_pt)); +this_pt_coords = all_coords(find(all_pts==this_pt),:); +this_pt_eeg = all_wake_data(:,find(all_pts==this_pt)); +%this_pt_delta_coh_abn(isnan(this_pt_delta_coh_abn)) = 0; + +this_pt_multiclass = 2*this_pt_soz+this_pt_eiz; + +this_pt_zscores = all_feat_zscores(find(all_pts==this_pt),:); + +this_pt_test = zeros(size(this_pt_res)); +this_pt_test([28,9,49])=1; + +figure(2);clf; +subplot(1,3,1) +imagesc(this_pt_zscores(28,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +subplot(1,3,2) +imagesc(this_pt_zscores(9,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +subplot(1,3,3) +imagesc(this_pt_zscores(49,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +% +figure(3); +subplot(3,1,1) +plot(this_pt_eeg(6001:12000,4)) +subplot(3,1,2) +plot(this_pt_eeg(6001:12000,23)) +subplot(3,1,3) +plot(this_pt_eeg(6001:12000,27)) + + +figure(4);clf; +final_elec_matrix = [this_pt_coords,this_pt_abn,ones(size(this_pt_res))]; +%final_elec_matrix = final_elec_matrix([2,23,50],:) +dlmwrite('render_elecs.node',final_elec_matrix,'delimiter',' ','precision',5) +BrainNet_MapCfg('BrainMesh_ICBM152_smoothed.nv','render_elecs.node','abnormality_render.mat') +%delete render_elecs.node + +%% patient-specific correlation between abnormality features +for s = 107:166 + pt_zscores = zscores_complete(find([all_class_labels(:,1)==s]),:); + + pt_zscores_univar = pt_zscores(:,1:5); + pt_zscores_bivar = pt_zscores(:,6:10); + + corrval2(s-106) = corr(pt_zscores_univar(:),pt_zscores_bivar(:),'Type','Pearson'); + + for f = 1:5 + corrval((s-106),f) = corr(pt_zscores(:,f),pt_zscores(:,(f+5)),'Type','Pearson'); + end + + corrval_outcome((s-106),1) = late_outcome(s-106); + corrval_lesion((s-106),1) = lesion_field(s-106); +end + +p1 = ranksum(corrval(corrval_outcome==1,1),corrval(corrval_outcome>1,1),'tail','right') +p2 = ranksum(corrval(corrval_outcome==1,2),corrval(corrval_outcome>1,2),'tail','right') +p3 = ranksum(corrval(corrval_outcome==1,3),corrval(corrval_outcome>1,3),'tail','right') +p4 = ranksum(corrval(corrval_outcome==1,4),corrval(corrval_outcome>1,4),'tail','right') +p5 = ranksum(corrval(corrval_outcome==1,5),corrval(corrval_outcome>1,5),'tail','right') + +p6 = ranksum(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1),'tail','right') +d6 = computeCohen_d(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1)) + +plot_matrix = padcat(corrval2(corrval_outcome==1)',corrval2(corrval_outcome>1)'); +figure(1);clf; +UnivarScatter(plot_matrix); +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylim([-0.3 0.7]) + + +%% *** Supplemental Figure 8 *** +% Resected vs non-resected epileptogenicity +a = 0; +for s = [107:166] + a = a+1; + pred_pt = all_pred3(all_pts==s); + mean_in(a,1) = mean(pred_pt(find(resect_binary(all_pts==s)))); + mean_out(a,1) = mean(pred_pt(find(~resect_binary(all_pts==s)))); + res_test(a,1) = mean(pred_pt(find(resect_binary(all_pts==s))))>mean(pred_pt(find(~resect_binary(all_pts==s)))) + +end + +mean_diff = mean_in-mean_out; + +figure(1);clf; +hold on +plot_matrix = padcat(mean_diff(late_outcome==1),mean_diff(late_outcome>1)); +UnivarScatter(plot_matrix); +xticks([1:2]) +xticklabels({'Engel 1','Engel 2+'}) +plot([0.5 2.5],[0 0],'k-.') +ylabel('Mean resected - mean non-resected epileptogenicity') + + + +%% *** Figure 6D *** +% compute AUPRC/DRS at a per-patient level +a = 0; +for s = [107:132,134:139,141:155,157:161,163:166] + s + a = a+1; + this_pt_label = all_class_labels(find(all_class_labels(:,1)==s),4); + this_pt_pred = all_pred_analysis(find(all_class_labels(:,1)==s)); + this_pt_pred2 = all_pred_analysis2(find(all_class_labels(:,1)==s)); + this_pt_pred3 = all_pred_analysis3(find(all_class_labels(:,1)==s)); + [X1,Y1,T,AUC1] = perfcurve(this_pt_label,this_pt_pred,1,'XCrit','TPR','YCrit','PPV'); + [X2,Y2,T,AUC2] = perfcurve(this_pt_label,this_pt_pred2,1,'XCrit','TPR','YCrit','PPV'); + [X3,Y3,T,AUC3] = perfcurve(this_pt_label,this_pt_pred3,1,'XCrit','TPR','YCrit','PPV'); + + auc_pt_outcome(a,1) = late_outcome(s-106); + pt_auc(a,1:3) = [AUC1,AUC2,AUC3]; +end + +good_inds = find(auc_pt_outcome==1); +poor_inds = find(auc_pt_outcome>1); + +plot_matrix = padcat(pt_auc(good_inds,3),pt_auc(poor_inds,3)); +figure(1);clf; +UnivarScatter(plot_matrix) + +ranksum(plot_matrix(:,1),plot_matrix(:,2)) +computeCohen_d(plot_matrix(:,1),plot_matrix(:,2)) +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylabel('Area under precision-recall curve') + +%% *** Figure 6C *** +% patient-level abnormality score all regions + +pt_uninvolved = []; +pt_eiz = []; +pt_soz = []; + +for s = 107:166 + + pt_uninvolved_inds = find([all_class_labels(:,2)==0].*[all_class_labels(:,1)==s]); + pt_eiz_inds = find([all_class_labels(:,3)==1].*[all_class_labels(:,1)==s]); + pt_soz_inds = find([all_class_labels(:,4)==1].*[all_class_labels(:,1)==s]); + + pt_uninvolved = [pt_uninvolved; nanmedian(all_pred_analysis(pt_uninvolved_inds))]; + pt_eiz = [pt_eiz; nanmedian(all_pred_analysis(pt_eiz_inds))]; + pt_soz = [pt_soz; nanmedian(all_pred_analysis(pt_soz_inds))]; +end + +plot_matrix = padcat(pt_uninvolved,pt_eiz,pt_soz); + +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','irritative zone','seizure onset zone'}) +ylabel('predicted epileptogenicity score') +ylim([0 0.8]) + +p1 = ranksum(pt_uninvolved,pt_eiz) +d1 = computeCohen_d(pt_uninvolved,pt_eiz) +p2 = ranksum(pt_eiz,pt_soz) +d2 = computeCohen_d(pt_eiz,pt_soz) +p3 = ranksum(pt_uninvolved,pt_soz) +d3 = computeCohen_d(pt_uninvolved,pt_soz) + +%% *** Figure 5B *** +% correlation length of diagnosis + +single_univariate_feat_2 = median(univariate_feats')'; + +pt_in_rz = []; +pt_out_rz = []; + +pt_in_rz2 = []; +pt_out_rz2 = []; + +for s = 107:166 + + w_inds = find([all_class_labels(:,1)==s]); + o_inds = find([all_class_labels(:,1)==s]); + + length_diagnosis = age_surgery - age_onset; + + pt_in_rz = [pt_in_rz; nanmedian(single_bivariate_feat_2(w_inds))]; + pt_out_rz = [pt_out_rz; nanmedian(single_bivariate_feat_2(o_inds))]; + + +end + +figure(1);clf; +hold on +plot(length_diagnosis, pt_out_rz,'b.','Color',[106 178 180]/255,'MarkerSize',12) +p = polyfit(length_diagnosis,pt_out_rz,1); +f = polyval(p,length_diagnosis); +plot(length_diagnosis,f,'b-','LineWidth',2) +[R,P] = corr(length_diagnosis, pt_out_rz,'type','Pearson','tail','right') + +%% calculate spectral density differences between dataset 2 and dataset 3 + +load extra_data_1 +extra1_wake = wake_clip; + +load extra_data_2 +extra2_wake = wake_clip; + +[pxx,f] = pwelch(extra1_wake,400,200,[0.5,2:80],200); +pxx_norm1 = pxx./sum(pxx); + +[pxx,f] = pwelch(extra2_wake,400,200,[0.5,2:80],200); +pxx_norm2 = pxx./sum(pxx); + +for freq = 1:5 + if freq==1 + feat_vals1 = median(pxx_norm1(1:4,:))'; + feat_vals2 = median(pxx_norm2(1:4,:))'; + elseif freq==2 + feat_vals1 = median(pxx_norm1(4:8,:))'; + feat_vals2 = median(pxx_norm2(4:8,:))'; + elseif freq==3 + feat_vals1 = median(pxx_norm1(8:13,:))'; + feat_vals2 = median(pxx_norm2(8:13,:))'; + elseif freq==4 + feat_vals1 = median(pxx_norm1(13:30,:))'; + feat_vals2 = median(pxx_norm2(13:30,:))'; + elseif freq==5 + feat_vals1 = median(pxx_norm1(30:80,:))'; + feat_vals2 = median(pxx_norm2(30:80,:))'; + end + + freq_feats1(freq,1:3431) = feat_vals1'; + freq_feats2(freq,1:3431) = feat_vals2'; + +end + +for freq = 1:5 + freq_corr(freq) = corr(freq_feats1(freq,:)',freq_feats2(freq,:)') +end + +%% + +load AED_recordings.mat +figure(1);clf +a = 0; + +for pt = [1,2,4] + + freq_data_pt = []; + + for c = 1:4 + data_clip = AED_recordings(pt).clip(c).data; + [pxx,f] = pwelch(data_clip,400,200,[0.5,2:80],200); + pxx_norm_AED = pxx./sum(pxx); + + % calculate power data + for freq = 1:5 + if freq==1 + feat_vals_AED = median(pxx_norm_AED(1:4,:))'; + elseif freq==2 + feat_vals_AED = median(pxx_norm_AED(4:8,:))'; + elseif freq==3 + feat_vals_AED = median(pxx_norm_AED(8:13,:))'; + elseif freq==4 + feat_vals_AED = median(pxx_norm_AED(13:30,:))'; + elseif freq==5 + feat_vals_AED = median(pxx_norm_AED(30:80,:))'; + end + + freq_data_pt(c,:,freq) = feat_vals_AED'; + + end + + end + + % analyze across clips + for freq = 1:5 + a = a+1; + subplot(3,5,a) + hold on + plot(freq_data_pt(2,:,freq),freq_data_pt(4,:,freq),'b.','MarkerSize',8) + lsline + [RHO,PVAL] = corr(freq_data_pt(2,:,freq)',freq_data_pt(4,:,freq)','type','Pearson'); + pvals(freq,pt) = PVAL; + rhovals(freq,pt) = RHO; + title(sprintf('rho = %.2f',RHO)) + hold off + end + +end + +%% do a connectivity calculation here +conn_band = {'delta','theta','alpha','beta','gamma'}; + +for pt = 5 + + % loop through subjects + for freq = 1:5 + + for c = 1:4 + + % extract subject data from all the data (columns are channels) + pt_data = AED_recordings(pt).clip(c).data; + + % call functional connectivity code (200 is Hz, 1 is window size in seconds) + [pt_adj] = compute_FC(pt_data, 200, 1, 'coh',conn_band{freq}); + + % assign median adjacency matrix into cell structure + AED_conn(pt).freq(freq).data(:,:,c) = pt_adj; + end + + end +end + +m = 1 + +%% +figure(1);clf; +a = 0; +for pt = [1,2,4] + for freq = 1:5 + matrix1 = -1.*log(AED_conn(pt).freq(freq).data(:,:,2)); + matrix2 = -1.*log(AED_conn(pt).freq(freq).data(:,:,4)); + + a = a+1; + subplot(3,5,a) + hold on + plot(matrix1(:),matrix2(:),'b.','MarkerSize',2) + lsline + xlim([0,6]) + ylim([0,6]) + [RHO,PVAL] = corr(matrix1(:),matrix2(:),'type','Pearson'); + pvals_conn(freq,pt) = PVAL; + rhovals_conn(freq,pt) = RHO; + title(sprintf('rho = %.2f',RHO)) + hold off + end +end + +%% correlate feature space + +load('zscores_complete_rev1.mat') +zscores_complete_rev1 = zscores_complete; + +load('zscores_complete_rev2.mat') +zscores_complete_rev2 = zscores_complete; + +for f = 1:10 +a = 0; + for s = 106:166 + a = a+1; + [dataset_r, dataset_p] = corr(zscores_complete_rev1((all_class_labels(:,1)==s),f),zscores_complete_rev2((all_class_labels(:,1)==s),f)); + dataset_feat_corr_r(a,f) = dataset_r; + dataset_feat_corr_p(a,f) = dataset_p; + end +end + +figure(1);clf; +UnivarScatter(dataset_feat_corr_r) +xticks([1:10]) +xticklabels({'delta','theta','alpha','beta','gamma','delta','theta','alpha','beta','gamma'}) +ylabel('Correlation') + +%% correlate predicted epileptogenicity + +load('all_pred3_rev1.mat') +all_pred3_rev1 = all_pred3; + +load('all_pred3_rev2.mat') +all_pred3_rev2 = all_pred3; + +a = 0; +for s = 107:166 + a = a+1; + [dataset_r, dataset_p] = corr(all_pred3_rev1(all_pts==s),all_pred3_rev2(all_pts==s)); + dataset_pred_corr_r(a,1) = dataset_r; + dataset_pred_corr_p(a,1) = dataset_p; +end + +figure(1);clf; +histogram(dataset_pred_corr_r, 'NumBins', 10) +xlabel('Pearson correlation') +ylabel('Number of subjects') +%histogram(dataset_pred_corr_p) + +%% +figure(1);clf; +hold on +plot(this_thresh,thresh_auc_abs,'b--o','LineWidth',1,'Color',color1,'MarkerSize',8,'MarkerEdgeColor',color1,'MarkerFaceColor',color1) +plot(this_thresh,thresh_auc_nonabs,'r--o','LineWidth',1,'Color',color2,'MarkerSize',8,'MarkerEdgeColor',color2,'MarkerFaceColor',color2) +legend('Absolute value','Directional features','Location','NorthWest') +xlabel('Percentile threshold') +ylabel('AUROC') + + diff --git a/atlas_synthesis.m b/iEEG_atlas_analysis.m similarity index 53% rename from atlas_synthesis.m rename to iEEG_atlas_analysis.m index cc7d433..7ef74b3 100644 --- a/atlas_synthesis.m +++ b/iEEG_atlas_analysis.m @@ -29,6 +29,7 @@ engel_scores = rmmissing(engel_scores); early_outcome(s) = floor(engel_scores(1)); % early outcome is 6 months late_outcome(s) = floor(engel_scores(end)); % late outcome is latest time point up to 24 months + late_outcome_raw(s) = engel_scores(end); end therapy_field = metadata{:,6}; @@ -42,7 +43,7 @@ therapy_field = strcmp(therapy_field,'Ablation'); % 1 ffor ablation load('data/MNI_atlas_new.mat') -load('data/HUP_atlas_May25.mat') +load('data/HUP_atlas_final.mat') all_pts = [Patient; (patient_no+106)]; % combine patient numbers for MNI & HUP @@ -67,6 +68,8 @@ color3 = [255, 248, 209]./255; color4 = [103 55 155]/255; color6 = [78 172 91]/255; +color7 = [0.9290, 0.6940, 0.1250]; +color8 = [0.8500, 0.3250, 0.0980]; % assign wake clips from MNI and HUP into a final data matrix all_wake_data = [Data_W,wake_clip]; @@ -113,13 +116,13 @@ spike_ind = [spike_24h>spike_thresh]; % define all abnormal channels -abnormal_ch = find([resected_ch+spike_ind+soz_ch]>0)+1772; +abnormal_ch = find([resected_ch+spike_ind+soz_ch+(HUP_outcome_all>1)]>0)+1772; % define all seizure onset indices soz_ch_inds = find(soz_ch)+1772; % define normal HUP channels -normal_HUP_ch = find([resected_ch+spike_ind+soz_ch]==0)+1772; +normal_HUP_ch = find([resected_ch+spike_ind+soz_ch+(HUP_outcome_all>1)]==0)+1772; % define normal MNI channels normal_MNI_ch = [1:1772]'; @@ -130,17 +133,6 @@ % define exclusive irritative zone EIZ_ch = find(spike_ind - spike_ind.*soz_ch)+1772; -% split EIZ into good vs poor -EIZ_good_ch = find((spike_ind - spike_ind.*soz_ch).*(HUP_outcome_all==1))+1772; -EIZ_poor_ch = find((spike_ind - spike_ind.*soz_ch).*(HUP_outcome_all>1))+1772; - -% split SOZ into good vs poor -soz_good_ch = find(soz_ch.*(HUP_outcome_all==1))+1772; -soz_poor_ch = find(soz_ch.*(HUP_outcome_all>1))+1772; - -% define non-EIZ non-SOZ poor outcome channels -poor_out_ch = find([resected_ch+spike_ind+soz_ch+(HUP_outcome_all==1)]==0)+1772; - %% count HUP and MNI in original atlas and then new atlas mni_roi_old = all_roi(1:1772); hup_roi_old = all_roi(1773:end); @@ -202,7 +194,7 @@ %% calculate power spectrum all_wake_data = [Data_W,wake_clip]; -[pxx,f] = pwelch(all_wake_data,200,100,[0.5:0.5:80],200); +[pxx,f] = pwelch(all_wake_data,400,200,[0.5:0.5:80],200); % original was 200/100 pxx_norm = pxx./sum(pxx); %% Analyze power spectral density between HUP/MNI, Normal/irritative zone/SOZ @@ -221,8 +213,8 @@ % (combining data from L and R hemispheres) roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_EIZ = intersect([EIZ_good_ch;EIZ_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_soz = intersect([soz_good_ch;soz_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_EIZ = intersect([EIZ_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_soz = intersect([soz_ch_inds],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); % median power spectral density across channels median_mni = median(pxx_norm(:,roi_mni),2); @@ -281,11 +273,11 @@ txt_b = '\beta'; txt_g = '\gamma'; txt_sig = '*'; - text(2,0.15,txt_d); - text(5.5,0.15,txt_t); - text(9.5,0.15,txt_a); - text(18,0.15,txt_b); - text(40,0.15,txt_g); + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); if pval_matrix_norm(i,1)<(0.05/100) text(2,0.12,txt_sig,'FontSize', 20); end @@ -322,19 +314,19 @@ plot([30,30],[0,0.4],'k-') plot(x_axis,median_normal) patch([x_axis fliplr(x_axis)], [prct_25_75_normal(1,:) fliplr(prct_25_75_normal(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') - plot(x_axis,median_EIZ) - patch([x_axis fliplr(x_axis)], [prct_25_75_EIZ(1,:) fliplr(prct_25_75_EIZ(2,:))], color6, 'FaceAlpha',0.5, 'EdgeColor','none') - try plot(x_axis,median_soz) + plot(x_axis,median_EIZ,'Color',color7) + patch([x_axis fliplr(x_axis)], [prct_25_75_EIZ(1,:) fliplr(prct_25_75_EIZ(2,:))], color7, 'FaceAlpha',0.5, 'EdgeColor','none') + try plot(x_axis,median_soz,'Color',color2) patch([x_axis fliplr(x_axis)], [prct_25_75_soz(1,:) fliplr(prct_25_75_soz(2,:))], color2, 'FaceAlpha',0.5, 'EdgeColor','none') catch anyerror end % add greek letters for frequency bands - text(2,0.15,txt_d); - text(5.5,0.15,txt_t); - text(9.5,0.15,txt_a); - text(18,0.15,txt_b); - text(40,0.15,txt_g); + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); if pval_matrix_EIZ(i,1)<(0.05/100) text(2,0.12,txt_sig,'FontSize', 20); end @@ -384,10 +376,10 @@ try plot(x_axis,median_soz) % power spectrum and entropy abnormality all_f = {'delta','theta','alpha','beta','gamma'}; % frequency bands -all_m = {'bandpower','entropy'}; % metrics +all_m = {'bandpower'};%,'entropy'}; % metrics c = 0; % c for column in z score matrix -for m = 1:2 % m for different metrics +for m = 1 % m for different metrics for f = 1:5 % f for frequency bands c = c+1; @@ -409,169 +401,8 @@ try plot(x_axis,median_soz) end end -%save('univariate_zscores.mat','univariate_zscores') - -%% mass univariate testing -% 21 ROI x 5 = 105 tests to correct for -freq_inds = [1, 8; % 0.5 - 4 Hz - 9, 16; % 4 - 8 Hz - 17, 26; % 8 - 13 Hz - 27, 60; % 13 - 30 Hz - 61, 160]; % 30 - 80 Hz - -for i = 1:20 - for j = 1:5 - roi_soz = intersect((find(soz_ch)+1772),[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_eiz = intersect((EIZ_ch),[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - soz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_soz))'; - eiz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_eiz))'; - HUP_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_hup))'; - MNI_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_mni))'; - [p,h] = ranksum(HUP_feat,MNI_feat); % HUP vs MNI - pval_matrix_norm(i,j) = p; - [p1,h] = ranksum([HUP_feat;MNI_feat],soz_feat); - pval_matrix_soz(i,j) = p1; - [p2,h] = ranksum([HUP_feat;MNI_feat],eiz_feat); - pval_matrix_eiz(i,j) = p2; - end -end - -figure(1);clf; -subplot(1,3,1) -imagesc((pval_matrix_norm.*100)<0.05) -title('HUP vs. MNI normal channels') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) -subplot(1,3,2) -imagesc((pval_matrix_eiz.*100)<0.05) -title('Irritative zone vs. composite atlas') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) -subplot(1,3,3) -imagesc((pval_matrix_soz.*100)<0.05) -title('Seizure onset zone vs. composite atlas') -xticks([1:5]) -xticklabels({'Delta','Theta','Alpha','Beta','Gamma'}) -yticks([1:20]) -set(gca,'TickLabelInterpreter', 'none'); -yticklabels(split_roi_name) - -%% wavelet entropy validation across all regions -for i = 1:60 - for j = 1:size(all_wake_data,2) - start_inds = (i-1)*200+1; - end_inds = 200*i; - all_wentropy(i,j) = wentropy(all_wake_data((start_inds:end_inds),j),'shannon'); - end -end - -all_mean_wentropy = log(-1*median(all_wentropy)); - -r1 = all_mean_wentropy([normal_HUP_ch;normal_MNI_ch]); -r2 = all_mean_wentropy(poor_out_ch); -r3 = all_mean_wentropy(EIZ_good_ch); -r4 = all_mean_wentropy(EIZ_poor_ch); -r5 = all_mean_wentropy(soz_good_ch); -r6 = all_mean_wentropy(soz_poor_ch); - -ranksum(all_mean_wentropy(normal_HUP_ch),all_mean_wentropy(normal_MNI_ch)) - -p1 = ranksum(r1,r3) -p2 = ranksum(r3,r5) -p3 = ranksum(r1,r5) -p4 = ranksum(r2,r4) -p5 = ranksum(r4,r6) -p6 = ranksum(r2,r6) -p7 = ranksum(r1,r2) -p8 = ranksum(r3,r4) -p9 = ranksum(r5,r6) - -figure(1);clf; -hold on -scatter(ones(length(r1),1),r1,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([0.75 1.25],[median(r1) median(r1)],'k-','LineWidth',2) -scatter(2*ones(length(r2),1),r2,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([1.75 2.25],[median(r2) median(r2)],'k-','LineWidth',2) -scatter(3*ones(length(r3),1),r3,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([2.75 3.25],[median(r3) median(r3)],'k-','LineWidth',2) -scatter(4*ones(length(r4),1),r4,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([3.75 4.25],[median(r4) median(r4)],'k-','LineWidth',2) -scatter(5*ones(length(r5),1),r5,'MarkerEdgeColor',color1,'MarkerFaceColor',color1,'jitter','on') -plot([4.75 5.25],[median(r5) median(r5)],'k-','LineWidth',2) -scatter(6*ones(length(r6),1),r6,'MarkerEdgeColor',color2,'MarkerFaceColor',color2,'jitter','on') -plot([5.75 6.25],[median(r6) median(r6)],'k-','LineWidth',2) -xlim([0.5, 6.5]) -ylim([2,22]) -ylabel('Wavelet entropy') -xticks([1:6]) -xticklabels({'Engel 1, non-SOZ, non-EIZ','Engel 2+, non-SOZ, non-EIZ', 'Engel 1 EIZ','Engel 2+ EIZ','Engel 1 SOZ','Engel 2+ SOZ'}) -hold off - -%% -figure(1);clf; -for i = 1:20 - % extract indices that correspond to a given ROI for HUP and MNI data - % (combining data from L and R hemispheres) - roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_EIZ = intersect([EIZ_good_ch;EIZ_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - roi_soz = intersect([soz_good_ch;soz_poor_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); - - entropy_normal = all_mean_wentropy([roi_mni;roi_hup]); - entropy_EIZ = all_mean_wentropy(roi_EIZ); - entropy_soz = all_mean_wentropy(roi_soz); - - pval_EIZ(i,1) = ranksum(entropy_normal,entropy_EIZ); - pval_soz(i,1) = ranksum(entropy_normal,entropy_soz); - pval_EIZ_vs_soz(i,1) = ranksum(entropy_EIZ,entropy_soz); - - % median power spectral density across channels - %entropy_mni = median(pxx_norm(:,roi_mni),2); - %median_hup = median(pxx_norm(:,roi_hup),2); - plot_cell{1,1} = entropy_normal; - plot_cell{1,2} = entropy_EIZ; - plot_cell{1,3} = entropy_soz; - - subplot(4,5,i) - hold on - violin(plot_cell,'xlabel',{'','',''},'facecolor',[color1;color6;color2],'mc',[],'medc','k');%,'edgecolor','k'); - legend('off') - txt_sig1 = '+'; - txt_sig2 = '*'; - txt_sig3 = '#'; - ylim([7 23]) - - if pval_EIZ(i,1) < (0.05./20) - plot([1,2],[20,20],'k-','LineWidth',2) - text(1.45,20.25,txt_sig2,'FontSize', 20); - end - if pval_soz(i,1) < (0.05./20) - plot([1,3],[21.75,21.75],'k-','LineWidth',2) - text(2,22,txt_sig2,'FontSize', 20); - end - if pval_EIZ_vs_soz(i,1) < (0.05./20) - plot([2,3],[20,20],'k-','LineWidth',2) - text(2.5,20.25,txt_sig2,'FontSize', 20); - end - - ylabel('Probability') - xlabel('- Log Entropy') - this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality - split_roi_name{i,1} = this_roi_name(1); - title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title - - hold off - -end +%save('univariate_zscores.mat','univariate_zscores') %% create atlas data - > edge-wise and patient-wise should be options @@ -582,7 +413,7 @@ try plot(x_axis,median_soz) conn_type = {'corr','coh'}; % loop through types of connectivity and frequency bands -for c = 1:2 +for c = 2 for f = 1:5 % loop through subjects @@ -590,7 +421,7 @@ try plot(x_axis,median_soz) s % extract subject data from all the data (columns are channels) - pt_data = sleep_clip;%all_wake_data(:,[all_pts==s]); + pt_data = all_wake_data(:,[all_pts==s]); % call functional connectivity code (200 is Hz, 1 is window size in seconds) [pt_adj] = compute_FC(pt_data, 200, 1, conn_type{c},conn_band{f}); @@ -634,13 +465,13 @@ try plot(x_axis,median_soz) %% functional cconnectivity atlas construction & z-score code conn_band = {'delta','theta','alpha','beta','gamma'}; -conn_type = {'corr','coh'}; +conn_type = {'coh'}; % which feature feat = 0; % loop through connectivity and frequency -for c = 1:2 +for c = 1 for f = 1:5 feat = feat+1; @@ -662,6 +493,9 @@ try plot(x_axis,median_soz) % call for HUP patients [conn_edge_hup, std_edge_hup, samples_edge_hup, ~, raw_atlas_edge_hup] = create_atlas_by_edge(all_pt_adj(107:166), pt_loc(107:166), all_abn(107:166), [1:40]', 1); + this_feat_conn(feat).mni = conn_edge_mni; + this_feat_conn(feat).hup = conn_edge_hup; + % loop through all patients for s = 1:166 % change this back a = a+1; @@ -681,6 +515,13 @@ try plot(x_axis,median_soz) end %% correlate connectivity between MNI and HUP +for feat = 1:5 + a = conn_band{feat}; + aa = 'coherence'; + + conn_edge_mni = this_feat_conn(feat).mni; + conn_edge_hup = this_feat_conn(feat).hup; + conn_edge_all = conn_edge_mni(:)+conn_edge_hup(:); conn_edge_mni(isnan(conn_edge_all)) = []; conn_edge_hup(isnan(conn_edge_all)) = []; @@ -689,22 +530,31 @@ try plot(x_axis,median_soz) std_edge_mni(isnan(std_edge_all)) = []; std_edge_hup(isnan(std_edge_all)) = []; -figure(1);clf; -plot(exp(conn_edge_mni),exp(conn_edge_hup),'ko') -[r,p] = corr(exp(conn_edge_mni)',exp(conn_edge_hup)') -xlabel('MNI delta coherence') -ylabel('HUP delta coherence') -title('Correlation between atlas edges, r = 0.43, p = 4.6e-19') -xlim([0 0.1]) -ylim([0 0.1]) +figure(1); +subplot(1,5,feat) +hold on +plot((conn_edge_mni),(conn_edge_hup),'k.','Color',[128 128 128]/255) +[r,p1] = corr((conn_edge_mni)',(conn_edge_hup)') +rvals(feat) = r; +xlabel('MNI connectivity') +ylabel('HUP connectivity') +p = polyfit((conn_edge_mni),(conn_edge_hup),1); +f = polyval(p,(conn_edge_mni)); +plot((conn_edge_mni),f,'k-','LineWidth',2) +title(sprintf('%s %s, r = %f, p = %f',a,aa,r,p1)) +hold off +end + +% figure(2);clf +% bar(rvals) %% process bivariate edge features into nodal features -for f = 1:10 +for f = 1:5 this_feat = []; abs_feat = []; for s = 1:166 - % take 80th percentile across all edges of each node - abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),80)']; + % take 75th percentile across all edges of each node + abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),75)']; end % assign into feature matrix @@ -713,386 +563,493 @@ try plot(x_axis,median_soz) % for the single bivariate feature we want the maximum absolute Z score % across all 10 bivariate features -single_bivariate_feat = max(abs_bivariate_feats')'; +[single_bivariate_feat] = median(abs_bivariate_feats')'; % assign previously calculated univariate z scores into bivariate scores -all_feat_zscores = [univariate_zscores, abs_bivariate_feats]; +all_feat_zscores = [abs(univariate_zscores), abs_bivariate_feats]; % incomplete channels are ones which are unlocalized in new atlas incomplete_channels = find(isnan(sum(all_feat_zscores'))); -single_bivariate_feat(incomplete_channels) = []; +%single_bivariate_feat(incomplete_channels) = []; + +univariate_feats = all_feat_zscores(:,1:5); +univariate_feats(incomplete_channels,:) = []; +%[single_univariate_feat, which_feat_uni] = max(abs(univariate_feats),[],2); + +%% mean Z-scores for MTL (regions 17-20) in SOZ vs IZ vs uninvolved at nodal level +MTL_uninvolved = find([all_possible_labels(:,2)==0].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); +MTL_EIZ = find([all_possible_labels(:,3)==1].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); +MTL_SOZ = find([all_possible_labels(:,4)==1].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); + +mean_uninv = single_bivariate_feat_2(MTL_uninvolved); +mean_EIZ = single_bivariate_feat_2(MTL_EIZ); +mean_SOZ = single_bivariate_feat_2(MTL_SOZ); + +plot_matrix = padcat(mean_uninv,mean_EIZ,mean_SOZ); +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','EIZ','SOZ'}) +ylim([0 10]) +p1 = ranksum(mean_EIZ,mean_uninv) +cd1 = computeCohen_d(mean_EIZ,mean_uninv) +p2 = ranksum(mean_SOZ,mean_uninv) +cd2 = computeCohen_d(mean_SOZ,mean_uninv) +p3 = ranksum(mean_EIZ,mean_SOZ) +cd3 = computeCohen_d(mean_SOZ,mean_EIZ) %% combine univariate and bivariate %all_feat_zscores = [univariate_zscores, all_bivariate_feats]; all_feat_zscores(isinf(all_feat_zscores)) = 0; incomplete_channels = find(isnan(sum(all_feat_zscores'))); +complete_channels = find(~isnan(sum(all_feat_zscores'))); all_feat_zscores2 = all_feat_zscores; all_feat_zscores2(incomplete_channels,:) = []; EIZ = [zeros(1772,1);[spike_ind - spike_ind.*soz_ch]]; resect_ch_all = [zeros(1772,1);resected_ch]; % 1: patient, 2: abnormal, 3: EIZ, 4: SOZ, 5: resect, 6:outcome -all_possible_labels = [all_pts,[EIZ+[zeros(1772,1);soz_ch]],EIZ,[zeros(1772,1);soz_ch],resect_ch_all,[ones(1772,1);HUP_outcome_all]]; -multi_class_all = [~all_possible_labels(:,2)+2.*all_possible_labels(:,3)+3.*all_possible_labels(:,4)]; +all_possible_labels = [all_pts,[EIZ+[zeros(1772,1);soz_ch]],EIZ,[zeros(1772,1);soz_ch],resect_ch_all,[ones(1772,1);HUP_outcome_all],[[zeros(1772,1);soz_ch].*resect_ch_all]]; all_possible_labels(incomplete_channels,:) = []; -multi_class = [~all_possible_labels(:,2)+2.*all_possible_labels(:,3)+3.*all_possible_labels(:,4)]; - % get only non EIZ normal versus SOZ (in good outcome patients) -non_normal_non_soz = [zeros(1772,1);[soz_ch==0]]; -non_normal_non_soz(incomplete_channels) = []; -eliminate_channels = find(all_possible_labels(:,3)==1);%find(non_normal_non_soz==1); +eliminate_channels = [find(all_possible_labels(:,3)==1);find((all_possible_labels(:,4)+all_possible_labels(:,5))==1)];%find(non_normal_non_soz==1); +retain_channels = find(all_possible_labels(:,3)==0); all_feat_zscores3 = all_feat_zscores2; all_feat_zscores3(eliminate_channels,:) = []; new_labels = all_possible_labels; new_labels(eliminate_channels,:) = []; -all_feats = all_feat_zscores2; -all_labels = new_labels(:,4); +all_feats = all_feat_zscores3; +all_labels = new_labels(:,7); % was column 5 before all_pred_1 = zeros(size(all_labels)); all_pred_2 = zeros(size(all_labels)); all_pred_3 = zeros(size(all_labels)); -[part] = make_xval_partition(length(all_labels), 10); -for p = 1:10 +num_parts = 10; + +[pt_part] = make_xval_partition(166, num_parts); + +% cross-validate based on patients +for p = 1:num_parts + indices = find(pt_part==p) + for i = 1:length(indices) + part(find(indices(i)==new_labels(:,1))) = p; + + part2(find(indices(i)==all_possible_labels(:,1))) = p; + end +end + +% Define classes, 0 = non seizure, 1 = seizure +ClassNames = [0,1]; + +% Define costs -> penalty here is 5 +cost.ClassNames = ClassNames; +cost.ClassificationCosts = [0 1; 1 0]; + +for p = 1:num_parts p X_train = all_feats(find(part~=p),:); Y_train = all_labels(find(part~=p)); X_test = all_feats(find(part==p),:); Y_test = all_labels(find(part==p)); - Mdl1 = TreeBagger(100,X_train(:,1:10),Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); + X_analysis = all_feat_zscores2(find(part2==p),:); + + Mdl1 = TreeBagger(100,X_train(:,1:5),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); - Mdl2 = TreeBagger(100,X_train(:,11:20),Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); + Mdl2 = TreeBagger(100,X_train(:,6:10),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); Mdl3 = TreeBagger(100,X_train,Y_train,'Method',... - 'classification','PredictorSelection','curvature','OOBPredictorImportance','on'); + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); - [Yfit,pihat1] = predict(Mdl1,X_test(:,1:10)); - [Yfit,pihat2] = predict(Mdl2,X_test(:,11:20)); + [Yfit,pihat1] = predict(Mdl1,X_test(:,1:5)); + [Yfit,pihat2] = predict(Mdl2,X_test(:,6:10)); [Yfit,pihat3] = predict(Mdl3,X_test); + [Yfit, pihat_analysis] = predict(Mdl3,X_analysis); + [Yfit, pihat_analysis2] = predict(Mdl1,X_analysis(:,1:5)); + [Yfit, pihat_analysis3] = predict(Mdl2,X_analysis(:,6:10)); [~, Y_pred] = max(pihat1,[],2); all_pred_1(find(part==p)) = pihat1(:,2); all_pred_2(find(part==p)) = pihat2(:,2); all_pred_3(find(part==p)) = pihat3(:,2); + all_pred_analysis(find(part2==p),1) = pihat_analysis(:,2); + all_pred_analysis2(find(part2==p),1) = pihat_analysis2(:,2); + all_pred_analysis3(find(part2==p),1) = pihat_analysis3(:,2); xfold_acc(p) = sum(Y_pred==Y_test)./length(Y_test); Y_pred = round(pihat1(:,2)); - + end - [X1,Y1,T,AUC1] = perfcurve(all_labels,all_pred_1,1) + [X1,Y1,T,AUC1] = perfcurve(all_labels,all_pred_1,1); + [X2,Y2,T,AUC2] = perfcurve(all_labels,all_pred_2,1); + [X3,Y3,T,AUC3] = perfcurve(all_labels,all_pred_3,1); - [X2,Y2,T,AUC2] = perfcurve(all_labels,all_pred_2,1) - - [X3,Y3,T,AUC3] = perfcurve(all_labels,all_pred_3,1) + AUC1 + AUC2 + AUC3 + + Mdl3 = TreeBagger(100,all_feats,all_labels,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); -% random forest feature importance +%% random forest feature importance imp = Mdl3.OOBPermutedPredictorDeltaError; figure(1);clf; bar(imp); title('Curvature Test'); ylabel('Predictor importance estimates'); -xticks([1:20]) +xticks([1:10]) figure(2);clf; hold on -plot(X1,Y1,'LineWidth',2)%4 -plot(X2,Y2,'LineWidth',2)%7 -plot(X3,Y3,'LineWidth',2)%10 -legend(sprintf('Univariate - AUC: %f2',AUC1),sprintf('Bivariate - AUC: %f2',AUC2),sprintf('All - AUC: %f2',AUC3),'Location','SouthEast')%86/85/91 +plot(X1,Y1,'LineWidth',2,'Color',color4)%4 +plot(X2,Y2,'LineWidth',2,'Color',color7)%7 +plot(X3,Y3,'LineWidth',2,'Color',color1)%10 +legend(sprintf('Univariate - AUC: %.2f',AUC1),sprintf('Bivariate - AUC: %.2f',AUC2),sprintf('All - AUC: %.2f',AUC3),'Location','SouthEast')%86/85/91 xlabel('False positive rate') ylabel('True positive rate') title('All abnormal vs. normal channel classification') - +%% correlation +[r,p] = corr(all_feat_zscores2); +figure(1);clf; +subplot(1,2,1) +imagesc(r-eye(10)) +caxis([-1 1]) +colormap(color_bar) +colorbar +subplot(1,2,2) +imagesc(p<(0.05./100)) %% need to map back onto all channels and render -all_pred = ones(5203,1); +all_pred = zeros(5203,1); all_pred(incomplete_channels) = 0; -all_pred(all_pred==1) = 1-xfold_pred3; +all_pred(complete_channels) = all_pred_analysis; + +all_pred2 = zeros(5203,1); +all_pred2(incomplete_channels) = 0; +all_pred2(complete_channels) = all_pred_analysis2; + +all_pred3 = zeros(5203,1); +all_pred3(incomplete_channels) = 0; +all_pred3(complete_channels) = all_pred_analysis3; %% render abnormalities soz_ch_all = [zeros(1772,1);soz_ch]; -this_pt = 111; %114 res - 45/50 EIZ - 52/64, non-inv - 20,30 % 120 - res - 2., EIZ - 23./64, non-inv - 33/50. -this_pt_abn = all_pred_3(find(all_pts==this_pt)); -this_pt_soz = multi_class_all(find(all_pts==this_pt)); +this_pt = 116;%116/122/154 / 154 %114 res - 45/50 EIZ - 52/64, non-inv - 20,30 % 120 - res - 2., EIZ - 23./64, non-inv - 33/50. +this_pt_abn = all_pred(find(all_pts==this_pt)); +this_pt_soz = soz_ch_all(find(all_pts==this_pt)); +this_pt_eiz = EIZ(find(all_pts==this_pt)); this_pt_res = resect_ch_all(find(all_pts==this_pt)); this_pt_coords = all_coords(find(all_pts==this_pt),:); this_pt_eeg = all_wake_data(:,find(all_pts==this_pt)); %this_pt_delta_coh_abn(isnan(this_pt_delta_coh_abn)) = 0; +this_pt_multiclass = 2*this_pt_soz+this_pt_eiz; + this_pt_zscores = all_feat_zscores(find(all_pts==this_pt),:); -figure(1);clf; -hold on -for i = 1:10 - plot([1:12000],[this_pt_eeg(:,4*i)+i*175],'k-') -end -hold off +this_pt_test = zeros(size(this_pt_res)); +this_pt_test([28,9,49])=1; -figure(1);clf; -imagesc(all_feat_zscores(find(all_pts==this_pt),:)) -colorbar -caxis([0,5]) +% figure(1);clf; +% hold on +% for i = 1:10 +% plot([1:2000],[this_pt_eeg(:,4*i)+i*175],'k-') +% end +% hold off + +% figure(1);clf; +% imagesc(all_feat_zscores(find(all_pts==this_pt),:)) +% colorbar +% caxis([0,5]) figure(2);clf; subplot(1,3,1) -imagesc(this_pt_zscores(2,:)') -caxis([0,4]) +imagesc(this_pt_zscores(28,:)') +caxis([0,3]) +colormap('jet') +colorbar colorbar subplot(1,3,2) -imagesc(this_pt_zscores(23,:)') -caxis([0,4]) +imagesc(this_pt_zscores(9,:)') +caxis([0,3]) +colormap('jet') +colorbar colorbar subplot(1,3,3) -imagesc(this_pt_zscores(50,:)') -caxis([0,4]) +imagesc(this_pt_zscores(49,:)') +caxis([0,3]) +colormap('jet') colorbar - +colorbar +% figure(3); subplot(3,1,1) -plot(this_pt_eeg(:,2)) +plot(this_pt_eeg(6001:12000,4)) subplot(3,1,2) -plot(this_pt_eeg(:,23)) +plot(this_pt_eeg(6001:12000,23)) subplot(3,1,3) -plot(this_pt_eeg(:,50)) +plot(this_pt_eeg(6001:12000,27)) figure(4);clf; -final_elec_matrix = [this_pt_coords,this_pt_abn,ones(size(this_pt_coords,1),1)]; +final_elec_matrix = [this_pt_coords,this_pt_abn,ones(size(this_pt_res))]; %final_elec_matrix = final_elec_matrix([2,23,50],:) dlmwrite('render_elecs.node',final_elec_matrix,'delimiter',' ','precision',5) BrainNet_MapCfg('BrainMesh_ICBM152_smoothed.nv','render_elecs.node','abnormality_render.mat') -delete render_elecs.node +%delete render_elecs.node + +%% patient-specific correlation between abnormality features +for s = 107:166 + pt_zscores = all_feat_zscores2(find([all_possible_labels(:,1)==s]),:); + + pt_zscores_univar = pt_zscores(:,1:5); + pt_zscores_bivar = pt_zscores(:,6:10); + + corrval2(s-106) = corr(pt_zscores_univar(:),pt_zscores_bivar(:),'Type','Spearman'); + + for f = 1:5 + corrval((s-106),f) = corr(pt_zscores(:,f),pt_zscores(:,(f+5)),'Type','Spearman'); + end + + corrval_outcome((s-106),1) = late_outcome(s-106); + corrval_lesion((s-106),1) = lesion_field(s-106); +end + +p1 = ranksum(corrval(corrval_outcome==1,1),corrval(corrval_outcome>1,1),'tail','right') +p2 = ranksum(corrval(corrval_outcome==1,2),corrval(corrval_outcome>1,2),'tail','right') +p3 = ranksum(corrval(corrval_outcome==1,3),corrval(corrval_outcome>1,3),'tail','right') +p4 = ranksum(corrval(corrval_outcome==1,4),corrval(corrval_outcome>1,4),'tail','right') +p5 = ranksum(corrval(corrval_outcome==1,5),corrval(corrval_outcome>1,5),'tail','right') + +p6 = ranksum(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1),'tail','right') +d6 = computeCohen_d(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1)) -%% -single_bivariate_feat_2 = single_bivariate_feat; -new_roi_2 = new_roi; -new_roi_2(incomplete_channels) = []; -new_roi_2(new_roi_2==17) = 18; -for ch = 1:4702 - if multi_class(ch)==1 +plot_matrix = padcat(corrval2(corrval_outcome==1)',corrval2(corrval_outcome>1)'); +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylim([-0.3 0.5]) +%ylabel('Correlation of beta power & coherence |Z|') + +%% compute AUPRC/DRS at a per-patient level +a = 0; +for s = [107:127,128:132,134:139,141:155,157:161,163:166] + s + a = a+1; + this_pt_label = all_possible_labels(find(all_possible_labels(:,1)==s),4); + this_pt_pred = all_pred_analysis(find(all_possible_labels(:,1)==s)); + this_pt_pred2 = all_pred_analysis2(find(all_possible_labels(:,1)==s)); + this_pt_pred3 = all_pred_analysis3(find(all_possible_labels(:,1)==s)); + [X1,Y1,T,AUC1] = perfcurve(this_pt_label,this_pt_pred,1,'XCrit','TPR','YCrit','PPV'); + [X2,Y2,T,AUC2] = perfcurve(this_pt_label,this_pt_pred2,1,'XCrit','TPR','YCrit','PPV'); + [X3,Y3,T,AUC3] = perfcurve(this_pt_label,this_pt_pred3,1,'XCrit','TPR','YCrit','PPV'); + + auc_pt_outcome(a,1) = late_outcome(s-106); + pt_auc(a,1:3) = [AUC1,AUC2,AUC3]; +end + +good_inds = find(auc_pt_outcome==1); +poor_inds = find(auc_pt_outcome>1); + +plot_matrix = padcat(pt_auc(good_inds,3),pt_auc(poor_inds,3)); +figure(1);clf; +UnivarScatter(plot_matrix) + +ranksum(plot_matrix(:,1),plot_matrix(:,2)) +computeCohen_d(plot_matrix(:,1),plot_matrix(:,2)) +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylabel('Area under precision-recall curve') +%title('Comparison between AUPRC, rank-sum p = 0.033') + +%% +[single_univariate_feat, which_feat_uni] = max(all_feat_zscores(:,1:5),[],2); +[single_bivariate_feat, which_feat_bi] = max(all_feat_zscores(:,6:10),[],2); +figure(1);clf; +plot(abs(single_univariate_feat),single_bivariate_feat,'ko') +%% all-node abnormality score all regions +for ch = 1:4717 + if all_possible_labels(ch,2)==0 grouping{ch,1} = 'uninvolved'; class_ind(ch,1) = 1; - elseif multi_class(ch)==2 + elseif all_possible_labels(ch,3)==1 grouping{ch,1} = 'EIZ'; class_ind(ch,1) = 2; - else + elseif all_possible_labels(ch,4)==1 grouping{ch,1} = 'SOZ'; class_ind(ch,1) = 3; end end -univariate_feats = all_feat_zscores(:,1:10); -univariate_feats(incomplete_channels,:) = []; +new_roi_2 = new_roi; +new_roi_2(incomplete_channels) = []; + +single_bivariate_feat_2 = single_bivariate_feat; +single_bivariate_feat_2(incomplete_channels) = []; + +for f = 1:5 + for r = [1:16,17:20] + + new_roi_2(new_roi_2==(2*r-1)) = 2*r; + + pred_uninvolved = all_feat_zscores2(find((class_ind==1).*(new_roi_2==2*r)),(f+5)); + pred_eiz = all_feat_zscores2(find((class_ind==2).*(new_roi_2==2*r)),(f+5)); + pred_soz = all_feat_zscores2(find((class_ind==3).*(new_roi_2==2*r)),(f+5)); + + median(pred_uninvolved) + median(pred_eiz) + median(pred_soz) + + plot_matrix = padcat(pred_uninvolved, pred_eiz,pred_soz); + + p1_region(f,r) = ranksum(pred_uninvolved,pred_eiz) + p2_region(f,r) = ranksum(pred_uninvolved,pred_soz) + p3_region(f,r) = ranksum(pred_eiz,pred_soz) + + cohen_d2(f,r) = computeCohen_d(pred_uninvolved,pred_soz) + cohen_d3(f,r) = computeCohen_d(pred_eiz,pred_soz) + + end +end -single_univariate_feat = prctile(abs(univariate_feats'),100)'; +figure(1);clf; +subplot(3,1,1) +imagesc([p3_region]) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +colorbar +subplot(3,1,2) +imagesc([p3_region]<0.05) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +subplot(3,1,3) +imagesc([cohen_d3]) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +caxis([0,1]) +colorbar -single_univariate_feat(isinf(single_univariate_feat)) = nanmedian(single_univariate_feat); -single_bivariate_feat_2(isinf(single_bivariate_feat_2)) = nanmedian(single_bivariate_feat_2); +%% patient-level abnormality score all regions -single_univariate_feat = single_univariate_feat(new_roi_2==18); -single_bivariate_feat_2 = single_bivariate_feat_2(new_roi_2==18); -grouping = grouping(new_roi_2==18); -class_ind = class_ind(new_roi_2==18); + pt_uninvolved = []; + pt_eiz = []; + pt_soz = []; + +for s = 107:166 + + pt_uninvolved_inds = find([all_possible_labels(:,2)==0].*[all_possible_labels(:,1)==s]);%.*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); + pt_eiz_inds = find([all_possible_labels(:,3)==1].*[all_possible_labels(:,1)==s]);%.*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); + pt_soz_inds = find([all_possible_labels(:,4)==1].*[all_possible_labels(:,1)==s]);%.*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); + + pt_uninvolved = [pt_uninvolved; nanmean(all_pred_analysis(pt_uninvolved_inds))]; + pt_eiz = [pt_eiz; nanmean(all_pred_analysis(pt_eiz_inds))]; + pt_soz = [pt_soz; nanmean(all_pred_analysis(pt_soz_inds))]; +end +plot_matrix = padcat(pt_uninvolved,pt_eiz,pt_soz); figure(1);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat_2,'Group',grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,grouping,'orientation','horizontal',... - 'label',{'','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat_2,grouping,'orientation','horizontal',... - 'label', {'','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) +plot_cell{1,1} = rmmissing(pt_uninvolved); +plot_cell{1,2} = rmmissing(pt_eiz); +plot_cell{1,3} = rmmissing(pt_soz); +violin(plot_cell,'xlabel',{'','',''},'facecolor',[color1;color7;color2],'mc',[],'medc','k');%,'edgecolor','k'); +xlim([0.5 3.5]) +ylim([0 1]) -d1 = computeCohen_d(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) -p1 = ranksum(single_univariate_feat(find(class_ind==2)),single_univariate_feat(find(class_ind==1))) -d2 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) -p2 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==2))) -d3 = computeCohen_d(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) -p3 = ranksum(single_univariate_feat(find(class_ind==3)),single_univariate_feat(find(class_ind==1))) - -d4 = computeCohen_d(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) -p4 = ranksum(single_bivariate_feat_2(find(class_ind==2)),single_bivariate_feat_2(find(class_ind==1))) -d5 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) -p5 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==2))) -d6 = computeCohen_d(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) -p6 = ranksum(single_bivariate_feat_2(find(class_ind==3)),single_bivariate_feat_2(find(class_ind==1))) - -[d1 d2 d3 d4 d5 d6] -[p1 p2 p3 p4 p5 p6] -%% do node abormality (max univariate+bivariate in plot) for following: -% 1. within/outside soz for lesional / non-lesional - -all_lesion = zeros(5203,1); -for i = 1:166 - if i>106 - all_lesion(all_pts==i) = lesion_field(i-106); - end -end - -within_soz_lesional = find(all_lesion.*soz_ch_all); -within_soz_nonlesional = find([all_lesion==0].*soz_ch_all); -outside_soz_lesional = find(all_lesion.*[soz_ch_all==0]); -outside_soz_nonlesional = find([all_lesion==0].*[soz_ch_all==0]); - -for ch = 1:5203 - if ismember(ch,within_soz_lesional) - lesional_grouping{ch,1} = 'within_soz_lesional'; - elseif ismember(ch,within_soz_nonlesional) - lesional_grouping{ch,1} = 'within_soz_nonlesional'; - elseif ismember(ch,outside_soz_lesional) - lesional_grouping{ch,1} = 'outside_soz_lesional'; - elseif ismember(ch,outside_soz_nonlesional) - lesional_grouping{ch,1} = 'outside_soz_nonlesional'; +figure(2);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','irritative zone','seizure onset zone'}) +ylabel('predicted epileptogenicity score') +ylim([0 0.8]) + +p1 = ranksum(pt_uninvolved,pt_eiz) +d1 = computeCohen_d(pt_uninvolved,pt_eiz) +p2 = ranksum(pt_eiz,pt_soz) +d2 = computeCohen_d(pt_eiz,pt_soz) +p3 = ranksum(pt_uninvolved,pt_soz) +d3 = computeCohen_d(pt_uninvolved,pt_soz) + +%% patient-level abnormality score res in/out good/poor outcome +good_rz = []; +good_out = []; + +poor_rz = []; +poor_out = []; + +for s = 107:166 + this_pt_outcome = late_outcome(s-106); + this_pt_ind_rz = find([all_possible_labels(:,1)==s].*[all_possible_labels(:,5)==1]); + this_pt_ind_out = find([all_possible_labels(:,1)==s].*[all_possible_labels(:,5)==0]); + + if this_pt_outcome==1 + good_rz = [good_rz; nanmedian(all_pred_analysis(this_pt_ind_rz))]; + good_out = [good_out; nanmedian(all_pred_analysis(this_pt_ind_out))]; + else + poor_rz = [poor_rz; nanmedian(all_pred_analysis(this_pt_ind_rz))]; + poor_out = [poor_out; nanmedian(all_pred_analysis(this_pt_ind_out))]; end end -univariate_feats = all_feat_zscores(:,1:10); -single_univariate_feat = prctile(abs(univariate_feats'),100)'; -single_bivariate_feat = prctile(abs_bivariate_feats',100)'; - -uni_within_soz_lesional = single_univariate_feat(within_soz_lesional); -uni_within_soz_nonlesional = single_univariate_feat(within_soz_nonlesional); -uni_outside_soz_lesional = single_univariate_feat(outside_soz_lesional); -uni_outside_soz_nonlesional = single_univariate_feat(outside_soz_nonlesional); -bi_within_soz_lesional = single_bivariate_feat(within_soz_lesional); -bi_within_soz_nonlesional = single_bivariate_feat(within_soz_nonlesional); -bi_outside_soz_lesional = single_bivariate_feat(outside_soz_lesional); -bi_outside_soz_nonlesional = single_bivariate_feat(outside_soz_nonlesional); - -p1 = ranksum(uni_within_soz_lesional,uni_within_soz_nonlesional); -p2 = ranksum(uni_within_soz_lesional,uni_outside_soz_lesional); -p3 = ranksum(uni_within_soz_nonlesional,uni_outside_soz_nonlesional); -p4 = ranksum(uni_outside_soz_lesional,uni_outside_soz_nonlesional); -p5 = ranksum(bi_within_soz_lesional,bi_within_soz_nonlesional); -p6 = ranksum(bi_within_soz_lesional,bi_outside_soz_lesional); -p7 = ranksum(bi_within_soz_nonlesional,bi_outside_soz_nonlesional); -p8 = ranksum(bi_outside_soz_lesional,bi_outside_soz_nonlesional); - -single_univariate_feat(incomplete_channels) = []; -single_bivariate_feat(incomplete_channels) = []; -lesional_grouping(incomplete_channels) = []; - +plot_matrix = padcat(good_rz,good_out,poor_rz,poor_out); figure(2);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',lesional_grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,lesional_grouping,'orientation','horizontal',... - 'label',{'','','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat,lesional_grouping,'orientation','horizontal',... - 'label', {'','','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) +UnivarScatter(plot_matrix) +%% patient-level abnormality score soz in/out lesional -[p1 p2 p3 p4 p5 p6 p7 p8] -%% 2. within/outside ablation zone for good / poor outcome -all_outcome = ones(5203,1); -all_abl = NaN(5203,1); -for i = 1:166 - if i>106 - all_outcome(all_pts==i) = late_outcome(i-106); - if therapy_field(i-106) - all_abl(all_pts==i) = resect_ch_all(all_pts==i); - end - end -end +single_univariate_feat_2 = median(univariate_feats')'; -within_abl_good = find([all_outcome==1].*[all_abl==1]); -outside_abl_good = find([all_outcome==1].*[all_abl==0]); -within_abl_poor = find([all_outcome>1].*[all_abl==1]); -outside_abl_poor = find([all_outcome>1].*[all_abl==0]); +pt_in_rz = []; +pt_out_rz = []; -incomplete_channels2 = incomplete_channels; +pt_in_rz2 = []; +pt_out_rz2 = []; -abl_grouping = cell(5203,1); +for s = 107:166 + + w_inds = find([all_possible_labels(:,1)==s])%;.*[all_possible_labels(:,4)==1]); + o_inds = find([all_possible_labels(:,1)==s])%;.*[all_possible_labels(:,4)==0]); + + length_diagnosis = age_surgery - age_onset; + + pt_in_rz = [pt_in_rz; nanmedian(single_bivariate_feat_2(w_inds))]; + pt_out_rz = [pt_out_rz; nanmedian(single_bivariate_feat_2(o_inds))]; + + pt_in_rz2 = [pt_in_rz2; nanmedian(single_univariate_feat_2(w_inds))]; + pt_out_rz2 = [pt_out_rz2; nanmedian(single_univariate_feat_2(o_inds))]; -for ch = 1:5203 - if ismember(ch,within_abl_good) - abl_grouping{ch,1} = 'within_abl_good'; - elseif ismember(ch,outside_abl_good) - abl_grouping{ch,1} = 'outside_abl_good'; - elseif ismember(ch,within_abl_poor) - abl_grouping{ch,1} = 'within_abl_poor'; - elseif ismember(ch,outside_abl_poor) - abl_grouping{ch,1} = 'outside_abl_poor'; - else - incomplete_channels2 = [incomplete_channels2, ch]; - end end -incomplete_channels2 = unique(incomplete_channels2); - -univariate_feats = all_feat_zscores(:,1:10); -single_univariate_feat = prctile(abs(univariate_feats'),100)'; -single_bivariate_feat = prctile(abs_bivariate_feats',100)'; - -uni_within_abl_good = single_univariate_feat(within_abl_good); -uni_outside_abl_good = single_univariate_feat(outside_abl_good); -uni_within_abl_poor = single_univariate_feat(within_abl_poor); -uni_outside_abl_poor = single_univariate_feat(outside_abl_poor); -bi_within_abl_good = single_bivariate_feat(within_abl_good); -bi_outside_abl_good = single_bivariate_feat(outside_abl_good); -bi_within_abl_poor = single_bivariate_feat(within_abl_poor); -bi_outside_abl_poor = single_bivariate_feat(outside_abl_poor); - -p1 = ranksum(uni_within_abl_good,uni_outside_abl_good); -p2 = ranksum(uni_within_abl_poor,uni_outside_abl_poor); -p3 = ranksum(uni_within_abl_good,uni_within_abl_poor); -p4 = ranksum(uni_outside_abl_good,uni_outside_abl_poor); -p5 = ranksum(bi_within_abl_good,bi_outside_abl_good); -p6 = ranksum(bi_within_abl_poor,bi_outside_abl_poor); -p7 = ranksum(bi_within_abl_good,bi_within_abl_poor); -p8 = ranksum(bi_outside_abl_good,bi_outside_abl_poor); - -single_univariate_feat(incomplete_channels2) = []; -single_bivariate_feat(incomplete_channels2) = []; -abl_grouping(incomplete_channels2) = []; - -figure(3);clf; -h = scatterhist(single_univariate_feat,single_bivariate_feat,'Group',abl_grouping,'Marker','.','MarkerSize',12) - -clr = get(h(1),'colororder'); -boxplot(h(2),single_univariate_feat,abl_grouping,'orientation','horizontal',... - 'label',{'','','',''},'color',clr); - xlim(h(2),[0,30]) -boxplot(h(3),single_bivariate_feat,abl_grouping,'orientation','horizontal',... - 'label', {'','','',''},'color',clr); - xlim(h(3),[0,30]) - set(h(2:3),'XTickLabel',''); -view(h(3),[270,90]); % Rotate the Y plot -axis(h(1),'auto'); -xlim([0 10]) -ylim([0 10]) +figure(1);clf; +hold on +plot(length_diagnosis, pt_out_rz,'b.','Color',[106 178 180]/255,'MarkerSize',12) +p = polyfit(length_diagnosis,pt_out_rz,1); +f = polyval(p,length_diagnosis); +plot(length_diagnosis,f,'b-','LineWidth',2) +[a2,b2] = corr(length_diagnosis, pt_out_rz,'type','Pearson','tail','right') +xlim([0 45]) +ylim([0.6 2.6]) -[p1 p2 p3 p4 p5 p6 p7 p8] +[a3,b3] = corr(length_diagnosis, pt_out_rz2,'type','Pearson','tail','right') %% Examine correlations between z scores [r,p] = corr(all_feat_zscores2) figure(1);clf; +r(p>0.0001) = NaN; imagesc(r) colorbar - diff --git a/iEEG_atlas_revision.m b/iEEG_atlas_revision.m new file mode 100644 index 0000000..914f4df --- /dev/null +++ b/iEEG_atlas_revision.m @@ -0,0 +1,1216 @@ + %% Set up workspace +clear all + +% set up path - change iEEG_atlas_path to where you download the repository +iEEG_atlas_path = '/Users/jbernabei/Documents/PhD_Research/atlas_project/iEEG_atlas_dev'; + +% other paths may stay the same +metadata = readtable('data/atlas_metadata_final.xlsx'); +custom_atlas = readtable('data/custom_atlas.xlsx'); + +% anatomical atlas information +region_list = zeros(1,90); % for the 90 AAL regions we will be using +region_names = cell(1,90); +fi = fopen("localization/AAL116_WM.txt"); +for j = 1:90 + label = split(fgetl(fi)); + region_list(j) = str2double(label{3}); + region_names{j} = label{2}; +end + +% extract metadata +age_onset = metadata{:,11}; +age_surgery = metadata{:,12}; +HUP_outcome = floor(metadata{:,3}); + +% extract outcome +for s = 1:length(age_surgery) + engel_scores = metadata{s,3:5}; + engel_scores = rmmissing(engel_scores); + early_outcome(s) = floor(engel_scores(1)); % early outcome is 6 months + late_outcome(s) = floor(engel_scores(end)); % late outcome is latest time point up to 24 months + late_outcome_raw(s) = engel_scores(end); +end + +therapy_field = metadata{:,6}; +implant_field = metadata{:,7}; +target_field = metadata{:,8}; +laterality_field = metadata{:,9}; +lesion_field = metadata{:,10}; +gender_field = metadata{:,13}; + +lesion_field = strcmp(lesion_field,'Lesional'); % 1 for lesional +therapy_field = strcmp(therapy_field,'Ablation'); % 1 ffor ablation + +load('data/MNI_atlas_new.mat') +load('data/HUP_atlas_final.mat') + +all_pts = [Patient; (patient_no+106)]; % combine patient numbers for MNI & HUP + +HUP_outcome_all = zeros(length(soz_ch),1); +for i = 1:60 + HUP_outcome_all(find(patient_no==i)) = late_outcome(i); % we use outcome as late outcome +end + +clear NodesLeft +clear NodesLeftInflated +clear NodesRightInflated +clear NodesRight +clear NodesRegionLeft +clear NodesRegionRight +clear Data_N2 +clear Data_N3 +clear Data_R + +% some colors for plotting, etc +color1 = [0, 0.4470, 0.7410]; +color2 = [0.6350, 0.0780, 0.1840]; +color3 = [255, 248, 209]./255; +color4 = [103 55 155]/255; +color6 = [78 172 91]/255; +color7 = [0.9290, 0.6940, 0.1250]; +color8 = [0.8500, 0.3250, 0.0980]; + +% assign wake clips from MNI and HUP into a final data matrix +all_wake_data = [Data_W,wake_clip]; + +%% perform localization based on AAL116 atlas (which has an internal WM mask also) +all_coords = [ChannelPosition;mni_coords]; + +try [coords, all_roi, NN_ind] = nifti_values(all_coords,'localization/AAL116_WM.nii'); +catch ME + for i = 1:size(all_coords,1) + if mod(i,100) + fprintf('%d\n',i) + end + try [mni_coords1, mni_roi, NN_ind] = nifti_values(all_coords(i,:),'localization/AAL116_WM.nii'); + catch ME + mni_roi = NaN; + end + all_roi(i,1) = mni_roi; + end +end + +clear NN_ind +clear coords + +% assign AAL regions into new, aggregate regions of custom atlas. (116 -> 40 regions) +region_matrix = custom_atlas{:,2:4}; +for c = 1:length(all_roi) + old_ind = find(region_list==all_roi(c)); + try [row,col,V] = find(region_matrix==old_ind); + try new_roi(c,1) = row; + catch anyerror + new_roi(c,1) = 0; + end + catch anyerror + new_roi(c,1) = 0; + end + +end +%% define normal MNI/HUP and abnormal channels (SOZ/high-spike/RZ) +% set threshold on spike counts +spike_thresh = 24; % this is empirical, 1 spike/hour + +% find indices for which spikes are greater than threshold +spike_ind = [spike_24h>spike_thresh]; + +% define all abnormal channels +abnormal_ch = find([spike_ind+soz_ch+resected_ch+(HUP_outcome_all>1)]>0)+1772; % w/o resected ch + +% define all seizure onset indices +soz_ch_inds = find(soz_ch)+1772; + +% define normal HUP channels +normal_HUP_ch = find([spike_ind+soz_ch+resected_ch+(HUP_outcome_all>1)]==0)+1772; % w/o resected ch + +% define normal MNI channels +normal_MNI_ch = [1:1772]'; + +% define all normal channels +all_normal_ch = [normal_MNI_ch;normal_HUP_ch]; + +% define exclusive irritative zone +EIZ_ch = find(spike_ind - spike_ind.*soz_ch)+1772; + +%% count HUP and MNI in original atlas and then new atlas +mni_roi_old = all_roi(1:1772); +hup_roi_old = all_roi(1773:end); + +for i = 1:45 + roi_count(i,1) = sum(mni_roi_old==region_list(2*i-1))+sum(mni_roi_old==region_list(2*i)); + roi_count(i,2) = sum(all_roi(normal_HUP_ch)==region_list(2*i-1))+sum(all_roi(normal_HUP_ch)==region_list(2*i)); + this_roi_name = split(string(region_names{2*i}),'_R'); % extract roi name without laterality + split_roi_name_old{i,1} = this_roi_name(1); +end + +for i = 1:20 + roi_count_new_normal(i,1) = length(intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); + roi_count_new_normal(i,2) = length(intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)])); + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); +end + +new_counts = [roi_count_new_normal(:,1),(roi_count_new_normal(:,2))]; + +% Supplemental figure +figure(1);clf; +subplot(1,2,1) +barh(roi_count) +set(gca,'TickLabelInterpreter', 'none'); +xlabel('Number of nodes per region') +title('AAL regions') +yticks([1:45]) +yticklabels(split_roi_name_old) +subplot(1,2,2) +barh(new_counts) +set(gca,'TickLabelInterpreter', 'none'); +title('Composite regions') +yticks([1:21]) +yticklabels(split_roi_name) +xlabel('Number of nodes per region') + +%% calculate power spectrum +all_wake_data = [Data_W,wake_clip]; + +[pxx,f] = pwelch(all_wake_data,400,200,[0.5:0.5:80],200); % original was 200/100 +pxx_norm = pxx./sum(pxx); +%pxx_norm = pxx; +%% Analyze power spectral density between HUP/MNI, Normal/irritative zone/SOZ (Figure 4) +% change this to be on a per-patient level rather than a per-node level + +freq_inds = [1, 8; % 0.5 - 4 Hz + 9, 16; % 4 - 8 Hz + 17, 26; % 8 - 13 Hz + 27, 60; % 13 - 30 Hz + 61, 160]; % 30 - 80 Hz + + +figure(1);clf % supplemental +figure(2);clf % supplemental +for i = 1:20 + % extract indices that correspond to a given ROI for HUP and MNI data + % (combining data from L and R hemispheres) + roi_mni = intersect(normal_MNI_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_hup = intersect(normal_HUP_ch,[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_EIZ = intersect([EIZ_ch],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + roi_soz = intersect([soz_ch_inds],[find(new_roi==(2*(i-1)+1));find(new_roi==2*i)]); + + % median power spectral density across channels + median_mni = median(pxx_norm(:,roi_mni),2); + median_hup = median(pxx_norm(:,roi_hup),2); + median_normal = median(pxx_norm(:,[roi_mni;roi_hup]),2); + median_EIZ = median(pxx_norm(:,roi_EIZ),2); + median_soz = median(pxx_norm(:,roi_soz),2); + + for j = 1:5 + HUP_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_hup))'; + MNI_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_mni))'; + normal_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],[roi_mni;roi_hup]))'; + EIZ_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_EIZ))'; + soz_feat = median(pxx_norm([freq_inds(j,1):freq_inds(j,2)],roi_soz))'; + + [p,h] = ranksum(HUP_feat,MNI_feat); % HUP vs MNI + pval_matrix_norm(i,j) = p; + + [p,h] = ranksum(normal_feat,EIZ_feat); % normal vs EIZ + pval_matrix_EIZ(i,j) = p; + + [p,h] = ranksum(normal_feat,soz_feat); % normal vs SOZ + pval_matrix_soz(i,j) = p; + + [p,h] = ranksum(EIZ_feat,soz_feat); % EIZ vs SOZ + pval_matrix_EIZ_vs_soz(i,j) = p; + + end + + % find upper and lower bounds (here 25th and 75th percentiles) + prct_25_75_mni = prctile(pxx_norm(:,roi_mni)',[25,75]); + prct_25_75_hup = prctile(pxx_norm(:,roi_hup)',[25,75]); + prct_25_75_normal = prctile(pxx_norm(:,[roi_mni;roi_hup])',[25,75]); + prct_25_75_EIZ = prctile(pxx_norm(:,roi_EIZ)',[25,75]); + prct_25_75_soz = prctile(pxx_norm(:,roi_soz)',[25,75]); + + x_axis = [0.5:0.5:80]; + + figure(1) + % 3 x 7 plot for 21 regions (aggregating across hemispheres) + subplot(4,5,i) + hold on + plot([4,4],[0,0.4],'k-') + plot([8,8],[0,0.4],'k-') + plot([13,13],[0,0.4],'k-') + plot([30,30],[0,0.4],'k-') + plot(x_axis,median_mni) + patch([x_axis fliplr(x_axis)], [prct_25_75_mni(1,:) fliplr(prct_25_75_mni(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') + plot(x_axis,median_hup) + patch([x_axis fliplr(x_axis)], [prct_25_75_hup(1,:) fliplr(prct_25_75_hup(2,:))], color6, 'FaceAlpha',0.5, 'EdgeColor','none') + + % add greek letters for frequency bands + txt_d = '\delta'; + txt_t = '\theta'; + txt_a = '\alpha'; + txt_b = '\beta'; + txt_g = '\gamma'; + txt_sig = '*'; + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); + if pval_matrix_norm(i,1)<(0.05/100) + text(2,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,2)<(0.05/100) + text(5.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,3)<(0.05/100) + text(9.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,4)<(0.05/100) + text(18,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_norm(i,5)<(0.05/100) + text(40,0.12,txt_sig,'FontSize', 20); + end + + % set axes, labels, subgraph titles + set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) + ylabel('Spectral Density') + xlabel('Frequency (Hz)') + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); + title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title + xlim([0.5,80]) % xlimit 0.5-80 Hz + ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) + hold off + + figure(2) + subplot(4,5,i) + hold on + plot([4,4],[0,0.4],'k-') + plot([8,8],[0,0.4],'k-') + plot([13,13],[0,0.4],'k-') + plot([30,30],[0,0.4],'k-') + plot(x_axis,median_normal) + patch([x_axis fliplr(x_axis)], [prct_25_75_normal(1,:) fliplr(prct_25_75_normal(2,:))], color1, 'FaceAlpha',0.5, 'EdgeColor','none') + plot(x_axis,median_EIZ,'Color',color7) + patch([x_axis fliplr(x_axis)], [prct_25_75_EIZ(1,:) fliplr(prct_25_75_EIZ(2,:))], color7, 'FaceAlpha',0.5, 'EdgeColor','none') + try plot(x_axis,median_soz,'Color',color2) + patch([x_axis fliplr(x_axis)], [prct_25_75_soz(1,:) fliplr(prct_25_75_soz(2,:))], color2, 'FaceAlpha',0.5, 'EdgeColor','none') + catch anyerror + end + + % add greek letters for frequency bands + text(2,0.15,txt_d,'FontSize', 12); + text(5.5,0.15,txt_t,'FontSize', 12); + text(9.5,0.15,txt_a,'FontSize', 12); + text(18,0.15,txt_b,'FontSize', 12); + text(40,0.15,txt_g,'FontSize', 12); + if pval_matrix_EIZ(i,1)<(0.05/100) + text(2,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,2)<(0.05/100) + text(5.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,3)<(0.05/100) + text(9.5,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,4)<(0.05/100) + text(18,0.12,txt_sig,'FontSize', 20); + end + if pval_matrix_EIZ(i,5)<(0.05/100) + text(40,0.12,txt_sig,'FontSize', 20); + end + + if pval_matrix_soz(i,1)<(0.05/100) + text(2,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,2)<(0.05/100) + text(5.5,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,3)<(0.05/100) + text(9.5,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,4)<(0.05/100) + text(18,0.17,txt_sig,'FontSize', 20); + end + if pval_matrix_soz(i,5)<(0.05/100) + text(40,0.17,txt_sig,'FontSize', 20); + end + + % set axes, labels, subgraph titles + set(gca, 'XScale','log', 'XTick',[0.5,4,8,13,30,80], 'XTickLabel',[0.5,4,8,13,30,80]) + ylabel('Spectral Density') + xlabel('Frequency (Hz)') + this_roi_name = split(string(custom_atlas{2*i,1}),'_R'); % extract roi name without laterality + split_roi_name{i,1} = this_roi_name(1); + title(sprintf('%s',this_roi_name(1)), 'Interpreter', 'none') % use that as title + xlim([0.5,80]) % xlimit 0.5-80 Hz + ylim([0 0.2]) % ylimit 0-0.2 (arbitrary, just to visualize density) + hold off + +end + +%% do paired boxplot of channel-wise and patient-wise distributions for +% one MTL and one cortical activity +activity_pt_1 = []; +activity_ch_1 = []; +activity_pt_2 = []; +activity_ch_2 = []; +activity_pt_3 = []; +activity_ch_3 = []; + +for s = 1:166 + chs_1 = find([all_pts==s].*([new_roi==17]+[new_roi==18])); + chs_2 = find([all_pts==s].*([new_roi==37]+[new_roi==38])); + chs_3 = find([all_pts==s].*([new_roi==7]+[new_roi==8])); + + pt1 = [mean(mean(pxx_norm(1:8,chs_1))),mean(mean(pxx_norm(9:16,chs_1))),mean(mean(pxx_norm(17:26,chs_1))),mean(mean(pxx_norm(27:60,chs_1))),mean(mean(pxx_norm(61:160,chs_1)))]; + ch1 = [mean(pxx_norm(1:8,chs_1))',mean(pxx_norm(9:16,chs_1))',mean(pxx_norm(17:26,chs_1))',mean(pxx_norm(27:60,chs_1))',mean(pxx_norm(61:160,chs_1))']; + pt2 = [mean(mean(pxx_norm(1:8,chs_2))),mean(mean(pxx_norm(9:16,chs_2))),mean(mean(pxx_norm(17:26,chs_2))),mean(mean(pxx_norm(27:60,chs_2))),mean(mean(pxx_norm(61:160,chs_2)))]; + ch2 = [mean(pxx_norm(1:8,chs_2))',mean(pxx_norm(9:16,chs_2))',mean(pxx_norm(17:26,chs_2))',mean(pxx_norm(27:60,chs_2))',mean(pxx_norm(61:160,chs_2))']; + pt3 = [mean(mean(pxx_norm(1:8,chs_3))),mean(mean(pxx_norm(9:16,chs_3))),mean(mean(pxx_norm(17:26,chs_3))),mean(mean(pxx_norm(27:60,chs_3))),mean(mean(pxx_norm(61:160,chs_3)))]; + ch3 = [mean(pxx_norm(1:8,chs_3))',mean(pxx_norm(9:16,chs_3))',mean(pxx_norm(17:26,chs_3))',mean(pxx_norm(27:60,chs_3))',mean(pxx_norm(61:160,chs_3))']; + + try activity_pt_1 = [activity_pt_1; pt1]; + catch + end + + try activity_ch_1 = [activity_ch_1; ch1]; + catch + end + + try activity_pt_2 = [activity_pt_2; pt2]; + catch + end + + try activity_ch_2 = [activity_ch_2; ch2]; + catch + end + + try activity_pt_3 = [activity_pt_3; pt3]; + catch + end + + try activity_ch_3 = [activity_ch_3; ch3]; + catch + end + +end + +colors = [color1; color2; color1; color2; color1; color2; color1; color2; color1; color2]; + +plot_data1 = padcat(activity_pt_1(:,1),activity_ch_1(:,1),activity_pt_1(:,2),activity_ch_1(:,2),activity_pt_1(:,3),activity_ch_1(:,3),activity_pt_1(:,4),activity_ch_1(:,4),activity_pt_1(:,5),activity_ch_1(:,5)); +plot_data2 = padcat(activity_pt_2(:,1),activity_ch_2(:,1),activity_pt_2(:,2),activity_ch_2(:,2),activity_pt_2(:,3),activity_ch_2(:,3),activity_pt_2(:,4),activity_ch_2(:,4),activity_pt_2(:,5),activity_ch_2(:,5)); +plot_data3 = padcat(activity_pt_3(:,1),activity_ch_3(:,1),activity_pt_3(:,2),activity_ch_3(:,2),activity_pt_3(:,3),activity_ch_3(:,3),activity_pt_3(:,4),activity_ch_3(:,4),activity_pt_3(:,5),activity_ch_3(:,5)); + +figure(1);clf; +subplot(1,3,1) +boxplot(plot_data1) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Normalized spectral density') + +subplot(1,3,2) +boxplot(plot_data2) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Normalized spectral density') + +subplot(1,3,3) +boxplot(plot_data3) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Normalized spectral density') +%% calculate univariate z-scores of neural activity +% power spectrum and entropy abnormality + +all_f = {'delta','theta','alpha','beta','gamma'}; % frequency bands +all_m = {'bandpower'};%,'entropy'}; % metrics + +c = 0; % c for column in z score matrix +for m = 1 % m for different metrics + for f = 1:5 % f for frequency bands + c = c+1; + + % first create the univariate atlas, taking only normal channels + % the univariate atlas is mean/median + stdev of a given feature + % across ROI + [mean_vec,median_vec,std_vec,feat_vals] = create_univariate_atlas(all_wake_data(:,find(all_normal_ch)),[1:40],new_roi(all_normal_ch),sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); + + % we call this function again on all of the data, but are only + % using it to extract all of the feature values on each channel + [~,~,~,feat_vals] = create_univariate_atlas(all_wake_data,[1:40],new_roi,sprintf('%s',all_m{m}),sprintf('%s',all_f{f})); + + % now we test the univariate atlas to generate z scores + [zscores, pt_zscores, pt_ind, roi_ind] = test_univariate_atlas(median_vec, std_vec, feat_vals, all_pts, [1:40], new_roi,'patient'); + + % store the z scores + univariate_zscores(:,c) = zscores; + all_pt_zscores(:,c) = pt_zscores; + + end +end + +%save('univariate_zscores.mat','univariate_zscores') + + +%% create atlas data - > edge-wise and patient-wise should be options +% for HUP data, aggregate electrodes to ignore +% these should be poor-outcome/RZ/SOZ/>1 spike per hour + +conn_band = {'delta','theta','alpha','beta','gamma'}; +conn_type = {'corr','coh'}; + +% loop through types of connectivity and frequency bands +for c = 2 + for f = 2:5 + + % loop through subjects + for s = 107:max(all_pts) + s + + % extract subject data from all the data (columns are channels) + pt_data = all_wake_data(:,[all_pts==s]); + + % call functional connectivity code (200 is Hz, 1 is window size in seconds) + [pt_adj] = compute_FC(pt_data, 200, 1, conn_type{c},conn_band{f}); + + % assign median adjacency matrix into cell structure + all_pt_adj{s} = pt_adj; + + end + + % save adjacency matrices + %save(sprintf('adj_matrices_rev2/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f}),'all_pt_adj') + + end +end + +%% assign atlas locations per patient into cell structure +for s = 1:max(all_pts) + pt_loc{s} = new_roi(all_pts==s)'; +end + +% indices of all abnormal channels +all_abnormal_ch = zeros(5203,1); +all_abnormal_ch(abnormal_ch) = 1; + +for s = 1:166 + + % check abnormal channels for HUP patients + if s>106 + + % abnormal channel indices + all_abn{s} = find(all_abnormal_ch(all_pts==s)); + if isempty(all_abn{s}) + all_abn{s} = []; + end + else + % if MNI patient, no abnormal channels by definition + all_abn{s} = []; + end +end + +%% functional cconnectivity atlas construction & z-score code + +conn_band = {'delta','theta','alpha','beta','gamma'}; +conn_type = {'coh'}; + +% which feature +feat = 0; + +% loop through connectivity and frequency +for c = 1 + for f = 1:5 + + feat = feat+1; + a = 0 + + load(sprintf('adj_matrices/all_adj_%s_%s_May25.mat',conn_type{c},conn_band{f})); + + if c==2 + for s = 1:166 + all_pt_adj{s} = log(all_pt_adj{s}); + end + end + + [conn_edge, std_edge, samples_edge, sem_edge, raw_atlas_edge] = create_atlas_by_edge(all_pt_adj, pt_loc, all_abn, [1:40]', 2); + + [median_conn, std_conn, num_samples, sem_conn, raw_atlas] = create_atlas(all_pt_adj, pt_loc, all_abn, [1:40]', 2); + + % call for MNI patients + [conn_edge_mni, std_edge_mni, samples_edge_mni, ~, raw_atlas_edge_mni] = create_atlas_by_edge(all_pt_adj(1:106), pt_loc(1:106), all_abn(1:106), [1:40]', 1); + + % call for HUP patients + [conn_edge_hup, std_edge_hup, samples_edge_hup, ~, raw_atlas_edge_hup] = create_atlas_by_edge(all_pt_adj(107:166), pt_loc(107:166), all_abn(107:166), [1:40]', 1); + + this_feat_conn(feat).mni = conn_edge_mni; + this_feat_conn(feat).hup = conn_edge_hup; + + all_raw_atlas_edge{f} = raw_atlas_edge; + all_raw_atlas{f} = raw_atlas; + + % loop through all patients + for s = 1:166 % change this back + a = a+1; + + % do nodal space + try [native_adj_scores, corr_val] = test_native_adj(all_pt_adj{s}, pt_loc{s}, conn_edge, std_edge, [1:40]'); + + bivariate_native(feat).subj(a).data = native_adj_scores; + + catch anyerror + bivariate_native(feat).subj(a).data = NaN(length(pt_loc{s})); + end + + end + end + +end + + +%% do paired boxplot of channel-wise and patient-wise distributions for +% one MTL and one cortical connectivity +conn_pt_1 = -1.*log([rmmissing(squeeze(all_raw_atlas{1}(36,40,:))),rmmissing(squeeze(all_raw_atlas{2}(36,40,:))),rmmissing(squeeze(all_raw_atlas{3}(36,40,:))),rmmissing(squeeze(all_raw_atlas{4}(36,40,:))),rmmissing(squeeze(all_raw_atlas{5}(36,40,:)))]) +conn_ch_1 = -1.*log([rmmissing(squeeze(all_raw_atlas_edge{1}(36,40,:))),rmmissing(squeeze(all_raw_atlas_edge{2}(36,40,:))),rmmissing(squeeze(all_raw_atlas_edge{3}(36,40,:))),rmmissing(squeeze(all_raw_atlas_edge{4}(36,40,:))),rmmissing(squeeze(all_raw_atlas_edge{5}(36,40,:)))]) + +conn_pt_2 = -1.*log([rmmissing(squeeze(all_raw_atlas{1}(5,37,:))),rmmissing(squeeze(all_raw_atlas{2}(5,37,:))),rmmissing(squeeze(all_raw_atlas{3}(5,37,:))),rmmissing(squeeze(all_raw_atlas{4}(5,37,:))),rmmissing(squeeze(all_raw_atlas{5}(5,37,:)))]) +conn_ch_2 = -1.*log([rmmissing(squeeze(all_raw_atlas_edge{1}(5,37,:))),rmmissing(squeeze(all_raw_atlas_edge{2}(5,37,:))),rmmissing(squeeze(all_raw_atlas_edge{3}(5,37,:))),rmmissing(squeeze(all_raw_atlas_edge{4}(5,37,:))),rmmissing(squeeze(all_raw_atlas_edge{5}(5,37,:)))]) + +conn_pt_3 = -1.*log([rmmissing(squeeze(all_raw_atlas{1}(4,8,:))),rmmissing(squeeze(all_raw_atlas{2}(4,8,:))),rmmissing(squeeze(all_raw_atlas{3}(4,8,:))),rmmissing(squeeze(all_raw_atlas{4}(4,8,:))),rmmissing(squeeze(all_raw_atlas{5}(4,8,:)))]) +conn_ch_3 = -1.*log([rmmissing(squeeze(all_raw_atlas_edge{1}(4,8,:))),rmmissing(squeeze(all_raw_atlas_edge{2}(4,8,:))),rmmissing(squeeze(all_raw_atlas_edge{3}(4,8,:))),rmmissing(squeeze(all_raw_atlas_edge{4}(4,8,:))),rmmissing(squeeze(all_raw_atlas_edge{5}(4,8,:)))]) + +plot_data1 = padcat(conn_pt_1(:,1),conn_ch_1(:,1),conn_pt_1(:,2),conn_ch_1(:,2),conn_pt_1(:,3),conn_ch_1(:,3),conn_pt_1(:,4),conn_ch_1(:,4),conn_pt_1(:,5),conn_ch_1(:,5)); +plot_data2 = padcat(conn_pt_2(:,1),conn_ch_2(:,1),conn_pt_2(:,2),conn_ch_2(:,2),conn_pt_2(:,3),conn_ch_1(:,3),conn_pt_2(:,4),conn_ch_2(:,4),conn_pt_2(:,5),conn_ch_2(:,5)); +plot_data3 = padcat(conn_pt_3(:,1),conn_ch_3(:,1),conn_pt_3(:,2),conn_ch_3(:,2),conn_pt_3(:,3),conn_ch_1(:,3),conn_pt_3(:,4),conn_ch_3(:,4),conn_pt_3(:,5),conn_ch_3(:,5)); + +figure(1);clf; +subplot(1,3,1) +boxplot(plot_data1) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Negative log coherence') + + +subplot(1,3,2) +boxplot(plot_data2) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Negative log coherence') + + +subplot(1,3,3) +boxplot(plot_data3) + +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),colors(j,:),'FaceAlpha',.5); +end + +xticks([1.5 3.5 5.5 7.5 9.5]) +xticklabels({'delta','theta','alpha','beta','gamma'}) +ylabel('Negative log coherence') + + +num_samples2 = samples_edge; +for i = 1:40 + num_samples2(i,i) = 0; +end + +figure(2);clf +subplot(2,2,1) +imagesc(num_samples) +colorbar +subplot(2,2,2) +imagesc(num_samples2) +colorbar +subplot(2,2,3) +histogram(num_samples(:)) +subplot(2,2,4) +histogram(num_samples2(:)) +%% process bivariate edge features into nodal features + +for f = 1:5 + this_feat = []; + abs_feat = []; + bivariate_feat = []; + for s = 1:166 + % take 75th percentile across all edges of each node + %abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),75)']; + abs_feat = [abs_feat;prctile(abs(bivariate_native(f).subj(s).data),100)']; % test with 50th percentile + bivariate_feat = [bivariate_feat;prctile(bivariate_native(f).subj(s).data,75)']; + end + + % assign into feature matrix + abs_bivariate_feats(:,f) = abs_feat; + non_abs_bivariate_feats(:,f) = bivariate_feat; +end + +% for the single bivariate feature we want the maximum absolute Z score +% across all 10 bivariate features +[single_bivariate_feat] = median(abs_bivariate_feats')'; + +% assign previously calculated univariate z scores into bivariate scores +all_feat_zscores = [abs(univariate_zscores), abs_bivariate_feats]; %% original +%all_feat_zscores = [univariate_zscores, non_abs_bivariate_feats]; + +% incomplete channels are ones which are unlocalized in new atlas +incomplete_channels = find(isnan(sum(all_feat_zscores'))); +%single_bivariate_feat(incomplete_channels) = []; + +univariate_feats = all_feat_zscores(:,1:5); +univariate_feats(incomplete_channels,:) = []; +%[single_univariate_feat, which_feat_uni] = max(abs(univariate_feats),[],2); + + +%% Figure 6A: combine univariate and bivariate +%all_feat_zscores = [univariate_zscores, all_bivariate_feats]; +all_feat_zscores(isinf(all_feat_zscores)) = 0; +incomplete_channels = find(isnan(sum(all_feat_zscores'))); +complete_channels = find(~isnan(sum(all_feat_zscores'))); +all_feat_zscores2 = all_feat_zscores; +all_feat_zscores2(incomplete_channels,:) = []; % remove channels that are incomplete +EIZ = [zeros(1772,1);[spike_ind - spike_ind.*soz_ch]]; +resect_ch_all = [zeros(1772,1);resected_ch]; + +% 1: patient, 2: abnormal, 3: EIZ, 4: SOZ, 5: resect, 6:outcome +all_possible_labels = [all_pts,[EIZ+[zeros(1772,1);soz_ch]],EIZ,[zeros(1772,1);soz_ch],resect_ch_all,[ones(1772,1);HUP_outcome_all],[[zeros(1772,1);soz_ch].*resect_ch_all]]; +all_possible_labels(incomplete_channels,:) = []; + +% get only non EIZ normal versus resected SOZ (in good outcome patients) +% here we eliminate EIZ channels, as well as SOZ +eliminate_channels = [find(all_possible_labels(:,3)==1);find((all_possible_labels(:,4)+all_possible_labels(:,5))==1)];%find(non_normal_non_soz==1); +%retain_channels = find(all_possible_labels(:,2)==0); +all_feat_zscores3 = all_feat_zscores2; +all_feat_zscores3(eliminate_channels,:) = []; +new_labels = all_possible_labels; +new_labels(eliminate_channels,:) = []; + +all_feats = all_feat_zscores3; +all_labels = new_labels(:,7); % was column 5 before +all_pred_1 = zeros(size(all_labels)); +all_pred_2 = zeros(size(all_labels)); +all_pred_3 = zeros(size(all_labels)); + +num_parts = 10; + +[pt_part] = make_xval_partition(166, num_parts); + +% cross-validate based on patients +for p = 1:num_parts + indices = find(pt_part==p) + for i = 1:length(indices) + part(find(indices(i)==new_labels(:,1))) = p; + + part2(find(indices(i)==all_possible_labels(:,1))) = p; + end +end + +% Define classes, 0 = non seizure, 1 = seizure +ClassNames = [0,1]; + +% Define costs -> penalty here is 5 +cost.ClassNames = ClassNames; +cost.ClassificationCosts = [0 1; 1 0]; + +for p = 1:num_parts + p + X_train = all_feats(find(part~=p),:); + Y_train = all_labels(find(part~=p)); + X_test = all_feats(find(part==p),:); + Y_test = all_labels(find(part==p)); + + X_analysis = all_feat_zscores2(find(part2==p),:); + + Mdl1 = TreeBagger(100,X_train(:,1:5),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + Mdl2 = TreeBagger(100,X_train(:,6:10),Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + Mdl3 = TreeBagger(100,X_train,Y_train,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + [Yfit,pihat1] = predict(Mdl1,X_test(:,1:5)); + [Yfit,pihat2] = predict(Mdl2,X_test(:,6:10)); + [Yfit,pihat3] = predict(Mdl3,X_test); + [Yfit, pihat_analysis] = predict(Mdl3,X_analysis); + [Yfit, pihat_analysis2] = predict(Mdl1,X_analysis(:,1:5)); + [Yfit, pihat_analysis3] = predict(Mdl2,X_analysis(:,6:10)); + + [~, Y_pred] = max(pihat1,[],2); + all_pred_1(find(part==p)) = pihat1(:,2); + all_pred_2(find(part==p)) = pihat2(:,2); + all_pred_3(find(part==p)) = pihat3(:,2); + all_pred_analysis(find(part2==p),1) = pihat_analysis(:,2); + all_pred_analysis2(find(part2==p),1) = pihat_analysis2(:,2); + all_pred_analysis3(find(part2==p),1) = pihat_analysis3(:,2); + + xfold_acc(p) = sum(Y_pred==Y_test)./length(Y_test); + + Y_pred = round(pihat1(:,2)); + + +end + + [X1,Y1,T,AUC1] = perfcurve(all_labels,all_pred_1,1); + [X2,Y2,T,AUC2] = perfcurve(all_labels,all_pred_2,1); + [X3,Y3,T,AUC3] = perfcurve(all_labels,all_pred_3,1); + + AUC1 + AUC2 + AUC3 + + Mdl3 = TreeBagger(100,all_feats,all_labels,'Method',... + 'classification','PredictorSelection','curvature','OOBPredictorImportance','on','cost',cost); + + +%% Figure 6B: random forest feature importance +imp = Mdl3.OOBPermutedPredictorDeltaError; + +figure(1);clf; +bar(imp); +title('Curvature Test'); +ylabel('Predictor importance estimates'); +xticks([1:10]) + +figure(2);clf; +hold on +plot(X1,Y1,'LineWidth',2,'Color',color4)%4 +plot(X2,Y2,'LineWidth',2,'Color',color7)%7 +plot(X3,Y3,'LineWidth',2,'Color',color1)%10 +legend(sprintf('Univariate - AUC: %.2f',AUC1),sprintf('Bivariate - AUC: %.2f',AUC2),sprintf('All - AUC: %.2f',AUC3),'Location','SouthEast')%86/85/91 +xlabel('False positive rate') +ylabel('True positive rate') +title('All abnormal vs. normal channel classification') + +%% Figure 5A: mean Z-scores for MTL (regions 17-20) in SOZ vs IZ vs uninvolved +new_roi_2 = new_roi; +new_roi_2(incomplete_channels) = []; + +single_bivariate_feat_2 = single_bivariate_feat; +single_bivariate_feat_2(incomplete_channels) = []; + +MTL_uninvolved = []; +MTL_EIZ = []; +MTL_SOZ = []; + +for pt = 1:166 + uninv_pt = find([all_possible_labels(:,1)==pt].*[all_possible_labels(:,2)==0].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + EIZ_pt = find([all_possible_labels(:,1)==pt].*[all_possible_labels(:,3)==1].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + SOZ_pt = find([all_possible_labels(:,1)==pt].*[all_possible_labels(:,4)==1].*[[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==18]+[new_roi_2==20]]); + + MTL_uninvolved = [MTL_uninvolved; mean(single_bivariate_feat_2(uninv_pt))]; + MTL_EIZ = [MTL_EIZ; mean(single_bivariate_feat_2(EIZ_pt))]; + MTL_SOZ = [MTL_SOZ; mean(single_bivariate_feat_2(SOZ_pt))]; +end + +mean_uninv = MTL_uninvolved; +mean_EIZ = MTL_EIZ; +mean_SOZ = MTL_SOZ; + + +% MTL_uninvolved = find([all_possible_labels(:,2)==0].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); +% MTL_EIZ = find([all_possible_labels(:,3)==1].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); +% MTL_SOZ = find([all_possible_labels(:,4)==1].*[[new_roi_2==18]+[new_roi_2==17]+[new_roi_2==19]+[new_roi_2==20]]); +% +% mean_uninv = single_bivariate_feat_2(MTL_uninvolved); +% mean_EIZ = single_bivariate_feat_2(MTL_EIZ); +% mean_SOZ = single_bivariate_feat_2(MTL_SOZ); + +plot_matrix = padcat(mean_uninv,mean_EIZ,mean_SOZ); +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','EIZ','SOZ'}) +ylim([0 5]) + +p1 = ranksum(mean_uninv,mean_EIZ,'tail','left') +cd1 = computeCohen_d(mean_uninv,mean_EIZ) +p2 = ranksum(mean_uninv,mean_SOZ,'tail','left') +cd2 = computeCohen_d(mean_uninv,mean_SOZ) +p3 = ranksum(mean_EIZ,mean_SOZ,'tail','left') +cd3 = computeCohen_d(mean_SOZ,mean_EIZ) + +%% need to map back onto all channels and render +all_pred = zeros(5203,1); +all_pred(incomplete_channels) = 0; +all_pred(complete_channels) = all_pred_analysis; + +all_pred2 = zeros(5203,1); +all_pred2(incomplete_channels) = 0; +all_pred2(complete_channels) = all_pred_analysis2; + +all_pred3 = zeros(5203,1); +all_pred3(incomplete_channels) = 0; +all_pred3(complete_channels) = all_pred_analysis3; + +%% render abnormalities +soz_ch_all = [zeros(1772,1);soz_ch]; +this_pt = 116;%116/122/154 / 154 %114 res - 45/50 EIZ - 52/64, non-inv - 20,30 % 120 - res - 2., EIZ - 23./64, non-inv - 33/50. +this_pt_abn = all_pred(find(all_pts==this_pt)); +this_pt_soz = soz_ch_all(find(all_pts==this_pt)); +this_pt_eiz = EIZ(find(all_pts==this_pt)); +this_pt_res = resect_ch_all(find(all_pts==this_pt)); +this_pt_coords = all_coords(find(all_pts==this_pt),:); +this_pt_eeg = all_wake_data(:,find(all_pts==this_pt)); +%this_pt_delta_coh_abn(isnan(this_pt_delta_coh_abn)) = 0; + +this_pt_multiclass = 2*this_pt_soz+this_pt_eiz; + +this_pt_zscores = all_feat_zscores(find(all_pts==this_pt),:); + +this_pt_test = zeros(size(this_pt_res)); +this_pt_test([28,9,49])=1; + +% figure(1);clf; +% hold on +% for i = 1:10 +% plot([1:2000],[this_pt_eeg(:,4*i)+i*175],'k-') +% end +% hold off + +% figure(1);clf; +% imagesc(all_feat_zscores(find(all_pts==this_pt),:)) +% colorbar +% caxis([0,5]) + +figure(2);clf; +subplot(1,3,1) +imagesc(this_pt_zscores(28,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +subplot(1,3,2) +imagesc(this_pt_zscores(9,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +subplot(1,3,3) +imagesc(this_pt_zscores(49,:)') +caxis([0,3]) +colormap('jet') +colorbar +colorbar +% +figure(3); +subplot(3,1,1) +plot(this_pt_eeg(6001:12000,4)) +subplot(3,1,2) +plot(this_pt_eeg(6001:12000,23)) +subplot(3,1,3) +plot(this_pt_eeg(6001:12000,27)) + + +figure(4);clf; +final_elec_matrix = [this_pt_coords,this_pt_abn,ones(size(this_pt_res))]; +%final_elec_matrix = final_elec_matrix([2,23,50],:) +dlmwrite('render_elecs.node',final_elec_matrix,'delimiter',' ','precision',5) +BrainNet_MapCfg('BrainMesh_ICBM152_smoothed.nv','render_elecs.node','abnormality_render.mat') +%delete render_elecs.node + +%% patient-specific correlation between abnormality features +for s = 107:166 + pt_zscores = all_feat_zscores2(find([all_possible_labels(:,1)==s]),:); + + pt_zscores_univar = pt_zscores(:,1:5); + pt_zscores_bivar = pt_zscores(:,6:10); + + corrval2(s-106) = corr(pt_zscores_univar(:),pt_zscores_bivar(:),'Type','Spearman'); + + for f = 1:5 + corrval((s-106),f) = corr(pt_zscores(:,f),pt_zscores(:,(f+5)),'Type','Spearman'); + end + + corrval_outcome((s-106),1) = late_outcome(s-106); + corrval_lesion((s-106),1) = lesion_field(s-106); +end + +p1 = ranksum(corrval(corrval_outcome==1,1),corrval(corrval_outcome>1,1),'tail','right') +p2 = ranksum(corrval(corrval_outcome==1,2),corrval(corrval_outcome>1,2),'tail','right') +p3 = ranksum(corrval(corrval_outcome==1,3),corrval(corrval_outcome>1,3),'tail','right') +p4 = ranksum(corrval(corrval_outcome==1,4),corrval(corrval_outcome>1,4),'tail','right') +p5 = ranksum(corrval(corrval_outcome==1,5),corrval(corrval_outcome>1,5),'tail','right') + +p6 = ranksum(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1),'tail','right') +d6 = computeCohen_d(corrval2(corrval_outcome==1),corrval2(corrval_outcome>1)) + +plot_matrix = padcat(corrval2(corrval_outcome==1)',corrval2(corrval_outcome>1)'); +figure(1);clf; +UnivarScatter(plot_matrix); +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylim([-0.3 0.5]) +ylabel('Correlation of beta power & coherence |Z|') + +%% Figure 5C: compute AUPRC/DRS at a per-patient level +a = 0; +for s = [107:127,128:132,134:139,141:155,157:161,163:166] + s + a = a+1; + this_pt_label = all_possible_labels(find(all_possible_labels(:,1)==s),4); + this_pt_pred = all_pred_analysis(find(all_possible_labels(:,1)==s)); + this_pt_pred2 = all_pred_analysis2(find(all_possible_labels(:,1)==s)); + this_pt_pred3 = all_pred_analysis3(find(all_possible_labels(:,1)==s)); + [X1,Y1,T,AUC1] = perfcurve(this_pt_label,this_pt_pred,1,'XCrit','TPR','YCrit','PPV'); + [X2,Y2,T,AUC2] = perfcurve(this_pt_label,this_pt_pred2,1,'XCrit','TPR','YCrit','PPV'); + [X3,Y3,T,AUC3] = perfcurve(this_pt_label,this_pt_pred3,1,'XCrit','TPR','YCrit','PPV'); + + auc_pt_outcome(a,1) = late_outcome(s-106); + pt_auc(a,1:3) = [AUC1,AUC2,AUC3]; +end + +good_inds = find(auc_pt_outcome==1); +poor_inds = find(auc_pt_outcome>1); + +plot_matrix = padcat(pt_auc(good_inds,1),pt_auc(poor_inds,1)); +figure(1);clf; +UnivarScatter(plot_matrix) + +ranksum(plot_matrix(:,1),plot_matrix(:,2)) +computeCohen_d(plot_matrix(:,1),plot_matrix(:,2)) +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylabel('Area under precision-recall curve') +title('Comparison between AUPRC, rank-sum p = 0.033') +%% Figure 6C: patient-level abnormality score all regions + +pt_uninvolved = []; +pt_eiz = []; +pt_soz = []; + +% loop through HUP patients +for s = 107:166 + + pt_uninvolved_inds = find([all_possible_labels(:,2)==0].*[all_possible_labels(:,1)==s]); + pt_eiz_inds = find([all_possible_labels(:,3)==1].*[all_possible_labels(:,1)==s]); + pt_soz_inds = find([all_possible_labels(:,4)==1].*[all_possible_labels(:,1)==s]); + + pt_uninvolved = [pt_uninvolved; nanmean(all_pred_analysis(pt_uninvolved_inds))]; + pt_eiz = [pt_eiz; nanmean(all_pred_analysis(pt_eiz_inds))]; + pt_soz = [pt_soz; nanmean(all_pred_analysis(pt_soz_inds))]; +end + +plot_matrix = padcat(pt_uninvolved,pt_eiz,pt_soz); + +figure(1);clf; +UnivarScatter(plot_matrix) +xticks([1:3]) +xticklabels({'uninvolved','irritative zone','seizure onset zone'}) +ylabel('predicted epileptogenicity score') +ylim([0 0.8]) + +p1 = ranksum(pt_uninvolved,pt_eiz) +d1 = computeCohen_d(pt_uninvolved,pt_eiz) +p2 = ranksum(pt_eiz,pt_soz) +d2 = computeCohen_d(pt_eiz,pt_soz) +p3 = ranksum(pt_uninvolved,pt_soz) +d3 = computeCohen_d(pt_uninvolved,pt_soz) + +%% Figure 5B: correlation length of diagnosis + +single_univariate_feat_2 = median(univariate_feats')'; + +pt_in_rz = []; +pt_out_rz = []; + +pt_in_rz2 = []; +pt_out_rz2 = []; + +for s = 107:166 + + w_inds = find([all_possible_labels(:,1)==s])%;.*[all_possible_labels(:,4)==1]); + o_inds = find([all_possible_labels(:,1)==s])%;.*[all_possible_labels(:,4)==0]); + + length_diagnosis = age_surgery - age_onset; + + % bivariate features + pt_in_rz = [pt_in_rz; nanmedian(single_bivariate_feat_2(w_inds))]; + pt_out_rz = [pt_out_rz; nanmedian(single_bivariate_feat_2(o_inds))]; + + % univariate features + pt_in_rz2 = [pt_in_rz2; nanmedian(single_univariate_feat_2(w_inds))]; + pt_out_rz2 = [pt_out_rz2; nanmedian(single_univariate_feat_2(o_inds))]; + +end + +figure(1);clf; +hold on +plot(length_diagnosis, pt_out_rz,'b.','Color',[106 178 180]/255,'MarkerSize',12) +p = polyfit(length_diagnosis,pt_out_rz,1); +f = polyval(p,length_diagnosis); +plot(length_diagnosis,f,'b-','LineWidth',2) +[a2,b2] = corr(length_diagnosis, pt_out_rz,'type','Pearson','tail','right') +xlim([0 45]) +ylim([0.6 2.6]) + +[a3,b3] = corr(length_diagnosis, pt_out_rz2,'type','Pearson','tail','right') + +%% calculate spectral density differences between dataset 2 and dataset 3 + +load extra_data_1 +extra1_wake = wake_clip; + +load extra_data_2 +extra2_wake = wake_clip; + +[pxx,f] = pwelch(extra1_wake,400,200,[0.5,2:80],200); % original was 200/100 +pxx_norm1 = pxx./sum(pxx); + +[pxx,f] = pwelch(extra2_wake,400,200,[0.5,2:80],200); % original was 200/100 +pxx_norm2 = pxx./sum(pxx); + +for freq = 1:5 + if freq==1 + feat_vals1 = median(pxx_norm1(1:4,:))'; + feat_vals2 = median(pxx_norm2(1:4,:))'; + elseif freq==2 + feat_vals1 = median(pxx_norm1(4:8,:))'; + feat_vals2 = median(pxx_norm2(4:8,:))'; + elseif freq==3 + feat_vals1 = median(pxx_norm1(8:13,:))'; + feat_vals2 = median(pxx_norm2(8:13,:))'; + elseif freq==4 + feat_vals1 = median(pxx_norm1(13:30,:))'; + feat_vals2 = median(pxx_norm2(13:30,:))'; + elseif freq==5 + feat_vals1 = median(pxx_norm1(30:80,:))'; + feat_vals2 = median(pxx_norm2(30:80,:))'; + end + + freq_feats1(freq,1:3431) = feat_vals1'; + freq_feats2(freq,1:3431) = feat_vals2'; + +end + +for freq = 1:5 + freq_corr(freq) = corr(freq_feats1(freq,:)',freq_feats2(freq,:)') +end + +%% check connectivity from clips 2 and 3 + +%% correlate power spectral density in three patients +% pts = [34, 38, 44, 49, 60]; +load AED_recordings.mat +figure(1);clf +a = 0; + +for pt = [1,2,4] + + freq_data_pt = []; + + for c = 1:4 + data_clip = AED_recordings(pt).clip(c).data; + [pxx,f] = pwelch(data_clip,400,200,[0.5,2:80],200); % original was 200/100 + pxx_norm_AED = pxx./sum(pxx); + + % calculate power data + for freq = 1:5 + if freq==1 + feat_vals_AED = median(pxx_norm_AED(1:4,:))'; + elseif freq==2 + feat_vals_AED = median(pxx_norm_AED(4:8,:))'; + elseif freq==3 + feat_vals_AED = median(pxx_norm_AED(8:13,:))'; + elseif freq==4 + feat_vals_AED = median(pxx_norm_AED(13:30,:))'; + elseif freq==5 + feat_vals_AED = median(pxx_norm_AED(30:80,:))'; + end + + freq_data_pt(c,:,freq) = feat_vals_AED'; + + end + + end + + % analyze across clips + for freq = 1:5 + a = a+1; + subplot(3,5,a) + plot(freq_data_pt(2,:,freq),freq_data_pt(4,:,freq),'k.') + [RHO,PVAL] = corr(freq_data_pt(2,:,freq)',freq_data_pt(4,:,freq)','type','Pearson'); + pvals(freq,pt) = PVAL; + rhovals(freq,pt) = RHO; + title(sprintf('p = %.1d, R = %.2f',PVAL,RHO)) + end + +end + +%% do a connectivity calculation here +conn_band = {'delta','theta','alpha','beta','gamma'}; + +for pt = 5 + + % loop through subjects + for freq = 1:5 + + for c = 1:4 + + % extract subject data from all the data (columns are channels) + pt_data = AED_recordings(pt).clip(c).data; + + % call functional connectivity code (200 is Hz, 1 is window size in seconds) + [pt_adj] = compute_FC(pt_data, 200, 1, 'coh',conn_band{freq}); + + % assign median adjacency matrix into cell structure + AED_conn(pt).freq(freq).data(:,:,c) = pt_adj; + end + + end +end + +m = 1 + +%% correlate connectivity in three patients +figure(1);clf; +a = 0; +for pt = [1,2,4] + for freq = 1:5 + matrix1 = -1.*log(AED_conn(pt).freq(freq).data(:,:,2)); + matrix2 = -1.*log(AED_conn(pt).freq(freq).data(:,:,4)); + + a = a+1; + subplot(3,5,a) + plot(matrix1(:),matrix2(:),'k.') + [RHO,PVAL] = corr(matrix1(:),matrix2(:),'type','Pearson'); + pvals_conn(freq,pt) = PVAL; + rhovals_conn(freq,pt) = RHO; + title(sprintf('p = %.1d, R = %.2f',PVAL,RHO)) + end +end + diff --git a/multivar_model.mat b/multivar_model.mat deleted file mode 100644 index 65191e3..0000000 Binary files a/multivar_model.mat and /dev/null differ diff --git a/new_atlas_render.jpg b/new_atlas_render.jpg index 3aa4610..7176f9f 100644 Binary files a/new_atlas_render.jpg and b/new_atlas_render.jpg differ diff --git a/random_colors.mat b/random_colors.mat deleted file mode 100644 index acf635d..0000000 Binary files a/random_colors.mat and /dev/null differ diff --git a/render_elecs.node b/render_elecs.node new file mode 100644 index 0000000..28f44b1 --- /dev/null +++ b/render_elecs.node @@ -0,0 +1,69 @@ +-49.817 3.1753 -45.93 0.01 1 +-56.177 -2.1784 -38.013 0.01 1 +-60.86 -8.5388 -30.02 0.01 1 +-62.594 -28.457 -7.3278 4.0323e-05 1 +-62.526 -26.572 2.0837 0.01 1 +-61.47 -22.42 12.458 0 1 +-57.608 -18.487 22.76 0.060161 1 +-52.598 -18.168 31.387 0.17 1 +-48.176 -19.7 44.123 0.080384 1 +-63.549 -39.248 -12.401 0 1 +-61.593 -37.573 -2.9488 0 1 +-60.868 -35.303 5.3732 0 1 +-59.198 -32.034 15.662 0.1042 1 +-56.34 -29.647 25.937 0.093333 1 +-50.874 -28.735 34.426 0.19744 1 +-45.397 -28.373 45.584 0.29 1 +-59.272 -50.089 -11.015 0.0575 1 +-57.761 -46.798 -1.7 0.030667 1 +-57.663 -43.034 9.7729 0.02 1 +-57.256 -38.536 21.185 0.25333 1 +-54.839 -36.014 30.687 0.15693 1 +-49.706 -36.152 40.533 0.31333 1 +-44.149 -37.093 51.194 0.66 1 +-51.744 -62.193 -7.9696 0 1 +-54.266 -58.149 1.3359 0.030159 1 +-53.041 -54.084 11.979 0.010094 1 +-53.094 -50.265 23.504 0.070159 1 +-53.27 -45.923 34.252 0.05837 1 +-48.26 -44.835 44.618 0.17443 1 +-47.597 -71.6 -3.6689 0.070049 1 +-49.664 -66.517 7.9975 0.090159 1 +-49.109 -61.485 18.476 0.10357 1 +-47.269 -58.432 28.142 0.065574 1 +-45.473 -54.956 39.81 0.14258 1 +-42.799 -52.188 50.126 0.20421 1 +-38.008 -49.324 56.559 0.10009 1 +-43.873 -77.327 -0.54637 0.15321 1 +-43.026 -73.354 11.569 0.10321 1 +-41.729 -69.5 21.04 0.10037 1 +-40.472 -64.888 29.712 0.1725 1 +-38.71 -60.206 41.61 0.04 1 +-38.085 -85.076 3.5616 0.07062 1 +-36.638 -81.708 15.105 0.069086 1 +-34.412 -77.226 24.224 0.073046 1 +-33.946 -71.368 31.507 0.091 1 +-36.096 -66.24 41.528 0.26405 1 +-34.413 -60.604 52.35 0.20667 1 +-26.721 -54.495 58.901 0.22676 1 +-30.418 -90.628 5.5725 0.216 1 +-29.241 -87.621 14.649 0.11037 1 +-27.159 -84.327 23.049 0.090784 1 +-24.816 -80.093 32.731 0.051731 1 +-25.034 -74.945 43.296 0.072833 1 +-23.41 -69.104 52.197 0.040161 1 +-18.05 -60.588 61.372 0.057538 1 +-43.147 -32.737 -29.755 0 1 +-50.223 -40.183 -28.166 0 1 +-56.778 -48.305 -23.919 0.010667 1 +-1.2227 -64.512 34.691 0.36604 1 +-2.6651 -67.771 44.61 0.24719 1 +-7.7806 -69.377 54.31 0.57 1 +-13.982 -70.724 -7.0693 0.023495 1 +-23.014 -73.572 -12.763 0.24 1 +-34.794 -77.193 -14.838 0.23792 1 +-14.199 -87 -11.846 0.071161 1 +-8.68 -89.218 -7.6094 0.030161 1 +-6.8114 -92.566 0.31663 0.12421 1 +-7.0915 -95.105 9.1392 0.1525 1 +-10.589 -94.642 17.507 0.17093 1 diff --git a/renderings/.DS_Store b/renderings/.DS_Store new file mode 100644 index 0000000..9c44efe Binary files /dev/null and b/renderings/.DS_Store differ diff --git a/renderings/HUP065_pred.jpg b/renderings/HUP065_pred.jpg new file mode 100644 index 0000000..5dfa2c6 Binary files /dev/null and b/renderings/HUP065_pred.jpg differ diff --git a/renderings/HUP065_soz_eiz.jpg b/renderings/HUP065_soz_eiz.jpg new file mode 100644 index 0000000..bc423ba Binary files /dev/null and b/renderings/HUP065_soz_eiz.jpg differ diff --git a/renderings/HUP068_pred.jpg b/renderings/HUP068_pred.jpg new file mode 100644 index 0000000..310222f Binary files /dev/null and b/renderings/HUP068_pred.jpg differ diff --git a/renderings/HUP070_pred.jpg b/renderings/HUP070_pred.jpg new file mode 100644 index 0000000..e2054fb Binary files /dev/null and b/renderings/HUP070_pred.jpg differ diff --git a/renderings/HUP074_pred.jpg b/renderings/HUP074_pred.jpg new file mode 100644 index 0000000..831d37c Binary files /dev/null and b/renderings/HUP074_pred.jpg differ diff --git a/renderings/HUP075_pred.jpg b/renderings/HUP075_pred.jpg new file mode 100644 index 0000000..f92aa11 Binary files /dev/null and b/renderings/HUP075_pred.jpg differ diff --git a/renderings/HUP075_soz_eiz.jpg b/renderings/HUP075_soz_eiz.jpg new file mode 100644 index 0000000..1c0551b Binary files /dev/null and b/renderings/HUP075_soz_eiz.jpg differ diff --git a/renderings/HUP083_pred.jpg b/renderings/HUP083_pred.jpg new file mode 100644 index 0000000..1269006 Binary files /dev/null and b/renderings/HUP083_pred.jpg differ diff --git a/renderings/HUP083_soz_eiz.jpg b/renderings/HUP083_soz_eiz.jpg new file mode 100644 index 0000000..4190375 Binary files /dev/null and b/renderings/HUP083_soz_eiz.jpg differ diff --git a/renderings/HUP097_pred.jpg b/renderings/HUP097_pred.jpg new file mode 100644 index 0000000..b95fb3e Binary files /dev/null and b/renderings/HUP097_pred.jpg differ diff --git a/renderings/HUP097_soz_eiz.jpg b/renderings/HUP097_soz_eiz.jpg new file mode 100644 index 0000000..a2aa114 Binary files /dev/null and b/renderings/HUP097_soz_eiz.jpg differ diff --git a/renderings/HUP114_pred.jpg b/renderings/HUP114_pred.jpg new file mode 100644 index 0000000..a31f6c2 Binary files /dev/null and b/renderings/HUP114_pred.jpg differ diff --git a/renderings/HUP123_pred.jpg b/renderings/HUP123_pred.jpg new file mode 100644 index 0000000..7deab2a Binary files /dev/null and b/renderings/HUP123_pred.jpg differ diff --git a/renderings/HUP138_pred.jpg b/renderings/HUP138_pred.jpg new file mode 100644 index 0000000..f64ee0c Binary files /dev/null and b/renderings/HUP138_pred.jpg differ diff --git a/renderings/HUP165_pred.jpg b/renderings/HUP165_pred.jpg new file mode 100644 index 0000000..3eaa722 Binary files /dev/null and b/renderings/HUP165_pred.jpg differ diff --git a/renderings/HUP165_soz_eiz.jpg b/renderings/HUP165_soz_eiz.jpg new file mode 100644 index 0000000..e869066 Binary files /dev/null and b/renderings/HUP165_soz_eiz.jpg differ diff --git a/renderings/HUP179_pred.jpg b/renderings/HUP179_pred.jpg new file mode 100644 index 0000000..62e2dcd Binary files /dev/null and b/renderings/HUP179_pred.jpg differ diff --git a/supplemental_analysis.m b/supplemental_analysis.m new file mode 100644 index 0000000..4765d68 --- /dev/null +++ b/supplemental_analysis.m @@ -0,0 +1,178 @@ +% supplmental figures + + +%% do renderings of original and new atlas +random_colors = rand([40,3]); + +color_matrix = ones(117,3); + +for i = 1:40 + this_region = custom_atlas{i,2:4}; + for r = 1:length(rmmissing(this_region)) + color_matrix(this_region(r),:) = random_colors(i,:); + end +end + +% writes a file that can be used in BrainNet viewer to visualize +dlmwrite('atlas_render_color.txt',color_matrix(1:90,:),'delimiter',' ','precision',5) + +figure(1);clf; +imagesc(reshape(random_colors,40,1,3)) +set(gca,'TickLabelInterpreter', 'none'); +yticks([1:40]) +yticklabels(custom_atlas{:,1}) + +%% correlate connectivity between MNI and HUP +for feat = 1:5 + a = conn_band{feat}; + aa = 'coherence'; + + conn_edge_mni = this_feat_conn(feat).mni; + conn_edge_hup = this_feat_conn(feat).hup; + +conn_edge_all = conn_edge_mni(:)+conn_edge_hup(:); +conn_edge_mni(isnan(conn_edge_all)) = []; +conn_edge_hup(isnan(conn_edge_all)) = []; + +std_edge_all = std_edge_mni(:)+std_edge_hup(:); +std_edge_mni(isnan(std_edge_all)) = []; +std_edge_hup(isnan(std_edge_all)) = []; + +figure(1); +subplot(1,5,feat) +hold on +plot((conn_edge_mni),(conn_edge_hup),'k.','Color',[128 128 128]/255) +[r,p1] = corr((conn_edge_mni)',(conn_edge_hup)') +rvals(feat) = r; +xlabel('MNI connectivity') +ylabel('HUP connectivity') +p = polyfit((conn_edge_mni),(conn_edge_hup),1); +f = polyval(p,(conn_edge_mni)); +plot((conn_edge_mni),f,'k-','LineWidth',2) +title(sprintf('%s %s, r = %f, p = %f',a,aa,r,p1)) +hold off +end + +% figure(2);clf +% bar(rvals) +%% all-node abnormality score all regions +for ch = 1:4717 + if all_possible_labels(ch,2)==0 + grouping{ch,1} = 'uninvolved'; + class_ind(ch,1) = 1; + elseif all_possible_labels(ch,3)==1 + grouping{ch,1} = 'EIZ'; + class_ind(ch,1) = 2; + elseif all_possible_labels(ch,4)==1 + grouping{ch,1} = 'SOZ'; + class_ind(ch,1) = 3; + end +end + +new_roi_2 = new_roi; +new_roi_2(incomplete_channels) = []; + +single_bivariate_feat_2 = single_bivariate_feat; +single_bivariate_feat_2(incomplete_channels) = []; + +for f = 1:5 + for r = [1:16,17:20] + + new_roi_2(new_roi_2==(2*r-1)) = 2*r; + + pred_uninvolved = all_feat_zscores2(find((class_ind==1).*(new_roi_2==2*r)),(f+5)); + pred_eiz = all_feat_zscores2(find((class_ind==2).*(new_roi_2==2*r)),(f+5)); + pred_soz = all_feat_zscores2(find((class_ind==3).*(new_roi_2==2*r)),(f+5)); + + median(pred_uninvolved) + median(pred_eiz) + median(pred_soz) + + plot_matrix = padcat(pred_uninvolved, pred_eiz,pred_soz); + + p1_region(f,r) = ranksum(pred_uninvolved,pred_eiz) + p2_region(f,r) = ranksum(pred_uninvolved,pred_soz) + p3_region(f,r) = ranksum(pred_eiz,pred_soz) + + cohen_d2(f,r) = computeCohen_d(pred_uninvolved,pred_soz) + cohen_d3(f,r) = computeCohen_d(pred_eiz,pred_soz) + + end +end + +figure(1);clf; +subplot(3,1,1) +imagesc([p3_region]) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +colorbar +subplot(3,1,2) +imagesc([p3_region]<0.05) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +subplot(3,1,3) +imagesc([cohen_d3]) +xticks([1:20]) +xticklabels(split_roi_name) +xtickangle(-45) +caxis([0,1]) +colorbar + +%% patient-level abnormality score res in/out good/poor outcome +good_rz = []; +good_out = []; + +poor_rz = []; +poor_out = []; + +for s = 107:166 + this_pt_outcome = late_outcome(s-106); + this_pt_ind_rz = find([all_possible_labels(:,1)==s].*[all_possible_labels(:,5)==1]); + this_pt_ind_out = find([all_possible_labels(:,1)==s].*[all_possible_labels(:,5)==0]); + + if this_pt_outcome==1 + good_rz = [good_rz; nanmedian(all_pred_analysis(this_pt_ind_rz))]; + good_out = [good_out; nanmedian(all_pred_analysis(this_pt_ind_out))]; + else + poor_rz = [poor_rz; nanmedian(all_pred_analysis(this_pt_ind_rz))]; + poor_out = [poor_out; nanmedian(all_pred_analysis(this_pt_ind_out))]; + end +end + +plot_matrix = padcat(good_rz,good_out,poor_rz,poor_out); +figure(2);clf; +UnivarScatter(plot_matrix) + +%% compute AUC at a per-patient level +a = 0; +for s = [107:127,128:132,134:139,141:155,157:161,163:166] + s + a = a+1; + this_pt_label = all_possible_labels(find(all_possible_labels(:,1)==s),4); + this_pt_pred = all_pred_analysis(find(all_possible_labels(:,1)==s)); + this_pt_pred2 = all_pred_analysis2(find(all_possible_labels(:,1)==s)); + this_pt_pred3 = all_pred_analysis3(find(all_possible_labels(:,1)==s)); + [X1,Y1,T,AUC1] = perfcurve(this_pt_label,this_pt_pred,1); + [X2,Y2,T,AUC2] = perfcurve(this_pt_label,this_pt_pred2,1); + [X3,Y3,T,AUC3] = perfcurve(this_pt_label,this_pt_pred3,1); + + auc_pt_outcome(a,1) = late_outcome(s-106); + pt_auc(a,1:3) = [AUC1,AUC2,AUC3]; +end + +good_inds = find(auc_pt_outcome==1); +poor_inds = find(auc_pt_outcome>1); + +plot_matrix = padcat(pt_auc(good_inds,1),pt_auc(poor_inds,1)); +figure(1);clf; +UnivarScatter(plot_matrix) + +ranksum(plot_matrix(:,1),plot_matrix(:,2)) +computeCohen_d(plot_matrix(:,1),plot_matrix(:,2)) +xticks([1,2]) +xticklabels({'Engel 1','Engel 2+'}) +ylabel('Area under precision-recall curve') +title('Comparison between AUPRC, rank-sum p = 0.033') + diff --git a/support_files/compute_FC.m b/support_files/compute_FC.m index 21ca1bb..b005594 100644 --- a/support_files/compute_FC.m +++ b/support_files/compute_FC.m @@ -34,33 +34,33 @@ pt_adj = zeros(size(iEEG_signals,2),size(iEEG_signals,2),num_windows); window_size = window_len.*sample_rate; - % filter in order to extract desired frequency band - if strcmp(band,'delta') - % bandpass filter 0.5-4 Hz - [b1,a1] = butter(4,[0.5 4]/(200/2)); - iEEG_signals = filter(b1, a1, iEEG_signals); - - elseif strcmp(band,'theta') - % bandpass filter 4-8 Hz - [b2,a2] = butter(4,[4 8]/(200/2)); - iEEG_signals = filter(b2, a2, iEEG_signals); - - elseif strcmp(band,'alpha') - % bandpass filter 8-13 Hz - [b3,a3] = butter(4,[8 13]/(200/2)); - iEEG_signals = filter(b3, a3, iEEG_signals); - - elseif strcmp(band,'beta') - % bandpass filter 13-30 Hz - [b4,a4] = butter(4,[13 30]/(200/2)); - iEEG_signals = filter(b4, a4, iEEG_signals); - - elseif strcmp(band,'gamma') - % bandpass filter 30-80 Hz - [b5,a5] = butter(4,[30 80]/(200/2)); - iEEG_signals = filter(b5, a5, iEEG_signals); - - end +% % filter in order to extract desired frequency band +% if strcmp(band,'delta') +% % bandpass filter 0.5-4 Hz +% [b1,a1] = butter(4,[0.5 4]/(200/2)); +% iEEG_signals = filter(b1, a1, iEEG_signals); +% +% elseif strcmp(band,'theta') +% % bandpass filter 4-8 Hz +% [b2,a2] = butter(4,[4 8]/(200/2)); +% iEEG_signals = filter(b2, a2, iEEG_signals); +% +% elseif strcmp(band,'alpha') +% % bandpass filter 8-13 Hz +% [b3,a3] = butter(4,[8 13]/(200/2)); +% iEEG_signals = filter(b3, a3, iEEG_signals); +% +% elseif strcmp(band,'beta') +% % bandpass filter 13-30 Hz +% [b4,a4] = butter(4,[13 30]/(200/2)); +% iEEG_signals = filter(b4, a4, iEEG_signals); +% +% elseif strcmp(band,'gamma') +% % bandpass filter 30-80 Hz +% [b5,a5] = butter(4,[30 80]/(200/2)); +% iEEG_signals = filter(b5, a5, iEEG_signals); +% +% end diff --git a/support_files/compute_distinguishability.m b/support_files/compute_distinguishability.m new file mode 100644 index 0000000..4762c75 --- /dev/null +++ b/support_files/compute_distinguishability.m @@ -0,0 +1,52 @@ +function [D_rs_vec] = compute_distinguishability(adj_matrices,resected_elecs,metric_type,WM_eliminate, patient_roi) + + % calculate number of patients + num_patients = length(adj_matrices); + + % loop through patients + for pt = 1:num_patients + + patient_WM_inds = find(patient_roi{pt}==9171); + + % extract resected electrode indices + pt_res_elecs = resected_elecs{pt}; + + % loop through frequency bands + for f = 1:5 + + % extract adjacency matrix + patient_adj = adj_matrices{pt}(f).data; + resected_elec_bool = zeros(size(patient_adj,1),1); + if WM_eliminate + patient_adj(patient_WM_inds,:) = []; + patient_adj(:,patient_WM_inds) = []; + % do resected elecs + + resected_elec_bool(resected_elecs{pt}) = 1; + resected_elec_bool(patient_WM_inds) = []; + pt_res_elecs = find(resected_elec_bool); + end + + if strcmp(metric_type,'node_strength') + pt_metric{pt}.data(:,f) = zscore(sum(patient_adj)); + elseif strcmp(metric_type,'betweenness_centrality') + pt_metric{pt}.data(:,f) = zscore(betweenness_wei(patient_adj)); + elseif strcmp(metric_type,'control_centrality') + pt_metric{pt}.data(:,f) = zscore(control_centrality(patient_adj)); + elseif strcmp(metric_type,'participation_coefficient') + [S,Q] = modularity_und(patient_adj,1); + pt_metric{pt}.data(:,f) = participation_coef(patient_adj,S); + end + + res_result = pt_metric{pt}.data(pt_res_elecs,f); + non_res_result = pt_metric{pt}.data(:,f); + non_res_result(pt_res_elecs) = []; + + mwu_result = mwwtest(res_result',non_res_result'); + + D_rs_vec(pt,f) = 1-max(mwu_result.U(2))./(length(res_result).*length(non_res_result)); + + end + end + +end \ No newline at end of file diff --git a/support_files/create_atlas_by_edge_rev.m b/support_files/create_atlas_by_edge_rev.m new file mode 100644 index 0000000..412e967 --- /dev/null +++ b/support_files/create_atlas_by_edge_rev.m @@ -0,0 +1,160 @@ +function [mean_conn, std_conn, num_samples, sem_conn, raw_atlas] = create_atlas_by_edge_rev(all_conn, all_roi, all_resect, region_list, threshold) +% [mean_conn, std_conn] = create_atlas_by_edge(all_conn, all_roi, all_resect, region_list) +% takes in an array of conectivity structs, an array of 3D mni coordinate +% arrays, an array of resected electrode vectors, and a vector containing +% all mni labels, and outputs two matrices representing the average and +% standard deviation of connections between all regions. +% +% Input: +% all_conn (cell): cell array containing patient connectivity structs in +% order +% all_roi (cell): cell array containing regions of interest corresponding +% to each electrode for each patient in order +% all_resect (cell): cell array containing patient resected electrode +% arrays in order +% region_list (double): array containing all region labels +% band (int): frequency band to be used +% threshold (int): minimum sample size required for edges to be +% incorporated into the atlas. Default value = 1 +% +% Output: +% mean_conn (double): (i,j) matrix of mean connectivity strengths between +% region_list(i) and region_list(j) +% std_conn (double): (i,j) matrix of standard deviations of connectivity +% strengths between region_list(i) and region_list(j) +% num_samples (double): a matrix where entry (i,j) is the number of +% samples used to calculate the edge weight between region_list(i) and region_list(j) +% +% John Bernabei and Ian Ong +% johnbe@seas.upenn.edu +% ianzyong@seas.upenn.edu +% 7/30/2020 + +% helper function +get_data = @(x) x(triu(true(size(x)),1)); + +% sets default threshold value if none is given +if ~exist('threshold','var'), threshold = 1; end + +% get number of patients +num_patients = length(all_conn); + +% get number of regions +num_regions = length(region_list); + +% initialize output arrays +mean_conn = NaN(num_regions); +std_conn = NaN(num_regions); +sem_conn = NaN(num_regions); +num_samples = NaN(num_regions); +raw_atlas = NaN(num_regions,num_regions,3576); % change this number + +% load in relevant matrices and resected electrodes for each patient +band_matrices = cell(1,num_patients); +resect_booleans = cell(1,num_patients); + +for p = 1:num_patients + % get resected electrodes for the patient + res_elec_inds = all_resect{p}; + % orient vector vertically + res_elec_size = size(res_elec_inds); + if res_elec_size(1) == 1 + res_elec_inds = res_elec_inds.'; + end + + % calculate logical with resected indices + try + resect_booleans{p} = cat(2,accumarray(res_elec_inds,1).',... + zeros(1,length(all_roi{p})-max(res_elec_inds))); + catch any_error + resect_booleans{p} = zeros(1,length(all_roi{p})); + end + + % extract band data + band_matrix{p} = all_conn{p}; +end + +% double loop through regions +for i = 1:num_regions % first region + + for j = 1:num_regions % second region + + % set up array to hold all edge values + edge_values_for_region_pair = cell(num_patients,1); + + % counter for number of data points for this region pair + sample_counter = 0; + + for p = 1:num_patients + + % get electrodes contained within first region + first_reg_elec = (all_roi{p} == region_list(i) & ~resect_booleans{p}); + + % extract rows corresponding to electrodes in the first region + current_matrix = band_matrix{p}; + first_reg_strengths = current_matrix(first_reg_elec,:); + + % get electrodes contained within second region + second_reg_elec = (all_roi{p} == region_list(j) & ~resect_booleans{p}); + + % extract connection strengths between the two regions + patient_strengths = first_reg_strengths(:,second_reg_elec); + + if i == j % loop case + % get only the upper triangular data since it's symmetric + + % I edited this to quantify self connections (within ROI) - John + same_roi_data = current_matrix(first_reg_elec, second_reg_elec); + patient_strengths = same_roi_data(:); %get_data(patient_strengths); + + else % assuming regions are spatially disjoint + patient_strengths = patient_strengths(:); + end + + % add number of edges to the counter + sample_counter = sample_counter + length(patient_strengths); + + % store electrode edges for this region pair in this patient + edge_values_for_region_pair{p} = patient_strengths; + + end + + % record mean, std, and sem if the number of edges exceeds the threshold + if sample_counter >= threshold + % take the average and standard deviation of all edge values for + % the region pair and place both in the corresponding output arrays + + + mean_conn(i,j) = mean(cell2mat(edge_values_for_region_pair)); + std_conn(i,j) = std(cell2mat(edge_values_for_region_pair)); + sem_conn(i,j) = std_conn(i,j)/sqrt(sample_counter); + num_edge_pairs = length(cell2mat(edge_values_for_region_pair)); + num_current_edges = sum(isnan(raw_atlas(i,j,:))); + raw_atlas(i,j,((num_current_edges+1):(num_current_edges+num_edge_pairs))) = cell2mat(edge_values_for_region_pair); + + end + + % record the number of edges for this region pair + num_samples(i,j) = sample_counter; + + if i==j + mean_conn(i,j) = NaN; + std_conn(i,j) = NaN; + sem_conn(i,j) = NaN; + raw_atlas(i,j,:) = NaN; + end + + end + +end + +raw_atlas(raw_atlas==0) = NaN; + +% symmetrize all output matrices +symmetrize = @(x) triu(x) + tril(x.',-1); +mean_conn = symmetrize(mean_conn); +std_conn = symmetrize(std_conn); +sem_conn = symmetrize(sem_conn); +num_samples = symmetrize(num_samples); + +end diff --git a/support_files/create_atlas_rev.m b/support_files/create_atlas_rev.m new file mode 100644 index 0000000..a5c099f --- /dev/null +++ b/support_files/create_atlas_rev.m @@ -0,0 +1,145 @@ +function [median_conn, std_conn, num_samples, sem_conn, raw_atlas] = create_atlas_rev(all_conn, all_roi, all_resect, region_list, threshold) +% [mean_conn, std_conn] = create_atlas(all_conn, all_roi, all_resect, region_list) +% takes in an array of conectivity structs, an array of 3D mni coordinate +% arrays, an array of resected electrode vectors, and a vector containing +% all mni labels, and outputs two matrices representing the average and +% standard deviation of connections between all regions. +% +% Input: +% all_conn (cell): cell array containing patient connectivity structs in +% order +% all_roi (cell): cell array containing regions of interest corresponding +% to each electrode for each patient in order +% all_resect (cell): cell array containing patient resected electrode +% arrays in order +% region_list (double): array containing all region labels +% band (int): frequency band to be used +% threshold (int): minimum sample size required for edges to be +% incorporated into the atlas. Default value = 1 +% +% Output: +% mean_conn (double): (i,j) matrix of mean connectivity strengths between +% region_list(i) and region_list(j) +% std_conn (double): (i,j) matrix of standard deviations of connectivity +% strengths between region_list(i) and region_list(j) +% num_samples (double): a matrix where entry (i,j) is the number of +% samples used to calculate the edge weight between region_list(i) and region_list(j) +% +% John Bernabei and Ian Ong +% johnbe@seas.upenn.edu +% ianzyong@seas.upenn.edu +% 7/5/2020 + +% sets default threshold value if none is given +if ~exist('threshold','var'), threshold = 1; end + +% get number of patients +num_patients = length(all_conn); +num_roi = length(region_list); + +% initialize output array to store connection strengths +median_conn = NaN(num_roi,num_roi,num_patients); + +for p = 1:num_patients + + % get electrode regions for the patient + patient_electrode_regions = all_roi{p}; + % get resected electrodes for the patient + res_elec_inds = all_resect{p}; + % orient vector vertically + res_elec_size = size(res_elec_inds); + if res_elec_size(1) == 1 + res_elec_inds = res_elec_inds.'; + end + + % calculate logical with resected indices + try + resect_boolean = cat(2,accumarray(res_elec_inds,1).',... + zeros(1,length(patient_electrode_regions)-max(res_elec_inds))); + catch any_error + resect_boolean = zeros(1,length(patient_electrode_regions)); + end + + % get connection strengths between each pair of electrodes + band_matrix = all_conn{p}; + + + % double loop through regions + for i = 1:num_roi % first region + + % get electrodes contained within first region + if isempty(resect_boolean) + p + end + first_reg_elec = (patient_electrode_regions == region_list(i) & ~resect_boolean); + + % extract rows corresponding to electrodes in the first region + first_reg_strengths = band_matrix(first_reg_elec,:); + +% % calculate connections within the region +% patient_strengths = first_reg_strengths(:,first_reg_elec); +% same_elec_mask = eye(size(patient_strengths)); +% patient_strengths(same_elec_mask==1) = NaN; +% patient_strength = nanmedian(patient_strengths); +% %nanmedian(patient_strengths(triu(true(size(patient_strengths)),1))); +% %nanmean(patient_strengths(triu(true(size(patient_strengths)),1))); +% +% % add to output array if the region contains any electrodes +% if ~isnan(patient_strength) +% median_conn(i,i,p) = patient_strength; +% else +% % skip calculations for this region pair if the first region +% % does not contain any electrodes +% continue +% end + + for j = 1:num_roi % second region (used to be (i+1):90) + + % get electrodes contained within second region + second_reg_elec = (patient_electrode_regions == region_list(j) & ~resect_boolean); + + % extract connection strengths between the two regions + patient_strengths = first_reg_strengths(:,second_reg_elec); + + % average connection strengths between the two regions + patient_strength = nanmedian(patient_strengths(:)); + + % add to output array + median_conn(i,j,p) = patient_strength; + + if i==j + median_conn(i,j,p) = NaN; + end + end + end +end + +raw_atlas = median_conn; + +% take the standard deviation element-wise to get the std matrix +std_conn = std(median_conn,0,3,'omitnan'); + +% get the number of samples available for each edge +num_samples = sum(~isnan(median_conn),3); + +% get the sem +sem_conn = std_conn./sqrt(num_samples); +sem_conn(isinf(sem_conn)) = NaN; + +% divide out the number of patients element-wise to get the median matrix +median_conn = median(median_conn,3,'omitnan'); + +% remove edge values and standard deviations for edges with sample sizes +% less than the threshold +median_conn(num_samples < threshold) = NaN; +std_conn(num_samples < threshold) = NaN; +sem_conn(num_samples < threshold) = NaN; + +% symmetrize all output matrices +% symmetrize = @(x) triu(x) + tril(x.',-1); +% mean_conn = symmetrize(mean_conn); +% std_conn = symmetrize(std_conn); +% sem_conn = symmetrize(sem_conn); +% num_samples = symmetrize(num_samples); + +end diff --git a/support_files/create_univariate_atlas.m b/support_files/create_univariate_atlas.m index 308d3e2..fd7001b 100644 --- a/support_files/create_univariate_atlas.m +++ b/support_files/create_univariate_atlas.m @@ -18,47 +18,47 @@ median_vec = NaN(length(region_list),1); std_vec = NaN(length(region_list),1); - if strcmp(method,'entropy') - - % filter in order to extract desired frequency band - if strcmp(freq_band,'delta') - % bandpass filter 0.5-4 Hz - [b1,a1] = butter(4,[0.5 4]/(200/2)); - filt_data = filter(b1, a1, raw_data); - - elseif strcmp(freq_band,'theta') - % bandpass filter 4-8 Hz - [b2,a2] = butter(4,[4 8]/(200/2)); - filt_data = filter(b2, a2, raw_data); - - elseif strcmp(freq_band,'alpha') - % bandpass filter 8-13 Hz - [b3,a3] = butter(4,[8 13]/(200/2)); - filt_data = filter(b3, a3, raw_data); - - elseif strcmp(freq_band,'beta') - % bandpass filter 13-30 Hz - [b4,a4] = butter(4,[13 30]/(200/2)); - filt_data = filter(b4, a4, raw_data); - - elseif strcmp(freq_band,'gamma') - % bandpass filter 30-80 Hz - [b5,a5] = butter(4,[30 80]/(200/2)); - filt_data = filter(b5, a5, raw_data); - - end + if strcmp(method,'entropy') +% +% % filter in order to extract desired frequency band +% if strcmp(freq_band,'delta') +% % bandpass filter 0.5-4 Hz +% [b1,a1] = butter(4,[0.5 4]/(200/2)); +% filt_data = filter(b1, a1, raw_data); +% +% elseif strcmp(freq_band,'theta') +% % bandpass filter 4-8 Hz +% [b2,a2] = butter(4,[4 8]/(200/2)); +% filt_data = filter(b2, a2, raw_data); +% +% elseif strcmp(freq_band,'alpha') +% % bandpass filter 8-13 Hz +% [b3,a3] = butter(4,[8 13]/(200/2)); +% filt_data = filter(b3, a3, raw_data); +% +% elseif strcmp(freq_band,'beta') +% % bandpass filter 13-30 Hz +% [b4,a4] = butter(4,[13 30]/(200/2)); +% filt_data = filter(b4, a4, raw_data); +% +% elseif strcmp(freq_band,'gamma') +% % bandpass filter 30-80 Hz +% [b5,a5] = butter(4,[30 80]/(200/2)); +% filt_data = filter(b5, a5, raw_data); +% +% end % loop through time and channels to create entropy matrix for i = 1:60 - for j = 1:size(filt_data,2) + for j = 1:size(raw_data,2) start_inds = (i-1)*200+1; end_inds = 200*i; - all_wentropy(i,j) = wentropy(filt_data((start_inds:end_inds),j),'shannon'); + all_wentropy(i,j) = wentropy(raw_data((start_inds:end_inds),j),'shannon'); end end % take median of each channel's entropy across time - feat_vals = median(real(log(all_wentropy)))'; + feat_vals = median(real(log(-1*all_wentropy)))'; elseif strcmp(method,'bandpower') diff --git a/support_files/create_univariate_atlas_rev.m b/support_files/create_univariate_atlas_rev.m new file mode 100644 index 0000000..a103891 --- /dev/null +++ b/support_files/create_univariate_atlas_rev.m @@ -0,0 +1,99 @@ +function [mean_vec,median_vec,std_vec,feat_vals] = create_univariate_atlas_rev(raw_data,region_list,channel_roi,method,freq_band) + + % raw_data is the underlying data, which should be 12000 x N_channels + % in this project, given that we are using 60 seconds of data + % downsampled to 200 Hz + + % region_list is a vector of all possible atlas ROI numerical labels + + % channel_roi is an N_channel x 1 vector of each channel's atlas ROI + + + % method is a string with the following possibilities: + % 'entropy' for shannon wavelet entropy + % 'bandpower' for spectral density + + + mean_vec = NaN(length(region_list),1); + median_vec = NaN(length(region_list),1); + std_vec = NaN(length(region_list),1); + + if strcmp(method,'entropy') +% +% % filter in order to extract desired frequency band +% if strcmp(freq_band,'delta') +% % bandpass filter 0.5-4 Hz +% [b1,a1] = butter(4,[0.5 4]/(200/2)); +% filt_data = filter(b1, a1, raw_data); +% +% elseif strcmp(freq_band,'theta') +% % bandpass filter 4-8 Hz +% [b2,a2] = butter(4,[4 8]/(200/2)); +% filt_data = filter(b2, a2, raw_data); +% +% elseif strcmp(freq_band,'alpha') +% % bandpass filter 8-13 Hz +% [b3,a3] = butter(4,[8 13]/(200/2)); +% filt_data = filter(b3, a3, raw_data); +% +% elseif strcmp(freq_band,'beta') +% % bandpass filter 13-30 Hz +% [b4,a4] = butter(4,[13 30]/(200/2)); +% filt_data = filter(b4, a4, raw_data); +% +% elseif strcmp(freq_band,'gamma') +% % bandpass filter 30-80 Hz +% [b5,a5] = butter(4,[30 80]/(200/2)); +% filt_data = filter(b5, a5, raw_data); +% +% end + + % loop through time and channels to create entropy matrix + for i = 1:60 + for j = 1:size(raw_data,2) + start_inds = (i-1)*200+1; + end_inds = 200*i; + all_wentropy(i,j) = wentropy(raw_data((start_inds:end_inds),j),'shannon'); + end + end + + % take median of each channel's entropy across time + feat_vals = median(real(log(-1*all_wentropy)))'; + + elseif strcmp(method,'bandpower') + + % we compute a normalized broadband spectral decomposition first + [pxx,f] = pwelch(raw_data,400,200,[0.5,2:80],200); + pxx_norm = real(normalize(pxx,'norm',1)); + + % extract spectral density in desired frequency band + if strcmp(freq_band,'delta') + % power in 1-4 Hz band; + feat_vals = median(pxx_norm(1:4,:))'; + elseif strcmp(freq_band,'theta') + feat_vals = median(pxx_norm(4:8,:))'; + elseif strcmp(freq_band,'alpha') + feat_vals = median(pxx_norm(8:13,:))'; + elseif strcmp(freq_band,'beta') + feat_vals = median(pxx_norm(13:30,:))'; + elseif strcmp(freq_band,'gamma') + feat_vals = median(pxx_norm(30:80,:))'; + end + + end + + % map to atlas + for r = 1:length(region_list) + % select region int and find which channels belong to it + this_region = region_list(r); + all_region_channels = find(channel_roi==this_region); + + % find feature vals in this region + region_feat_vals = feat_vals(all_region_channels); + + % assign mean, median, std + mean_vec(r) = nanmean(region_feat_vals); + median_vec(r) = nanmedian(region_feat_vals); + std_vec(r) = nanstd(region_feat_vals); + end +end \ No newline at end of file diff --git a/support_files/test_univariate_atlas_rev.m b/support_files/test_univariate_atlas_rev.m new file mode 100644 index 0000000..b9be32d --- /dev/null +++ b/support_files/test_univariate_atlas_rev.m @@ -0,0 +1,46 @@ +function [zscores, pt_zscores, pt_ind,roi_ind] = test_univariate_atlas_rev(mean_vec, std_vec, feat_vals, subject_vec, region_list, channel_roi,type) + % mean_vec: N_roi x 1 vector of means of each ROI + % std_vec: N_roi x 1 vector of std of each ROI + % feat_vals: N_channel x 1 feature values of each channel we want to test + % region_list: N_roi x 1 vector of ROI integer indices + % channel ROI: N_channel x 1 vctor of ROI integer indices + + + zscores = NaN(size(feat_vals)); + pt_zscores = []; + pt_ind = []; + roi_ind = []; + + for r = 1:length(region_list) + % select region int and find which channels belong to it + this_region = region_list(r); + all_region_channels = find(channel_roi==this_region); + + % find feature vals in this region + region_feat_vals = feat_vals(all_region_channels); + + % get region mean and std val + region_mean_val = mean_vec(r); + region_std_val = std_vec(r); + + zscores(all_region_channels,1) = (region_feat_vals - region_mean_val)./region_std_val; + + end + if strcmp(type,'patient') + for i = 1:max(subject_vec) + clear region_abnormality + this_pt_zscores = zscores(find(subject_vec==i)); + this_pt_roi = channel_roi(find(subject_vec==i)); + this_pt_unique_roi = unique(this_pt_roi); + for r = 1:length(this_pt_unique_roi) + this_unique_roi = this_pt_unique_roi(r); + region_abnormality(r,1) = max(this_pt_zscores(this_pt_roi==this_unique_roi)); + pt_ind = [pt_ind;i]; + roi_ind = [roi_ind;r]; + end + pt_zscores = [pt_zscores;region_abnormality]; + end + end + + +end \ No newline at end of file diff --git a/univariate_zscores.mat b/univariate_zscores.mat deleted file mode 100644 index 635b587..0000000 Binary files a/univariate_zscores.mat and /dev/null differ