Authored by Anthony Boyer

Re-work of spikedetection so it can deal with inconsistent cropped stimulations.

@@ -2,110 +2,142 @@ function ImaGIN_SpikesDetection(S) @@ -2,110 +2,142 @@ function ImaGIN_SpikesDetection(S)
2 % This function extracts and concatenates the baselines from cropped, qualitatively controlled, stimulation files. 2 % This function extracts and concatenates the baselines from cropped, qualitatively controlled, stimulation files.
3 % It runs Delphos' spike detection on extracted baselines after bipolar montage computation. 3 % It runs Delphos' spike detection on extracted baselines after bipolar montage computation.
4 % It then computes corrected spike rates for each bipolar recording channel. 4 % It then computes corrected spike rates for each bipolar recording channel.
5 -% Manipulating the data of hundreds of stimulation files brings performance (RAM) issues...  
6 -% ...which are addressed by zeroing the baselines during the extraction loop, at the expense of readability. 5 +% /!\ Cropped files CAN have a variable number of channels and thus the monopolar <-> bipolar mapping has to be established for each cropped file
  6 +% and taken into account during baseline concatenation for consistency. /!\
7 7
8 - % Vars parsing 8 + %% Vars parsing
9 stims_files = S.stims_files; 9 stims_files = S.stims_files;
10 - path_out = S.path_out; 10 + path_out = S.path_out;
  11 + channels = struct;
11 12
12 - % Establish the channel mapping from monopolar to bipolar using a single stimulation file  
13 - mono_stimulation_file = fullfile(path_out,'single_stimulation');  
14 - mono_stimulation_obj = clone(spm_eeg_load(stims_files{1}),mono_stimulation_file);  
15 - monopolar_sensors = sensors(mono_stimulation_obj,'EEG');  
16 - Sbp.Fname = mono_stimulation_file;  
17 - bp_stimulation_obj = ImaGIN_BipolarMontage(Sbp);  
18 - clear channels_map;  
19 - for cc=1:nchannels(bp_stimulation_obj)  
20 - bipolar_chan = chanlabels(bp_stimulation_obj,cc);  
21 - mono_chans = regexp(bipolar_chan,'-','split');  
22 - channels_map(cc).bp_chan = bipolar_chan;  
23 - channels_map(cc).bp_idx = cc;  
24 - channels_map(cc).mono_chan1 = mono_chans{1}{1};  
25 - channels_map(cc).mono_chan2 = mono_chans{1}{2};  
26 - % Get the monopolar channels idxes.  
27 - chan1_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{1}));  
28 - chan2_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{2}));  
29 - % Test and correct for duplicated labels in the monopolar object  
30 - if numel(chan1_idx) == 0  
31 - error('Could not find the monopolar channel from the bipolar montage.');  
32 - elseif numel(chan1_idx) == 1 % If 1 match, as expected, do nothing  
33 - %pass  
34 - else  
35 - error('2 or more channels with the same label')  
36 - end  
37 - % Same logic as for chan1_idx  
38 - if numel(chan2_idx) == 0  
39 - error('Could not find the monopolar channel from the bipolar montage.');  
40 - elseif numel(chan2_idx) == 1  
41 - %pass  
42 - else  
43 - error('2 or more channels with the same label')  
44 - end  
45 - channels_map(cc).mono_idx1 = chan1_idx;  
46 - channels_map(cc).mono_idx2 = chan2_idx;  
47 - end  
48 - clear mono_stimulation_obj bp_stimulation_obj monopolar_sensors;  
49 -  
50 - % Extract baselines: get the pre-stimulation period for each stimulation and concatenate them for each channel.  
51 - concat_baselines = []; % Store timeseries  
52 - concat_zeroed = []; % Store booleans later used for the correction of spike rates 13 + %% Extract data from monopolar stimulation files
53 for ff=1:numel(stims_files) 14 for ff=1:numel(stims_files)
54 - clear D;  
55 - D = spm_eeg_load(stims_files{ff});  
56 - % Get the list of monopolar bad channels 15 + % Establish the mappping for the stimulation
  16 + stim_mapping = monopolar_to_bipolar_mapping(stims_files{ff},path_out);
  17 + % Load the stimulation for data extraction
  18 + D = spm_eeg_load(stims_files{ff});
  19 + % Get the samping frequency
  20 + sampling(ff) = fsample(D);
  21 + % Get the coordinates
  22 + sensors_info = sensors(D,'EEG');
  23 + % Get the list of bad monopolar channels
