From 38fbed44902713df424ac7860158e208bbe74e09 Mon Sep 17 00:00:00 2001 From: Harvey Huang Date: Tue, 23 Apr 2024 16:35:08 -0500 Subject: [PATCH 1/3] Changes to be consistent with upstream - Added personalDataPath to gitignore - removed duplicate files that persisted through perge --- .gitignore | 2 +- preprocNSDMef.m | 175 ---------------- script01_analyzePreprocScalarBB.m | 337 ------------------------------ script01_preprocNSDMef.m | 6 +- 4 files changed, 4 insertions(+), 516 deletions(-) delete mode 100644 preprocNSDMef.m delete mode 100644 script01_analyzePreprocScalarBB.m diff --git a/.gitignore b/.gitignore index fd389bd..308478b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ output/ -setPathsNSD.m +personalDataPath.m .DS_Store */.DS_Store diff --git a/preprocNSDMef.m b/preprocNSDMef.m deleted file mode 100644 index cfcfca0..0000000 --- a/preprocNSDMef.m +++ /dev/null @@ -1,175 +0,0 @@ -%% This is the preprocessing script for the clean NSD analysis. -% Requires each subject to have annotated channels, events, and electrodes SOZ - -%% Paths, load tables - -% cd here -setPathsNSD; -addpath('functions'); - -sub = 'X'; -[tasks, runs] = subjectConfig(sub); - -% electrodes path (for SOZ) -elecPath = fullfile(localRoot, sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg', sprintf('sub-%s_ses-ieeg01_electrodes.tsv', sub)); -elecs = ieeg_readtableRmHyphens(elecPath); - -% load file to indicate which tasks were performed in which sessions (to use as a possible random effect later) -sessions = readtable(fullfile('data', 'metadata', sprintf('%s_sessions.txt', sub)), 'FileType', 'text', 'Delimiter', '\t'); - -% All possible runs: col1 = mefPath, col2 = channels path, col3 = events path -dataPathStems = {}; -taskruns = []; -for ii = 1:length(runs) - for jj = 1:length(tasks) - dataPathStems = [dataPathStems; - {fullfile(localRoot, sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg', sprintf('sub-%s_ses-ieeg01_task-NSDspecial%02d_run-%02d', sub, tasks(jj), runs(ii)))}]; - taskruns = [taskruns; tasks(jj), runs(ii)]; % keep track of which NSD task and run - end -end - -% convert taskruns to table, this way 3rd column later will be for bad channels -taskruns = array2table(taskruns, 'VariableNames', {'tasknumber', 'run'}); -taskruns.badChs = repmat({'n/a'}, height(taskruns), 1); - -% contants -srate = 1200; - -%% Loop through NSD01 - NSD10 and all runs individually, concatenate output at the end - -Mephys = []; % ephys data matrix -Meye = []; % eyetracker data matrix -eventsST = []; % stim/test events - -for ii = 1:length(dataPathStems) - %% Load channels, events, mef data - - fprintf('Loading %s\n', dataPathStems{ii}); - - % load channels and events - try - channels = ieeg_readtableRmHyphens(sprintf('%s_channels.tsv', dataPathStems{ii})); - events = readtable(sprintf('%s_events.tsv', dataPathStems{ii}), 'FileType', 'text', 'Delimiter', '\t', 'TreatAsEmpty', 'n/a'); - catch - warning('Error loading channels/events file for %02d, skipping', ii); % e.g., no 2nd run for a certain NSD session - continue - end - - % force status_description to be text (so that it can be searched) in case all nans - events.status_description = string(events.status_description); - - % load mef data - [metadata, data] = readMef3(sprintf('%s_ieeg.mefd', dataPathStems{ii})); - - %% Reorganize and clean up events so that adjacent rest and stim events are combined - - % First event should be rest. Sometimes this pdiode onset is not acquired bc pdiode was on during presequence gray - if strcmp(events.trial_type(1), 'stim') - warning('Appending artificial rest trial to start of events'); - T1 = cell2table({events.onset(1) - 2, 2, events.sample_start(1) - 2*srate, 'rest', '', 'good', 'n/a'}, ... - 'VariableNames', events.Properties.VariableNames); - events = [T1; events]; - end - - % reassemble events table to be stim/test only. Add pre (rest) status and durations - idxST = find(ismember(events.trial_type, {'stim', 'test'})); - assert(all(strcmp(events.trial_type(idxST-1), 'rest')), 'Error: trials preceding some stim/test trials are not rest'); - eventsSTCurr = events(idxST, :); - eventsSTCurr.pre_duration = events.duration(idxST-1, :); - eventsSTCurr.pre_status = events.status(idxST-1, :); - eventsSTCurr.pre_status_description = events.status_description(idxST-1, :); - - % delete events that aren't good currently or in pre-rest, as well as high interictal - eventsSTCurr(~strcmp(eventsSTCurr.status, 'good') | ~strcmp(eventsSTCurr.pre_status, 'good'), :) = []; - eventsSTCurr(contains(eventsSTCurr.status_description, 'Categorical-level/High') | contains(eventsSTCurr.pre_status_description, 'Categorical-level/High'), :) = []; - - % add random effects time markers: task (NSD1-10), run (1, 2...), session (ordered blocks of time) - eventsSTCurr.tasknumber = repmat(taskruns.tasknumber(ii), height(eventsSTCurr), 1); - eventsSTCurr.run = repmat(taskruns.run(ii), height(eventsSTCurr), 1); - sess = sessions.session(sessions.task == taskruns.tasknumber(ii) & sessions.run == taskruns.run(ii)); - eventsSTCurr.session = repmat(sess, height(eventsSTCurr), 1); - - %% Preprocess and convert to trial struction - - % Highpass the SEEG channels - data(strcmp(channels.type, 'SEEG'), :) = ieeg_highpass(data(strcmp(channels.type, 'SEEG'), :)', srate)'; - - % Convert data to trial structure - samprange = -2*srate:2*srate-1; - - % remove any events missing the full sample range (e.g. due to crash) - eventsSTCurr(eventsSTCurr.sample_start+samprange(end)+1 > size(data, 2), :) = []; - - tt = samprange/srate; - M = zeros(size(data, 1), length(samprange), height(eventsSTCurr)); - for jj = 1:height(eventsSTCurr) - M(:, :, jj) = data(:, (eventsSTCurr.sample_start(jj)+samprange(1)+1) : (eventsSTCurr.sample_start(jj)+samprange(end)+1)); % convert to 1 indexing - end - - % check the DI1 channel to make sure that the first sample of 1 occurs at t=0 - %di1 = squeeze(M(strcmp(channels.name, 'DigitalInput1'), :, :)); - %figure; plot(tt, di1); xlim([-0.01, 0.01]); - - % Separate ephys from eyetracker data - MephysCurr = M(strcmp(channels.type, 'SEEG'), :, :); - if ii > 1 - assert(sum(strcmp(channels.type, 'SEEG')) == height(channelsEphys), 'Error: Mismatch between current number of ephys channels and prev number of ephys channels'); - end - channelsEphys = channels(strcmp(channels.type, 'SEEG'), :); - MeyeCurr = M(startsWith(channels.name, 'Eyetracker'), :, :); assert(size(MeyeCurr, 1) == 14, 'Error: Should have 14 eyetracker channels'); - eyetrackerLabels = channels.name(startsWith(channels.name, 'Eyetracker')); - - % Find all bad channels (including SOZ) - elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); % first set SOZ channels to bad - channelsEphys.status(ismember(channelsEphys.name, elecsSOZ)) = {'bad'}; - badChs = union(ieeg_getChsNR(channelsEphys, elecs), find(~strcmp(channelsEphys.status, 'good'))); % double check to see if any NR channels have been missed - - % CAR by lowest variance, with all trials considered the same group - opts = struct(); - opts.vartype = 'var'; - [MephysCurr, chsUsed] = ccep_CARVariance(tt, MephysCurr, srate, badChs, [], opts); - - %% Concatenate data, add info on good channels at current run - - Mephys = cat(3, Mephys, MephysCurr); - Meye = cat(3, Meye, MeyeCurr); - eventsST = [eventsST; eventsSTCurr]; - taskruns.badChs{ii} = badChs; - -end - -% remove taskruns entries that didn't happen -taskruns(strcmp(taskruns.badChs, 'n/a'), :) = []; - -%% Save all outputs - -% use most recent channelsephys (all are same height bc assertion), but remove status column. bad channels are indicated in taskruns variable -if any(ismember(channelsEphys.Properties.VariableNames, {'status', 'status_description'})) - channelsEphys.status = []; - channelsEphys.status_description = []; -end - -outdir = fullfile('output', 'fromMef', sub); -mkdir(outdir); -save(fullfile(outdir, sprintf('preprocessed_%s.mat', sub)), 'tt', 'srate', 'Mephys', 'Meye', 'eventsST', 'channelsEphys', 'taskruns', 'eyetrackerLabels', '-v7.3'); - -writetable(eventsST, fullfile(outdir, sprintf('eventsSTPreproc_%s.tsv', sub)), 'FileType', 'text', 'Delimiter', '\t'); - -% Plot and save a heatmap of ERPs, for channels good in at least 5 runs -% nan-out trials for channels in bad runs before averaging and plotting -MephysPlot = Mephys; -for ii = 1:height(taskruns) - MephysPlot(taskruns.badChs{ii}, :, eventsST.tasknumber == taskruns.tasknumber(ii) & eventsST.run == taskruns.run(ii)) = nan; -end -evoked = mean(MephysPlot, 3, 'omitnan'); % average across trials for each channel (only for runs where channels are good) -[goodChsOnN, badChsOnN] = goodChsNRuns(taskruns.badChs, height(channelsEphys), 5); % must be good in 5 runs to be a good channel to plot -evoked(badChsOnN, :) = nan; -figure('Position', [200, 200, 600, 1200]); -imagescNaN(tt, 1:height(channelsEphys), evoked, 'CLim', [-50, 50]); -xline(0, 'Color', 'r'); -xlim([-0.8, 0.8]); -set(gca, 'YTick', 1:2:height(channelsEphys), 'YTickLabels', channelsEphys.name(1:2:end)); -title('ERPs'); -xlabel('Time (s)'); ylabel('Channels'); -colorbar; -saveas(gcf, fullfile(outdir, sprintf('ERPsummary_%s.png', sub))); diff --git a/script01_analyzePreprocScalarBB.m b/script01_analyzePreprocScalarBB.m deleted file mode 100644 index bd7258d..0000000 --- a/script01_analyzePreprocScalarBB.m +++ /dev/null @@ -1,337 +0,0 @@ -%% Script to perform first-step ERP, broadband analysis on preprocessed Mef data -% First, load preprocessed data and all giftis - -% cd here -setPathsNSD; -addpath('functions'); -cmkjm = ieeg_kjmloccolormap(); - -sub = 'X'; - -% Determine which task/runs are bad to ignore, if any -fprintf('analyzePreprocScalarBB01.m on sub-%s\n', sub); -[~, ~, badtasks] = subjectConfig(sub); - -% Load electrodes -elecPath = fullfile(localRoot, sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg', sprintf('sub-%s_ses-ieeg01_electrodes.tsv', sub)); -elecs = ieeg_readtableRmHyphens(elecPath); - -% load preprocessed mef data from local copy -outdir = fullfile('output', 'fromMef', sub); -tic; load(fullfile(localRoot, sprintf('sub-%s', sub), 'NSDpreproc', sprintf('preprocessed_%s.mat', sub))); toc; - -% sort electrodes -elecs = ieeg_sortElectrodes(elecs, channelsEphys); - -% extract hemisphere information from electrodes -hemis = unique(cellfun(@(s) {s(1)}, lower(elecs.name))); -disp('Hemispheres:'); -disp(hemis); - -% Remove bad trials from events and taskruns table based on task/runs to ignore -badTrs = zeros(height(eventsST), 1, 'logical'); -for ii = 1:size(badtasks, 1) - badTrs(eventsST.tasknumber == badtasks(ii, 1) & eventsST.run == badtasks(ii, 2)) = true; - taskruns(taskruns.tasknumber == badtasks(ii, 1) & taskruns.run == badtasks(ii, 2), :) = []; -end -fprintf('Removing %d trials from bad tasks\n', sum(badTrs)); -eventsST(badTrs, :) = []; -Mephys(:, :, badTrs) = []; -Meye(:, :, badTrs) = []; - -% Make relevant output directories -mkdir(fullfile(outdir, 'ERP')); -mkdir(fullfile(outdir, 'broadband')); -mkdir(fullfile(outdir, 'giftiLabeled')); - -% Load giftis (pial, inflated, and white (used to transform to inflated)) -giis = cellfun(@(hemi) gifti(fullfile(localRoot, sprintf('sub-%s', sub), sprintf('pial.%s.surf.gii', hemi))), upper(hemis), 'UniformOutput', false); -giiInfs = cellfun(@(hemi) gifti(fullfile(localRoot, sprintf('sub-%s', sub), sprintf('inflated.%s.surf.gii', hemi))), upper(hemis), 'UniformOutput', false); -whites = cellfun(@(hemi) gifti(fullfile(localRoot, sprintf('sub-%s', sub), sprintf('white.%s.surf.gii', hemi))), upper(hemis), 'UniformOutput', false); -sulcs = cellfun(@(hemi) read_curv(fullfile(localRoot, sprintf('sub-%s', sub), sprintf('%sh.sulc', hemi))), lower(hemis), 'UniformOutput', false); -for ii = 1:length(hemis) - assert(length(giis{ii}.faces) == length(giiInfs{ii}.faces) && length(giis{ii}.faces) == length(whites{ii}.faces), 'Error: face mismatch between %dth gifti/giiInf/white', ii); -end - -% Create bipolar-rereferenced data, for broadband analysis. Load manually-written lead-segment metadata -segmentedLeads = readtable(fullfile('data', 'metadata', sprintf('%s_segmentedLeads.txt', sub)), 'FileType', 'text', 'Delimiter', '\t'); -seg5 = {}; seg6 = {}; -if ~isnumeric(segmentedLeads.seg5), seg5 = strtrim(split(upper(segmentedLeads.seg5), {' ', ','})); seg5(cellfun(@isempty, seg5)) = []; end -if ~isnumeric(segmentedLeads.seg6), seg6 = strtrim(split(upper(segmentedLeads.seg6), {' ', ','})); seg6(cellfun(@isempty, seg6)) = []; end -[~, bipolarNames, bipolarChans] = ieeg_bipolarSEEG(Mephys(:, :, 1)', channelsEphys.name, [], seg5, seg6); % initialize to get number of bipolar channels -MephysBip = nan(length(bipolarNames), length(tt), height(eventsST)); -for ii = 1:height(eventsST) - MephysBip(:, :, ii) = ieeg_bipolarSEEG(Mephys(:, :, ii)', channelsEphys.name, [], seg5, seg6, false)'; % don't consider bad channels here -end - -% Determine run-specific bipolar bad channels -for ii = 1:height(taskruns) - [~, ~, ~, badChs] = ieeg_bipolarSEEG(Mephys(:, :, 1)', channelsEphys.name, taskruns.badChs{ii}, seg5, seg6, false); - taskruns.badChsBip{ii} = badChs; -end - -% get interpolated bipolar pair coordinates -elecsBip = ieeg_Els2BipLocs(elecs, bipolarNames); - -% Make masks indicating for each channel with trials are bad (i.e., which task-runs are bad). % can add to this later other bad trials -badTrsMask = zeros(height(channelsEphys), height(eventsST), 'logical'); -badTrsMaskBip = zeros(length(bipolarNames), height(eventsST), 'logical'); -for ii = 1:height(taskruns) - badTrsMask(taskruns.badChs{ii}, eventsST.tasknumber == taskruns.tasknumber(ii) & eventsST.run == taskruns.run(ii)) = true; - badTrsMaskBip(taskruns.badChsBip{ii}, eventsST.tasknumber == taskruns.tasknumber(ii) & eventsST.run == taskruns.run(ii)) = true; -end - -%% Calculate inflated pial positions for elecs and bipolar-interpolated elecs, add to elecs and elecsBip tables - -fprintf('Plotting labeled pial and inflated pials\n'); - -% Calculate inflated bipolar electrode positions -elecs.xInf = nan(height(elecs), 1); -elecs.yInf = nan(height(elecs), 1); -elecs.zInf = nan(height(elecs), 1); -elecsBip.xInf = nan(height(elecsBip), 1); -elecsBip.yInf = nan(height(elecsBip), 1); -elecsBip.zInf = nan(height(elecsBip), 1); - -% check each hemisphere separately -for ii = 1:length(hemis) - - % Calculate inflated positions for real electrodes - xyzInf = getXyzsInf(elecs, hemis{ii}, whites{ii}, giiInfs{ii}, 4, false); % must be within 4 mm of Gray/white boundary - - % also check 2 mm of pial surface, add any additional discovered in case of thick GM - xyzInfPial = getXyzsInf(elecs, hemis{ii}, giis{ii}, giiInfs{ii}, 2, false); - extraElecsPial = isnan(xyzInf(:, 1)) & ~isnan(xyzInfPial(:, 1)); - xyzInf(extraElecsPial, :) = xyzInfPial(extraElecsPial, :); - - idxes = ~isnan(xyzInf(:, 1)); % inflated positions found for current hemisphere - elecs.xInf(idxes) = xyzInf(idxes, 1); elecs.yInf(idxes) = xyzInf(idxes, 2); elecs.zInf(idxes) = xyzInf(idxes, 3); - - % Calculate inflated positions for bipolar-interpolated electrode pairs - xyzInf = getXyzsInf(elecsBip, hemis{ii}, whites{ii}, giiInfs{ii}, 4, false); % must be within 4 mm of Gray/white boundary - - % also check 2 mm of pial surface, add any additional discovered in case of thick GM - xyzInfPial = getXyzsInf(elecsBip, hemis{ii}, giis{ii}, giiInfs{ii}, 2, false); - extraElecsPial = isnan(xyzInf(:, 1)) & ~isnan(xyzInfPial(:, 1)); - xyzInf(extraElecsPial, :) = xyzInfPial(extraElecsPial, :); - - idxes = ~isnan(xyzInf(:, 1)); % inflated positions found for current hemisphere - elecsBip.xInf(idxes) = xyzInf(idxes, 1); elecsBip.yInf(idxes) = xyzInf(idxes, 2); elecsBip.zInf(idxes) = xyzInf(idxes, 3); -end - -% find which channels are good in at least one run -goodChsOn1 = goodChsNRuns(taskruns.badChs, height(elecs), 1); -goodChsBipOn1 = goodChsNRuns(taskruns.badChsBip, height(elecsBip), 1); - -% Plot electrode positions to gifti for reference (good elecs only) -plotElectrodesToGifti(elecs(goodChsOn1, :), elecsBip(goodChsBipOn1, :), giis, giiInfs, sulcs, hemis, fullfile(outdir, 'giftiLabeled'), 'off'); - -%% Calculate PSDs and average (scalar) broadband power for CAR electrodes - -fprintf('CAR electrodes\n'); - -% include channels that are good on 5+ runs -[goodChsOn5, badChsOn5] = goodChsNRuns(taskruns.badChs, height(elecs), 5); - -% Calculate PSD for each trial -intervalStim = [0, 0.5]; intervalRest = [-0.5, 0]; % stim = first 500ms after stimulus, rest = last 500 ms before stimulus -pwelchWin = hann(srate/4); -psdStim = nan(height(channelsEphys), srate/2+1, height(eventsST)); -psdRest = nan(height(channelsEphys), srate/2+1, height(eventsST)); -for ii = 1:height(channelsEphys) - [psdStim(ii, :, :), f] = pwelch(squeeze(Mephys(ii, tt >= intervalStim(1) & tt < intervalStim(2), :)), pwelchWin, length(pwelchWin)/2, srate, srate); - psdRest(ii, :, :) = pwelch(squeeze(Mephys(ii, tt >= intervalRest(1) & tt < intervalRest(2), :)), pwelchWin, length(pwelchWin)/2, srate, srate); -end - -% Plot average stationary PSD for all channels and save -toplot = nan(height(channelsEphys), length(f)); -for ii = 1:height(channelsEphys) - toplot(ii, :) = mean(log10(psdStim(ii, :, ~badTrsMask(ii, :))), 3, 'omitnan') - mean(log10(psdRest(ii, :, ~badTrsMask(ii, :))), 3, 'omitnan'); -end -toplot(badChsOn5, :) = nan; -figure('Position', [200, 200, 600, 1200], 'Visible', 'off'); -imagescNaN(f, 1:height(channelsEphys), toplot, 'CLim', [-0.2, 0.2]); colormap(cmkjm); -xlim([0, 200]); -set(gca, 'YTick', 1:2:height(channelsEphys), 'YTickLabels', channelsEphys.name(1:2:end)); -xlabel('Frequency (Hz)'); -ylabel('channels'); -colorbar; -saveas(gcf, fullfile(outdir, sprintf('PSDsummary_%s.png', sub))); - -% Calculate average broadband power for each trial (geomean across frequency bins as in Yuasa et al., 2023) -BBpowerStim = squeeze(geomean(psdStim(:, f >= 70 & f <= 115 | f >= 125 & f <= 170, :), 2)); -BBpowerRest = squeeze(geomean(psdRest(:, f >= 70 & f <= 115 | f >= 125 & f <= 170, :), 2)); - -%% Plot Rsq of BB power stim vs rest as activations to brain - -% Calculate rsq for each channel, across trials that are good for each channel -rsqBB = nan(height(channelsEphys), 2); -for ii = 1:height(channelsEphys) - [rsq, p] = mnl_rsq(BBpowerStim(ii, ~badTrsMask(ii, :)), BBpowerRest(ii, ~badTrsMask(ii, :))); - rsqBB(ii, :) = [rsq, p]; -end -% set bad channels (based on N runs) to nan -rsqBB(badChsOn5, :) = nan; - -% Add to bipolar electrodes table and write as output -elecs.BBrsq = rsqBB(:, 1); -elecs.BBrsqP = rsqBB(:, 2); -writetable(elecs, fullfile(outdir, 'elecsBB1.tsv'), 'FileType', 'text', 'Delimiter', '\t'); - -% Plot broadband Rsq activations on gifti and inflated gifti -plotActivationsToGiftis([], elecs, elecs.BBrsq, giis, giiInfs, sulcs, hemis, 0.5, fullfile(outdir, 'broadband'), 'BBrsq', 'off'); - -% save a colorbar -fCbar = plotCtxWtsAdjCbar(0.5, 'SigFraction', 0.05, 'SizeJump', 2); -saveas(fCbar, fullfile(outdir, 'broadband', 'cbar_BBrsq.png')); - -%% Calculate noise ceiling SNR for each CAR electrode using the special100 images - -% find indices of the special100 images -special100Idxes = getSpecial100Idxes(eventsST.stim_file); - -% copy BBpowerStim, and put nans in bad chxtrs. These are removed in estimateNCSNR function -BBpowerStimNan = BBpowerStim; -BBpowerStimNan(badTrsMask) = nan; - -rng('default'); - -NCSNRs = nan(height(channelsEphys), 2); % cols = NCSNR, p -for ch = 1:height(channelsEphys) - - fprintf('.'); - if any(ch == badChsOn5), continue; end - - % pull out broadband power, normalize by mean rest power and convert to log power for normality - BBpowerRestMean = geomean(BBpowerRest(ch, ~badTrsMask(ch, :))); - BBlogpowerStimSpecial100 = cell(100, 1); - for ss = 1:100 - BBlogpowerStimSpecial100{ss} = log10(BBpowerStimNan(ch, special100Idxes{ss, 2}) / BBpowerRestMean); - end - - % calculate SNR and pvalue - [NCSNR, p] = estimateNCSNR(BBlogpowerStimSpecial100); - NCSNRs(ch, 1) = NCSNR; NCSNRs(ch, 2) = p; - -end -fprintf('\n'); - -elecs.BBncsnr = NCSNRs(:, 1); -elecs.BBncsnrP = NCSNRs(:, 2); -writetable(elecs, fullfile(outdir, 'elecsBB2.tsv'), 'FileType', 'text', 'Delimiter', '\t'); - -NCSNRsSig = NCSNRs(:, 1); -NCSNRsSig(NCSNRs(:, 2) > 0.05) = 0; -plotActivationsToGiftis([], elecs, NCSNRsSig, giis, giiInfs, sulcs, hemis, 0.6, fullfile(outdir, 'broadband'), 'BBncsnr', 'off'); - -% save a colorbar -fCbar = plotCtxWtsAdjCbar(0.6, 'SigFraction', 0.05, 'SizeJump', 2, 'InsigText', 'Insig.'); -saveas(fCbar, fullfile(outdir, 'broadband', 'cbar_BBncsnr.png')); - - - -%% BIPOLAR -%% Scalar average broadband power, calculated for bipolar-reref electrodes -% Add BB rsq and p-value to bipolar electrodes table, save, and also plot mean BB rsq on giftis and inflated giftis - -fprintf('Bipolar electrodes analysis\n'); - -% include bipolar channels that are good on 5+ runs -[goodChsBipOn5, badChsBipOn5] = goodChsNRuns(taskruns.badChsBip, height(elecsBip), 5); - -% Calculate PSD for each trial -intervalStim = [0, 0.5]; intervalRest = [-0.5, 0]; % stim = first 500ms after stimulus, rest = last 500 ms before stimulus -pwelchWin = hann(srate/4); -psdStimBip = nan(length(bipolarNames), srate/2+1, height(eventsST)); -psdRestBip = nan(length(bipolarNames), srate/2+1, height(eventsST)); -for ii = 1:length(bipolarNames) - [psdStimBip(ii, :, :), f] = pwelch(squeeze(MephysBip(ii, tt >= intervalStim(1) & tt < intervalStim(2), :)), pwelchWin, length(pwelchWin)/2, srate, srate); - psdRestBip(ii, :, :) = pwelch(squeeze(MephysBip(ii, tt >= intervalRest(1) & tt < intervalRest(2), :)), pwelchWin, length(pwelchWin)/2, srate, srate); -end - -% Plot average stationary PSD for all bipolar channels and save -toplot = nan(length(bipolarNames), length(f)); -for ii = 1:length(bipolarNames) - toplot(ii, :) = mean(log10(psdStimBip(ii, :, ~badTrsMaskBip(ii, :))), 3, 'omitnan') - mean(log10(psdRestBip(ii, :, ~badTrsMaskBip(ii, :))), 3, 'omitnan'); -end -toplot(badChsBipOn5, :) = nan; -figure('Position', [200, 200, 600, 1200], 'Visible', 'off'); -imagescNaN(f, 1:length(bipolarNames), toplot, 'CLim', [-0.2, 0.2]); colormap(cmkjm); -xlim([0, 200]); -set(gca, 'YTick', 1:length(bipolarNames), 'YTickLabels', bipolarNames); -xlabel('Frequency (Hz)'); -ylabel('channels'); -colorbar; -saveas(gcf, fullfile(outdir, sprintf('PSDsummaryBipolar_%s.png', sub))); - -% Calculate average broadband power for each trial (geomean across frequency bins as in Yuasa et al., 2023) -BBpowerStimBip = squeeze(geomean(psdStimBip(:, f >= 70 & f <= 115 | f >= 125 & f <= 170, :), 2)); -BBpowerRestBip = squeeze(geomean(psdRestBip(:, f >= 70 & f <= 115 | f >= 125 & f <= 170, :), 2)); - -% consider normalizing by sessions. But then again, R^2 should still work in the absence of adjusting for any other effects (univariate) -% % Normalize stim and rest power by mean rest power in each session, to account for large impedance differences over time -% sessions = unique(eventsST.session); -% for ii = 1:length(sessions) -% sessionBaseline = geomean(BBpowerRestBip(:, eventsST.session == sessions(ii)), 2); -% temp(:, eventsST.session == sessions(ii)) = temp(:, eventsST.session == sessions(ii)) ./ sessionBaseline; -% end - -%% Calculate rsq for each channel, across trials good for those channels - -rsqBBBip = nan(length(bipolarNames), 2); -for ii = 1:length(bipolarNames) - [rsq, p] = mnl_rsq(BBpowerStimBip(ii, ~badTrsMaskBip(ii, :)), BBpowerRestBip(ii, ~badTrsMaskBip(ii, :))); - rsqBBBip(ii, :) = [rsq, p]; -end -% set bad channels (in any run) to nan -rsqBBBip(badChsBipOn5, :) = nan; - -% Add to bipolar electrodes table and write as output -elecsBip.BBrsq = rsqBBBip(:, 1); -elecsBip.BBrsqP = rsqBBBip(:, 2); -writetable(elecsBip, fullfile(outdir, 'elecsBipBB1.tsv'), 'FileType', 'text', 'Delimiter', '\t'); - -% Plot broadband Rsq activations on gifti and inflated gifti -plotActivationsToGiftis(elecs, elecsBip, elecsBip.BBrsq, giis, giiInfs, sulcs, hemis, 0.5, fullfile(outdir, 'broadband'), 'BBrsqBip', 'off'); - -%% Calculate noise ceiling SNR for each bipolar electrode using the special100 images - -% find indices of the special100 images -special100Idxes = getSpecial100Idxes(eventsST.stim_file); - -% copy BBpowerStimBip, and put nans in bad chxtrs. These nans are removed in estimateNCSNR function -BBpowerStimBipNan = BBpowerStimBip; -BBpowerStimBipNan(badTrsMaskBip) = nan; - -rng('default'); - -NCSNRsBip = nan(length(bipolarNames), 2); % cols = NCSNR, p -for ch = 1:length(bipolarNames) - - fprintf('.'); - if any(ch == badChsBipOn5), continue; end - - % pull out broadband power, normalize by mean rest power and convert to log power for normality - BBpowerRestMean = geomean(BBpowerRestBip(ch, ~badTrsMaskBip(ch, :))); - BBlogpowerStimSpecial100 = cell(100, 1); - for ss = 1:100 - BBlogpowerStimSpecial100{ss} = log10(BBpowerStimBipNan(ch, special100Idxes{ss, 2}) / BBpowerRestMean); - end - - % calculate SNR and pvalue - [NCSNR, p] = estimateNCSNR(BBlogpowerStimSpecial100); - NCSNRsBip(ch, 1) = NCSNR; NCSNRsBip(ch, 2) = p; - -end -fprintf('\n'); - -elecsBip.BBncsnr = NCSNRsBip(:, 1); -elecsBip.BBncsnrP = NCSNRsBip(:, 2); -writetable(elecsBip, fullfile(outdir, 'elecsBipBB2.tsv'), 'FileType', 'text', 'Delimiter', '\t'); - -NCSNRsBipSig = NCSNRsBip(:, 1); -NCSNRsBipSig(NCSNRsBip(:, 2) > 0.05) = 0; -plotActivationsToGiftis(elecs, elecsBip, NCSNRsBipSig, giis, giiInfs, sulcs, hemis, 0.6, fullfile(outdir, 'broadband'), 'BBncsnrBip', 'off'); - - - diff --git a/script01_preprocNSDMef.m b/script01_preprocNSDMef.m index d92b542..f39f549 100644 --- a/script01_preprocNSDMef.m +++ b/script01_preprocNSDMef.m @@ -9,10 +9,10 @@ localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) addpath('functions'); -subjects = {'01','02','03','04','05','06','07'}; +% subject to preprocess +ss = 1; +sub_label = sprintf('%02d', ss); -ss = 7; -sub_label = subjects{ss}; ses_label = 'ieeg01'; % list of task labels that should be preprocessed and concatenated task_labels = {'NSDspecial01','NSDspecial02','NSDspecial03','NSDspecial04','NSDspecial05',... From b96a5f6d7ddc02799187c98408c5604845521260 Mon Sep 17 00:00:00 2001 From: Harvey Huang Date: Tue, 23 Apr 2024 16:35:46 -0500 Subject: [PATCH 2/3] Update .gitignore to ignore data directory --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 308478b..46184aa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +data/ output/ personalDataPath.m .DS_Store From d083c7277555da4578a683b57eeaec4c2858d18d Mon Sep 17 00:00:00 2001 From: Harvey Huang Date: Thu, 25 Apr 2024 10:05:37 -0500 Subject: [PATCH 3/3] Coded added to script02_writeMNI.m, plus other minor changes script02_writeMNI.m - added 2 blocks of code: writes MNI152 positions to a new electrodes.tsv (using SPM forward deformation fields), and writes MNI305 positions to a new electrodes.tsv using fsaverage gitignore - added ignores for *.asv and personalDataPath.m script01_preprocNSDMef - elecPath -> now elecsPath --- .gitignore | 1 + script01_preprocNSDMef.m | 10 +++---- script02_writeMNI.m | 58 +++++++++++++++++++++++++++++++--------- 3 files changed, 52 insertions(+), 17 deletions(-) diff --git a/.gitignore b/.gitignore index 46184aa..337276f 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ output/ personalDataPath.m .DS_Store */.DS_Store +*.asv \ No newline at end of file diff --git a/script01_preprocNSDMef.m b/script01_preprocNSDMef.m index f39f549..553b28e 100644 --- a/script01_preprocNSDMef.m +++ b/script01_preprocNSDMef.m @@ -10,7 +10,7 @@ addpath('functions'); % subject to preprocess -ss = 1; +ss = 17; sub_label = sprintf('%02d', ss); ses_label = 'ieeg01'; @@ -20,8 +20,8 @@ task_list = readtable(fullfile(localDataPath.input,['sub-' sub_label],['ses-' ses_label],['sub-' sub_label '_ses-' ses_label '_scans.tsv']), 'FileType', 'text', 'Delimiter', '\t', 'TreatAsEmpty', 'n/a'); % electrodes path (to exclude electrodes on the SOZ) -elecPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); -elecs = ieeg_readtableRmHyphens(elecPath); +elecsPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); +elecs = ieeg_readtableRmHyphens(elecsPath); % For all possible runs, set data_info data_info = []; @@ -253,8 +253,8 @@ %% render and plot noise ceiling SNR -elecPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); -elecs = ieeg_readtableRmHyphens(elecPath); +elecsPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); +elecs = ieeg_readtableRmHyphens(elecsPath); name = all_channels.name; all_channels_table = table(name); diff --git a/script02_writeMNI.m b/script02_writeMNI.m index 203345e..49240d2 100644 --- a/script02_writeMNI.m +++ b/script02_writeMNI.m @@ -1,24 +1,58 @@ +%% This script calculates MNI152 and MNI305 positions each subject and saves them -clear all localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) addpath('functions'); -subjects = {'01','02','03','04','05','06'}; +ss = 17; +sub_label = sprintf('%02d', ss); -ss = 4; -sub_label = subjects{ss}; ses_label = 'ieeg01'; - outdir = fullfile(localDataPath.output,'derivatives','preproc_car',['sub-' sub_label]); -% load(fullfile(outdir, ['sub-' sub_label '_desc-preprocCAR_ieeg.mat']), 'tt', 'srate', 'Mdata', 'eventsST', 'all_channels'); -load(fullfile(outdir, ['sub-' sub_label '_desc-preprocCARBB_ieeg.mat']), 'tt', 'srate', 'Mbb', 'eventsST', 'all_channels'); -readtable(fullfile(outdir, ['sub-' sub_label '_desc-preprocCAR_events.tsv']), 'FileType', 'text', 'Delimiter', '\t', 'TreatAsEmpty', 'n/a'); +% load electrodes +elecsPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); +elecs = ieeg_readtableRmHyphens(elecsPath); +elecmatrix = [elecs.x, elecs.y, elecs.z]; + +%% Get and save MNI152 positions to electrodesMni152.tsv (volumetric, SPM12) + +% locate forward deformation field from SPM. There are variabilities in session name, so we use dir to find a matching one +niiPath = dir(fullfile(localDataPath.input, 'sourcedata', 'spm_forward_deformation_fields', sprintf('sub-%s_ses-*_T1w_acpc.nii', sub_label))); +assert(length(niiPath) == 1, 'Error: did not find exactly one match in sourcedata T1w MRI'); % check for only one unique match +niiPath = fullfile(niiPath.folder, niiPath.name); + +% create a location in derivatives to save the transformed electrode images +rootdirMni = fullfile(localDataPath.input, 'derivatives', 'MNI152_electrode_transformations', sprintf('sub-%s', sub_label)); +mkdir(rootdirMni); + +% calculate MNI152 coordinates for electrodes +xyzMni152 = ieeg_getXyzMni(elecmatrix, niiPath, rootdirMni); + +% save as separate MNI 152 electrodes table +elecsMni152Path = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_space-' 'MNI152NLin2009' '_electrodes.tsv']); +elecsMni152 = elecs; +elecsMni152.x = xyzMni152(:, 1); elecsMni152.y = xyzMni152(:, 2); elecsMni152.z = xyzMni152(:, 3); +writetable(elecsMni152, elecsMni152Path, 'FileType', 'text', 'Delimiter', '\t'); + +fprintf('Saved to %s\n', elecsMni152Path); + +%% Get and save MNI305 positions (through fsaverage) + +% FS dir of current subject +FSdir = fullfile(localDataPath.input, 'sourcedata', 'freesurfer', sprintf('sub-%s', sub_label)); +FSsubjectsdir = fullfile(FSdir, '..'); -%% need to get MNI positions first +% calculate MNI305 coordinates for electrodes +[xyzMni305, vertIdxFsavg, minDists, surfUsed] = ieeg_mni305ThroughFsSphere(elecmatrix, elecs.hemisphere, FSdir, FSsubjectsdir, 'closest', 5); -% mni_coords = ieeg_mni305ThroughFsSphere(elecmatrix,hemi,FSdir,FSsubjectsdir) +% save as separate MNI 305 electrodes table +elecsMni305Path = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_space-' 'MNI305' '_electrodes.tsv']); +elecsMni305 = elecs; +elecsMni305.x = xyzMni305(:, 1); elecsMni305.y = xyzMni305(:, 2); elecsMni305.z = xyzMni305(:, 3); +elecsMni305.vertex_fsaverage = vertIdxFsavg; % also add a column to indicate vertex on fsavg, so we can easily get position for inflated brain +writetable(elecsMni305, elecsMni305Path, 'FileType', 'text', 'Delimiter', '\t'); +fprintf('Saved to %s\n', elecsMni305Path); %% Normalize bb power per run @@ -66,8 +100,8 @@ %% render and plot noise ceiling SNR -elecPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); -elecs = ieeg_readtableRmHyphens(elecPath); +elecsPath = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_electrodes.tsv']); +elecs = ieeg_readtableRmHyphens(elecsPath); name = all_channels.name; all_channels_table = table(name);