Merge pull request #47 from dorahermes/master

update to readme and correction to area 46
MultimodalNeuroimagingLab · Mar 22, 2022 · d1fe4e9 · d1fe4e9
2 parents f457172 + 48d1c71
commit d1fe4e9
Show file tree

Hide file tree

Showing 10 changed files with 214 additions and 327 deletions.
diff --git a/README.md b/README.md
@@ -1,6 +1,31 @@
 # mnl_ccepBids
-Repository to publicly share code to analyze/visualize CCEP data across development in BIDS
 
+This repository contains the scripts and functions necessary to reproduce the analyses and figures in this manuscript:
+
+- The development of efficient communication in the human connectome. D. van Blooijs, J.F. van der Aar, G.J.M. Huiskamp, G. Castegnaro, M. Demuru, W.J.E.M. Zweiphenning, P. van Eijsden, K. J. Miller, F.S.S. Leijten, D. Hermes bioRxiv 2022.03.14.484297; doi: https://doi.org/10.1101/2022.03.14.484297
+
+
+# Generating the figures
+Scripts to process the data and detect N1 responses:
+scripts/ccep01_averageCCEPs.m
+scripts/ccep02_loadN1.m
+
+Scripts to make the figure panels:
+makeFig1A_plotMNI.m
+makeFig1B_ExampleResponse.m
+makeFig2_plotResponsesAge.m
+makeFig3and4_plotN1_meanacrossage.m
+
+
+To make all functions work, an m-file called personalDataPath.m should be stored in the root dir. This file should have the following content:
+```
+function localDataPath = personalDataPath()
+% function that contains local data path, is ignored in .gitignore
+localDataPath.input = '/my/path/to/load/data/';
+localDataPath.output = '/my/path/to/save/data/';
+addpath('/my/path/to/fieldtrip')
+ft_defaults
+```
 
 # Dependencies
   Fieldtrip: http://www.fieldtriptoolbox.org/

diff --git a/functions/ccep_categorizeAnatomicalRegions.m b/functions/ccep_categorizeAnatomicalRegions.m
@@ -9,9 +9,9 @@
 
 % central areas
 % G_postcentral G_precentral S_central
-roi{2} = {'28','29','46'};
+roi{2} = {'28','29','45'};
 roi_name{2} = 'central';
-roi_central = [28 29 46];
+roi_central = [28 29 45];
 
 % parietal areas:
 % G_pariet_inf-Angular, G_pariet_inf-Supramar, G_parietal_sup
@@ -21,7 +21,7 @@
 
 % frontal areas:
 % G_front_inf-Triangul, G_front_middle, G_front_inf-Opercular
-roi{4} = {'14','15','12'}; % maybe add 16: G_front_sup
+roi{4} = {'14','15','12'}; 
 roi_name{4} = 'frontal';
 roi_frontal = [14 15 12]; 
 

diff --git a/scripts/makeFig1A_plotMNI.m b/scripts/makeFig1A_plotMNI.m
@@ -209,10 +209,10 @@
 ieeg_elAdd(els(ismember(all_hemi,'L') & ismember(allmni_labels,roi_parietal),:),[0 .5 0],15)
 ieeg_viewLight(v_d(1),v_d(2))
 
-% figureName = fullfile(myDataPath.output,'derivatives','render','leftMNIpial');
+figureName = fullfile(myDataPath.output,'derivatives','render','leftMNIpial');
 
-% set(gcf,'PaperPositionMode','auto')
-% print('-dpng','-r300',figureName)
+set(gcf,'PaperPositionMode','auto')
+print('-dpng','-r300',figureName)
 
 
 %% Plot figure with right pial with electrodes in mni space
@@ -222,7 +222,7 @@
 gr.faces = Rmnipial_face+1;
 gr.vertices = Rmnipial_vert;
 gr = gifti(gr);