57 all_channels_idxes = indchantype(D,'ALL'); 24 all_channels_idxes = indchantype(D,'ALL');
58 good_eeg_channels_idxes = indchantype(D,'EEG','GOOD'); 25 good_eeg_channels_idxes = indchantype(D,'EEG','GOOD');
59 mono_bad_channel_idxes = setdiff(all_channels_idxes,good_eeg_channels_idxes); 26 mono_bad_channel_idxes = setdiff(all_channels_idxes,good_eeg_channels_idxes);
60 - % Use the monopolar-bipolar channel mapping to identify monopolar channels later coupled with bad channels 27 + % Use the mapping to identify monopolar channels later coupled with bad channels
61 bp_bad_channels_idxes = []; 28 bp_bad_channels_idxes = [];
62 for cc=1:numel(mono_bad_channel_idxes) 29 for cc=1:numel(mono_bad_channel_idxes)
63 - bad_channel_1_map = find(mono_bad_channel_idxes(cc) == [channels_map.mono_idx1]);  
64 - bad_channel_2_map = find(mono_bad_channel_idxes(cc) == [channels_map.mono_idx2]);  
65 - circular_bp_bad_channels_idxes = unique([channels_map(bad_channel_1_map).mono_idx1,...  
66 - channels_map(bad_channel_1_map).mono_idx2,...  
67 - channels_map(bad_channel_2_map).mono_idx1,...  
68 - channels_map(bad_channel_2_map).mono_idx2]); 30 + bad_channel_1_map = find(mono_bad_channel_idxes(cc) == [stim_mapping.mono_idx1]);
  31 + bad_channel_2_map = find(mono_bad_channel_idxes(cc) == [stim_mapping.mono_idx2]);
  32 + circular_bp_bad_channels_idxes = unique([stim_mapping(bad_channel_1_map).mono_idx1,...
  33 + stim_mapping(bad_channel_1_map).mono_idx2,...
  34 + stim_mapping(bad_channel_2_map).mono_idx1,...
  35 + stim_mapping(bad_channel_2_map).mono_idx2]);
