-
Notifications
You must be signed in to change notification settings - Fork 3
/
NET_ICU_Pipeline.m
179 lines (144 loc) · 8.23 KB
/
NET_ICU_Pipeline.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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
%% NET_ICU Pipeline
% This code is still in progress and not a final version
% Resonsible authors:
% Miriam Han
% Charlotte Maschke
%% Define Parameters
% wPLI and dPLI parameters
frequencies = ["alpha"]; % This can be ["alpha" "theta" "delta"]
%frequencies = ["alpha"]; % This can be ["alpha" "theta" "delta"]
window_size = 10; % This is in seconds and will be how we chunk the whole dataset
number_surrogate = 20; % Number of surrogate wPLI to create / # of permutations
p_value = 0.05; % the p value to make our test on
step_size = window_size;
% HUB parameters
% threshold can be either an float between 0 and 1 for an absolute threshold or "MSG" standing for minimally
% spanning graph. If "MSG" is used, you need to define the range of thrsholds
%threshold = "SCT"; % Smallest connected Threshold OR a value between 0 to 1
threshold_range = 0.9:-0.01:0.01; % used ONLY if MSG More connected to less connected
threshold = 0.05;
%% Load clean EEG data set
% select the folder with the
waitfor(msgbox('Select the folder which contains the data in BIDS format.'));
datafolder=uigetdir(path);
waitfor(msgbox('Select the Saving directory'));
resultsfolder = uigetdir(path);
cd(datafolder)
cd("eeg")
files = dir('*.set*');
for f = 1:numel(files)
%% Load data and get information about the state
recording = load_set(files(f).name,pwd);
sampling_rate = recording.sampling_rate;
info = split(files(f).name,'_');
ID = info{1}(5:end);
task = info{2}(6:end);
hemispheres = ["Left", "Right","Whole"];
disp("load complete: " + ID + '_' + task)
%% Spectrogram
spectopo_prp = spectopo_prp_struct;
disp('spectopo_prp load complete')
outdir_spectrogram = fullfile(resultsfolder, ID,"Whole",'Spectrogram');
spectrogram = spectrogram_function(recording, spectopo_prp, ID, task, outdir_spectrogram);
%% Topographic Maps of Alpha and Theta Power
outdir_topographicmap = fullfile(resultsfolder, ID, "Whole",'Topographic Maps');
topographic_map = topographic_map_function(recording, spectopo_prp, ID, task, outdir_topographicmap);
%% wPLI analysis for alpha and theta frequencies
for fr = 1:length(frequencies)
frequency = frequencies(fr);
outdir = fullfile(resultsfolder , ID);
mkdir(fullfile(outdir,'fc_data'));
participant_out_path_wpli = strcat(fullfile(outdir,'fc_data'),filesep,'wpli_',frequency,'_',ID,'_',task,'.mat');
participant_out_path_dpli = strcat(fullfile(outdir,'fc_data'),filesep,'dpli_',frequency,'_',ID,'_',task,'.mat');
if frequency == "alpha"
frequency_band = [7 13]; % This is in Hz
elseif frequency == "theta"
frequency_band = [4 8]; % This is in Hz
elseif frequency == "delta"
frequency_band = [1 4]; % This is in Hz
end
%% Calculate Results for whole brain:
% Calculate the wpli
disp(strcat("Participant: ", ID , "_wPLI"));
result_wpli = na_wpli(recording, frequency_band, window_size, step_size, number_surrogate, p_value);
% Calculate the dpli
disp(strcat("Participant: ", ID , "_dPLI"));
result_dpli = na_dpli_corrected(recording, frequency_band, window_size, step_size, number_surrogate, p_value);
%% Save wPLI and dPLI for later Matlab and Python use
% save dPLI and wPLI for later in Matlab
save(participant_out_path_wpli,'result_wpli')
save(participant_out_path_dpli,'result_dpli')
%save wPLI and dPLI in an easier format for python
% add a 'py' at the end for python
participant_out_path_wpli = strcat(fullfile(outdir,'fc_data'),filesep,'wpli_',frequency,'_',ID,'_',task,'py.mat');
participant_out_path_dpli = strcat(fullfile(outdir,'fc_data'),filesep,'dpli_',frequency,'_',ID,'_',task,'py.mat');
participant_channel_path_wpli = strcat(fullfile(outdir,'fc_data'),filesep,'channel_wpli_',frequency,'_',ID,'_',task,'py.mat');
participant_channel_path_dpli = strcat(fullfile(outdir,'fc_data'),filesep,'channel_dpli_',frequency,'_',ID,'_',task,'py.mat');
% extract channels and data directly and save
channels = struct2cell(result_wpli.metadata.channels_location);
data = result_wpli.data.wpli;
save(participant_out_path_wpli,'data')
save(participant_channel_path_wpli,'channels')
% extract channels and data directly and save
channels = struct2cell(result_dpli.metadata.channels_location);
data = result_dpli.data.dpli;
save(participant_out_path_dpli,'data')
save(participant_channel_path_dpli,'channels')
%% Loop over Hemispheres and save images
for h = 1:length(hemispheres)
hemisphere = hemispheres(h);
% only keep the channels in the hemispere and reorder them
% according the indicated file
pattern_file = "biapt_egi129_" + hemisphere + ".csv";
%% reorder and save wPLI and dPLI
% reorder average wpli
data = result_wpli.data.avg_wpli;
channels = result_wpli.metadata.channels_location;
[ro_wpli, ro_w_channels, ro_w_regions] = filter_and_reorder_channels(data,channels,pattern_file);
% reorder average dpli
data = result_dpli.data.avg_dpli;
channels = result_dpli.metadata.channels_location;
[ro_dpli, ro_d_channels, ro_d_regions] = filter_and_reorder_channels(data,channels,pattern_file);
% define new outdir
outdir = fullfile(resultsfolder , ID, hemisphere);
mkdir(fullfile(outdir,'wPLI'));
mkdir(fullfile(outdir,'dPLI'));
plot_wPLI(ro_wpli, ID, frequency, task, hemisphere, fullfile(outdir,'wPLI'),ro_w_regions)
plot_dPLI(ro_dpli, ID, frequency, task, hemisphere, fullfile(outdir,'dPLI'),ro_d_regions)
if hemisphere == "Whole"
% initiate empty hub structre to fill in hub for every
% timepoint
timesteps = size(result_wpli.data.wpli);
timesteps = timesteps(1);
% create empty structure with hub weights (same nr of channesl as FC matrix)
result_hub_weights = zeros(timesteps, length(ro_d_channels));
% if we have the Smallest connected threshold we calculate it
% on the averaged dpli and apply it then to every time window
if string(threshold) == "SCT"
[finalthreshold] = find_smallest_connected_threshold(ro_wpli, threshold_range);
[b_wpli] = binarize_matrix(threshold_matrix(ro_wpli, finalthreshold));
elseif isfloat(threshold)
finalthreshold = threshold;
else
disp("Wrong Threshold selected. Can be either SCT or a value between 0 and 1 ");
end
for w=1:timesteps
%% calculate HUB at every window
disp(strcat("Participant: ", ID , "_HUB at window" , string(w) ));
data = squeeze(result_wpli.data.wpli(w,:,:));
channels = result_wpli.metadata.channels_location;
[ro_wpli, ro_w_channels, ro_w_regions] = filter_and_reorder_channels(data,channels,pattern_file);
[b_wpli] = binarize_matrix(threshold_matrix(ro_wpli, finalthreshold));
% here we are using only the degree and not the betweeness centrality
[~, hub_weights] = unnorm_binary_hub_location(b_wpli, ro_w_channels, 1.0, 0.0);
result_hub_weights(w,:) = hub_weights;
% do not normalize hub to z-score
% hub_norm_weights = (hub_weights - mean(hub_weights)) / std(hub_weights);
end
hub_avg = mean(result_hub_weights);
mkdir(fullfile(outdir,'HUB'));
plot_hub(hub_avg, ID, frequency, task, hemisphere, fullfile(outdir,'HUB'), ro_w_channels, finalthreshold)
end
end
end
end