-tH = ieeg_RenderGifti(gr); %#ok<NASGU>
+tH = ieeg_RenderGifti(gr); 
 
 % make sure electrodes pop out
 a_offset = .5*max(abs(allmni_coords(:,1)))*[cosd(v_d(1)-90)*cosd(v_d(2)) sind(v_d(1)-90)*cosd(v_d(2)) sind(v_d(2))];
@@ -236,10 +236,10 @@
 ieeg_elAdd(els(ismember(all_hemi,'R') & ismember(allmni_labels,roi_parietal),:),[0 .5 0],15)
 ieeg_viewLight(v_d(1),v_d(2))
 
-figureName = fullfile(myDataPath.output,'derivatives','render','rightMNIpial'); %#ok<NASGU>
+figureName = fullfile(myDataPath.output,'derivatives','render','rightMNIpial'); 
 
 set(gcf,'PaperPositionMode','auto')
-% print('-dpng','-r300',figureName)
+print('-dpng','-r300',figureName)
 
 
 %% Plot left inflated brain surface with electrodes in mni space

diff --git a/scripts/makeFig2_plotResponsesAge.m b/scripts/makeFig2_plotResponsesAge.m
@@ -21,8 +21,9 @@
     error('Answer to previous question is not recognized.')
 end
 
+myDataPath = setLocalDataPath(1);
+
 if select_amplitude==0 
-    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')
@@ -32,9 +33,8 @@
         disp('Run first script ccep02_loadN1.m')
     end
 elseif select_amplitude==8 % only 8 mA
-    myDataPath = setLocalDataPath(1);
     if exist(fullfile(myDataPath.output,'derivatives','av_ccep','n1Latencies_8ma.mat'),'file')
-        % if the n1Latencies_V1.mat was saved after ccep02_loadN1, load the n1Latencies structure here
+        % ?if the n1Latencies_8ma.mat was saved after ccep02_loadN1?, load the n1Latencies structure here
         load(fullfile(myDataPath.output,'derivatives','av_ccep','n1Latencies_8ma.mat'),'n1Latencies8ma')
 
         n1Latencies = n1Latencies8ma;
@@ -254,17 +254,17 @@
     end
 end
 
-% if select_amplitude==0
-%     figureName = fullfile(myDataPath.output,'derivatives','age',...
-%         ['AllSortAge_tmax' int2str(ttmax*1000)]);
-% elseif select_amplitude==8
-%     figureName = fullfile(myDataPath.output,'derivatives','age',...
-%         ['AllSortAge_tmax' int2str(ttmax*1000), '_8mA']);
-% end
+if select_amplitude==0
+    figureName = fullfile(myDataPath.output,'derivatives','age',...
+        ['AllSortAge_tmax' int2str(ttmax*1000)]);
+elseif select_amplitude==8
+    figureName = fullfile(myDataPath.output,'derivatives','age',...
+        ['AllSortAge_tmax' int2str(ttmax*1000), '_8mA']);
+end
 
-% set(gcf,'PaperPositionMode','auto')
-% print('-dpng','-r300',figureName)
-% print('-depsc','-r300',figureName)
+set(gcf,'PaperPositionMode','auto')
+print('-dpng','-r300',figureName)
+print('-depsc','-r300',figureName)
 
 
 %% figure with colormap
@@ -275,6 +275,7 @@
 axis off
 figureName = fullfile(myDataPath.output,'derivatives','age',...
     ['AllSortAge_tmax' int2str(ttmax*1000) '_cm']);
+
 % set(gcf,'PaperPositionMode','auto')
 % print('-dpng','-r300',figureName)
 % print('-depsc','-r300',figureName)

diff --git a/scripts/makeFig3_plotN1_meanacrossage.m → scripts/makeFig3and4_plotN1_meanacrossage.m b/scripts/makeFig3_plotN1_meanacrossage.m → scripts/makeFig3and4_plotN1_meanacrossage.m
@@ -18,18 +18,16 @@
 else
     error('Answer to previous question is not recognized.')
 end
+myDataPath = setLocalDataPath(1);
 
 if select_amplitude==0 