69 bp_bad_channels_idxes = unique([bp_bad_channels_idxes circular_bp_bad_channels_idxes]); 36 bp_bad_channels_idxes = unique([bp_bad_channels_idxes circular_bp_bad_channels_idxes]);
70 end 37 end
71 % Concatenate all identified bad channels 38 % Concatenate all identified bad channels
72 all_bad_channels_idxes = unique([mono_bad_channel_idxes bp_bad_channels_idxes]); 39 all_bad_channels_idxes = unique([mono_bad_channel_idxes bp_bad_channels_idxes]);
73 - % Extract baseline of the stimulation and zeroed all bad channels 40 + % Extract baseline of the stimulation and zero all bad channels
74 time_axis = time(D); 41 time_axis = time(D);
75 onset = timeonset(D); 42 onset = timeonset(D);
76 - negative_time_idxes = find(time_axis > onset+1 & time_axis < -1); % Arbitrary time margins of 1 s around prestimulus period  
77 -  
78 - disp(size(D));  
79 - disp(stims_files{ff});  
80 - 43 + negative_time_idxes = find(time_axis > onset+1 & time_axis < -1); % Arbitrary time margins of 1 s around prestimulus period
81 file_baselines = D(:,negative_time_idxes,1); 44 file_baselines = D(:,negative_time_idxes,1);
82 - file_baselines(all_bad_channels_idxes,:,1) = zeros(numel(all_bad_channels_idxes),numel(negative_time_idxes));  
83 - concat_baselines = [concat_baselines file_baselines]; 45 + file_baselines(all_bad_channels_idxes,:,1) = zeros(numel(all_bad_channels_idxes),numel(negative_time_idxes));
84 % Generate the corresponding boolean matrix to identify good (1) versus bad (0) samples 46 % Generate the corresponding boolean matrix to identify good (1) versus bad (0) samples
85 file_zeroed = ones(size(file_baselines)); 47 file_zeroed = ones(size(file_baselines));
86 file_zeroed(all_bad_channels_idxes,:,1) = zeros(numel(all_bad_channels_idxes),numel(negative_time_idxes)); 48 file_zeroed(all_bad_channels_idxes,:,1) = zeros(numel(all_bad_channels_idxes),numel(negative_time_idxes));
87 - concat_zeroed = [concat_zeroed file_zeroed]; 49 + % Update the channels struct
  50 + eeg_chans_idxes = indchantype(D,'EEG');
  51 + for cc=1:numel(eeg_chans_idxes)
  52 + channel_name = chanlabels(D,eeg_chans_idxes(cc));
  53 + channels.(channel_name{:}).timeserie{ff,:} = file_baselines(eeg_chans_idxes(cc),:,:);
  54 + channels.(channel_name{:}).zeroed{ff,:} = file_zeroed(eeg_chans_idxes(cc),:,:);
  55 + channels.(channel_name{:}).coordinates{ff,:} = sensors_info.elecpos(find(strcmp(channel_name,sensors_info.label)),:);
  56 + end
  57 + % Tidy the channels struct with empty inputs for consistency
  58 + temp_montage = fieldnames(channels);
  59 + for nn=1:numel(temp_montage)
  60 + if sum(strcmp(temp_montage{nn},chanlabels(D))) == 0
  61 + channels.(temp_montage{nn}).timeserie{ff} = [];
  62 + channels.(temp_montage{nn}).zeroed{ff} = [];
  63 + channels.(temp_montage{nn}).coordinates{ff} = [];
  64 + end
  65 + end
