Skip to content

Commit

Permalink
Edits made during revisions for JNM
Browse files Browse the repository at this point in the history 
Renamed:
- figS1_compareCARvCARLA > figS3_compareCARvCARLA
- figS3_simCCEPComponentsGlobalSig.m > figS5_simCCEPComponentsGlobalSig.m

Added:
- figS1a_*, figS1b_*, figS2a_*, figS2b_* code to generate figures for the new supplementary figures S1 and S2.
- CAR_subset.m, which is a copy of ccep_CARVariance.m from mnl_ieegBasics. This performs CARLA-like re-referencing but does not do optimization. Requires the user to define a percentile threshold.

Modifications:
- CARLA: Minor bugfix to require minimum CAR size to be 2 when there are very few channels (i.e., 10)
- main.m: Added in sections pertaining to new supplementary figures S1 and S2
- README: updated doi to point to arXiv preprint.
  • Loading branch information
@hharveygit
hharveygit committed Feb 12, 2024
1 parent 414f67e commit a630235
Show file tree
Hide file tree
Showing 12 changed files with 1,157 additions and 8 deletions.
2 changes: 1 addition & 1 deletion README.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
If this code is used in a publication, please cite the manuscript:
"CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, GA Worrell, KJ Miller, and D Hermes.
A preprint is available currently at doi: ___
A preprint is available currently at https://doi.org/10.48550/arXiv.2310.00185.

Correspondence:
H Huang: [email protected]; D Hermes: [email protected]
Expand Down
171 changes: 171 additions & 0 deletions figS1a_simCCEPTotalChs.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
%% This script is derived from simCCEPLoop, except an outer loop varies the number of total channels
% We simulate 30 reps for each combination of total channels and responsiveness and test CARLA's accuracy
%
% 2024/02/12
%
% If this code is used in a publication, please cite the manuscript:
% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
% GA Worrell, KJ Miller, and D Hermes.
%
% CARLA manuscript package.
% Copyright (C) 2023 Harvey Huang
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
%% Control parameters of simulated data

srate = 4800;
tt = -0.5:1/srate:1-1/srate;

ntrs = 12; % number of trials
reps = 30; % number of reps for each condition

cmSens = [1, 165/255, 0]; % use orange as color for the sensitive optimum

Aglobal = 0;

% Total number of channels to simulate for each condition
nChsToTest = [25, 20, 15, 10];


%% Outer loop through number of total channels

rng(nChsToTest(1)); % Set seed depending on what number of channels we want, for reproducibility

for nn = 1:length(nChsToTest)

nchs = nChsToTest(nn);
fprintf('Number of channels = %d\n', nchs);

%% Inner loop through number of responsive channels
for nresp = 0:0.2*nchs:0.8*nchs % 0 to 80% responsiveness, by 20% intervals

%% Create data

fprintf('%d of %d channels responsive\n', nresp, nchs);

outdir = fullfile('output', 'simLoopNChs', sprintf('nchs%d-%d', nchs, nresp));
mkdir(outdir);

chsResp = 1:nresp; % first x channels responsive, for ease. order doesn't matter. will sort and color them red

%%
for rr = 1:reps % multiple repetitions at each significance
fprintf('.');

% i) artifact only
V0 = zeros(length(tt), nchs);
Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel
artifact = sin(2*pi*600*tt)';
artifact(tt < 0 | tt > 0.002) = 0;
V0 = V0 + artifact*Aart';

% A) Add the evoked potentials
A = 100;
V1 = V0;
sig = genRandSig(tt, length(chsResp), A);
V1(:, chsResp) = V0(:, chsResp) + sig;

% B) Option to add a global noise
if Aglobal
sigCommon = genRandSig(tt, 1, Aglobal);
else
sigCommon = 0;
end
V2 = V1 + sigCommon;

% C) Add common noise to all channels at each trial
V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial
phLN = rand(ntrs, 3)*2*pi; % LN phases
LN = zeros(length(tt), ntrs);
for ii = 1:ntrs
LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
end

% brown noise shared across channels
BN = cumsum(0.4*randn(2*length(tt), ntrs)); % variable brown noise coefficient now
BN = ieeg_highpass(BN, srate, true);
BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);

noiseCommon = LN + BN;
V3 = V3 + shiftdim(noiseCommon, -1);