-    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.mat')
     end
 elseif select_amplitude==8 % only 8 mA
-    myDataPath = setLocalDataPath(1);
-
     if exist(fullfile(myDataPath.output,'derivatives','av_ccep','n1Latencies_8ma.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_8ma.mat'),'n1Latencies8ma')
@@ -261,18 +259,16 @@
 end
 
 if select_amplitude == 0
-
     figureName = fullfile(myDataPath.output,'derivatives','age',...
         'AgeVsLatency_N1_meanacrossage');
-
 elseif select_amplitude == 8
     figureName = fullfile(myDataPath.output,'derivatives','age',...
         'AgeVsLatency_N1_meanacrossage_8mA');
-
 end
-% set(gcf,'PaperPositionMode','auto')
-% print('-dpng','-r300',figureName)
-% print('-depsc','-r300',figureName)
+
+set(gcf,'PaperPositionMode','auto')
+print('-dpng','-r300',figureName)
+print('-depsc','-r300',figureName)
 
 
 %% display in command window the cod and delta for each subplot
@@ -292,7 +288,7 @@
 end
 
 
-%% find average velocities
+%% find average latencies
 clc
 
 delta_all = [];
@@ -356,169 +352,3 @@
 betterlinearfit(betterlinearfit<=0) = NaN;
 mean(betterlinearfit,'omitnan')
 