88 end 66 end
89 -  
90 - % Generate the baseline object.  
91 - baseline_file = fullfile(path_out,'Baseline');  
92 - baseline_obj = clone(D,baseline_file,[size(concat_baselines,1) size(concat_baselines,2) 1],3);  
93 - baseline_obj(:) = concat_baselines(:);  
94 - baseline_obj = events(baseline_obj,[],{''});  
95 - baseline_obj = timeonset(baseline_obj,0);  
96 - clear baseline_obj;  
97 67
98 - % Compute bipolar montage on baseline object to get the bipolar baseline object  
99 - S.Fname = baseline_file;  
100 - bp_baseline_obj = ImaGIN_BipolarMontage(S); 68 + %% Test for extracted baseline sizes and fill-in missing channels
  69 + global_montage = fieldnames(channels);
  70 + for ff=1:numel(stims_files)
  71 + for nn=1:numel(global_montage)
  72 + series_length(nn) = numel(channels.(global_montage{nn}).timeserie{ff,:});
  73 + end
  74 + empty_series = find(series_length == 0);
  75 + nonempty_series = find(series_length);
  76 + valid_length = series_length(nonempty_series);
  77 + if all(valid_length == valid_length(1))
  78 + for ee=1:numel(empty_series)
  79 + channels.(global_montage{empty_series(ee)}).timeserie{ff,:} = zeros(1,valid_length(1));
  80 + channels.(global_montage{empty_series(ee)}).zeroed{ff,:} = zeros(1,valid_length(1));
  81 + channels.(global_montage{empty_series(ee)}).coordinates{ff,:} = NaN(1,3);
  82 + end
  83 + else
  84 + error('Inconsistent length of channels.');
  85 + end
  86 + end
  87 +
  88 + %% Create the new monopolar baseline object
  89 + for nn=1:numel(global_montage)
  90 + channel{nn} = global_montage{nn};
  91 + timeseries(nn,:) = [channels.(global_montage{nn}).timeserie{:}];
  92 + zeroed(nn,:) = [channels.(global_montage{nn}).zeroed{:}];
  93 + coordinates(nn,:) = channels.(global_montage{nn}).coordinates{1,:};
  94 + end
  95 + % Create the new sensor
  96 + new_sensors_info.balance.current = 'none';
  97 + new_sensors_info.chanpos = coordinates;
  98 + new_sensors_info.chantype = cell(1,numel(global_montage)); new_sensors_info.chantype(:) = {'eeg'};
  99 + new_sensors_info.chanunit = cell(1,numel(global_montage)); new_sensors_info.chanunit(:) = {'V'};
  100 + new_sensors_info.elecpos = coordinates;
  101 + new_sensors_info.label = channel;
  102 + new_sensors_info.type = 'ctf';
  103 + new_sensors_info.unit = 'mm';
  104 + % Create the new files
  105 + file_header = fullfile(path_out,'monopolar_baseline.mat');
  106 + file_data = fullfile(path_out,'monopolar_baseline.dat');
  107 + if exist(file_header, 'file') == 2
  108 + delete(file_header);
  109 + end
  110 + if exist(file_data, 'file') == 2
  111 + delete(file_data);
  112 + end
  113 + % Create the new object
  114 + D = meeg(size(timeseries,1),size(timeseries,2),1);
  115 + D = fname(D,'monopolar_baseline');
  116 + D = path(D,path_out);
  117 + if all(sampling == sampling(1))
  118 + D = fsample(D,sampling(1));
  119 + else
  120 + error('Inconsistent sampling frequency across cropped files.');
  121 + end
  122 + D = sensors(D,'EEG',new_sensors_info)
  123 + D = blank(D) ;
  124 + D = link(D,file_data);
  125 + D(:,:,1) = timeseries;
  126 + save(D);
101 127
  128 + %% Compute bipolar montage on baseline object to get the bipolar baseline object
  129 + Sbp.Fname = file_header;
  130 + Sbp.FileOut = fullfile(path_out,'bipolar_baseline');
  131 + bp_baseline_obj = ImaGIN_BipolarMontage(Sbp);
  132 + global_mapping = monopolar_to_bipolar_mapping(Sbp.Fname,path_out);
  133 +
102 % Estimate bipolar boolean matrix from monopolar boolean matrix 134 % Estimate bipolar boolean matrix from monopolar boolean matrix
103 bp_concat_zeroed = zeros(size(bp_baseline_obj)); 135 bp_concat_zeroed = zeros(size(bp_baseline_obj));
104 - for cc=1:numel(channels_map)  
105 - both_good = (concat_zeroed(channels_map(cc).mono_idx1,:,1) == 1 & concat_zeroed(channels_map(cc).mono_idx2,:,1) == 1); 136 + for cc=1:numel(global_mapping)
  137 + both_good = (zeroed(global_mapping(cc).mono_idx1,:,1) & zeroed(global_mapping(cc).mono_idx2,:,1));