% D) add random brown noise
noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
for ii = 1:nchs
noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
end
noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
V4 = V3 + noiseRand;

%% Apply CARLA and plot outputs

[Vout, CAR, stats] = CARLA(tt, V4, srate, true); % get the sensitive output

% number of channels used for the CAR
nCAR = length(stats.chsUsed);
[~, nCARglob] = max(mean(stats.zMinMean, 2)); % number of channels at global maximum

% Plot average zmin across trials
figure('Position', [200, 200, 400, 300], 'Visible', 'off'); hold on
errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k-o');
plot(nCARglob, mean(stats.zMinMean(nCARglob, :), 2), 'b*'); % global max as blue
if nCARglob ~= nCAR; plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); end
yline(0, 'Color', 'k');
saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'png');
saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'svg');

% Sort and plot channels by increasing covariance, draw line at cutoff
V4MeanSorted = mean(V4(stats.order, :, :), 3);
respBool = antifind(chsResp, nchs);
respBool = respBool(stats.order); % logical array of where responsive channels are
cm = zeros(nchs, 3);
cm(respBool, 1) = 1; % make red
figure('Position', [200, 200, 250, 600], 'Visible', 'off');
yspace = 80;
ys = ieeg_plotTrials(tt, V4MeanSorted', yspace, [], cm, 'LineWidth', 1);
yline(ys(nCARglob)-yspace/2, 'Color', 'b', 'LineWidth', 1.5);
if nCARglob ~= nCAR; yline(ys(nCAR)-yspace/2, 'Color', cmSens, 'LineWidth', 1.5); end
xlim([-0.1, 0.5]); set(gca, 'xtick', [0, 0.5]);
xlabel('Time (s)'); ylabel('Channels');
saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'png');
saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'svg');

close all;

% Accuracy values. positive means responsive/excluded from CAR
% We keep these variables as named here, but note that FN and FP are now renamed RCM and NCM in the manuscript.
TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff)
TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff
FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most
FP = sum(find(~respBool) > nCAR); % NR channels incorrectly excluded from CAR

% same for the global threshold
TPglob = sum(find(respBool) > nCARglob);
TNglob = sum(find(~respBool) <= nCARglob);
FNglob = sum(find(respBool) <= nCARglob);
FPglob = sum(find(~respBool) > nCARglob);

fid = fopen(fullfile(outdir, sprintf('accuracy_rep%d.txt', rr)), 'w');
fprintf(fid, 'TP\t%d\nTN\t%d\nFN\t%d\nFP\t%d\n', TP, TN, FN, FP);
fprintf(fid, 'TPglob\t%d\nTNglob\t%d\nFNglob\t%d\nFPglob\t%d', TPglob, TNglob, FNglob, FPglob);
fclose(fid);

end

fprintf('\n');

end
end
137 changes: 137 additions & 0 deletions figS1b_summarizeSimCCEPTotalChs.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
%% This script takes outputs from figS1a_simCCEPTotalChs, and summarizes the accuracy for each total channel size and level of responsiveness
% Generates outputs for Figure S1A
%
% 2024/02/12
%
% If this code is used in a publication, please cite the manuscript:
% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
% GA Worrell, KJ Miller, and D Hermes.
%
% CARLA manuscript package.
% Copyright (C) 2023 Harvey Huang
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
%% Configure parameters and directories

parCm = brighten(parula(15), -0.5); % black is the 50 ch condition, yellow is lowest number of channels (10)
parCm(1, :) = [0, 0, 0];
parCm(4, :) = [0.4940 0.1840 0.5560];

% Configure the name of the folders to load accuracy scores from.
rootdir = 'simLoopNChs';
rootdir50 = 'simLoopNoGlob'; % We do not regenerate the 50 Chs condition, just pull from previous outputs

nChsToTest = [50, 25, 20, 15, 10]; % total number of channels
fracResp = 0:0.2:0.8; % fraction of all channels responsive
nReps = 30;

%% Read and store accuracy values at each number of total channels, then save to file

% each field is in order of forloop in simCCEP_testNChs: nchs x responsiveness x nReps
accuracy = struct();
accuracy.TP = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.TN = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.FN = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.FP = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.TPglob = nan(length(nChsToTest), length(fracResp), nReps); % for the global maxes
accuracy.TNglob = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.FNglob = nan(length(nChsToTest), length(fracResp), nReps);
accuracy.FPglob = nan(length(nChsToTest), length(fracResp), nReps);

