diff --git a/functions/getXyzsInf.m b/functions/getXyzsInf.m deleted file mode 100644 index 075c7f7..0000000 --- a/functions/getXyzsInf.m +++ /dev/null @@ -1,65 +0,0 @@ -%% 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/script01_preprocNSDMef.m b/script01_preprocNSDMef.m index 553b28e..e12b73c 100644 --- a/script01_preprocNSDMef.m +++ b/script01_preprocNSDMef.m @@ -10,7 +10,7 @@ addpath('functions'); % subject to preprocess -ss = 17; +ss = 3; sub_label = sprintf('%02d', ss); ses_label = 'ieeg01'; @@ -52,10 +52,10 @@ % 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 +all_channels.status = good_channel_bool; % 1 is good, zero is bad % SOZ channels elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); % first set SOZ channels to bad -all_channels.soz = ismember(all_channels.name, elecsSOZ); +all_channels.soz = ismember(all_channels.name, elecsSOZ); % 1 is SOZ, 0 is no SOZ %% Loop through NSD01 - NSD10 and all runs individually, concatenate output at the end diff --git a/script01b_preprocBipolarBroadband.m b/script01b_preprocBipolarBroadband.m new file mode 100644 index 0000000..4c9bdbb --- /dev/null +++ b/script01b_preprocBipolarBroadband.m @@ -0,0 +1,63 @@ +%% Load CAR trial data, bipolar-rereference, and calculate broadband to save + +localDataPath = personalDataPath(); + +ss = 1; +% for ss = 13:17 +% if ss == 4, continue; end + +sub = sprintf('sub-%02d', ss); + +%% 1) Load CAR data, bipolar-rereference + +fprintf('Loading CAR evoked potentials for %s\n', sub); +load(fullfile(localDataPath.input, 'derivatives', 'preproc_car', sub, sprintf('%s_desc-preprocCAR_ieeg.mat', sub))); + +% load segmented data to get segs for this subject +participants = readtable(fullfile(localDataPath.input, 'participants.tsv'), 'Filetype', 'text', 'Delimiter', '\t'); +seg5 = participants.seg5{strcmp(participants.participant_id, sub)}; +seg5 = strip(split(seg5, ',')); +seg6 = participants.seg6{strcmp(participants.participant_id, sub)}; +seg6 = strip(split(seg6, ',')); +if strcmp(seg5, 'n/a'), seg5 = {}; end +if strcmp(seg6, 'n/a'), seg6 = {}; end + +% Get SEEG channels +ephysIdxes = strcmp(all_channels.type, 'SEEG'); +MdataEphys = Mdata(ephysIdxes, :, :); +nEphys = sum(ephysIdxes); +chNamesEphys = all_channels.name(ephysIdxes); + +% First run twice, to get how many bipolar channels, and to get bipolar bad channel numbers, then bipolar SOZ channels +fprintf('Performing bipolar rereference per trial\n'); +[~, bipolarNames, bipolarChans, badChansBip] = ieeg_bipolarSEEG(MdataEphys(:, :, 1)', chNamesEphys, find(all_channels.status(ephysIdxes) == 0), seg5, seg6); +[~, ~, ~, sozBip] = ieeg_bipolarSEEG(MdataEphys(:, :, 1)', chNamesEphys, find(all_channels.soz(ephysIdxes) == 1), seg5, seg6, false); + +% Get all bipolar channels +MdataEphysBip = nan(length(bipolarNames), length(tt), height(eventsST)); +for ii = 1:height(eventsST) + MdataEphysBip(:, :, ii) = ieeg_bipolarSEEG(MdataEphys(:, :, ii)', chNamesEphys, [], seg5, seg6, false)'; +end + +%% 2) Calculate broadband on trial data and save +fprintf('Filtering for broadband per channel\n'); +bands = [70, 80; 80, 90; 90, 100; 100, 110; 130, 140; 140, 150; 150, 160; 160, 170]; % 10 hz bins, avoiding 60 and 120, matches NSDieegPrep CAR bands +MbbBip = nan(size(MdataEphysBip)); +for ii = 1:length(bipolarNames) + fprintf('.'); + MbbBip(ii, :, :) = ieeg_getHilbert(squeeze(MdataEphysBip(ii, :, :)), bands, srate, 'power'); +end +fprintf('\n'); + +fprintf('Saving bipolar broadband data...'); +MbbBip = single(MbbBip); % single for less storage +all_channels_bipolar = struct(); +all_channels_bipolar.name = bipolarNames; +all_channels_bipolar.originalChannels = bipolarChans; +all_channels_bipolar.badChannels = badChansBip; +all_channels_bipolar.soz = sozBip; +save(fullfile(localDataPath.output, 'derivatives', 'pipeline', sub, sprintf('%s_desc-preprocBipolarBB_ieeg.mat', sub)), ... + 'MbbBip', 'tt', 'srate', 'all_channels_bipolar', 'eventsST', '-v7.3'); +fprintf('.\n'); + +%end diff --git a/script02_writeMNI.m b/script02_writeMNI.m index 49240d2..a6a5307 100644 --- a/script02_writeMNI.m +++ b/script02_writeMNI.m @@ -1,6 +1,7 @@ %% This script calculates MNI152 and MNI305 positions each subject and saves them -localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) +%localDataPath = setLocalDataPath(1); % runs local PersonalDataPath (gitignored) +localDataPath = personalDataPath(); addpath('functions'); ss = 17; @@ -14,7 +15,11 @@ elecs = ieeg_readtableRmHyphens(elecsPath); elecmatrix = [elecs.x, elecs.y, elecs.z]; -%% Get and save MNI152 positions to electrodesMni152.tsv (volumetric, SPM12) +% FS dir at root level, then of subject +FSRootdir = fullfile(localDataPath.input, 'sourcedata', 'freesurfer'); +FSdir = fullfile(FSRootdir, sprintf('sub-%s', sub_label)); + +%% 1) 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))); @@ -26,33 +31,112 @@ mkdir(rootdirMni); % calculate MNI152 coordinates for electrodes -xyzMni152 = ieeg_getXyzMni(elecmatrix, niiPath, rootdirMni); +xyzsMni152 = 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); +elecsMni152.x = xyzsMni152(:, 1); elecsMni152.y = xyzsBipMni152(:, 2); elecsMni152.z = xyzsBipMni152(:, 3); writetable(elecsMni152, elecsMni152Path, 'FileType', 'text', 'Delimiter', '\t'); fprintf('Saved to %s\n', elecsMni152Path); -%% Get and save MNI305 positions (through fsaverage) +%% 2) Get and save MNI305 positions (through fsaverage) + +% calculate MNI305 coordinates for electrodes +[xyzsMni305, vertIdxFsavg, minDists, surfUsed] = ieeg_mni305ThroughFsSphere(elecmatrix, elecs.hemisphere, FSdir, FSRootdir, 'closest', 5); + +% save as separate MNI 305 electrodes table +elecsBipMni305Path = fullfile(localDataPath.input, ['sub-' sub_label], ['ses-' ses_label], 'ieeg', ['sub-' sub_label '_ses-' ses_label '_space-' 'MNI305' '_electrodes.tsv']); +elecsBipMni305 = elecs; +elecsBipMni305.x = xyzsMni305(:, 1); elecsBipMni305.y = xyzsMni305(:, 2); elecsBipMni305.z = xyzsMni305(:, 3); +elecsBipMni305.vertex_fsaverage = vertIdxFsavg; % also add a column to indicate vertex on fsavg, so we can easily get position for inflated brain +writetable(elecsBipMni305, elecsBipMni305Path, 'FileType', 'text', 'Delimiter', '\t'); + +fprintf('Saved to %s\n', elecsBipMni305Path); + +%% 3) Calculate bipolar electrodes in native space, save + +% note: we decide not to sort electrodes by channels before bipolar deriv, to avoid many n/a rows. +% Channels.tsv contains non-ephys channels anyway, which we would need to remove before sorting and before matching to elecs when loading + +% load segmented data to get segs for this subject +participants = readtable(fullfile(localDataPath.input, 'participants.tsv'), 'Filetype', 'text', 'Delimiter', '\t'); +seg5 = participants.seg5{strcmp(participants.participant_id, sprintf('sub-%s', sub_label))}; +seg5 = strip(split(seg5, ',')); +seg6 = participants.seg6{strcmp(participants.participant_id, sprintf('sub-%s', sub_label))}; +seg6 = strip(split(seg6, ',')); +if strcmp(seg5, 'n/a'), seg5 = {}; end +if strcmp(seg6, 'n/a'), seg6 = {}; end + +% create dummy data to get bipolar names/channels from +Mdummy = zeros(1, height(elecs)); % note warnings here are ok, because electrodes.tsv doesn't contain all channels +[~, bipolarNames, bipolarChans] = ieeg_bipolarSEEG(Mdummy, elecs.name, [], seg5, seg6); + +% get hemisphere info, checking that bipolar pair have the same hemi (or is outside of brain '') +hemiBip = elecs.hemisphere(bipolarChans); +assert(all(strcmp(hemiBip(:, 1), hemiBip(:, 2)) | any(strcmp(hemiBip, ''), 2)), ... + 'Error: hemisphere mismatch found between electrodes in bipolar pairs'); +hemiBip = hemiBip(:, 1); % keep just one col + +% get bipolar electrode positions +xyzsBip = nan(length(bipolarNames), 3); +for ii = 1:length(bipolarNames) + xyzsBip(ii, :) = mean(elecmatrix(bipolarChans(ii, :), :)); % center of the pair xyzs +end + +% get interpolated destrieux position (requires vistasoft +[labs, labs_val] = ieeg_getLabelXyzDestrieux(xyzsBip, FSdir, 3); + +% keep seizure_zone columns for each original electrode +sozs = elecs.seizure_zone; +sozs(cellfun(@isempty, sozs)) = {'n/a'}; +soz_labels = sozs(bipolarChans); + +% assemble table and save +elecsBip = table(bipolarNames, xyzsBip(:, 1), xyzsBip(:, 2), xyzsBip(:, 3), hemiBip, soz_labels(:, 1), soz_labels(:, 2), labs_val, labs, ... + 'VariableNames', {'name', 'x', 'y', 'z', 'hemisphere', 'seizure_zone_1', 'seizure_zone_2', 'Destrieux_label', 'Destrieux_label_text'}); +bipolarDir = fullfile(localDataPath.output, 'derivatives', 'pipeline', sprintf('sub-%s', sub_label)); +mkdir(bipolarDir); +writetable(elecsBip, fullfile(bipolarDir, sprintf('sub-%s_desc-bipolar_electrodes.tsv', sub_label)), 'Filetype', 'text', 'Delimiter', '\t'); -% FS dir of current subject -FSdir = fullfile(localDataPath.input, 'sourcedata', 'freesurfer', sprintf('sub-%s', sub_label)); -FSsubjectsdir = fullfile(FSdir, '..'); +%% 4) Bipolar - Get and save MNI152 positions (run 3 first) + +% 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. It will be hard to tell which were for bipolar, which were CAR elecs. Is this important? +rootdirMni = fullfile(localDataPath.input, 'derivatives', 'MNI152_electrode_transformations', sprintf('sub-%s', sub_label)); +mkdir(rootdirMni); + +% calculate MNI152 coordinates for bipolar electrodes +xyzsBipMni152 = ieeg_getXyzMni(xyzsBip, niiPath, rootdirMni); + +% save as separate MNI 152 electrodes table +elecsBipMni152Path = fullfile(bipolarDir, sprintf('sub-%s_space-MNI152NLin2009_desc-bipolar_electrodes.tsv', sub_label)); +elecsBipMni152 = elecsBip; +elecsBipMni152.x = xyzsBipMni152(:, 1); elecsBipMni152.y = xyzsBipMni152(:, 2); elecsBipMni152.z = xyzsBipMni152(:, 3); +writetable(elecsBipMni152, elecsBipMni152Path, 'FileType', 'text', 'Delimiter', '\t'); + +fprintf('Saved to %s\n', elecsBipMni152Path); + +%% 5) Bipolar - Get and save MNI305 positions (run 3 first) % calculate MNI305 coordinates for electrodes -[xyzMni305, vertIdxFsavg, minDists, surfUsed] = ieeg_mni305ThroughFsSphere(elecmatrix, elecs.hemisphere, FSdir, FSsubjectsdir, 'closest', 5); +[xyzsBipMni305, vertIdxFsavg] = ieeg_mni305ThroughFsSphere(xyzsBip, elecsBip.hemisphere, FSdir, FSRootdir, 'closest', 5); % 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'); +elecsBipMni305Path = fullfile(bipolarDir, sprintf('sub-%s_space-MNI305_desc-bipolar_electrodes.tsv', sub_label)); +elecsBipMni305 = elecsBip; +elecsBipMni305.x = xyzsBipMni305(:, 1); elecsBipMni305.y = xyzsBipMni305(:, 2); elecsBipMni305.z = xyzsBipMni305(:, 3); +elecsBipMni305.vertex_fsaverage = vertIdxFsavg; % also add a column to indicate vertex on fsavg, so we can easily get position for inflated brain +writetable(elecsBipMni305, elecsBipMni305Path, 'FileType', 'text', 'Delimiter', '\t'); + +fprintf('Saved to %s\n', elecsBipMni305Path); -fprintf('Saved to %s\n', elecsMni305Path); +return %% Normalize bb power per run