106 bp_concat_zeroed(cc,find(both_good)) = ones(1,numel(find(both_good))); 138 bp_concat_zeroed(cc,find(both_good)) = ones(1,numel(find(both_good)));
107 end 139 end
108 - 140 +
109 % Run Delphos' spike detection on baselines 141 % Run Delphos' spike detection on baselines
110 fs = bp_baseline_obj.fsample; 142 fs = bp_baseline_obj.fsample;
111 results = Delphos_detector(bp_baseline_obj,chanlabels(bp_baseline_obj),'SEEG',fs,{'Spk'},[],[],40,[]); 143 results = Delphos_detector(bp_baseline_obj,chanlabels(bp_baseline_obj),'SEEG',fs,{'Spk'},[],[],40,[]);
@@ -117,9 +149,9 @@ function ImaGIN_SpikesDetection(S) @@ -117,9 +149,9 @@ function ImaGIN_SpikesDetection(S)
117 zeroed_samples = numel(find(~bp_concat_zeroed(cc,:,1))); 149 zeroed_samples = numel(find(~bp_concat_zeroed(cc,:,1)));
118 total_times(cc) = (total_samples - zeroed_samples) /fs/60; % Total recording time in minutes 150 total_times(cc) = (total_samples - zeroed_samples) /fs/60; % Total recording time in minutes
119 end 151 end
120 - channels = [results.labels]'; 152 + chans = [results.labels]';
121 rates = round(results.n_Spk./total_times',1); % Round to nearest decimal 153 rates = round(results.n_Spk./total_times',1); % Round to nearest decimal
122 - rates_table = table(channels,rates); 154 + rates_table = table(chans,rates);
123 155
124 % Store spike features in a table 156 % Store spike features in a table
125 spks_positions = [results.markers.position]'; 157 spks_positions = [results.markers.position]';
@@ -139,10 +171,49 @@ function ImaGIN_SpikesDetection(S) @@ -139,10 +171,49 @@ function ImaGIN_SpikesDetection(S)
139 171
140 % Store rate table in an individual files for easy extraction 172 % Store rate table in an individual files for easy extraction
141 rate_table_file = 'SPK_bBaseline_rates.mat'; 173 rate_table_file = 'SPK_bBaseline_rates.mat';
142 - spike_rates = [channels,num2cell(rates)]; 174 + spike_rates = [chans,num2cell(rates)];
143 save(fullfile(path_out,rate_table_file), 'spike_rates'); 175 save(fullfile(path_out,rate_table_file), 'spike_rates');
144 176
145 % Ending step 177 % Ending step
146 set_final_status('OK'); 178 set_final_status('OK');
147 disp('Done'); 179 disp('Done');
148 end 180 end
  181 +
  182 +function channels_map = monopolar_to_bipolar_mapping(monopolar_file,path_out)
  183 + mono_stimulation_obj = spm_eeg_load(monopolar_file);
  184 + monopolar_sensors = sensors(mono_stimulation_obj,'EEG');
  185 + [mono_path,mono_name,mono_ext] = fileparts(monopolar_file);
  186 + Sbp.Fname = monopolar_file;
  187 + Sbp.FileOut = fullfile(path_out,['b' mono_name]);
  188 + bp_stimulation_obj = ImaGIN_BipolarMontage(Sbp);
  189 + for cc=1:nchannels(bp_stimulation_obj)
  190 + bipolar_chan = chanlabels(bp_stimulation_obj,cc);
  191 + mono_chans = regexp(bipolar_chan,'-','split');
  192 + channels_map(cc).bp_chan = bipolar_chan;
  193 + channels_map(cc).bp_idx = cc;
  194 + channels_map(cc).mono_chan1 = mono_chans{1}{1};
  195 + channels_map(cc).mono_chan2 = mono_chans{1}{2};
  196 + % Get the monopolar channels idxes.
  197 + chan1_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{1}));
  198 + chan2_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{2}));
  199 + % Test and correct for duplicated labels in the monopolar object
  200 + if numel(chan1_idx) == 0
  201 + error('Could not find the monopolar channel from the bipolar montage.');
  202 + elseif numel(chan1_idx) == 1 % If 1 match, as expected, do nothing
  203 + %pass
  204 + else
  205 + error('2 or more channels with the same label')
  206 + end
  207 + % Same logic as for chan1_idx
  208 + if numel(chan2_idx) == 0
  209 + error('Could not find the monopolar channel from the bipolar montage.');
  210 + elseif numel(chan2_idx) == 1
  211 + %pass
  212 + else
  213 + error('2 or more channels with the same label')
  214 + end
  215 + channels_map(cc).mono_idx1 = chan1_idx;
  216 + channels_map(cc).mono_idx2 = chan2_idx;
  217 + end
  218 + delete([Sbp.FileOut '.mat']);delete([Sbp.FileOut '.dat']);delete([Sbp.FileOut '_log.txt']);
  219 +end