Updates to preprocessing, mostly involving bipolar reref

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

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

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