-%% mean per age instead of first mean latency per subject --> this is not up to date anymore
-
-
-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
-    min_age = min([out{outInd}.sub(:).age]);
-    max_age = max([out{outInd}.sub(:).age]);
-    my_output = NaN(max_age-min_age,4);
-
-    % get latencies per age
-    counter = 1;
-    for kk = min_age:max_age
-        subInd = find([out{outInd}.sub(:).age] == kk);
-
-        my_output(counter,1) = kk; % age
-        this_age_latencies = horzcat(out{outInd}.sub(subInd).latencies);
-        my_output(counter,2) = mean(this_age_latencies,'omitnan');
-        this_nr_latencies = length(this_age_latencies(~isnan(this_age_latencies)));
-        my_output(counter,3) = std(this_age_latencies,'omitnan')./sqrt(this_nr_latencies);
-        my_output(counter,4) = this_nr_latencies;
-        clear this_nr_latencies this_age_latencies
-        counter = counter+1;
-    end
-
-    % remove ages with no subjects (NaN)
-    x_vals = my_output(~isnan(my_output(:,2)),1); % age for ~isnan
-    y_vals = 1000* my_output(~isnan(my_output(:,2)),2);
-
-    %     % average per 2 years of age
-    %     x_vals_temp = 2*floor(.5*x_vals);
-    %     x_vals_new = unique(x_vals_temp);
-    %     y_vals_new = zeros(size(x_vals_new));
-    %     for kk = 1:length(x_vals_new)
-    %         y_vals_new(kk) = mean(y_vals(ismember(x_vals_temp,x_vals_new(kk))));
-    %     end
-    %     x_vals = x_vals_new; y_vals = y_vals_new;
-
-    % Test fitting a first order polynomial (leave 1 out cross validation)
-    % y  =  w1*age + w2
-    cross_val_linear = NaN(length(y_vals),4);
-    % size latency (ms) X prediction (ms) X p1 (slope) X p2 (intercept) of left out
-    sub_counter = 0;
-    for kk = 1:length(y_vals)
-        sub_counter = sub_counter+1;
-        % leave out kk
-        theseSubsTrain = ~ismember(1:length(y_vals),kk)';
-        P = polyfit(x_vals(theseSubsTrain),y_vals(theseSubsTrain),1);
-        cross_val_linear(sub_counter,3:4) = P;
-        cross_val_linear(sub_counter,1) = y_vals(kk); % kk (left out) actual
-        cross_val_linear(sub_counter,2) = P(1)*x_vals(kk)+P(2); % kk (left out) prediction
-    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(y_vals),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:length(y_vals)
-        sub_counter = sub_counter+1;
-        % leave out kk
-        theseSubsTrain = ~ismember(1:length(y_vals),kk)';
-        P = polyfit(x_vals(theseSubsTrain),y_vals(theseSubsTrain),2);
-        cross_val_second(sub_counter,3:5) = P;
-        cross_val_second(sub_counter,1) = y_vals(kk);
-        cross_val_second(sub_counter,2) = P(1)*x_vals(kk).^2+P(2)*x_vals(kk)+P(3);
-    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(y_vals),6);
-    % 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:length(y_vals)
-        sub_counter = sub_counter+1;
-        % leave out kk
-        theseSubsTrain = ~ismember(1:length(y_vals),kk)';
-
-        % use our own function:
-        x = x_vals(theseSubsTrain);
-        y = y_vals(theseSubsTrain);
-        [pp] = lsqnonlin(@(pp) ccep_fitpiecewiselinear(pp,y,x),...
-            [60 -1 1 20],[20 -Inf -Inf 10],[40 0 Inf 30],my_options);
-
-        x_fit = x_vals(kk);
-        y_fit = (pp(1) + pp(2)*min(pp(4),x_fit) + pp(3)*max(x_fit-pp(4),0));
-        % --> intercept = pp(1)
-        % --> tipping point = pp(4)
-        % --> slope before tipping point = pp(2)
-        % --> slope after tipping point = pp(3)
-
-        cross_val_piecewiselin(sub_counter,1) = y_vals(kk);
-        cross_val_piecewiselin(sub_counter,2) = y_fit;
-        cross_val_piecewiselin(sub_counter,3:6) = pp;
-    end
-    cod_out(outInd,3) = calccod(cross_val_piecewiselin(:,2),cross_val_piecewiselin(:,1),1);
-    cod_out(outInd,4) = length(y_vals);
-
-    subplot(4,4,outInd),hold on
-
-    %     % add this if you want to see every single subject and the effect of
-    %     % averaging within an age group
-    %     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)])
-
-    if cod_out(outInd,1)>cod_out(outInd,3)
-        % PLOT LINEAR FIT WITH CI
-        x_age = 1:1:51;
-        % get my crossval y distribution
-        y_n1LatCross = cross_val_linear(:,3)*x_age + cross_val_linear(:,4);
-        % PLOT SECOND ORDER POLYNOMIAL FIT WITH CI
-        % 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))])
-
-    else
-        % PLOT PIECEWISE LINEAR FIT (from our own function)
-        x_age = 1:1:51;
-        y_n1LatCross = NaN(size(cross_val_piecewiselin,1),size(x_age,2));
-        for rr = 1:size(cross_val_piecewiselin,1)% run across all crossvalidations to get confidence intervals
-            y_n1LatCross(rr,:) = (cross_val_piecewiselin(rr,3) + cross_val_piecewiselin(rr,4)*min(cross_val_piecewiselin(rr,6),x_age) + ...
-                cross_val_piecewiselin(rr,5)*max(x_age-cross_val_piecewiselin(rr,6),0));
-        end
-        % 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)],[1 .7 0],'EdgeColor',[1 .7 0])
-
-        % put COD in title
-        title(['COD=' int2str(cod_out(outInd,3))])
-
-    end
-
-    plot(x_vals,y_vals,'k.','MarkerSize',6)
-
-    xlim([0 60]), ylim([0 80])
-
-    %     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_meanacrossage');
diff --git a/scripts/makeSubFig3_only8maSubs.m → scripts/makeSubFig1_only8maSubs.m b/scripts/makeSubFig3_only8maSubs.m → scripts/makeSubFig1_only8maSubs.m
diff --git a/scripts/makeSubFig4_nrResponsesAge.m → scripts/makeSubFig2_nrResponsesAge.m b/scripts/makeSubFig4_nrResponsesAge.m → scripts/makeSubFig2_nrResponsesAge.m