% First store the 50 ch condition by pulling from previous output folder
for nR = 1:length(fracResp)
fprintf('.');
nResp = round(fracResp(nR)*50);
indir = fullfile('output', rootdir50, sprintf('nchs50-%d', nResp));
for rep = 1:nReps

T_acc = readtable(fullfile(indir, sprintf('accuracy_50-%d_rep%d.txt', nResp, rep)), 'FileType', 'text', 'Delimiter', '\t');
accuracy.TP(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
accuracy.TN(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
accuracy.FN(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
accuracy.FP(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
accuracy.TPglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob'));
accuracy.TNglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob'));
accuracy.FNglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob'));
accuracy.FPglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob'));

end
end

% now load for each total number of channels from rootdir
for nC = 2:length(nChsToTest)
for nR = 1:length(fracResp)
fprintf('.');
nChs = nChsToTest(nC);
indir = fullfile('output', rootdir, sprintf('nchs%d-%d', nChs, round(fracResp(nR)*nChs)));
for rep = 1:nReps
T_acc = readtable(fullfile(indir, sprintf('accuracy_rep%d.txt', rep)), 'FileType', 'text', 'Delimiter', '\t');
accuracy.TP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
accuracy.TN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
accuracy.FN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
accuracy.FP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
accuracy.TPglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob'));
accuracy.TNglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob'));
accuracy.FNglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob'));
accuracy.FPglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob'));
end
end
end
fprintf('\n');
save(fullfile('output', rootdir, 'accuracies_all.mat'), 'accuracy');

%% Plot FN, FP vs resp traces for each noise level, for first-peak method, as bars

load(fullfile('output', rootdir, 'accuracies_all.mat'), 'accuracy');

% displacements for the bars
disps = [-0.04, -0.02, 0, 0.02, 0.04];

% Based on more sensitive, first-peak optimum
figure('Position', [200, 200, 700, 800]);
subplot(2, 1, 1); % FN = responsive channels included in the CAR
hold on
for nc = 1:size(accuracy.FN, 1) % starting at 2 excludes the 50
FN = squeeze(accuracy.FN(nc, :, :));
FN = FN ./ (fracResp'*nChsToTest(nc)); % normalize by number of responsive channels
FN(isnan(FN)) = 0; % when there are 0 responsive channels, change these nans to 0
iqr = prctile(FN, [25, 75], 2);
med = median(FN, 2);
bar(fracResp + disps(nc), med, 0.08, 'FaceColor', parCm(nc*3-2, :));
errorbar(fracResp + disps(nc), med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineStyle', 'none', 'LineWidth', 1, 'CapSize', 8);
end
hold off
xlim([-0.05, 0.85]); ylim([-0.05, 1.05]);
set(gca, 'xtick', 0:0.2:0.8, 'xticklabels', 0:20:80, 'ytick', 0:0.2:1, 'yticklabels', 0:20:100, 'box', 'off');
ylabel('% Responsive Channels Missed');

subplot(2, 1, 2);
hold on
for nc = 1:size(accuracy.FN, 1)
FP = squeeze(accuracy.FP(nc, :, :));
FP = FP ./ ((1-fracResp)'*nChsToTest(nc)); % normalize by number of non-responsive channels
iqr = prctile(FP, [25, 75], 2);
med = median(FP, 2);
bar(fracResp + disps(nc), med, 0.08, 'FaceColor', parCm(nc*3-2, :));
errorbar(fracResp + disps(nc), med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineStyle', 'none', 'LineWidth', 1, 'CapSize', 8);
end
hold off
xlim([-0.05, 0.85]); ylim([-0.05, 1.05]);
set(gca, 'xtick', 0:0.2:0.8, 'xticklabels', 0:20:80, 'ytick', 0:0.2:1, 'yticklabels', 0:20:100, 'box', 'off');
xlabel('Fraction of Channels with Response');
ylabel('% Non-responsive Channels Missed');
saveas(gcf, fullfile('output', rootdir, 'FNFP_bars_acrossNChs'), 'png');
saveas(gcf, fullfile('output', rootdir, 'FNFP_bars_acrossNChs'), 'svg');
Loading

0 comments on commit a630235

Please sign in to comment.