diff --git a/data/metadata/sub-X_sessions.txt b/data/metadata/sub-X_sessions.txt deleted file mode 100644 index c9aa0aa..0000000 --- a/data/metadata/sub-X_sessions.txt +++ /dev/null @@ -1,11 +0,0 @@ -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/functions/ieeg_nsdParseEvents.m b/functions/ieeg_nsdParseEvents.m new file mode 100644 index 0000000..80e4029 --- /dev/null +++ b/functions/ieeg_nsdParseEvents.m @@ -0,0 +1,70 @@ + +function [events_status,nsd_idx,shared_idx,nsd_repeats] = ieeg_nsdParseEvents(all_events) +% +% function [events_status,nsd_idx,shared_idx,nsd_repeats] = ieeg_nsdParseEvents(all_events) +% +% Input +% all_events: BIDS events table +% +% Output +% events_status 2 columns, +% column 1: 1 if the event is bad +% column 2: 1 for low interictal, 2 for high interictal +% +% nsd_idx: 1530x1 column vector containing the nsd index numbers found +% from the original file names for all shared images +% ('shared0814_nsd59080_prepped1.png' --> 59080). NaN: probe images. + +% shared_idx: 1530x1 column vector containing the shared index numbers +% found from the original file names for all shared images +% ('shared0814_nsd59080_prepped1.png' --> 814). NaN: probe images. + +% nsd repeats: 0: nonrepeats, 1-6: repeats, NaN: probe +% +% +% DH, MM, Multimodal Neuroimaging Lab 2023 + + +% events status (0:good, 1:bad) and interictal (0:no, 1:low, 2:high) +events_status = zeros(height(all_events),2); + +% events status (0:good, 1:bad) +events_status(ismember(all_events.status,'bad'),1) = 1; + +% find interictal label (0:no, 1:low, 2:high) from events_table status_description +for ii_evts = 1:height(all_events) + if ~isempty(all_events.status_description{ii_evts}) + this_set_of_descriptions = strip(split(all_events.status_description{ii_evts},',')); + if sum(ismember(this_set_of_descriptions,'/Interictal-findings/Attribute/Categorical/Categorical-level/Low'))>0 + events_status(ii_evts,2) = 1; + elseif sum(ismember(this_set_of_descriptions,'/Interictal-findings/Attribute/Categorical/Categorical-level/High'))>0 + events_status(ii_evts,2) = 2; + end + end +end + +% get NSD idx +nsd_idx = NaN(height(all_events),1); +shared_idx = NaN(height(all_events),1); +stim_image_file = all_events.stim_file; +temp = split(stim_image_file, '\'); +temp = temp(:,end); % last split is filename +for ii_evts = 1:height(all_events) + if ~isequal(temp{ii_evts}(1),'z') % note a probe trial + temp_2 = extractBetween(temp{ii_evts},'nsd','_prepped'); + nsd_idx(ii_evts) = str2double(temp_2{1}); + temp_3 = extractBetween(temp{ii_evts},'shared','_nsd'); + shared_idx(ii_evts) = str2double(temp_3{1}); + end +end +clear temp_2 temp_3 + +% code repeats, 0: nonrepeats, 1-6: repeats, NaN: probe +nsd_repeats = NaN(height(all_events),1); +for ii_evts = 1:height(all_events) + if length(find(nsd_idx == nsd_idx(ii_evts)))>1 % repeat + nsd_repeats(ii_evts) = length(find(nsd_idx(1:ii_evts) == nsd_idx(ii_evts))); % how many times has it occured before ii_evts + elseif length(find(nsd_idx == nsd_idx(ii_evts)))==1 % nonrepeat + nsd_repeats(ii_evts) = 0; + end +end diff --git a/script01_preprocNSDMef.m b/script01_preprocNSDMef.m new file mode 100644 index 0000000..d92b542 --- /dev/null +++ b/script01_preprocNSDMef.m @@ -0,0 +1,350 @@ +%% This is the preprocessing script for the clean NSD analysis. +% Requires each subject to have annotated channels, events, and electrodes +% with SOZ annotated + +%% Set paths, get filenames and tables for one subject +% generates data_info with information for all runs +% generates all_channels with good channels across runs marked + +localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) +addpath('functions'); + +subjects = {'01','02','03','04','05','06','07'}; + +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',... + 'NSDspecial06','NSDspecial07','NSDspecial08','NSDspecial09','NSDspecial10'}; +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); + +% For all possible runs, set data_info +data_info = []; +set_counter = 0; +for ii = 1:height(task_list.filename) + % check if task is an accepted task label + if ismember(extractBetween(task_list.filename{ii},'task-','_run-'),task_labels) + set_counter = set_counter+1; + % get core of the iEEG filename + coreName = extractBetween(task_list.filename{ii},'ieeg/','_ieeg.mefd'); + % data iEEG filename + data_info(set_counter).filename = [coreName{1} '_ieeg.mefd']; + % events name + data_info(set_counter).eventsname = [coreName{1} '_events.tsv']; + % channels name + data_info(set_counter).channelsname = [coreName{1} '_channels.tsv']; + % track good channels throughout all runs + channels_table = ieeg_readtableRmHyphens(fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', data_info(set_counter).channelsname)); + if ii==1 + good_channel_bool = ismember(channels_table.type,{'SEEG','ECOG'}) & ismember(channels_table.status,'good'); + else + good_channel_bool = ismember(channels_table.type,{'SEEG','ECOG'}) & ismember(channels_table.status,'good') & good_channel_bool==1; + end + data_info(set_counter).general = task_list(ii,:); + end +end + +% Channel info across runs +all_channels.name = channels_table.name; +all_channels.type = channels_table.type; +all_channels.status = good_channel_bool; % 1 is good, zero is bad, -1 is SOZ +% SOZ channels +elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); % first set SOZ channels to bad +all_channels.soz = ismember(all_channels.name, elecsSOZ); + + +%% Loop through NSD01 - NSD10 and all runs individually, concatenate output at the end + +Mdata = []; % data matrix +Mbb = []; % data matrix for broadband filtered data +eventsST = []; % events + +for ii = 1:length(data_info) + fprintf('Loading %s\n', data_info(ii).filename); + + %%% ------------------------------------------------------------%%% + %%% Load events + try + events = readtable(fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ... + data_info(ii).eventsname), 'FileType', 'text', 'Delimiter', '\t', 'TreatAsEmpty', 'n/a'); + catch + warning('Error loading events file for %02d, skipping', ii); % e.g., no 2nd run for a certain NSD session + continue + end + + % make sure status description is a cell array to concatenate + if ~iscell(events.status_description) + events.status_description = cellstr(string(events.status_description)); + end + + %%% ------------------------------------------------------------%%% + %%% Load mef data + [metadata, data] = readMef3(fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', data_info(ii).filename)); + srate = metadata.time_series_metadata.section_2.sampling_frequency; + + %%% ------------------------------------------------------------%%% + %%% Reorganize and clean up events so that adjacent rest and stim events are combined + samprange = -2*srate:2*srate-1; % -2:+2 seconds + tt = samprange/srate; + + % 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. + idxST = find(ismember(events.trial_type, {'stim', 'test'})); + eventsSTCurr = events(idxST, :); % only keep stim and test events + % add pre (rest) status and durations + assert(all(strcmp(events.trial_type(idxST-1), 'rest')), 'Error: trials preceding some stim/test trials are not rest'); + eventsSTCurr.pre_duration = events.duration(idxST-1, :); + eventsSTCurr.pre_status = events.status(idxST-1, :); + eventsSTCurr.pre_status_description = events.status_description(idxST-1, :); + + % add task (NSD1-10), run (1 2...), session (ordered blocks of time) info + eventsSTCurr.tasknumber = repmat(data_info(ii).general.nsd_special_nr, height(eventsSTCurr), 1); + eventsSTCurr.run = repmat(data_info(ii).general.run, height(eventsSTCurr), 1); + eventsSTCurr.recording_set = repmat(data_info(ii).general.recording_set, height(eventsSTCurr), 1); + + % remove any events missing the full sample range (e.g. due to crash) + eventsSTCurr(eventsSTCurr.sample_start+samprange(end)+1 > size(data, 2), :) = []; + + %%% ------------------------------------------------------------%%% + %%% Preprocess: highpass and CAR by lowest variance + + % highpass the SEEG channels + data(strcmp(all_channels.type, 'SEEG'), :) = ieeg_highpass(data(strcmp(all_channels.type, 'SEEG'), :)', srate)'; + + % re-reference: CAR by lowest variance (across data from first to last events) + opts = struct(); + opts.vartype = 'var'; + opts.pctThresh = 25; % leave at default of 25% + opts.winResp = [eventsSTCurr.onset(1) eventsSTCurr.onset(end)]; + tt_all = (1:size(data,2))/srate; + badChs = find(all_channels.status==0 | all_channels.soz==1); % bad channels or soz + [car_tmp, chsUsed] = ccep_CARVariance(tt_all, data, srate, badChs, [], opts); + % only add CAR-reref data for sEEG channels + data_car = data; + data_car(strcmp(all_channels.type, 'SEEG'), :) = car_tmp(strcmp(all_channels.type, 'SEEG'), :); + clear car_tmp + + %%% ------------------------------------------------------------%%% + %%% Epoch data + M = zeros(size(data_car, 1), length(samprange), height(eventsSTCurr)); + for jj = 1:height(eventsSTCurr) + M(:, :, jj) = data_car(:, (round(eventsSTCurr.onset(jj)*srate)+samprange(1)+1) : (round(eventsSTCurr.onset(jj)*srate)+samprange(end)+1)); % convert to 1 indexing + end + + % after epoching - check the DI1 channel to make sure that the first sample of 1 occurs at t=0 + %di1 = squeeze(M(strcmp(all_channels.name, 'DigitalInput1'), :, :)); + %figure; plot(tt, di1); xlim([-0.01, 0.01]); + + %%% ------------------------------------------------------------%%% + %%% Extract broadband and epoch data + bands = [70 80;80 90; 90 100; 100 110; 130 140; 140 150; 150 160; 160 170]; % 10 hz bins, avoiding 60 and 120 + bb_power = data_car; % initizalize with all channels + bb_power(strcmp(all_channels.type, 'SEEG'), :) = ieeg_getHilbert(data_car(strcmp(all_channels.type, 'SEEG'), :)', bands, srate, 'power')'; + % epoch broadband + BB = zeros(size(data_car, 1), length(samprange), height(eventsSTCurr)); + for jj = 1:height(eventsSTCurr) + BB(:,:,jj) = bb_power(:, (round(eventsSTCurr.onset(jj)*srate)+samprange(1)+1) : (round(eventsSTCurr.onset(jj)*srate)+samprange(end)+1)); % convert to 1 indexing + end + + %%% ------------------------------------------------------------%%% + %%% Concatenate data and events + Mdata = cat(3, Mdata, M); + Mbb = cat(3, Mbb, BB); + eventsST = [eventsST; eventsSTCurr]; + +end + +% save as single to save space +Mbb = single(Mbb); +Mdata = single(Mdata); + +% Save all outputs +outdir = fullfile(localDataPath.output,'derivatives','preproc_car',['sub-' sub_label]); +if ~exist(outdir,'dir') + sprintf(['making output directory: ' outdir]) + mkdir(outdir); +end +save(fullfile(outdir, ['sub-' sub_label '_desc-preprocCAR_ieeg.mat']), 'tt', 'srate', 'Mdata', 'eventsST', 'all_channels', '-v7.3'); +save(fullfile(outdir, ['sub-' sub_label '_desc-preprocCARBB_ieeg.mat']), 'tt', 'srate', 'Mbb', 'eventsST', 'all_channels', '-v7.3'); +writetable(eventsST, fullfile(outdir, ['sub-' sub_label '_desc-preprocCAR_events.tsv']), 'FileType', 'text', 'Delimiter', '\t'); + + +%% Load BB output to check + +clear all +localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) +addpath('functions'); + +subjects = {'01','02','03','04','05','06'}; + +ss = 5; +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'); + +%% Normalize bb power per run + +% Initialize normalized log power of BB +Mbb_norm = log10(Mbb); + +% Indicate the interval for baseline, used in normalization +norm_int = find(tt>-.2 & tt<0); + +% Normalize per run +for run_idx = 1:max(eventsST.tasknumber) + this_run = find(eventsST.tasknumber==run_idx); % out of 1500 + + % find pre-stim events with 'good' status + trials_norm = find(ismember(eventsST.pre_status,'good') & eventsST.tasknumber==run_idx); + + Mbb_norm(:,:,this_run) = minus(Mbb_norm(:,:,this_run),mean(Mbb_norm(:,norm_int,trials_norm),[2 3],'omitnan')); +end + +%% Find repeated images, calculate SNR + +eventsST.status_description = cellstr(string(eventsST.status_description)); +[events_status,nsd_idx,shared_idx,nsd_repeats] = ieeg_nsdParseEvents(eventsST); + +all_chan_snr = NaN(size(Mbb_norm,1),1); +t_avg = tt>0.1 & tt<.5; + +for el_nr = 1:size(Mbb_norm,1) + + if ismember(all_channels.type(el_nr),'SEEG') && all_channels.status(el_nr)==1 + bb_strength = squeeze(mean(Mbb_norm(el_nr,t_avg==1,:),2)); + + all_repeats = find(nsd_repeats>0); + shared_idx_repeats = unique(shared_idx(all_repeats)); % 100 images + repeats_bb_strength = cell(length(shared_idx_repeats),1); + for kk = 1:length(shared_idx_repeats) + these_trials = find(shared_idx==shared_idx_repeats(kk)); % for this repeat, find the correct 6 trial numbers out of the 1500 and get the image and the data + repeats_bb_strength{kk} = bb_strength(these_trials); + end + + [NCSNR, p, NCSNRNull] = estimateNCSNR(repeats_bb_strength, 1000); + all_chan_snr(el_nr) = NCSNR; + end +end + +%% imagesc broadband and render prelim + +good_channel_nrs = find(all_channels.status==1); + +figure +imagesc(tt,1:length(good_channel_nrs),mean(Mbb_norm(good_channel_nrs,:,:),3),[-.2 .2]) +set(gca,'YTick',1:length(good_channel_nrs),'YTickLabels',all_channels.name(good_channel_nrs)) + +%% 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); + +name = all_channels.name; +all_channels_table = table(name); +elecs = ieeg_sortElectrodes(elecs, all_channels_table, 0); + +% load pial and inflated giftis +gL = gifti(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],['white.L.surf.gii'])); +gR = gifti(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],['white.R.surf.gii'])); +gL_infl = gifti(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],['inflated.L.surf.gii'])); +gR_infl = gifti(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],['inflated.R.surf.gii'])); + +% snap electrodes to surface and then move to inflated +xyz_inflated = ieeg_snap2inflated(elecs,gR,gL,gR_infl,gL_infl,4); + +% render with Wang labeling +Wang_ROI_Names = {... + 'V1v' 'V1d' 'V2v' 'V2d' 'V3v' 'V3d' 'hV4' 'VO1' 'VO2' 'PHC1' 'PHC2' ... + 'TO2' 'TO1' 'LO2' 'LO1' 'V3B' 'V3A' 'IPS0' 'IPS1' 'IPS2' 'IPS3' 'IPS4' ... + 'IPS5' 'SPL1' 'FEF'}; + +for hh = 1:2 + + if hh==1 + hemi = 'l'; + g = gL_infl; + views_plot = {[40,-30],[-45,-10],[-90,20]}; + elseif hh==2 + hemi = 'r'; + g = gR_infl; + views_plot = {[-40,-30],[45,-10],[90,20]}; + end + + % surface labels + surface_labels = MRIread(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],'surf',... + [hemi 'h.wang15_mplbl.mgz'])); + vert_label = surface_labels.vol(:); + + % sulcal labels + sulcal_labels = read_curv(fullfile(localDataPath.input,'sourcedata','freesurfer',['sub-' sub_label],'surf',... + [hemi 'h.sulc'])); + + % load the color map + cmap = make_WangColormap(); + + electrodes_thisHemi = find(ismember(elecs.hemisphere,upper(hemi))); + + % make a plot with electrode dots + for vv = 1:length(views_plot) + v_d = [views_plot{vv}(1),views_plot{vv}(2)]; + + % get the inflated coordinates + els = xyz_inflated; + % calculate popout so we can better read labels and see amplitude + a_offset = .1*max(abs(els(:,1)))*[cosd(v_d(1)-90)*cosd(v_d(2)) sind(v_d(1)-90)*cosd(v_d(2)) sind(v_d(2))]; + els_pop = els+repmat(a_offset,size(els,1),1); + + % no labels + figure + tH = ieeg_RenderGiftiLabels(g,vert_label,cmap,Wang_ROI_Names,sulcal_labels); + ieeg_elAdd(els(electrodes_thisHemi,:),[.99 .99 .99],25) % add electrode positions + ieeg_elAdd(els(electrodes_thisHemi,:),[.1 .1 .1],15) % add electrode positions + ieeg_viewLight(v_d(1),v_d(2)) % change viewing angle + set(gcf,'PaperPositionMode','auto') + print('-dpng','-r300',fullfile(localDataPath.input,'derivatives','render',['sub-' sub_label],... + ['inflated_dots_sub-' sub_label '_WangAreas_v' int2str(v_d(1)) '_' int2str(v_d(2)) '_' hemi])) + + % with labels + figure + tH = ieeg_RenderGiftiLabels(g,vert_label,cmap,Wang_ROI_Names,sulcal_labels); + ieeg_elAdd(els_pop(electrodes_thisHemi,:),[.1 .1 .1],15) % add electrode positions + ieeg_label(els_pop(electrodes_thisHemi,:),20,6,elecs.name(electrodes_thisHemi)) % add electrode names + ieeg_viewLight(v_d(1),v_d(2)) % change viewing angle + set(gcf,'PaperPositionMode','auto') + print('-dpng','-r300',fullfile(localDataPath.input,'derivatives','render',['sub-' sub_label],... + ['inflated_labels_sub-' sub_label '_WangAreas_v' int2str(v_d(1)) '_' int2str(v_d(2)) '_' hemi])) + +% % with activity +% figure +% tH = ieeg_RenderGiftiLabels(g,vert_label,cmap,Wang_ROI_Names,sulcal_labels); +% all_chan_snr_plot = all_chan_snr; +% all_chan_snr_plot(all_chan_snr_plot<.2) = 0; +% ieeg_elAdd_sizable(els_pop(electrodes_thisHemi,:),all_chan_snr_plot(electrodes_thisHemi),.8,40) % add electrode positions +% ieeg_viewLight(v_d(1),v_d(2)) % change viewing angle +% set(gcf,'PaperPositionMode','auto') +% print('-dpng','-r300',fullfile(localDataPath.input,'derivatives','render',['sub-' sub_label],... +% ['NCSNR_sub-' sub_label '_WangAreas_v' int2str(v_d(1)) '_' int2str(v_d(2)) '_' hemi])) + end + close all + +end + + + diff --git a/script02_writeMNI.m b/script02_writeMNI.m new file mode 100644 index 0000000..203345e --- /dev/null +++ b/script02_writeMNI.m @@ -0,0 +1,83 @@ + +clear all +localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) +addpath('functions'); + +subjects = {'01','02','03','04','05','06'}; + +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'); + +%% need to get MNI positions first + +% mni_coords = ieeg_mni305ThroughFsSphere(elecmatrix,hemi,FSdir,FSsubjectsdir) + + +%% Normalize bb power per run + +% Initialize normalized log power of BB +Mbb_norm = log10(Mbb); + +% Indicate the interval for baseline, used in normalization +norm_int = find(tt>-.2 & tt<0); + +% Normalize per run +for run_idx = 1:max(eventsST.tasknumber) + this_run = find(eventsST.tasknumber==run_idx); % out of 1500 + + % find pre-stim events with 'good' status + trials_norm = find(ismember(eventsST.pre_status,'good') & eventsST.tasknumber==run_idx); + + Mbb_norm(:,:,this_run) = minus(Mbb_norm(:,:,this_run),mean(Mbb_norm(:,norm_int,trials_norm),[2 3],'omitnan')); +end + +%% Find repeated images, calculate SNR + +eventsST.status_description = cellstr(string(eventsST.status_description)); +[events_status,nsd_idx,shared_idx,nsd_repeats] = ieeg_nsdParseEvents(eventsST); + +all_chan_snr = NaN(size(Mbb_norm,1),1); +t_avg = tt>0.1 & tt<.5; + +for el_nr = 1:size(Mbb_norm,1) + + if ismember(all_channels.type(el_nr),'SEEG') && all_channels.status(el_nr)==1 + bb_strength = squeeze(mean(Mbb_norm(el_nr,t_avg==1,:),2)); + + all_repeats = find(nsd_repeats>0); + shared_idx_repeats = unique(shared_idx(all_repeats)); % 100 images + repeats_bb_strength = cell(length(shared_idx_repeats),1); + for kk = 1:length(shared_idx_repeats) + these_trials = find(shared_idx==shared_idx_repeats(kk)); % for this repeat, find the correct 6 trial numbers out of the 1500 and get the image and the data + repeats_bb_strength{kk} = bb_strength(these_trials); + end + + [NCSNR, p, NCSNRNull] = estimateNCSNR(repeats_bb_strength, 1000); + all_chan_snr(el_nr) = NCSNR; + end +end + +%% 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); + +name = all_channels.name; +all_channels_table = table(name); +elecs = ieeg_sortElectrodes(elecs, all_channels_table, 0); + +% load pial and inflated giftis +gL = gifti(fullfile(localDataPath.input,'derivatives','freesurfer',['sub-' sub_label],['white.L.surf.gii'])); +gR = gifti(fullfile(localDataPath.input,'derivatives','freesurfer',['sub-' sub_label],['white.R.surf.gii'])); +gL_infl = gifti(fullfile(localDataPath.input,'derivatives','freesurfer',['sub-' sub_label],['inflated.L.surf.gii'])); +gR_infl = gifti(fullfile(localDataPath.input,'derivatives','freesurfer',['sub-' sub_label],['inflated.R.surf.gii'])); + +% snap electrodes to surface and then move to inflated +xyz_inflated = ieeg_snap2inflated(elecs,gR,gL,gR_infl,gL_infl,4); \ No newline at end of file diff --git a/script03_analyzePreprocScalarBB.m b/script03_analyzePreprocScalarBB.m new file mode 100644 index 0000000..bd7258d --- /dev/null +++ b/script03_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'); + + + diff --git a/setLocalDataPath.m b/setLocalDataPath.m index 41296ad..4b9b0c3 100644 --- a/setLocalDataPath.m +++ b/setLocalDataPath.m @@ -1,6 +1,6 @@ function localDataPath = setLocalDataPath(varargin) % function LocalDataPath = setLocalDataPath(varargin) -% Return the path to the root CCEP directory and add paths in this repo +% Return the path to the root directory and add paths in this repo % % input: % personalDataPath: optional, set to 1 if adding personalDataPath @@ -21,13 +21,13 @@ if isempty(varargin) rootPath = which('setLocalDataPath'); - ccepRepoPath = fileparts(rootPath); + RepoPath = fileparts(rootPath); % add path to functions - addpath(genpath(ccepRepoPath)); + addpath(genpath(RepoPath)); % add localDataPath default - localDataPath = fullfile(ccepRepoPath,'data'); + localDataPath = fullfile(RepoPath,'data'); elseif ~isempty(varargin) % add path to data if varargin{1}==1 && exist('personalDataPath','file') @@ -46,8 +46,8 @@ % add path to functions rootPath = which('setLocalDataPath'); - ccepRepoPath = fileparts(rootPath); - addpath(genpath(ccepRepoPath)); + RepoPath = fileparts(rootPath); + addpath(genpath(RepoPath)); end return