diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fd389bd --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +output/ +setPathsNSD.m +.DS_Store +*/.DS_Store diff --git a/README_preprocessing_notes_HH.txt b/README_preprocessing_notes_HH.txt new file mode 100644 index 0000000..1a5d02c --- /dev/null +++ b/README_preprocessing_notes_HH.txt @@ -0,0 +1,30 @@ +Follow order of code in main.m + +1. preprocNSDMef + - subjectConfig: which tasks and runs (copies) to load? Which out of these are bad to skip? + - metadata/sub_sessions.txt -- loaded and added as column in events to use as possible predictor + - load mef data + - create events structure that only has stim and test trials. Rest are "pre_duration" + - Keep only trials that are good DURING trial and in rest before trial + - Also remove trials with high interictals during or in rest before trial + - random effects predictors: task number run, session + - ieeg_highpass data + - convert data to match events, on interval [-2, 2)s. Remove any trials missing this full range + - re-referencing by CARVariance, with 'var' option and bottom 25%ile. badChs are SOZ and "bad" chs + - concatenate all trials across all runs + - badChs unique to each run stored in taskruns variable + - save data as .mat file + - save eventsST as text file + - currently, EKG channel is not saved + +2. analyzePreprocScalarBB + - Remove all trials corresponding to bad tasks/runs + - bipolar-reference electrodes, using manually written lead segment data in metadata/. + - for each run, find which bipolar electrodes are bad. Create masks (chs x trs) that equal 1 when the channel is bad in the run in which the trial is located + - for electrode localization, I use getXyzsInf to inflate electrodes. Checks 4mm of white, and also 2mm of pial, which is assigned if none found in proximity to white. + - aCAR electrode analysis + - PSDs: hann(srate/4), pwelch on 0-0.5s, normalized by -0.5-0s. + - scalar broadband power = geometric mean of PSD from 70 to 170 Hz in 1 Hz bins, ignoring 115-125 range + - rsqs are calculated using power. Keep only channels with at least 5 good runs. Saved as elecsBB1.tsv -> do this in log power as well + - Noise ceiling SNR and permutation p-value calculated with estimateNCSNR on log-transformed scalar broadband power. Saved as elecsBB2.tsv + - The same thing as above is done for the bipolar-referenced electrodes \ No newline at end of file diff --git a/data/metadata/sub-X_segmentedLeads.txt b/data/metadata/sub-X_segmentedLeads.txt new file mode 100644 index 0000000..2d98430 --- /dev/null +++ b/data/metadata/sub-X_segmentedLeads.txt @@ -0,0 +1,2 @@ +seg5 seg6 +LB LA, LC \ No newline at end of file diff --git a/data/metadata/sub-X_sessions.txt b/data/metadata/sub-X_sessions.txt new file mode 100644 index 0000000..c9aa0aa --- /dev/null +++ b/data/metadata/sub-X_sessions.txt @@ -0,0 +1,11 @@ +task run session notes +1 1 1 +2 1 1 +3 1 2 +4 1 2 +5 1 3 +6 1 3 notes about something important in this run +7 1 3 +8 1 4 +9 1 4 +10 1 4 \ No newline at end of file diff --git a/data/stims/special100Names.txt b/data/stims/special100Names.txt new file mode 100644 index 0000000..4df3b66 --- /dev/null +++ b/data/stims/special100Names.txt @@ -0,0 +1,100 @@ +shared0004_nsd03078 +shared0008_nsd03172 +shared0022_nsd03914 +shared0030_nsd04424 +shared0033_nsd04668 +shared0052_nsd05584 +shared0064_nsd06522 +shared0069_nsd06802 +shared0073_nsd07208 +shared0137_nsd11636 +shared0139_nsd11726 +shared0140_nsd11797 +shared0145_nsd11943 +shared0157_nsd12923 +shared0159_nsd13139 +shared0163_nsd13614 +shared0186_nsd15365 +shared0194_nsd16345 +shared0197_nsd16617 +shared0211_nsd17375 +shared0234_nsd19437 +shared0267_nsd21951 +shared0287_nsd22880 +shared0300_nsd23994 +shared0307_nsd24531 +shared0310_nsd24641 +shared0318_nsd25112 +shared0326_nsd25372 +shared0334_nsd25703 +shared0350_nsd26459 +shared0358_nsd27243 +shared0362_nsd27436 +shared0369_nsd28069 +shared0378_nsd28350 +shared0382_nsd28596 +shared0404_nsd30396 +shared0405_nsd30408 +shared0425_nsd31660 +shared0463_nsd34830 +shared0474_nsd35987 +shared0487_nsd36975 +shared0488_nsd36979 +shared0491_nsd37225 +shared0498_nsd38023 +shared0507_nsd38487 +shared0520_nsd39370 +shared0530_nsd40549 +shared0535_nsd40847 +shared0568_nsd42698 +shared0570_nsd42852 +shared0579_nsd43225 +shared0588_nsd43676 +shared0589_nsd43690 +shared0591_nsd43819 +shared0610_nsd44845 +shared0614_nsd45214 +shared0616_nsd45596 +shared0623_nsd45982 +shared0634_nsd46373 +shared0646_nsd47071 +shared0650_nsd47294 +shared0689_nsd50115 +shared0694_nsd50756 +shared0695_nsd50812 +shared0700_nsd51078 +shared0727_nsd53053 +shared0730_nsd53156 +shared0733_nsd53371 +shared0745_nsd54258 +shared0746_nsd54362 +shared0754_nsd54914 +shared0764_nsd55679 +shared0768_nsd55969 +shared0786_nsd56785 +shared0789_nsd56912 +shared0790_nsd56949 +shared0797_nsd57554 +shared0811_nsd59024 +shared0825_nsd59995 +shared0853_nsd61874 +shared0857_nsd62210 +shared0869_nsd63183 +shared0876_nsd63932 +shared0882_nsd64499 +shared0896_nsd65254 +shared0905_nsd65800 +shared0910_nsd66005 +shared0925_nsd66977 +shared0936_nsd67830 +shared0941_nsd68340 +shared0944_nsd68742 +shared0948_nsd68898 +shared0960_nsd69814 +shared0962_nsd69840 +shared0968_nsd70194 +shared0969_nsd70233 +shared0974_nsd70506 +shared0986_nsd71411 +shared0991_nsd72016 +shared0999_nsd72720 diff --git a/functions/estimateNCSNR.m b/functions/estimateNCSNR.m new file mode 100644 index 0000000..31c17dd --- /dev/null +++ b/functions/estimateNCSNR.m @@ -0,0 +1,90 @@ +%% This function estimates the NCSNR of input signal strengths, as the SD of signal divided by the SD of noise +% Derived from estimateNoiseCeiling.m +% +% Variance due to noise is estimated from trials of the same stimulus, averaged across all stimuli. The remaining +% variance is explained by signal, at the single trial level. If numAvg is given, will return the noise ceiling +% corresponding to averaging numAvg trials together (noise variance is reduced numAvg-fold, assuming Gaussian). +% +% NCSNR = estimateNoiseCeiling(strength); +% [NCSNR, p] = estimateNoiseCeiling(strength, nperm); +% strength = nx1 cell, measure of scalar signal strength for each trial with n unique stimuli and variable trials per stimulus. +% Examples include average bandpower or a coefficient fit for the trial on some template (e.g. fmri beta) +% nPerm = (optional) positive integer, number of permutations to use to construct permutation null model. Default = 1000 +% +% Returns: +% NCSNR = double, SNR as (SD of Signal)/(SD of noise). +% +% Based on: +% Allen, E.J., St-Yves, G., Wu, Y. et al. A massive 7T fMRI dataset to bridge cognitive neuroscience and artificial intelligence. +% Nat Neurosci 25, 116?126 (2022). https://doi.org/10.1038/s41593-021-00962-x +% +% For derivations and detailed explanation, please see noise ceiling estimation in Methods section of article. +% - Some more info is in "Functional data (NSD)/Noise ceiling" in NSD data manual, https://cvnlab.slite.com/p/channel/CPyFRAyDYpxdkPK6YbB5R1 +% +% Demo: +% NCSNRtest = zeros(10000, 1); +% for ii = 1:10000 +% xx = mat2cell(rand*10*randn(600, 1) + rand*10, 6*ones(100, 1)); +% NCSNRtest(ii) = estimateNCSNR(xx); +% end +% figure; histogram(NCSNRtest, 30); +% +% HH 2022/05 +% +function [NCSNR, p, NCSNRNull] = estimateNCSNR(strength, nPerm) + + + % Formatting: vertical cell array and vertical array in each cell, remove nan entries + strength = strength(:); + strength = cellfun(@(x) x(:), strength, 'UniformOutput', false); + strength = cellfun(@(x) x(~isnan(x)), strength, 'UniformOutput', false); + + % number of entries in each cell + cellLengths = cellfun(@length, strength); + assert(sum(cellLengths >= 3) > 0.5*length(cellLengths), 'Error: More than half of conditions have <3 trials. Noise SD is poorly estimated'); + + % Zscore strength so SD = var = 1 + meanStrength = mean(cell2mat(strength)); + SDStrength = std(cell2mat(strength), 0); % using N-1 unbiased estimator + Zstr = cellfun(@(x) (x - meanStrength)/SDStrength, strength, 'UniformOutput', false); + + % Calculate NCSNR + NCSNR = calcNCSNR(Zstr); + + if nargout < 2, return; end % no p-value requested + + % Calculate null distribution and p-value + if nargin < 2 || isempty(nPerm), nPerm = 1000; end + assert(nPerm > 0, 'Error: nPerm must be a positive integer'); + NCSNRNull = nan(nPerm, 1); + ZstrAll = cell2mat(Zstr); + n = length(ZstrAll); % total number of strength entries + for ii = 1:nPerm + Perm = randperm(n); + ZstrNull = mat2cell(ZstrAll(Perm), cellLengths); + NCSNRNull(ii) = calcNCSNR(ZstrNull); + end + p = sum(NCSNRNull >= NCSNR) / nPerm; + +end + +% Helper function so as to call it directly with normalized ZstrNull +function NCSNR = calcNCSNR(Zstr) + + cellLengths = cellfun(@length, Zstr); + + % residual variance across trials of each stimuli, averaged across stimuli. + % Reason to avg stimuli with >= 3 trials bc Bessel correction is still biased for SD (underestimates noise SD here). 3 matches sample size in NSD fMRI + varNoise = mean(cellfun(@(x) var(x, 0), Zstr(cellLengths >= 3))); + sdNoise = sqrt(varNoise); + + varSig = 1 - varNoise; % variability explained by signal. Check this distribution using demo. White noise is expected to have varSig ~0 on avg + + if varSig <= 0 + sdSig = 0; % half-rectification of signal SD (signal variance can be negative like R^2) + else + sdSig = sqrt(varSig); % sqrt of variance explained by signal + end + NCSNR = sdSig/sdNoise; % noise ceiling signal-to-noise ratio + +end \ No newline at end of file diff --git a/functions/getSpecial100Idxes.m b/functions/getSpecial100Idxes.m new file mode 100644 index 0000000..55f29ea --- /dev/null +++ b/functions/getSpecial100Idxes.m @@ -0,0 +1,27 @@ +%% Finds trial indices corresponding to each special 100 image, by path name +% Returns a 100 x 2 cell array. First column = name of special100 file, second column = trial indices +% +function special100Idxes = getSpecial100Idxes(stimFiles) + + if exist('data/stims/special100Names.txt', 'file') % read stim names directly if they've been written + stimNames = readcell('data/stims/special100Names.txt', 'Delimiter', ','); + assert(length(stimNames) == 100, 'Error: special100Names did not contain 100 entries'); + else % get names from image set and save as output + % Get names of stimuli + stimPaths = glob('data/stims/special100/*.png'); + assert(length(stimPaths) == 100, 'Error: Did not find 100 special100 files'); + stimNames = extractBetween(stimPaths, 'special100/', '.png'); + writecell(stimNames, 'data/stims/special100Names.txt'); + end + + % find indices of each special100 image + idxes = cell(100, 1); + for ii = 1:100 + idxes{ii} = find(contains(stimFiles, stimNames{ii})); + assert(~isempty(idxes{ii}), 'Error: No trial index found for special100 image %d\n', ii); + end + + % Add name and indices + special100Idxes = [stimNames, idxes]; + +end \ No newline at end of file diff --git a/functions/getXyzsInf.m b/functions/getXyzsInf.m new file mode 100644 index 0000000..075c7f7 --- /dev/null +++ b/functions/getXyzsInf.m @@ -0,0 +1,65 @@ +%% This function returns inflated gifti XYZ positions from an electrodes.tsv table +% Adapted from DH ieeg_snap2inflate.m to perform separately for each hemisphere +% +% xyzsInf = getXyzsInf(electrodes, hemi, gii, giiInf); +% xyzsInf = getXyzsInf(electrodes, hemi, gii, giiInf, distThresh); +% electrodes = nx_ table, read from an electrodes.tsv file. Must contain column "hemisphere" (str, 'R' or 'L'). +% Must also contain columns "Destrieux_label" (numerical) and "Destrieux_label_text" if is given as true. +% hemi = char, 'r' or 'l', corresponding to the hemisphere of the pial and inflated giftis +% gii = gifti object, pial surface of one hemisphere +% giiInf = gifti object, inflated surface of one hemisphere +% distThresh = double (optional), Only electrodes within mm of a pial vertex are kept in output. +% Default = 5. +% checks = logical (optional), whether to more strictly check that each electrode is cortical and not WM/thalamic/outside brain. Default == false +% +% Returns: +% xyzsInf = nx3 double, XYZ positions of gray matter electrodes in inflated gifti space. Electrodes in +% white matter, contralateral hemisphere, or beyond are returned as [NaN, NaN, NaN] rows. +% +% HH 2021 +% +function xyzsInf = getXyzsInf(electrodes, hemi, gii, giiInf, distThresh, checks) + + if nargin < 6, checks = false; end % to perform checks based on Destrieux atlas names + if nargin < 5 || isempty(distThresh), distThresh = 5; end % maximum distance allowed between valid electrode and gifti vertex + + hemi = lower(hemi); + assert(ismember(hemi, {'r', 'l'}), "hemi must be given as 'r' or 'l'"); + assert(size(gii.vertices, 1) == size(giiInf.vertices, 1), 'gii and giiInf do not have the same number of vertices'); + + xyzs = [electrodes.x, electrodes.y, electrodes.z]; + hemiLabs = lower(electrodes.hemisphere); + + if checks + anatLabs = electrodes.Destrieux_label; % numerical labels + anatText = electrodes.Destrieux_label_text; % char labels + else + anatLabs = zeros(height(electrodes), 1); + anatText = cell(height(electrodes), 1); + end + + xyzsInf = nan(size(xyzs)); + for ii = 1:size(xyzs, 1) + + if ~strcmp(hemiLabs(ii), hemi) % wrong hemisphere, continue + continue + end + + if checks && (isnan(anatLabs(ii)) || ismember(anatLabs(ii), [0, 41, 49, 10])) % non existent, outside brain, WM, or L/R thalamic + continue + end + + if checks && (~startsWith(anatText{ii}, sprintf('%sh', hemi))) % must start with lh or rh to be a cortical site + continue + end + + dists = vecnorm(gii.vertices - xyzs(ii, :), 2, 2); % distance from each vertex in gifti to current electrode + [minDist, idx] = min(dists); + + if minDist <= distThresh + xyzsInf(ii, :) = giiInf.vertices(idx, :); + end + + end + +end \ No newline at end of file diff --git a/functions/goodChsNRuns.m b/functions/goodChsNRuns.m new file mode 100644 index 0000000..ec40121 --- /dev/null +++ b/functions/goodChsNRuns.m @@ -0,0 +1,33 @@ +% This function returns which channels are good across at least n runs, where the input is bad chs for each run +% +% goodChs = goodChsNRuns(badChs, numChs, n); +% [goodChs, badChs] = goodChsNRuns(badChs, numChs, n); +% +% INPUTS +% badChsRuns = 1xn numerical cell array. bad channels in each run +% numChs = num. Total number of channels +% n = num. Minimum number of runs for a channel to be good in +% +% OUTPUTS +% goodChsOnN = 1xp num. Indices of good channels across >= n runs +% badChsOnN = 1xq num. Indices of bad channels (complement of goodChs on 1:numChs) +% +function [goodChsOnN, badChsOnN] = goodChsNRuns(badChsRuns, numChs, n) + + numRuns = length(badChsRuns); + assert(n >= 1 && n <= numRuns, 'Error: n must be specified between 1 and %d runs', numRuns); + + % chs x runs. 1 indicates a bad channel, 0 indicates a good channel + badChsLogical = zeros(numChs, numRuns, 'logical'); + for ii = 1:numRuns + badChsLogical(badChsRuns{ii}, ii) = true; + end + + % find which rows have at least n 0s + goodChsLogical = sum(~badChsLogical, 2) >= n; + + % convert to numerical indices, and also return complement as bad chs + goodChsOnN = find(goodChsLogical); + badChsOnN = find(~goodChsLogical); + +end \ No newline at end of file diff --git a/functions/mnl_rsq.m b/functions/mnl_rsq.m new file mode 100644 index 0000000..726217b --- /dev/null +++ b/functions/mnl_rsq.m @@ -0,0 +1,31 @@ +%% Computes the R-squared value for two one-dimensional distributions +% Taken by HH from legacy DH/KJM code, 04/2021 +function [rsq, p] = mnl_rsq(x1, x2) + + x1=double(x1(:)); + x2=double(x2(:)); + + sum1=sum(x1); + sum2=sum(x2); + n1=length(x1); + n2=length(x2); + sumsqu1=sum(x1.*x1); + sumsqu2=sum(x2.*x2); + + G=((sum1+sum2)^2)/(n1+n2); + + % /var(X,Y) + rsq = (sum1^2/n1 + sum2^2/n2 - G)/(sumsqu1 + sumsqu2 - G); + if mean(x2) > mean(x1) % inverse sign if second vector larger than 1st vector + rsq = -rsq; + end + + %% Calculate p-value from F-test + + SStot = (length(x1) + length(x2))*var([x1; x2], 1); % total sum of squares (var(_, 1) normalized by N) + SSres = length(x1)*var(x1, 1) + length(x2)*var(x2, 1); % residual sum of squares using respective means + F = (SStot - SSres)/(SSres/(n1+n2-2)); % F statistic using 1 df numerator and N-2 df denominator + + p = 1 - fcdf(F, 1, n1+n2-2); % p-value + +end \ No newline at end of file diff --git a/functions/plotActivationsToGiftis.m b/functions/plotActivationsToGiftis.m new file mode 100644 index 0000000..931d82e --- /dev/null +++ b/functions/plotActivationsToGiftis.m @@ -0,0 +1,68 @@ +%% Plot activations on gifti and inflated gifti, at lateral, medial, and inferior angles +% +% INPUTS +% elecs = electrodes table of real electrodes, or give empty if you don't wish to plot +% elecsAct = "electrodes table" format for the positions activations are at (i.e., bipolar-coordinate electrodes table) +% giis = pial giftis, matching hemis +% giiInfs = inflated pial giftis, matching hemis +% sulcs = sulcus objects, matching hemis +% hemis = cell array of hemispheres in subject +% activations = vector of activations that matches electrodes +% wm = weight max +% outdir = output directory to save plots +% activationStr = string of the type of activation being plotted (e.g., BBrsq) +% visi = Visibility of plots as they appear +% +function plotActivationsToGiftis(elecs, elecsAct, activations, giis, giiInfs, sulcs, hemis, wm, outdir, activationStr, visi) + + % View giftis and inflated gifti from 3 different angle. Configure theta and phi for each + angleText = {'lateral', 'medial', 'inferior'}; + + for ii = 1:length(hemis) + + % theta and phi, configured for each hemisphere + ths = [choose(strcmpi(hemis{ii}, 'r'), 90, -90), choose(strcmpi(hemis{ii}, 'r'), -90, 90), choose(strcmpi(hemis{ii}, 'r'), 90, -90)]; + phis = [0, 0, -90]; + + elecsActHemi = elecsAct(strcmpi(hemis{ii}, elecsAct.hemisphere), :); + activationsHemi = activations(strcmpi(hemis{ii}, elecsAct.hemisphere)); + + % normal gifti + figure('Position', [200, 200, 600, 400], 'Visible', visi); ieeg_RenderGifti(giis{ii}); alpha 0.3 + hold on + + % isolate the electrodes and activation values in current hemi and plot + if ~isempty(elecs) + elecsHemi = elecs(strcmpi(hemis{ii}, elecs.hemisphere), :); + plot3(elecsHemi.x, elecsHemi.y, elecsHemi.z, 'o', 'MarkerSize', 4, 'Color', 'k', 'MarkerFaceColor', 'w'); + end + + % plot weighted activations on gifti + plotCortexWeightsAdjustable([elecsActHemi.x, elecsActHemi.y, elecsActHemi.z], activationsHemi, 'SigFraction', 0.05, 'MaxWt', wm, 'SizeJump', 2); + hold off + + % rotate camera to each of 3 different views + for jj = 1:3 + ieeg_viewLight(ths(jj), phis(jj)); + saveas(gcf, fullfile(outdir, sprintf('gifti%s_%s_%s.png', upper(hemis{ii}), activationStr, angleText{jj}))); + end + close(gcf); + + + % for inflated gifti, a new figure needs to be made for each angle, because opaque and els_popout + for jj = 1:3 + xyzInf = els_popout([elecsActHemi.xInf, elecsActHemi.yInf, elecsActHemi.zInf], ths(jj), phis(jj), 0.02); + figure('Position', [200, 200, 600, 400], 'Visible', visi); ieeg_RenderGifti(giiInfs{ii}, sulcs{ii}); + hold on + + % don't plot actual electrodes themselves on inflated brain (too busy) + + plotCortexWeightsAdjustable(xyzInf, activationsHemi, 'SigFraction', 0.05, 'MaxWt', wm, 'SizeJump', 2); + hold off + ieeg_viewLight(ths(jj), phis(jj)); + saveas(gcf, fullfile(outdir, sprintf('giftiInf%s_%s_%s.png', upper(hemis{ii}), activationStr, angleText{jj}))); + close(gcf); + end + + end +end \ No newline at end of file diff --git a/main.m b/main.m new file mode 100644 index 0000000..d8486c0 --- /dev/null +++ b/main.m @@ -0,0 +1,40 @@ +%% Outline of which scripts to run for formal NSD manuscript analysis + +% subjectConfig.m configures subject-specific inf + +% i: Preprocessing mef data +preprocNSDMef + +%% Part 1: stationary maps + +% Perform analyses using scalar broadband values +script01_analyzePreprocScalarBB + +% % Merge across subjects +% script02_mniAcrossSubjects +% +% %% Part 2: correlate NSD fMRI with broadbabnd +% +% % Save beta weights for each of the 8 fMRI subjects for the shared1000 images +% saveFmriShared1000 +% +% % Save mean z-scored beta weights and save together for all subjects +% script03_meanShared1000Betas +% +% % calculates the average vertex-wide activation for each subject and across subjects +% script03b_fmriExpected.m +% +% % correlation maps between NSD subjects and fMRI +% script04_correlateBBFmri +% +% % expected correlations between maximum fMRI vertex and rest of cortex, for each iEEG electrode +% script05_correlateFmriVert2Cort +% +% %% Abstract for VSS +% +% % Calculate stats for all electrodes BBsnr -- how many significant, anatomical locations +% scriptAbstract1_statsBBncsnrAllSubs.m +% +% % For the different fMRI patterns, look at which NSD shared1000 stimuli are preferred +% scriptAbstract2_patternsAnalysis.m + diff --git a/preprocNSDMef.m b/preprocNSDMef.m new file mode 100644 index 0000000..cfcfca0 --- /dev/null +++ b/preprocNSDMef.m @@ -0,0 +1,175 @@ +%% 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 new file mode 100644 index 0000000..bd7258d --- /dev/null +++ b/script01_analyzePreprocScalarBB.m @@ -0,0 +1,337 @@ +%% 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'); + + +