-
Notifications
You must be signed in to change notification settings - Fork 2
/
pcc_fig3b_singleTrials.m
162 lines (132 loc) · 7.83 KB
/
pcc_fig3b_singleTrials.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
clearvars, clc, close all
%% Plot pannels with CCEPs coming from single-pair
addpath(genpath(pwd));
% set local path to your BIDS directory:
myPath = setLocalDataPath(1);
localDataPath = myPath.input;
% load the metadata across subjects
all_subjects = {'01','02','03','04','05','06','07','08'};
all_hemi = {'r','r','r','l','r','l','l','r'};
all_runs = {'01','01','01','01','01','01','01','01'};
out = []; % prepare to concatinate outputs across subjects
for ss = 1:length(all_subjects)
bids_sub = all_subjects{ss};
bids_ses = 'ieeg01';
bids_task = 'ccep';
bids_run = all_runs{ss};
% set var names: mefd, events, channels and electrodes
fileName = fullfile(localDataPath,['sub-' bids_sub],['ses-' bids_ses],'ieeg',...
['sub-' bids_sub '_ses-' bids_ses '_task-' bids_task '_run-' bids_run '_ieeg.mefd']);
events_tsv_name = fullfile(localDataPath,['sub-' bids_sub],['ses-' bids_ses],'ieeg',...
['sub-' bids_sub '_ses-' bids_ses '_task-' bids_task '_run-' bids_run '_events.tsv']);
channels_tsv_name = fullfile(localDataPath,['sub-' bids_sub],['ses-' bids_ses],'ieeg',...
['sub-' bids_sub '_ses-' bids_ses '_task-' bids_task '_run-' bids_run '_channels.tsv']);
electrodes_tsv_name = fullfile(localDataPath,['sub-' bids_sub],['ses-' bids_ses],'ieeg',...
['sub-' bids_sub '_ses-' bids_ses '_electrodes.tsv']);
% load metadata, events, channels, and electrodes
[metadata] = readMef3(fileName); % metadata table
events_table = readtable(events_tsv_name,'FileType','text','Delimiter','\t','TreatAsEmpty',{'N/A','n/a'}); % events table
channels_table = readtable(channels_tsv_name,'FileType','text','Delimiter','\t','TreatAsEmpty',{'N/A','n/a'}); % channels table
electrodes_table= readtable(electrodes_tsv_name,'FileType','text','Delimiter','\t','TreatAsEmpty',{'N/A','n/a'}); % electrodes table
% ------------------------------------------------------------------- %
% get necessary parameters from data
srate = metadata.time_series_metadata.section_2.sampling_frequency; % sampling frequency
nr_channels = length(metadata.time_series_channels); % number of channels
% list of channel names
channel_names = cell(nr_channels,1);
for kk = 1:nr_channels
channel_names{kk} = metadata.time_series_channels(kk).name;
end
% find good sEEG/ECoG channels
good_channels = find(ismember(channels_table.type,{'SEEG'}) & ismember(channels_table.status,'good'));
% load event data, clip out other stim amplitudes
events_table_clipped = bids_clipEvents(events_table,'electrical_stimulation_current', {'4.0 mA', '6.0 mA'}); % keep only events with stim current == 4.0 or 6.0 mA
% see which stim pairs are in areas of interest
areas_interest = [17 18 11106 11107 11108 11109 11110 11123 10 53 54 12106 12107 12108 12109 12110 12123 49]; % Destrieux_label
% get areas fo each channel name and each ccep stim pair name
channel_names = channels_table.name;
channel_areas = zeros(size(channel_names));
ccep_stim_names = unique(events_table_clipped.electrical_stimulation_site);
ccep_stim_areas = zeros(length(ccep_stim_names),2);
% loop through channel names to get channel areas
for kk = 1:length(channel_names)
% which number has this channel in the electrodes.tsv file?
thisElPos = find(ismember(electrodes_table.name,channel_names{kk}));
if ismember(electrodes_table.Destrieux_label(thisElPos),areas_interest)
channel_areas(kk) = areas_interest(ismember(areas_interest,electrodes_table.Destrieux_label(thisElPos)));
end
end
% get areas for each stim pair name
for kk = 1:length(ccep_stim_names)
if sum(ismember(ccep_stim_names{kk},'-'))==1
ccep_stim_areas(kk,1) = channel_areas(ismember(channel_names,extractBefore(ccep_stim_names{kk},'-')));
ccep_stim_areas(kk,2) = channel_areas(ismember(channel_names,extractAfter(ccep_stim_names{kk},'-')));
else % assume the second - if there is a - in the channel name
dash_in_name = find(ismember(ccep_stim_names{kk},'-'));
el1_name = ccep_stim_names{kk}(1:dash_in_name(2)-1);
el2_name = ccep_stim_names{kk}(dash_in_name(2)+1:end);
ccep_stim_areas(kk,1) = channel_areas(ismember(channel_names,el1_name));
ccep_stim_areas(kk,2) = channel_areas(ismember(channel_names,el2_name));
end
end
% save filename, channel areas and good channels
out(ss).fileName = fileName;
out(ss).channel_names = channel_names;
out(ss).channel_areas = channel_areas;
out(ss).good_channels = good_channels;
% only keep events where we stimulated the limbic network
out(ss).ccep_stim_areas_limbic = ccep_stim_areas(sum(ccep_stim_areas,2) > 0,:);
out(ss).ccep_stim_names_limbic = ccep_stim_names(sum(ccep_stim_areas,2)>0,:);
out(ss).events_table_clipped = bids_clipEvents(events_table_clipped,'electrical_stimulation_site', ccep_stim_names(sum(ccep_stim_areas,2)>0,:)); % keep only events with stim current == 4.0 or 6.0 mA
end
%% Plot single trials for Sub-01
ss = 1;
el_stim = {'RC1-RC2','RC2-RC3','RC3-RC4'};
el_record = 'RY2';
figure('Position',[0 0 550 800]), hold on; % prepare to plot
for kk = 1:length(el_stim)
events_table_thisSite = bids_clipEvents(out(ss).events_table_clipped,'electrical_stimulation_site', el_stim{kk});
% Settings for CRP
baseline_t = [-0.5 -0.05];
t_win_crp = [0.015 1];
% use adjusted CAR
[average_ccep, average_ccep_names, tt, srate, crp_out, single_trials] = ...
ccepPCC_loadAverageSubset(out(ss).fileName, events_table_thisSite, ...
out(ss).good_channels, out(ss).channel_areas, baseline_t, t_win_crp, 1, 1);
chan_ind = find(ismember(out(ss).channel_names,el_record)); % get channel index
plot_responses = squeeze(single_trials(chan_ind,:,:));
% plot single and average responses
subplot(3,1,kk), hold on
plot(tt(tt >0.010 & tt <2),zeros(size(tt(tt >0.010 & tt <2))),'k:') % plot zero line
plot(tt,plot_responses,'Color',[.7 .7 .7 .3]) % single responses in gray
plot(tt,mean(plot_responses),'k','LineWidth',1.5) % mean in thicker black line
% set axis and labels
xlabel('Time(s)'), xlim([-.2 2]);
ylabel('Voltage (uV)'),ylim([-500 1000]);
end
%% --------------------------------------------------
%% Plot single trials for Sub-08
ss = 8;
el_stim = {'RC1-RC2'}; % stim site with 55 trials
el_record = {'RZ1','RPO1','RQ1'}; % all PCC sites
events_table_thisSite = bids_clipEvents(out(ss).events_table_clipped,'electrical_stimulation_site', el_stim);
% Settings for CRP
baseline_t = [-0.5 -0.05];
t_win_crp = [0.015 1];
% load and preprocess data
[average_ccep, average_ccep_names, tt, srate, crp_out, single_trials] = ...
ccepPCC_loadAverageSubset(out(ss).fileName, events_table_thisSite, ...
out(ss).good_channels, out(ss).channel_areas, baseline_t, t_win_crp, 1, 1);
figure('Position',[0 0 550 800]), hold on;
for kk = 1:length(el_record)
chan_ind = find(ismember(out(ss).channel_names,el_record{kk})); % get channel index
plot_responses = squeeze(single_trials(chan_ind,:,:));
% plot single and average responses
subplot(3,1,kk), hold on
plot(tt(tt >0.010 & tt <2),zeros(size(tt(tt >0.010 & tt <2))),'k:') % plot zero line
plot(tt,plot_responses,'Color',[.7 .7 .7 .3]) % single responses in gray
plot(tt,mean(plot_responses),'k','LineWidth',1.5) % mean in thicker black line
% set axis and labels
xlabel('Time(s)'), xlim([-.2 2]);
ylabel('Voltage (uV)'),ylim([-500 1000]);
end