From ad596e154d7d80bae0553712d604e2ae9487e927 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Mon, 21 Nov 2022 20:01:18 -0600 Subject: [PATCH 01/10] Update script names in readme --- README.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 9b5480c..c1519a6 100644 --- a/README.md +++ b/README.md @@ -14,7 +14,8 @@ Scripts to process the data, detect N1 responses and add tract information: Scripts to make the figure panels: - makeFig1A_plotMNI.m -- makeFig1BandC_plotResponsesAge.m +- makeFig1B_subjectResponses.m +- makeFig1C_heatmapResponsesAge.m - makeFig2and3_transmissionAcrossAge.m From 08154d8a0c749fb351719bcebebe879233f0fd97 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Wed, 23 Nov 2022 10:44:05 -0600 Subject: [PATCH 02/10] Remove to avoid accidental use of this function --- functions/read_tsv.m | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 functions/read_tsv.m diff --git a/functions/read_tsv.m b/functions/read_tsv.m deleted file mode 100644 index fde9161..0000000 --- a/functions/read_tsv.m +++ /dev/null @@ -1,3 +0,0 @@ -function tsv = read_tsv(filename) -fprintf('reading %s\n', filename); -tsv = readtable(filename, 'Delimiter', 'tab', 'FileType', 'text', 'ReadVariableNames', true); From e9f5d71489cb1d09407ebbe23c799c02a450d362 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Tue, 24 Jan 2023 15:33:35 -0600 Subject: [PATCH 03/10] Integrated supplement fig 2 into script for supplement 3 to 5 --- scripts/makeSupFig2_noN1s.m | 322 ------------------ ...responses.m => makeSupFig2to5_responses.m} | 96 ++++-- 2 files changed, 77 insertions(+), 341 deletions(-) delete mode 100644 scripts/makeSupFig2_noN1s.m rename scripts/{makeSupFig3to5_responses.m => makeSupFig2to5_responses.m} (77%) diff --git a/scripts/makeSupFig2_noN1s.m b/scripts/makeSupFig2_noN1s.m deleted file mode 100644 index 1de28c3..0000000 --- a/scripts/makeSupFig2_noN1s.m +++ /dev/null @@ -1,322 +0,0 @@ -% -% This script produces supplementary figures 2. It also checks how common -% N1s are detected on connections that are not governed by major fiber -% pathways according to Yeh et all., 2018 -% -% Dora Hermes, Dorien van Blooijs, Max van den Boom, 2022 - -clear -close all -clc - -%% -% Load the ccepData and ccepAverages from the derivatives - -myDataPath = setLocalDataPath(1); -if exist(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepData_V2.mat'), 'file') - load(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepData_V2.mat'), 'ccepData') -else - disp('Run scripts ccep02_aggregateToStruct.m and ccep03_addtracts.m first') -end -if exist(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepAverages.mat'), 'file') - load(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepAverages.mat')) -else - disp('Run script ccep04_averageConnections.m first') -end - -stimStimElec_excludeDist = 18; % the distance between the stimulated electrodes (in mm) above which N1s are excluded, 0 = not excluding -respStimElec_excludeDist = 13; % the distance between a stimulated and response electrode (in mm) within which N1s are excluded, 0 = not excluding - -% load tracts and their corresponding end-point ROIs -rois = ccep_categorizeAnatomicalRegions(); - - -%% -% Define the connection to be displayed and prepare the data -conn_matrix = {[2 1 0], [3 1 0], [3 2 0], [1 1 0]; ... - [2 1 1], [3 1 1], [3 2 1], [1 1 1]; ... - }; - -% get N1s for tested connections: -out = cell(size(conn_matrix)); -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - iTr = conn_matrix{iRow, iCol}(1); - iSubTr = conn_matrix{iRow, iCol}(2); - iDir = conn_matrix{iRow, iCol}(3); - - % extract the connection name (includes tract-subtract and direction) - strName = rois(iTr).tract_name; - subDir = split(rois(iTr).sub_tract(iSubTr).name, '-'); - strName = [strName, ' - ',[subDir{iDir + 1}, ' -> ', subDir{~iDir + 1}]]; - - % - disp(['Running - ', strName]); - out{iRow, iCol}.name = strName; - - - % - % latencies, number of N1s and ratio of N1s - % - - % extract the latencies, number of N1s and ratio of N1s between the end-point ROIs for a specific (sub-)tract and direction - metrics = ccep_N1sBetweenRegions( ccepData, ... - rois(iTr).sub_tract(iSubTr).(['roi', num2str(iDir + 1)]), ... - rois(iTr).sub_tract(iSubTr).(['roi', num2str(~iDir + 1)]), ... - stimStimElec_excludeDist, respStimElec_excludeDist); - - % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); - for iSubj = 1:length(metrics) - subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; - subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; - subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); - end - - out{iRow, iCol}.subjectsN1Values = subjectsN1Values; - - end -end - -%% -%% Generate supplementary figure 2 that displays the ratio of #N1s per #channels for each of the connections between the end-point areas -%% - -p_all = []; -r_all = []; -n_all = []; - -figure('position', [0 0 1200 600]) -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subjectsN1Values = out{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s - subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; - % plot(subjectsN1Values(~isnan(subjectsN1Values(:, 2)), 1), subjectsN1Values(~isnan(subjectsN1Values(:, 2)), 5), 'k.', 'MarkerSize', 10); - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 6); - - plot(subjectsN1Values(total_possible_N1>0, 1),number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0), 'k.', 'MarkerSize', 10); - - [r,p] = corr(subjectsN1Values(total_possible_N1>0, 1),number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),'Type','Spearman'); - - r_all = [r_all r]; - p_all = [p_all p]; - n_all = [n_all outInd]; - - title(strrep(out{iRow, iCol}.name, '_', '\_')); - - hold off; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); - -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; - xlim([0 60]); ylim([0 1]); - if iRow == size(conn_matrix, 1), xlabel('age'); end - if iCol == 1, ylabel('ratio N1s'); end - - % - text(40, 0.9, ['\rho=', num2str(r_all(outInd), 2) ]); - text(40, 0.8, ['P_f_d_r=', num2str(p_all_fdr(outInd), 2)]); - end -end - - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2_AgeVsRatioN1s'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - - - - - - - - - - - - - -%% get N1s for non-present connection -% Per Yeh et al., 2018) these connections should not be present -% Precentral <-> fusiform /// 29 G_precentral <-> 21 G_oc-temp_lat-fusifor -% Postcentral <-> inferior temporal gyrus /// 28 G_postcentral <-> 37 G_temporal_inf -% Supramarginal <-> fusiform /// 26 G_pariet_inf-Supramar <-> 21 G_oc-temp_lat-fusifor - -nonrois(1).roi1 = [29]; -nonrois(1).roi2 = [21]; -nonrois(1).name = 'precentral-fusiform'; -nonrois(2).roi1 = [28]; -nonrois(2).roi2 = [37]; -nonrois(2).name = 'postcentral-inferiortemporal'; -nonrois(3).roi1 = [26]; -nonrois(3).roi2 = [21]; -nonrois(3).name = 'supramarginal-fusiform'; - -dir_rois = [0 1]; - -nonout = cell(length(nonrois),2); -for iRow = 1:length(nonrois) - for iCol = 1:2 - - iDir = dir_rois(iCol); - - % extract the connection name (includes tract-subtract and direction) - subDir = split(nonrois(iRow).name, '-'); - strName = [subDir{iDir + 1}, ' -> ', subDir{~iDir + 1}]; - - % - disp(['Running - ', strName]); - nonout{iRow, iCol}.name = strName; - - % - % latencies, number of N1s and ratio of N1s - % - - % extract the latencies, number of N1s and ratio of N1s between the end-point ROIs for a specific (sub-)tract and direction - metrics = ccep_N1sBetweenRegions( ccepData, ... - nonrois(iRow).(['roi', num2str(iDir + 1)]), ... - nonrois(iRow).(['roi', num2str(~iDir + 1)]), ... - stimStimElec_excludeDist, respStimElec_excludeDist); - - % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); - for iSubj = 1:length(metrics) - subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; - subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; - subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); - end - - nonout{iRow, iCol}.subjectsN1Values = subjectsN1Values; - - end -end -clear subjectsN1Values - - -%% distribution plot of the percentage of N1s in a connection - -p_all = NaN(size(conn_matrix, 1),size(conn_matrix, 2)); - -figure('position', [0 0 200 300]) -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - % outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subjectsN1Values = out{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s - subplot(size(conn_matrix, 1), 1, iRow); hold on; - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 7); - - % boxplot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - distributionPlot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),... - 'xValues',iCol,'histOpt',2,'addSpread',1) - - [h,p,ci,stats] = ttest(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - - p_all(iRow,iCol) = p; - - hold off; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); -set(gca,'XTick',[1:4]) - -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - subplot(size(conn_matrix, 1), 1, iRow); hold on; - if p_all_fdr(iRow,iCol)<0.05 - plot(iCol,-0.1,'r*') - end - ylim([-0.15 1.02]) - xlim([0.5 size(conn_matrix, 2)+.5]) - end -end - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2b_RatioN1sConn'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - - -%% distribution plot of the percentage of N1s in a connection - -p_all = NaN(length(nonrois),2); - -figure('position', [0 0 160 300]) -for iRow = 1:length(nonrois) - for iCol = 1:2 - subjectsN1Values = nonout{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s -% subplot(2,length(nonrois), iCol); hold on; - subplot(2,1, iCol); hold on; - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 7); - -% boxplot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - distributionPlot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),... - 'xValues',iRow,'histOpt',2,'addSpread',1) - - [h,p,ci,stats] = ttest(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - - p_all(iRow,iCol) = p; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); - -for iRow = 1:length(nonrois) - for iCol = 1:2 - subplot(2,1, iCol); hold on; - if p_all_fdr(iRow,iCol)<0.05 - plot(iCol,-0.1,'r*') - end - xlim([0.5 3.5]) - ylim([-0.15 1.02]) - end -end - -set(gca,'XTick',[1:3]) - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2c_RatioN1sNoconn'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - diff --git a/scripts/makeSupFig3to5_responses.m b/scripts/makeSupFig2to5_responses.m similarity index 77% rename from scripts/makeSupFig3to5_responses.m rename to scripts/makeSupFig2to5_responses.m index d1dfb4f..97a2475 100644 --- a/scripts/makeSupFig3to5_responses.m +++ b/scripts/makeSupFig2to5_responses.m @@ -40,7 +40,7 @@ all_varlat_p = []; all_meanvarlat_p = []; -% all_ratioN1s_p = []; +all_ratioN1s_p = []; all_n1Widths_p = []; % @@ -71,13 +71,25 @@ stimStimElec_excludeDist, respStimElec_excludeDist); % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); + % x + subjectsN1Values = NaN(length(metrics), 8); for iSubj = 1:length(metrics) subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); + subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); % mean in latencies + subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); % variance in latencies + subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); % variance in (latencies * 1000) + subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; % number or response electrodes on ROI + subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; % number or stim electrodes on ROI + subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); % number of N1s + + % relative number of N1s + total_potential_num_N1s = metrics(iSubj).numElecRespROI * metrics(iSubj).numElecStimROI; + if total_potential_num_N1s == 0 + subjectsN1Values(iSubj, 8) = nan; + else + subjectsN1Values(iSubj, 8) = length(metrics(iSubj).latencies) / total_potential_num_N1s; + end + end out{iRow, iCol}.subjectsN1Values = subjectsN1Values; @@ -95,14 +107,12 @@ out{iRow, iCol}.meanvarlat_p = meanvarlat_p; all_meanvarlat_p(end + 1, :) = [iRow, iCol, meanvarlat_p]; - - % Figure S2 moved to separate script, can later be placed back here -% % calculate and store the r and p for the age vs ratio of N1s -% [ratioN1s_r, ratioN1s_p] = corr(subjectsN1Values(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 5)), 1), ... -% subjectsN1Values(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 5)), 5), 'Type', 'Spearman'); -% out{iRow, iCol}.ratioN1s_r = ratioN1s_r; -% out{iRow, iCol}.ratioN1s_p = ratioN1s_p; -% all_ratioN1s_p(end + 1, :) = [iRow, iCol, ratioN1s_p]; + % calculate and store the r and p for the age vs ratio of N1s + [ratioN1s_r, ratioN1s_p] = corr(subjectsN1Values(~isnan(subjectsN1Values(:, 1)) & ~isnan(subjectsN1Values(:, 8)), 1), ... + subjectsN1Values(~isnan(subjectsN1Values(:, 1)) & ~isnan(subjectsN1Values(:, 8)), 8), 'Type', 'Spearman'); + out{iRow, iCol}.ratioN1s_r = ratioN1s_r; + out{iRow, iCol}.ratioN1s_p = ratioN1s_p; + all_ratioN1s_p(end + 1, :) = [iRow, iCol, ratioN1s_p]; % @@ -127,17 +137,56 @@ % calculate the FDR corrected P-values [~, ~, ~, all_varlat_p(:, 4)] = fdr_bh(all_varlat_p(:, 3), 0.05, 'pdep'); [~, ~, ~, all_meanvarlat_p(:, 4)] = fdr_bh(all_meanvarlat_p(:, 3), 0.05, 'pdep'); -% [~, ~, ~, all_ratioN1s_p(:, 4)] = fdr_bh(all_ratioN1s_p(:, 3), 0.05, 'pdep'); +[~, ~, ~, all_ratioN1s_p(:, 4)] = fdr_bh(all_ratioN1s_p(:, 3), 0.05, 'pdep'); [~, ~, ~, all_n1Widths_p(:, 4)] = fdr_bh(all_n1Widths_p(:, 3), 0.05, 'pdep'); for iP = 1:size(all_varlat_p, 1) out{all_varlat_p(iP, 1), all_varlat_p(iP, 2)}.varlat_p_fdr = all_varlat_p(iP, 4); out{all_meanvarlat_p(iP, 1), all_meanvarlat_p(iP, 2)}.meanvarlat_p_fdr = all_meanvarlat_p(iP, 4); -% out{all_ratioN1s_p(iP, 1), all_ratioN1s_p(iP, 2)}.ratioN1s_p_fdr = all_ratioN1s_p(iP, 4); + out{all_ratioN1s_p(iP, 1), all_ratioN1s_p(iP, 2)}.ratioN1s_p_fdr = all_ratioN1s_p(iP, 4); out{all_n1Widths_p(iP, 1), all_n1Widths_p(iP, 2)}.n1Widths_p_fdr = all_n1Widths_p(iP, 4); end +%% +% Generate supplementary figure 2 that displays the ratio of #N1s per #channels for each of the connections between the end-point areas + +figure('position', [0 0 1200 600]) +for iRow = 1:size(conn_matrix, 1) + for iCol = 1:size(conn_matrix, 2) + outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; + subjectsN1Values = out{iRow, iCol}.subjectsN1Values; + + % age vs relative number of N1s + subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; + plot(subjectsN1Values(:, 1), subjectsN1Values(:, 8), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 8))); + + % + title(strrep(out{iRow, iCol}.name, '_', '\_')); + xlim([0 60]); ylim([0 1]); + if iRow == size(conn_matrix, 1), xlabel('Age (years)'); end + if iCol == 1, ylabel('Relative number of N1s'); end + + % + text(38, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.ratioN1s_r, 2)), ' (n=', num2str(n, 2), ')']); + text(38, 0.8, ['P_f_d_r=', num2str(round(out{iRow, iCol}.ratioN1s_p_fdr, 2))]); + + hold off; + end +end + + +figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2_AgeVsRatioN1s'); +set(gcf,'PaperPositionMode', 'auto') +print('-dpng', '-r300', figureName) +print('-depsc', '-r300', figureName) + + + + %% % Generate supplementary figure 3 that displays the variance in latencies for each of the connections between the end-point areas @@ -150,6 +199,9 @@ % age vs variance latencies subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1Values(:, 1), 1000 * subjectsN1Values(:, 3), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 3))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -158,7 +210,7 @@ if iCol == 1, ylabel('Varience in N1 latency (ms)'); end % - text(40, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.varlat_r, 2))]); + text(40, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.varlat_r, 2)), ' (n=', num2str(n, 2), ')']); text(40, 0.8, ['P_f_d_r=', num2str(round(out{iRow, iCol}.varlat_p_fdr, 2))]); hold off; @@ -185,6 +237,9 @@ % mean latency (* 1000) vs latency variance (base values already multiplied before by 1000) subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1Values(:, 2) * 1000, subjectsN1Values(:, 4), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 4))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -194,7 +249,7 @@ if iCol == 1, ylabel('Varience in N1 latency (ms)'); end % - text(15, max(subjectsN1Values(:, 4) * .95), ['\rho=', num2str(round(out{iRow, iCol}.meanvarlat_r, 2))]); + text(15, max(subjectsN1Values(:, 4) * .95), ['\rho=', num2str(round(out{iRow, iCol}.meanvarlat_r, 2)), ' (n=', num2str(n, 2), ')']); if out{iRow, iCol}.meanvarlat_p_fdr < .001 text(15, max(subjectsN1Values(:, 4) * .85), 'P_f_d_r < 0.001'); elseif out{iRow, iCol}.meanvarlat_p_fdr < .01 @@ -238,6 +293,9 @@ % mean latency (* 1000) vs FWGM (* 1000) subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1WidthValues(:, 2) * 1000, subjectsN1WidthValues(:, 3) * 1000, 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1WidthValues(:, 2)) & ~isnan(subjectsN1WidthValues(:, 3))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -247,7 +305,7 @@ ylim([5, (max(subjectsN1WidthValues(:, 3) * 1000)) + 2]) % - text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .95), ['\rho=', num2str(round(out{iRow, iCol}.n1Widths_r, 2))]); + text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .95), ['\rho=', num2str(round(out{iRow, iCol}.n1Widths_r, 2)), ' (n=', num2str(n, 2), ')']); if out{iRow, iCol}.n1Widths_p_fdr < .001 text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .85), 'P_f_d_r < 0.001'); elseif out{iRow, iCol}.n1Widths_p_fdr < .01 From b446ad1e9b4d24ba4fcbe68f8f09f5abf3c96ed2 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Tue, 24 Jan 2023 15:35:11 -0600 Subject: [PATCH 04/10] Added sample size display in fig 2-5, and integrated supplement fig 2 into script for supplement 3 to 5 --- scripts/makeSupFig2_noN1s.m | 322 ------------------ ...responses.m => makeSupFig2to5_responses.m} | 96 ++++-- 2 files changed, 77 insertions(+), 341 deletions(-) delete mode 100644 scripts/makeSupFig2_noN1s.m rename scripts/{makeSupFig3to5_responses.m => makeSupFig2to5_responses.m} (77%) diff --git a/scripts/makeSupFig2_noN1s.m b/scripts/makeSupFig2_noN1s.m deleted file mode 100644 index 1de28c3..0000000 --- a/scripts/makeSupFig2_noN1s.m +++ /dev/null @@ -1,322 +0,0 @@ -% -% This script produces supplementary figures 2. It also checks how common -% N1s are detected on connections that are not governed by major fiber -% pathways according to Yeh et all., 2018 -% -% Dora Hermes, Dorien van Blooijs, Max van den Boom, 2022 - -clear -close all -clc - -%% -% Load the ccepData and ccepAverages from the derivatives - -myDataPath = setLocalDataPath(1); -if exist(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepData_V2.mat'), 'file') - load(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepData_V2.mat'), 'ccepData') -else - disp('Run scripts ccep02_aggregateToStruct.m and ccep03_addtracts.m first') -end -if exist(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepAverages.mat'), 'file') - load(fullfile(myDataPath.output, 'derivatives', 'av_ccep', 'ccepAverages.mat')) -else - disp('Run script ccep04_averageConnections.m first') -end - -stimStimElec_excludeDist = 18; % the distance between the stimulated electrodes (in mm) above which N1s are excluded, 0 = not excluding -respStimElec_excludeDist = 13; % the distance between a stimulated and response electrode (in mm) within which N1s are excluded, 0 = not excluding - -% load tracts and their corresponding end-point ROIs -rois = ccep_categorizeAnatomicalRegions(); - - -%% -% Define the connection to be displayed and prepare the data -conn_matrix = {[2 1 0], [3 1 0], [3 2 0], [1 1 0]; ... - [2 1 1], [3 1 1], [3 2 1], [1 1 1]; ... - }; - -% get N1s for tested connections: -out = cell(size(conn_matrix)); -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - iTr = conn_matrix{iRow, iCol}(1); - iSubTr = conn_matrix{iRow, iCol}(2); - iDir = conn_matrix{iRow, iCol}(3); - - % extract the connection name (includes tract-subtract and direction) - strName = rois(iTr).tract_name; - subDir = split(rois(iTr).sub_tract(iSubTr).name, '-'); - strName = [strName, ' - ',[subDir{iDir + 1}, ' -> ', subDir{~iDir + 1}]]; - - % - disp(['Running - ', strName]); - out{iRow, iCol}.name = strName; - - - % - % latencies, number of N1s and ratio of N1s - % - - % extract the latencies, number of N1s and ratio of N1s between the end-point ROIs for a specific (sub-)tract and direction - metrics = ccep_N1sBetweenRegions( ccepData, ... - rois(iTr).sub_tract(iSubTr).(['roi', num2str(iDir + 1)]), ... - rois(iTr).sub_tract(iSubTr).(['roi', num2str(~iDir + 1)]), ... - stimStimElec_excludeDist, respStimElec_excludeDist); - - % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); - for iSubj = 1:length(metrics) - subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; - subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; - subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); - end - - out{iRow, iCol}.subjectsN1Values = subjectsN1Values; - - end -end - -%% -%% Generate supplementary figure 2 that displays the ratio of #N1s per #channels for each of the connections between the end-point areas -%% - -p_all = []; -r_all = []; -n_all = []; - -figure('position', [0 0 1200 600]) -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subjectsN1Values = out{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s - subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; - % plot(subjectsN1Values(~isnan(subjectsN1Values(:, 2)), 1), subjectsN1Values(~isnan(subjectsN1Values(:, 2)), 5), 'k.', 'MarkerSize', 10); - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 6); - - plot(subjectsN1Values(total_possible_N1>0, 1),number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0), 'k.', 'MarkerSize', 10); - - [r,p] = corr(subjectsN1Values(total_possible_N1>0, 1),number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),'Type','Spearman'); - - r_all = [r_all r]; - p_all = [p_all p]; - n_all = [n_all outInd]; - - title(strrep(out{iRow, iCol}.name, '_', '\_')); - - hold off; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); - -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; - xlim([0 60]); ylim([0 1]); - if iRow == size(conn_matrix, 1), xlabel('age'); end - if iCol == 1, ylabel('ratio N1s'); end - - % - text(40, 0.9, ['\rho=', num2str(r_all(outInd), 2) ]); - text(40, 0.8, ['P_f_d_r=', num2str(p_all_fdr(outInd), 2)]); - end -end - - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2_AgeVsRatioN1s'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - - - - - - - - - - - - - -%% get N1s for non-present connection -% Per Yeh et al., 2018) these connections should not be present -% Precentral <-> fusiform /// 29 G_precentral <-> 21 G_oc-temp_lat-fusifor -% Postcentral <-> inferior temporal gyrus /// 28 G_postcentral <-> 37 G_temporal_inf -% Supramarginal <-> fusiform /// 26 G_pariet_inf-Supramar <-> 21 G_oc-temp_lat-fusifor - -nonrois(1).roi1 = [29]; -nonrois(1).roi2 = [21]; -nonrois(1).name = 'precentral-fusiform'; -nonrois(2).roi1 = [28]; -nonrois(2).roi2 = [37]; -nonrois(2).name = 'postcentral-inferiortemporal'; -nonrois(3).roi1 = [26]; -nonrois(3).roi2 = [21]; -nonrois(3).name = 'supramarginal-fusiform'; - -dir_rois = [0 1]; - -nonout = cell(length(nonrois),2); -for iRow = 1:length(nonrois) - for iCol = 1:2 - - iDir = dir_rois(iCol); - - % extract the connection name (includes tract-subtract and direction) - subDir = split(nonrois(iRow).name, '-'); - strName = [subDir{iDir + 1}, ' -> ', subDir{~iDir + 1}]; - - % - disp(['Running - ', strName]); - nonout{iRow, iCol}.name = strName; - - % - % latencies, number of N1s and ratio of N1s - % - - % extract the latencies, number of N1s and ratio of N1s between the end-point ROIs for a specific (sub-)tract and direction - metrics = ccep_N1sBetweenRegions( ccepData, ... - nonrois(iRow).(['roi', num2str(iDir + 1)]), ... - nonrois(iRow).(['roi', num2str(~iDir + 1)]), ... - stimStimElec_excludeDist, respStimElec_excludeDist); - - % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); - for iSubj = 1:length(metrics) - subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; - subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; - subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); - end - - nonout{iRow, iCol}.subjectsN1Values = subjectsN1Values; - - end -end -clear subjectsN1Values - - -%% distribution plot of the percentage of N1s in a connection - -p_all = NaN(size(conn_matrix, 1),size(conn_matrix, 2)); - -figure('position', [0 0 200 300]) -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - % outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; - subjectsN1Values = out{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s - subplot(size(conn_matrix, 1), 1, iRow); hold on; - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 7); - - % boxplot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - distributionPlot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),... - 'xValues',iCol,'histOpt',2,'addSpread',1) - - [h,p,ci,stats] = ttest(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - - p_all(iRow,iCol) = p; - - hold off; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); -set(gca,'XTick',[1:4]) - -for iRow = 1:size(conn_matrix, 1) - for iCol = 1:size(conn_matrix, 2) - subplot(size(conn_matrix, 1), 1, iRow); hold on; - if p_all_fdr(iRow,iCol)<0.05 - plot(iCol,-0.1,'r*') - end - ylim([-0.15 1.02]) - xlim([0.5 size(conn_matrix, 2)+.5]) - end -end - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2b_RatioN1sConn'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - - -%% distribution plot of the percentage of N1s in a connection - -p_all = NaN(length(nonrois),2); - -figure('position', [0 0 160 300]) -for iRow = 1:length(nonrois) - for iCol = 1:2 - subjectsN1Values = nonout{iRow, iCol}.subjectsN1Values; - - % Plot age vs ratio of N1s -% subplot(2,length(nonrois), iCol); hold on; - subplot(2,1, iCol); hold on; - - % total possible N1s: - total_possible_N1 = subjectsN1Values(:, 5) .* subjectsN1Values(:, 6); - - % number of N1s: - number_N1 = subjectsN1Values(:, 7); - -% boxplot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - distributionPlot(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0),... - 'xValues',iRow,'histOpt',2,'addSpread',1) - - [h,p,ci,stats] = ttest(number_N1(total_possible_N1>0)./total_possible_N1(total_possible_N1>0)); - - p_all(iRow,iCol) = p; - end -end - -[~,~,~,p_all_fdr] = fdr_bh(p_all, 0.05, 'pdep'); - -for iRow = 1:length(nonrois) - for iCol = 1:2 - subplot(2,1, iCol); hold on; - if p_all_fdr(iRow,iCol)<0.05 - plot(iCol,-0.1,'r*') - end - xlim([0.5 3.5]) - ylim([-0.15 1.02]) - end -end - -set(gca,'XTick',[1:3]) - -figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2c_RatioN1sNoconn'); -set(gcf,'PaperPositionMode', 'auto') -print('-dpng', '-r300', figureName) -print('-depsc', '-r300', figureName) - - diff --git a/scripts/makeSupFig3to5_responses.m b/scripts/makeSupFig2to5_responses.m similarity index 77% rename from scripts/makeSupFig3to5_responses.m rename to scripts/makeSupFig2to5_responses.m index d1dfb4f..97a2475 100644 --- a/scripts/makeSupFig3to5_responses.m +++ b/scripts/makeSupFig2to5_responses.m @@ -40,7 +40,7 @@ all_varlat_p = []; all_meanvarlat_p = []; -% all_ratioN1s_p = []; +all_ratioN1s_p = []; all_n1Widths_p = []; % @@ -71,13 +71,25 @@ stimStimElec_excludeDist, respStimElec_excludeDist); % retrieve metrics per subject, output format: - % x - subjectsN1Values = NaN(length(metrics), 4); + % x + subjectsN1Values = NaN(length(metrics), 8); for iSubj = 1:length(metrics) subjectsN1Values(iSubj, 1) = metrics(iSubj).age; - subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); - subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); + subjectsN1Values(iSubj, 2) = mean(metrics(iSubj).latencies, 'omitnan'); % mean in latencies + subjectsN1Values(iSubj, 3) = var(metrics(iSubj).latencies, 'omitnan'); % variance in latencies + subjectsN1Values(iSubj, 4) = var(1000 * metrics(iSubj).latencies, 'omitnan'); % variance in (latencies * 1000) + subjectsN1Values(iSubj, 5) = metrics(iSubj).numElecRespROI; % number or response electrodes on ROI + subjectsN1Values(iSubj, 6) = metrics(iSubj).numElecStimROI; % number or stim electrodes on ROI + subjectsN1Values(iSubj, 7) = length(metrics(iSubj).latencies); % number of N1s + + % relative number of N1s + total_potential_num_N1s = metrics(iSubj).numElecRespROI * metrics(iSubj).numElecStimROI; + if total_potential_num_N1s == 0 + subjectsN1Values(iSubj, 8) = nan; + else + subjectsN1Values(iSubj, 8) = length(metrics(iSubj).latencies) / total_potential_num_N1s; + end + end out{iRow, iCol}.subjectsN1Values = subjectsN1Values; @@ -95,14 +107,12 @@ out{iRow, iCol}.meanvarlat_p = meanvarlat_p; all_meanvarlat_p(end + 1, :) = [iRow, iCol, meanvarlat_p]; - - % Figure S2 moved to separate script, can later be placed back here -% % calculate and store the r and p for the age vs ratio of N1s -% [ratioN1s_r, ratioN1s_p] = corr(subjectsN1Values(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 5)), 1), ... -% subjectsN1Values(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 5)), 5), 'Type', 'Spearman'); -% out{iRow, iCol}.ratioN1s_r = ratioN1s_r; -% out{iRow, iCol}.ratioN1s_p = ratioN1s_p; -% all_ratioN1s_p(end + 1, :) = [iRow, iCol, ratioN1s_p]; + % calculate and store the r and p for the age vs ratio of N1s + [ratioN1s_r, ratioN1s_p] = corr(subjectsN1Values(~isnan(subjectsN1Values(:, 1)) & ~isnan(subjectsN1Values(:, 8)), 1), ... + subjectsN1Values(~isnan(subjectsN1Values(:, 1)) & ~isnan(subjectsN1Values(:, 8)), 8), 'Type', 'Spearman'); + out{iRow, iCol}.ratioN1s_r = ratioN1s_r; + out{iRow, iCol}.ratioN1s_p = ratioN1s_p; + all_ratioN1s_p(end + 1, :) = [iRow, iCol, ratioN1s_p]; % @@ -127,17 +137,56 @@ % calculate the FDR corrected P-values [~, ~, ~, all_varlat_p(:, 4)] = fdr_bh(all_varlat_p(:, 3), 0.05, 'pdep'); [~, ~, ~, all_meanvarlat_p(:, 4)] = fdr_bh(all_meanvarlat_p(:, 3), 0.05, 'pdep'); -% [~, ~, ~, all_ratioN1s_p(:, 4)] = fdr_bh(all_ratioN1s_p(:, 3), 0.05, 'pdep'); +[~, ~, ~, all_ratioN1s_p(:, 4)] = fdr_bh(all_ratioN1s_p(:, 3), 0.05, 'pdep'); [~, ~, ~, all_n1Widths_p(:, 4)] = fdr_bh(all_n1Widths_p(:, 3), 0.05, 'pdep'); for iP = 1:size(all_varlat_p, 1) out{all_varlat_p(iP, 1), all_varlat_p(iP, 2)}.varlat_p_fdr = all_varlat_p(iP, 4); out{all_meanvarlat_p(iP, 1), all_meanvarlat_p(iP, 2)}.meanvarlat_p_fdr = all_meanvarlat_p(iP, 4); -% out{all_ratioN1s_p(iP, 1), all_ratioN1s_p(iP, 2)}.ratioN1s_p_fdr = all_ratioN1s_p(iP, 4); + out{all_ratioN1s_p(iP, 1), all_ratioN1s_p(iP, 2)}.ratioN1s_p_fdr = all_ratioN1s_p(iP, 4); out{all_n1Widths_p(iP, 1), all_n1Widths_p(iP, 2)}.n1Widths_p_fdr = all_n1Widths_p(iP, 4); end +%% +% Generate supplementary figure 2 that displays the ratio of #N1s per #channels for each of the connections between the end-point areas + +figure('position', [0 0 1200 600]) +for iRow = 1:size(conn_matrix, 1) + for iCol = 1:size(conn_matrix, 2) + outInd = (iRow - 1) * size(conn_matrix, 2) + iCol; + subjectsN1Values = out{iRow, iCol}.subjectsN1Values; + + % age vs relative number of N1s + subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; + plot(subjectsN1Values(:, 1), subjectsN1Values(:, 8), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 8))); + + % + title(strrep(out{iRow, iCol}.name, '_', '\_')); + xlim([0 60]); ylim([0 1]); + if iRow == size(conn_matrix, 1), xlabel('Age (years)'); end + if iCol == 1, ylabel('Relative number of N1s'); end + + % + text(38, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.ratioN1s_r, 2)), ' (n=', num2str(n, 2), ')']); + text(38, 0.8, ['P_f_d_r=', num2str(round(out{iRow, iCol}.ratioN1s_p_fdr, 2))]); + + hold off; + end +end + + +figureName = fullfile(myDataPath.output, 'derivatives', 'age', 'SupFigS2_AgeVsRatioN1s'); +set(gcf,'PaperPositionMode', 'auto') +print('-dpng', '-r300', figureName) +print('-depsc', '-r300', figureName) + + + + %% % Generate supplementary figure 3 that displays the variance in latencies for each of the connections between the end-point areas @@ -150,6 +199,9 @@ % age vs variance latencies subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1Values(:, 1), 1000 * subjectsN1Values(:, 3), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 3))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -158,7 +210,7 @@ if iCol == 1, ylabel('Varience in N1 latency (ms)'); end % - text(40, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.varlat_r, 2))]); + text(40, 0.9, ['\rho=', num2str(round(out{iRow, iCol}.varlat_r, 2)), ' (n=', num2str(n, 2), ')']); text(40, 0.8, ['P_f_d_r=', num2str(round(out{iRow, iCol}.varlat_p_fdr, 2))]); hold off; @@ -185,6 +237,9 @@ % mean latency (* 1000) vs latency variance (base values already multiplied before by 1000) subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1Values(:, 2) * 1000, subjectsN1Values(:, 4), 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1Values(:, 2)) & ~isnan(subjectsN1Values(:, 4))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -194,7 +249,7 @@ if iCol == 1, ylabel('Varience in N1 latency (ms)'); end % - text(15, max(subjectsN1Values(:, 4) * .95), ['\rho=', num2str(round(out{iRow, iCol}.meanvarlat_r, 2))]); + text(15, max(subjectsN1Values(:, 4) * .95), ['\rho=', num2str(round(out{iRow, iCol}.meanvarlat_r, 2)), ' (n=', num2str(n, 2), ')']); if out{iRow, iCol}.meanvarlat_p_fdr < .001 text(15, max(subjectsN1Values(:, 4) * .85), 'P_f_d_r < 0.001'); elseif out{iRow, iCol}.meanvarlat_p_fdr < .01 @@ -238,6 +293,9 @@ % mean latency (* 1000) vs FWGM (* 1000) subplot(size(conn_matrix, 1), size(conn_matrix, 2), outInd); hold on; plot(subjectsN1WidthValues(:, 2) * 1000, subjectsN1WidthValues(:, 3) * 1000, 'k.', 'MarkerSize', 10); + + % determine the number of samples + n = sum(~isnan(subjectsN1WidthValues(:, 2)) & ~isnan(subjectsN1WidthValues(:, 3))); % title(strrep(out{iRow, iCol}.name, '_', '\_')); @@ -247,7 +305,7 @@ ylim([5, (max(subjectsN1WidthValues(:, 3) * 1000)) + 2]) % - text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .95), ['\rho=', num2str(round(out{iRow, iCol}.n1Widths_r, 2))]); + text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .95), ['\rho=', num2str(round(out{iRow, iCol}.n1Widths_r, 2)), ' (n=', num2str(n, 2), ')']); if out{iRow, iCol}.n1Widths_p_fdr < .001 text(75, max(subjectsN1WidthValues(:, 3) * 1000 * .85), 'P_f_d_r < 0.001'); elseif out{iRow, iCol}.n1Widths_p_fdr < .01 From bf2dc3b8cadfc8f4a1fd7ce1cc83bb233d97224e Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Tue, 24 Jan 2023 15:44:55 -0600 Subject: [PATCH 05/10] Rename to match figures --- ...heatmapResponsesAge.m => makeFig1Cand3A_heatmapResponsesAge.m} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/{makeFig1C_heatmapResponsesAge.m => makeFig1Cand3A_heatmapResponsesAge.m} (100%) diff --git a/scripts/makeFig1C_heatmapResponsesAge.m b/scripts/makeFig1Cand3A_heatmapResponsesAge.m similarity index 100% rename from scripts/makeFig1C_heatmapResponsesAge.m rename to scripts/makeFig1Cand3A_heatmapResponsesAge.m From 9cf3b9987a895a5b99cbd70d7999f2abf4838529 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Tue, 24 Jan 2023 17:16:24 -0600 Subject: [PATCH 06/10] Added sample size display for fig 1, 2 and 3 --- scripts/makeFig1Cand3A_heatmapResponsesAge.m | 6 ++-- scripts/makeFig2and3_transmissionAcrossAge.m | 33 +++++++++++--------- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/scripts/makeFig1Cand3A_heatmapResponsesAge.m b/scripts/makeFig1Cand3A_heatmapResponsesAge.m index cf0ac2e..c9aab6c 100644 --- a/scripts/makeFig1Cand3A_heatmapResponsesAge.m +++ b/scripts/makeFig1Cand3A_heatmapResponsesAge.m @@ -1,6 +1,6 @@ % % -% This code is used to create Figure 1C - the normalized CCEPs of all patients in order of age, +% This code is used to the plot for Figure 1C and Figure 3A - the normalized CCEPs of all patients in order of age, % % @@ -71,7 +71,7 @@ %% -% Figure 1C - Plots for each stimulated and responding region with normalized CCEPs + N1 sorted by age +% Plots for each stimulated and responding region with normalized CCEPs + N1 sorted by age ttmin = 0.010; ttmax = .100; @@ -101,7 +101,7 @@ set(gca,'XTick',20:20:80,'YTick',[]) axis tight - strSign = [' (p\_fdr = ', num2str(all_p_fdr{iTr}{iSubTr}(iDir + 1)), ')']; + strSign = [' (p\_fdr = ', num2str(all_p_fdr{iTr}{iSubTr}(iDir + 1)), ', n=', num2str(length(sortedCCEPs{iTr}{iSubTr}{iDir + 1}.age)), ')']; if all_p_fdr{iTr}{iSubTr}(iDir + 1) < .05, strSign = [strSign, ' *']; end title([rois(iTr).tract_name, ' - ', strSubTitle, strSign]); diff --git a/scripts/makeFig2and3_transmissionAcrossAge.m b/scripts/makeFig2and3_transmissionAcrossAge.m index eb4dcb3..5f9e34a 100644 --- a/scripts/makeFig2and3_transmissionAcrossAge.m +++ b/scripts/makeFig2and3_transmissionAcrossAge.m @@ -1,5 +1,5 @@ % -% This script produces Figure 3 of the article +% This script produces the plots for Figure 2 and 3 (B/C) of the article % @@ -127,6 +127,11 @@ n1LatencyMeans(exclMeans) = []; n1SpeedMeans(exclMeans) = []; + % store the number of age-samples (age can be averages over more than one subject) + out{iTr}{iSubTr}{iDir + 1}.numAges = length(ages); + + + % initialize statistics variables lat_cross_linear = NaN(length(n1LatencyMeans), 4); % x spd_cross_linear = NaN(length(n1LatencyMeans), 4); % x lat_cross_second = NaN(length(n1LatencyMeans), 5); % x @@ -441,15 +446,15 @@ % Display in command window the cod and delta for each subplot % this info is displayed in Figure 3 as well. -long_tracts = {'Y1065_TPAT - Parietal -> Temporal','Y1065_TPAT - Temporal -> Parietal',... - 'Y1065_AF - Frontal -> Temporal','Y1065_AF - Temporal -> Frontal',... - 'Y1065_SLF2 - Frontal -> Parietal','Y1065_SLF2 - Parietal -> Frontal',... - 'Y1065_SLF2 - Frontal -> Central','Y1065_SLF2 - Central -> Frontal'}; +long_tracts = { 'Y1065_TPAT - Parietal -> Temporal','Y1065_TPAT - Temporal -> Parietal',... + 'Y1065_AF - Frontal -> Temporal','Y1065_AF - Temporal -> Frontal',... + 'Y1065_SLF2 - Frontal -> Parietal','Y1065_SLF2 - Parietal -> Frontal',... + 'Y1065_SLF2 - Frontal -> Central','Y1065_SLF2 - Central -> Frontal'}; long_cod_lat = NaN(1,length(long_tracts)); long_cod_spd = NaN(1,length(long_tracts)); U_tracts = {'Y842_U - PreCentral -> PostCentral','Y842_U - PostCentral -> PreCentral',... - 'Y842_U - Frontal -> Frontal','Y842_U - Parietal -> Parietal'}; + 'Y842_U - Frontal -> Frontal','Y842_U - Parietal -> Parietal'}; U_cod_lat = NaN(1,length(U_tracts)); U_cod_spd = NaN(1,length(U_tracts)); @@ -460,13 +465,13 @@ for iDir = [false true] if strcmp(out{iTr}{iSubTr}{iDir + 1}.lat_fit, 'linear') - fprintf(' - %s: best fit is linear, with COD = %2.0f, and delta %1.2f \n', ... - out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.lat_cod, out{iTr}{iSubTr}{iDir + 1}.lat_delta) + fprintf(' - %s: n=%2.0f, best fit is linear, with COD = %2.0f, and delta %1.2f \n', ... + out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.numAges, out{iTr}{iSubTr}{iDir + 1}.lat_cod, out{iTr}{iSubTr}{iDir + 1}.lat_delta) elseif strcmp(out{iTr}{iSubTr}{iDir + 1}.lat_fit, 'second') - fprintf(' - %s: best fit is second, with COD = %2.0f, and delta: age0-10 = %1.2f, age10-20 = %1.2f, age20-30 = %1.2f, age30-40 = %1.2f, age40-50 = %1.2f \n', ... - out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.lat_cod, out{iTr}{iSubTr}{iDir + 1}.lat_delta); + fprintf(' - %s: n=%2.0f, best fit is second, with COD = %2.0f, and delta: age0-10 = %1.2f, age10-20 = %1.2f, age20-30 = %1.2f, age30-40 = %1.2f, age40-50 = %1.2f \n', ... + out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.numAges, out{iTr}{iSubTr}{iDir + 1}.lat_cod, out{iTr}{iSubTr}{iDir + 1}.lat_delta); end @@ -488,13 +493,13 @@ for iDir = [false true] if strcmp(out{iTr}{iSubTr}{iDir + 1}.spd_fit, 'linear') - fprintf(' - %s: best fit is linear, with COD = %2.0f, and delta %1.2f \n', ... - out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.spd_cod, out{iTr}{iSubTr}{iDir + 1}.spd_delta) + fprintf(' - %s: n=%2.0f, best fit is linear, with COD = %2.0f, and delta %1.2f \n', ... + out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.numAges, out{iTr}{iSubTr}{iDir + 1}.spd_cod, out{iTr}{iSubTr}{iDir + 1}.spd_delta) elseif strcmp(out{iTr}{iSubTr}{iDir + 1}.spd_fit, 'second') - fprintf(' - %s: best fit is second, with COD = %2.0f, and delta: age0-10 = %1.2f, age10-20 = %1.2f, age20-30 = %1.2f, age30-40 = %1.2f, age40-50 = %1.2f \n', ... - out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.spd_cod, out{iTr}{iSubTr}{iDir + 1}.spd_delta); + fprintf(' - %s: n=%2.0f, best fit is second, with COD = %2.0f, and delta: age0-10 = %1.2f, age10-20 = %1.2f, age20-30 = %1.2f, age30-40 = %1.2f, age40-50 = %1.2f \n', ... + out{iTr}{iSubTr}{iDir + 1}.name, out{iTr}{iSubTr}{iDir + 1}.numAges, out{iTr}{iSubTr}{iDir + 1}.spd_cod, out{iTr}{iSubTr}{iDir + 1}.spd_delta); end % get averages From 245976b82710fc6158a4d4a11f5040fd970ba742 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Tue, 24 Jan 2023 17:41:09 -0600 Subject: [PATCH 07/10] Removed script irrelevant to article --- scripts/ccepExtra01_plotN1.m | 379 ----------------------------------- 1 file changed, 379 deletions(-) delete mode 100644 scripts/ccepExtra01_plotN1.m diff --git a/scripts/ccepExtra01_plotN1.m b/scripts/ccepExtra01_plotN1.m deleted file mode 100644 index 92045bc..0000000 --- a/scripts/ccepExtra01_plotN1.m +++ /dev/null @@ -1,379 +0,0 @@ -%% ccep03_plotN1 -% this code: -% 1. loads the latencies from ccep02_loadN1.m -% 2. displays the age distribution -% 3. displays a scatter plot for latencies of connections from one region -% to another for different ages. -% 4. 3 repeated for all combinations of regions. Now a first and second -% order polynomial are fitted on the data. -% 5. multilinear regression is used to predict age from latencies. In the -% figure, the first and second order polynomial are fitted and the R^2 -% are displayed for both. -% 6. displays an adjacency matrix with all destrieux regions stimulated, -% showing CCEPs in other destrieux regoins -% 7. tests leave one out cross validation of a first and second order -% polynomial and plots the COD of the first order polynomial fitting. - -% none of those figures are used in the paper. - -% Dora Hermes, Dorien van Blooijs, 2020 - - -%% load the n1Latencies from the derivatives -clear -clc -close all - -myDataPath = setLocalDataPath(1); -if exist(fullfile(myDataPath.output,'derivatives','av_ccep','n1Latencies_V1.mat'),'file') - - % if the n1Latencies_V1.mat was saved after ccep02_loadN1, load the n1Latencies structure here - load(fullfile(myDataPath.output,'derivatives','av_ccep','n1Latencies_V1.mat'),'n1Latencies') -else - disp('Run first ccep02_loadN1.m') -end - -%% age distribution --> not used in paper - -all_ages = zeros(length(n1Latencies),1); -for kk = 1:length(n1Latencies) - all_ages(kk) = n1Latencies(kk).age; -end - -figure, -histogram(all_ages,'BinMethod','integers') -title('Age distribution') -xlabel('Age (years)') -ylabel('# subjects') - -%% connections from one region to another --> not used in paper -% shows scatter plot with latencies per patient with connection from a -% specific region to another region. In the title, the pearson correlation -% for latency (ms) per age is displayed. - -region_start = input('Choose roi where connections start [temporal, frontal, parietal, occipital, central]: ','s'); -region_end = input('Choose roi where connections end [temporal, frontal, parietal, occipital, central]: ','s'); - -% region_start = 'occipital'; -% region_end = 'occipital'; - -out = ccep_connectRegions(n1Latencies,region_start,region_end); - -% initialize output: age, mean and variance in latency per subject -my_output = NaN(length(out.sub),3); - -% get variable per subject -for kk = 1:length(out.sub) - my_output(kk,1) = out.sub(kk).age; - my_output(kk,2) = mean(out.sub(kk).latencies,'omitnan'); - my_output(kk,3) = var(out.sub(kk).latencies,'omitnan'); -end - -figure -plot(my_output(:,1),1000*my_output(:,2),'k.','MarkerSize',10) -xlabel('age (years)'),ylabel('mean dT (ms)') -[r,p] = corr(my_output(~isnan(my_output(:,2)),1),my_output(~isnan(my_output(:,2)),2),'Type','Pearson'); -title(['r=' num2str(r,3) ' p=' num2str(p,3)]) - -%% scatter plot with CCEP latency from one region to another related to age --> not used in paper -% we fitted a second order polynomial and first order polynomial. The r- -% and p-value of the pearson correlation is shown in the title of each -% subplot. -% subplots with all combinations - -clear out -regions_connect = {'temporal','central','parietal','frontal'}; -conn_matrix = [1 1; 1 2; 1 3; 1 4; 2 1; 2 2; 2 3; 2 4; 3 1; 3 2; 3 3; 3 4; 4 1; 4 2; 4 3; 4 4]; -out = cell(1,size(conn_matrix,1)); - -for outInd = 1:size(conn_matrix,1) - region_start = regions_connect{conn_matrix(outInd,1)}; - region_end = regions_connect{conn_matrix(outInd,2)}; - temp = ccep_connectRegions(n1Latencies,region_start,region_end); - out{outInd} = temp; - out{outInd}.name = {region_start,region_end}; -end - - -figure('position',[0 0 700 700]) -for outInd = 1:size(conn_matrix,1) - - % initialize output: age, mean and variance in latency per subject - my_output = NaN(length(out{outInd}.sub),3); - - % get variable per subject - for kk = 1:length(out{outInd}.sub) - my_output(kk,1) = out{outInd}.sub(kk).age; - my_output(kk,2) = mean(out{outInd}.sub(kk).latencies,'omitnan'); - my_output(kk,3) = var(out{outInd}.sub(kk).latencies,'omitnan'); - end - - % age vs mean CCEP - subplot(4,4,outInd),hold on - plot(my_output(:,1),1000*my_output(:,2),'k.','MarkerSize',10) - xlabel('age (years)'),ylabel('mean dT (ms)') - [r,p] = corr(my_output(~isnan(my_output(:,2)),1),my_output(~isnan(my_output(:,2)),2),'Type','Pearson'); -% title([out(outInd).name ' to ' out(outInd).name ', r=' num2str(r,3) ' p=' num2str(p,3)]) - title(['r=' num2str(r,3) ' p=' num2str(p,3)]) - xlim([0 max(all_ages)+1])%, ylim([0 100]) - - % Yeatman et al., fit a second order polynomial: - % y = w1* age^2 * w2*age + w3 - [P] = polyfit(my_output(~isnan(my_output(:,2)),1),1000*my_output(~isnan(my_output(:,2)),2),2); - x_age = 1:1:max(all_ages); - y_fit = P(1)*x_age.^2 + P(2)*x_age + P(3); - plot(x_age,y_fit,'r') - - % Let's fit a first order polynomial: - % y = w1*age + w2 - [P,S] = polyfit(my_output(~isnan(my_output(:,2)),1),1000*my_output(~isnan(my_output(:,2)),2),1); - x_age = 1:1:max(all_ages); - y_fit = P(1)*x_age + P(2); - plot(x_age,y_fit,'r') -end - -if ~exist(fullfile(myDataPath.output,'derivatives','age'),'dir') - mkdir(fullfile(myDataPath.output,'derivatives','age')); -end - -% figureName = fullfile(myDataPath.output,'derivatives','age','AgeVsLatency_N1'); - -% set(gcf,'PaperPositionMode','auto') -% print('-dpng','-r300',figureName) -% print('-depsc','-r300',figureName) - - -%% Predict age from mean/var --> not used in paper -% apply multilinear regression - -figure('position',[0 0 700 700]) -for outInd = 1:size(conn_matrix,1) - - % initialize output: age, mean and variance in latency per subject - my_output = NaN(length(out{outInd}.sub),3); - - % get variable per subject - for kk = 1:length(out{outInd}.sub) - my_output(kk,1) = out{outInd}.sub(kk).age; - my_output(kk,2) = mean(out{outInd}.sub(kk).latencies,'omitnan'); - my_output(kk,3) = std(out{outInd}.sub(kk).latencies,'omitnan'); - end - - % age vs mean CCEP - subplot(4,4,outInd),hold on - plot(1000*my_output(:,2),my_output(:,1),'k.','MarkerSize',10) - ylabel('age (years)'),xlabel('mean dT (ms)') -% [r,p] = corr(my_output(~isnan(my_output(:,2)),2),my_output(~isnan(my_output(:,2)),1),'Type','Pearson'); - stats2 = regstats(my_output(~isnan(my_output(:,2)),1),[my_output(~isnan(my_output(:,2)),2) my_output(~isnan(my_output(:,2)),3)]); % linear regression: age vs [mean latencies, SD latencies] - stats1 = regstats(my_output(~isnan(my_output(:,2)),1),my_output(~isnan(my_output(:,2)),2)); % linear regression: age vs mean latencies -% title([out(outInd).name ' to ' out(outInd).name ', r=' num2str(r,3) ' p=' num2str(p,3)]) - title(['R_2=' num2str(stats2.adjrsquare,3) 'R_1=' num2str(stats1.rsquare,3)]) - xlim([0 max(all_ages)])%, ylim([0 100]) - - % Yeatman et al., fit a second order polynomial: - % y = w1* age^2 * w2*age + w3 - [P,~] = polyfit(1000*my_output(~isnan(my_output(:,2)),2),my_output(~isnan(my_output(:,2)),1),2); - x_latency = 20:1:max(all_ages); - y_fit = P(1)*x_latency.^2 + P(2)*x_latency + P(3); - plot(x_latency,y_fit,'r') - - % Let's fit a first order polynomial: - % y = w1*age + w2 - [P,S] = polyfit(1000*my_output(~isnan(my_output(:,2)),2),my_output(~isnan(my_output(:,2)),1),1); - x_latency = 2:1:60; - y_fit = P(1)*x_latency + P(2); - plot(x_latency,y_fit,'r') -end - -if ~exist(fullfile(myDataPath.output,'derivatives','age'),'dir') - mkdir(fullfile(myDataPath.output,'derivatives','age')); -end - -% figureName = fullfile(myDataPath.output,'derivatives','age','LatencyVSage_N1'); - -% set(gcf,'PaperPositionMode','auto') -% print('-dpng','-r300',figureName) -% print('-depsc','-r300',figureName) - - -%% overview of the number of connections from one region to another --> not used in paper -% adjacency matrix of subject 1 and all subjects together - -DestAmatall = zeros(size(n1Latencies,2),75,75); -for kk = 1:length(n1Latencies) % loop subjects - % initialize connections as empty - DestAmat = zeros(75,75); % roi_start --> roi_end - - for ll = 1:length(n1Latencies(kk).run) % loop runs - % run through all first stimulated channels (nr 1 in pair) - for chPair = 1:size(n1Latencies(kk).run(ll).n1_peak_sample,2) - - % find region label of chPair - chPairDest = str2double(n1Latencies(kk).run(ll).average_ccep_DestrieuxNr(chPair,:)); - - % find N1 responders - chSig = find(~isnan(n1Latencies(kk).run(ll).n1_peak_sample(:,chPair))==1); - % find region label of N1 responders - chSigDest = str2double(n1Latencies(kk).run(ll).channel_DestrieuxNr(chSig)); - - - DestAmat(chSigDest(~isnan(chSigDest) & chSigDest ~= 0),chPairDest(~isnan(chPairDest) & chPairDest ~= 0)) = ... - DestAmat(chSigDest(~isnan(chSigDest) & chSigDest ~= 0),chPairDest(~isnan(chPairDest) & chPairDest ~= 0))+1; - end - end - - DestAmatall(kk,:,:) = DestAmat; -end - -DestAmatsum = squeeze(sum(DestAmatall,1)); - -% visualize all existing connections per patient and in total - -subj = 1; - -figure, -subplot(1,2,1), -imagesc(squeeze(DestAmatall(subj,:,:)),[0 10]) -title(sprintf('Connections between regions in subject %1.0f',subj)) -xlabel('Responding Destrieux region') -ylabel('Stimulated Destrieux region') - -subplot(1,2,2), -imagesc(DestAmatsum,[0 50]) -title('Connections between regions in all subjects') -xlabel('Responding Destrieux region') -ylabel('Stimulated Destrieux region') - -%% Test fitting a first and second order polynomial with leave 1 out cross validation --> not used in paper -% - -nsubs = length(out{outInd}.sub); -cod_out = zeros(size(conn_matrix,1),2); -figure('position',[0 0 700 600]) -for outInd = 1:size(conn_matrix,1) - - % initialize output: age, mean and variance in latency per subject - my_output = NaN(nsubs,3); - - % get variable per subject - for kk = 1:nsubs - my_output(kk,1) = out{outInd}.sub(kk).age; - my_output(kk,2) = mean(out{outInd}.sub(kk).latencies,'omitnan'); - my_output(kk,3) = std(out{outInd}.sub(kk).latencies,'omitnan')./sqrt(length(out{outInd}.sub(kk).latencies)); - end - - % Test fitting a first order polynomial (leave 1 out cross validation) - % y = w1*age + w2 - cross_val_linear = NaN(length(find(~isnan(my_output(:,2)))),4); - % size latency (ms) X prediction (ms) X p1 (slope) X p2 (intercept) of left out - sub_counter = 0; - for kk = 1:nsubs - if ~isnan(my_output(kk,2)) - sub_counter = sub_counter+1; - % leave out kk - theseSubsTrain = ~isnan(my_output(:,2)) & ~ismember(1:nsubs,kk)'; - P = polyfit(my_output(theseSubsTrain,1),1000*my_output(theseSubsTrain,2),1); - cross_val_linear(sub_counter,3:4) = P; - cross_val_linear(sub_counter,1) = 1000*my_output(kk,2); % kk (left out) actual - cross_val_linear(sub_counter,2) = P(1)*my_output(kk,1)+P(2); % kk (left out) prediction - end - end - cod_out(outInd,1) = calccod(cross_val_linear(:,2),cross_val_linear(:,1),1); - - % Like Yeatman et al., for DTI fit a second order polynomial: - cross_val_second = NaN(length(find(~isnan(my_output(:,2)))),5); - % size latency (ms) X prediction (ms) X p1 (age^2) X p2 (age) X p3 (intercept) of left out - sub_counter = 0; - for kk = 1:nsubs - if ~isnan(my_output(kk,2)) - sub_counter = sub_counter+1; - % leave out kk - theseSubsTrain = ~isnan(my_output(:,2)) & ~ismember(1:nsubs,kk)'; - P = polyfit(my_output(theseSubsTrain,1),1000*my_output(theseSubsTrain,2),2); - cross_val_second(sub_counter,3:5) = P; - cross_val_second(sub_counter,1) = 1000*my_output(kk,2); - cross_val_second(sub_counter,2) = P(1)*my_output(kk,1).^2+P(2)*my_output(kk,1)+P(3); - end - end - cod_out(outInd,2) = calccod(cross_val_second(:,2),cross_val_second(:,1),1); - - % Fit with a piecewise linear model: - cross_val_piecewiselin = NaN(length(find(~isnan(my_output(:,2)))),5); - % size latency (ms) X prediction (ms) X p1 (age^2) X p2 (age) X p3 (intercept) of left out - sub_counter = 0; - my_options = optimoptions(@lsqnonlin,'Display','off','Algorithm','trust-region-reflective'); - for kk = 1:nsubs - if ~isnan(my_output(kk,2)) - sub_counter = sub_counter+1; - % leave out kk - theseSubsTrain = ~isnan(my_output(:,2)) & ~ismember(1:nsubs,kk)'; - -% % use the slmengine tool from here: -% % John D'Errico (2020). SLM - Shape Language Modeling (https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling), MATLAB Central File Exchange. Retrieved July 14, 2020. -% slm = slmengine(my_output(theseSubsTrain,1),1000*my_output(theseSubsTrain,2),'degree',1,'plot','off','knots',3,'interiorknots','free'); -% % predicted value at left out x -% yhat = slmeval(my_output(kk,1),slm,0); -% cross_val_piecewiselin(sub_counter,1) = 1000*my_output(kk,2); -% cross_val_piecewiselin(sub_counter,2) = yhat; - - % use our own function: - x = my_output(theseSubsTrain,1); - y = 1000*my_output(theseSubsTrain,2); - [pp] = lsqnonlin(@(pp) ccep_fitpiecewiselinear(pp,y,x),... - [40 -1 0 20],[0 -Inf -Inf 10],[40 0 Inf 30],my_options); - x_fit = my_output(kk,1); - y_fit = (pp(1) + pp(2)*min(pp(4),x_fit) + pp(3)*max(pp(4),x_fit)); - cross_val_piecewiselin(sub_counter,1) = 1000*my_output(kk,2); - cross_val_piecewiselin(sub_counter,2) = y_fit; - end - end - cod_out(outInd,3) = calccod(cross_val_piecewiselin(:,2),cross_val_piecewiselin(:,1),1); - - subplot(4,4,outInd),hold on - - % figure,hold on - for kk = 1:nsubs - % plot histogram per subject in the background - if ~isnan(my_output(kk,2)) - distributionPlot(1000*out{outInd}.sub(kk).latencies','xValues',out{outInd}.sub(kk).age,... - 'color',[.6 .6 .6],'showMM',0,'histOpt',2) - end -% % plot mean+sterr per subject -% plot([my_output(kk,1) my_output(kk,1)],[1000*(my_output(kk,2)-my_output(kk,3)) 1000*(my_output(kk,2)+my_output(kk,3))],... -% 'k','LineWidth',1) - end - % plot mean per subject in a dot - plot(my_output(:,1),1000*my_output(:,2),'ko','MarkerSize',6) - - [r,p] = corr(my_output(~isnan(my_output(:,2)),1),my_output(~isnan(my_output(:,2)),2),'Type','Pearson'); -% title([out(outInd).name ' to ' out(outInd).name ', r=' num2str(r,3) ' p=' num2str(p,3)]) - - % plot confidence intervals for linear fit - x_age = 1:1:max(all_ages); - % get my crossval y distribution - y_n1LatCross = cross_val_linear(:,3)*x_age + cross_val_linear(:,4); - % y_n1LatCross = cross_val_second(:,3)*x_age.^2 + cross_val_second(:,4)*x_age + cross_val_second(:,5); - - % get 95% confidence intervals - low_ci = quantile(y_n1LatCross,.025,1); - up_ci = quantile(y_n1LatCross,.975,1); - fill([x_age x_age(end:-1:1)],[low_ci up_ci(end:-1:1)],[0 .7 1],'EdgeColor',[0 .7 1]) - - % put COD in title - title(['COD=' int2str(cod_out(outInd,1)) ' p=' num2str(p,2)]) - - xlim([0 60]), ylim([0 100]) - -% xlabel('age (years)'),ylabel('mean dT (ms)') - set(gca,'XTick',10:10:50,'YTick',20:20:100,'FontName','Arial','FontSize',12) - -end - -if ~exist(fullfile(myDataPath.output,'derivatives','age'),'dir') - mkdir(fullfile(myDataPath.output,'derivatives','age')); -end -figureName = fullfile(myDataPath.output,'derivatives','age','AgeVsLatency_N1'); - -% set(gcf,'PaperPositionMode','auto') -% print('-dpng','-r300',figureName) -% print('-depsc','-r300',figureName) From 8840147c1adadc04c67cd3497a623c81d66e787d Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Fri, 27 Jan 2023 10:49:37 -0600 Subject: [PATCH 08/10] Updated correlation symbil --- scripts/makeSupFig1_only8maSubs.m | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/scripts/makeSupFig1_only8maSubs.m b/scripts/makeSupFig1_only8maSubs.m index bd7f04a..84d59b4 100644 --- a/scripts/makeSupFig1_only8maSubs.m +++ b/scripts/makeSupFig1_only8maSubs.m @@ -1,7 +1,8 @@ % -% This script outputs supplementary figure 1 - relation between age (years) and latency (ms) at different stimulation currents +% This script produces supplementary figure 1 - relation between age (years) and latency (ms) at different stimulation currents % -% Dorien van Blooijs, Max van den Boom, 2022 +% Max van den Boom, Dorien van Blooijs, 2023 + clear close all @@ -24,7 +25,7 @@ %% -% ... +% Extract the different currents from the dataset all_unique_currents =[]; all_currents = struct(); @@ -104,7 +105,7 @@ all_currents.(['mA', num2str(subject_unique_currents(iCurrent))]) = []; end - % for this sbuject and current, add the age and latency + % for this subject and current, add the age and latency ind = length(all_currents.(['mA', num2str(subject_unique_currents(iCurrent))])); all_currents.(['mA', num2str(subject_unique_currents(iCurrent))])(ind + 1).age = ccepData(iSubj).age; all_currents.(['mA', num2str(subject_unique_currents(iCurrent))])(ind + 1).latency = subject_currents.(['mA', num2str(subject_unique_currents(iCurrent))]); @@ -178,15 +179,15 @@ % plot the points plot(ages, latencies * 1000, '.', 'Color', curColor, 'MarkerSize', 16, 'DisplayName', [num2str(all_unique_currents(iCurrent)) ' mA (n=', num2str(length(latencies)), ')']); + % check if there are enough points if length(latencies) >= minDataPoints - [P, S] = polyfit(ages(~isnan(ages) & ~isnan(latencies)), ... - latencies(~isnan(ages) & ~isnan(latencies)) * 1000, 1); + [P, S] = polyfit(ages(~isnan(ages) & ~isnan(latencies)), latencies(~isnan(ages) & ~isnan(latencies)) * 1000, 1); [y_fit, ~] = polyval(P, ages, S); - % Plot polyfit throught data points + % Plot fit throught data points plot(ages, y_fit, 'Color', curColor, 'LineWidth', 2, 'DisplayName', ... - ['r=' num2str(all_unique_currents_r(iCurrent), 3) ' p=' num2str(all_unique_currents_p_fdr(iCurrent), 3)]); + ['\rho=' num2str(all_unique_currents_r(iCurrent), 3) ' p=' num2str(all_unique_currents_p_fdr(iCurrent), 3)]); end end From e3cc79bd361c9e6c6ff02c66b23a715b844bafa0 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Fri, 27 Jan 2023 10:50:25 -0600 Subject: [PATCH 09/10] Make boxplot more visible in violins --- scripts/makeSupFig6and7_exploreSOZ.m | 55 ++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 15 deletions(-) diff --git a/scripts/makeSupFig6and7_exploreSOZ.m b/scripts/makeSupFig6and7_exploreSOZ.m index f36eb8c..06212dc 100644 --- a/scripts/makeSupFig6and7_exploreSOZ.m +++ b/scripts/makeSupFig6and7_exploreSOZ.m @@ -259,16 +259,24 @@ % Generate supplement 6 - Violin plot responses in or outside SOZ respall = cell(1); +SEs = cell(1); names = cell(1); count = 1; for iSubj = 1:size(n1Latencies, 2) + % check whether there are N1 latencies in electrodes that were on SOZs and electrodes that were not on SOZs if ~isempty(n1Latencies(iSubj).respSOZlatencies) && ~isempty(n1Latencies(iSubj).respNonSOZlatencies) + + % add N1 latencies for electrodes that were on SOZs to be plotted respall{count} = n1Latencies(iSubj).respSOZlatencies * 1000; + SEs{count} = std(respall{count}) / sqrt(length(respall{count})); names{count} = cellstr(repmat([n1Latencies(iSubj).id ' soz'],size(respall{count}, 2), 1)); count = count + 1; + + % add N1 latencies for electrodes that were on SOZs to be plotted respall{count} = n1Latencies(iSubj).respNonSOZlatencies * 1000; + SEs{count} = std(respall{count}) / sqrt(length(respall{count})); names{count} = cellstr(repmat([n1Latencies(iSubj).id ' nsoz'],size(respall{count}, 2), 1)); count = count + 1; @@ -282,23 +290,30 @@ end +% open the figure and draw the violins close all cmap = colormap('parula'); +close(figure(1)) ymax = ceil(max(horzcat(respall{:}))); temp = horzcat(respall{:}); ymin = floor(min(temp(temp > 0))); -h = figure('position',[0 0 2000 1200]); -vs = violinplot(horzcat(respall{:}), vertcat(names{:}), 'ViolinColor', cmap(1, :), 'Width', 0.3); -for iSubj = 1:3:size(vs, 2) - vs(iSubj).ViolinPlot.FaceColor = cmap(128, :); - vs(iSubj).ScatterPlot.MarkerFaceColor = cmap(128, :); - vs(iSubj).ScatterPlot.SizeData = 70; - vs(iSubj).MedianPlot.SizeData = 70; +h = figure('position',[0 0 2400 1200]); +vs = violinplot(horzcat(respall{:}), vertcat(names{:}), 'ViolinColor', cmap(1, :), 'BoxColor', [0.4 0.4 0.6], 'Width', 0.4, 'BoxWidth', .08, 'ShowData', false); + +% for each subject/violin plot pair (non-SOZ violin & SOZ violin) +for iViolin = 1:3:size(vs, 2) + vs(iViolin).ViolinPlot.FaceColor = cmap(128, :); + vs(iViolin).ScatterPlot.MarkerFaceColor = cmap(128, :); + vs(iViolin).ScatterPlot.SizeData = 70; + vs(iViolin).MedianPlot.SizeData = 70; + vs(iViolin).ViolinPlot.FaceColor = cmap(128, :); end hold on count = 1; for iSubj = 1:size(n1Latencies, 2) + + % check whether there are N1 latencies in electrodes that were on SOZs and electrodes that were not on SOZs if ~isempty(n1Latencies(iSubj).respSOZ_pFDR) && ~isnan(n1Latencies(iSubj).respSOZ_pFDR) if n1Latencies(iSubj).respSOZ_pFDR < .05 text(count + 0.3, 1.03 * ymax, '*', 'FontSize', 24) @@ -306,6 +321,7 @@ end count = count + 3; end + end hold off @@ -316,7 +332,7 @@ ax.XTickLabelRotation = 0; ax.XTick = 1.5:3:size(vs, 2); ax.XTickLabel = 1:size([n1Latencies(:).respSOZ_p], 2); -ax.YLabel.String = 'Latency (ms)'; +ax.YLabel.String = 'N1 Peak Latency (ms)'; ax.XLabel.String = 'Subjects'; ax.Title.String = 'Latency of measured electrodes in or outside SOZ'; lgd = legend([vs(1).ViolinPlot, vs(2).ViolinPlot], 'Responses measured on non-SOZ', 'Responses measured on SOZ'); @@ -340,8 +356,9 @@ count = 1; for iSubj = 1:size(n1Latencies, 2) - % check if there are + % check whether there are N1 latencies in electrodes that were on SOZs and electrodes that were not on SOZs if ~isempty(n1Latencies(iSubj).stimSOZlatencies) && ~isempty(n1Latencies(iSubj).stimNonSOZlatencies) + stimall{count} = n1Latencies(iSubj).stimSOZlatencies * 1000; names{count} = cellstr(repmat([n1Latencies(iSubj).id ' soz'],size(stimall{count}, 2), 1)); count = count + 1; @@ -359,22 +376,29 @@ end +% open the figure and draw the violins close all cmap = colormap('parula'); +close(figure(1)) ymax = ceil(max(horzcat(stimall{:}))); temp = horzcat(stimall{:}); ymin = floor(min(temp(temp > 0))); h = figure('position',[0 0 2000 1200]); -vs = violinplot(horzcat(stimall{:}), vertcat(names{:}), 'ViolinColor', cmap(1, :), 'Width', 0.3); -for iSubj = 1:3:size(vs, 2) - vs(iSubj).ViolinPlot.FaceColor = cmap(128, :); - vs(iSubj).ScatterPlot.MarkerFaceColor = cmap(128, :); - vs(iSubj).ScatterPlot.SizeData = 70; +vs = violinplot(horzcat(stimall{:}), vertcat(names{:}), 'ViolinColor', cmap(1, :), 'BoxColor', [0.4 0.4 0.6], 'Width', 0.4, 'BoxWidth', .08, 'ShowData', false); + +% for each subject/violin plot pair (non-SOZ violin & SOZ violin) +for iViolin = 1:3:size(vs, 2) + vs(iViolin).ViolinPlot.FaceColor = cmap(128, :); + vs(iViolin).ScatterPlot.MarkerFaceColor = cmap(128, :); + vs(iViolin).ScatterPlot.SizeData = 70; + end hold on count = 1; for iSubj = 1:size(n1Latencies, 2) + + % check whether there are N1 latencies in electrodes that were on SOZs and electrodes that were not on SOZs if ~isempty(n1Latencies(iSubj).stimSOZ_pFDR) && ~isnan(n1Latencies(iSubj).stimSOZ_pFDR) if n1Latencies(iSubj).stimSOZ_pFDR < .05 text(count + 0.3, 1.03 * ymax, '*', 'FontSize', 24) @@ -382,6 +406,7 @@ end count = count + 3; end + end hold off @@ -392,7 +417,7 @@ ax.XTickLabelRotation = 0; ax.XTick = 1.5:3:size(vs, 2); ax.XTickLabel = 1:size([n1Latencies(:).stimSOZ_p], 2); -ax.YLabel.String = 'Latency (ms)'; +ax.YLabel.String = 'N1 Peak Latency (ms)'; ax.XLabel.String = 'Subjects'; ax.Title.String = 'Latency of stim-pairs when stimulating in or outside SOZ'; lgd = legend([vs(1).ViolinPlot, vs(2).ViolinPlot], 'Responses when not stimulating SOZ', 'Responses when stimulating SOZ'); From 127d726dfdf1202cc6e124d3d0f03eb06565da71 Mon Sep 17 00:00:00 2001 From: MaxvandenBoom Date: Fri, 27 Jan 2023 10:51:34 -0600 Subject: [PATCH 10/10] Small improvement to make whiskers more visible --- external/Violinplot-Matlab/Violin.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/Violinplot-Matlab/Violin.m b/external/Violinplot-Matlab/Violin.m index e1c74ca..560f78b 100644 --- a/external/Violinplot-Matlab/Violin.m +++ b/external/Violinplot-Matlab/Violin.m @@ -173,7 +173,7 @@ hiwhisker = quartiles(3) + 1.5*IQR; hiwhisker = min(hiwhisker, max(data(data < hiwhisker))); if ~isempty(lowhisker) && ~isempty(hiwhisker) - obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker]); + obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker], 'LineWidth', 1.2); end obj.MedianPlot = scatter(pos, quartiles(2), [], [1 1 1], 'filled');