forked from PINE-Lab/HAPPE
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathHAPPE_v3.m
1039 lines (970 loc) · 52.8 KB
/
HAPPE_v3.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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% HAPPE Version 3.0
%
% Developed at Northeastern University's PINE Lab
%
% For a detailed description of the pipeline and user options, please see
% the following manuscripts:
% Gabard-Durnam, et al., 2018 - HAPPE v1.
% Lopez, et al. (in press) - HAPPILEE
% Monachino, et al. (in press) - HAPPE+ER
% ...
%
% Contributors to HAPPE:
% Laurel Joy Gabard-Durnam ([email protected])
% Adriana S. Mendez Leal ([email protected])
% Carol L. Wilkinson ([email protected])
% April R. Levin ([email protected])
% ----
% Alexa D. Monachino ([email protected])
% Kelsie L. Lopez ([email protected])
%
% HAPPE includes code that is dependent on the following third-party
% software. Please reference this third-party software in any manuscripts
% making use of HAPPE as below:
% EEGLAB: A Delorme & S Makeig (2004) EEGLAB: an open source toolbox for
% analysis of single-trial EEG dynamics. Journal of Neuroscience
% Methods 134:9-21
% Cleanline by Tim Mullen (as an EEGLAB plug-in): Mullen, T. (2012).
% NITRC: CleanLine: Tool/Resource Info. Available online at:
% http://www.nitrc.org/projects/cleanline.
% FASTER segment-level channel interpolation code: Nolan, H., Whelan, R.,
% & Reilly, R.B. (2010). FASTER: Fully Automated Statistical Thresholding
% for EEG artifact Rejection. Journal of Neuroscience Methods, 192,
% 152-162.
% ICLabel EEGLAB plugin (for muscIL option only) - Pion-Tonachini et al.
% (2019). ICLabel: An automated electroencephalographic independent
% component classifier, dataset, and website. Neuroimage, 198, 181–197.
% REST EEGLAB plugin (for REST re-reference only): Li Dong*, Fali Li,
% Qiang Liu, Xin Wen, Yongxiu Lai, Peng Xu and Dezhong Yao*. MATLAB
% Toolboxes for Reference Electrode Standardization Technique (REST)
% of Scalp EEG. Frontiers in Neuroscience, 2017:11(601).
%
% HAPPE includes code adapted from third-party software. See the relevant
% function for further documentation:
% butterFilt.m adapted from ERPLAB (Lopez-Calderon & Luck, 2014)
% adaptedREST.m adapted from REST (Yao, 2001), (Dong et al., 2019)
%
% Any code that is not part of the third-party dependencies is released
% under the GNU General Public License version 3.
%
% Copyright 2018, 2021 Alexa Monachino, Kelsie Lopez, Laurel Gabard-Durnam
%
% This program is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License (version 3) as
% published by the Free Software Foundation.
%
% This software is being distributed with the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See GNU General
% Public License for more details.
%
% In no event shall Boston Children’s Hospital (BCH), the BCH Division of
% Developmental Medicine, the Laboratories of Cognitive Neuroscience (LCN),
% Northeastern University, the Center for Cognitive and Brain Health (CCBH),
% the Department of Psychology at Northeastern University, the Plasticity
% in Neurodevelopment (PINE) Lab, or software contributors be liable to any
% party for direct, indirect, special, incidental, or consequential damages,
% including lost profits, arising out of the use of this software and its
% documentation, even if any of the above listed entities and/or contributors
% have been advised of the possibility of such damage. Software and
% documentation is provided “as is.” The listed entities and/or software
% contributors are under no obligation to provide maintenance, support,
% updates, enhancements, or modifications.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% CLEAR THE WORKSPACE
clear ;
%% SET FOLDERS FOR HAPPE AND EEGLAB PATHS
fprintf('Preparing HAPPE...\n') ;
% SET HAPPE AND EEGLAB PATHS USING THE RUNNING SCRIPT
happeDir = fileparts(which(mfilename('fullpath'))) ;
eeglabDir = [happeDir filesep 'Packages' filesep 'eeglab2022.0'] ;
% ADD HAPPE AND REQUIRED FOLDERS TO PATH
addpath([happeDir filesep 'acquisition_layout_information'], ...
[happeDir filesep 'scripts'], ...
[happeDir filesep 'scripts' filesep 'UI_scripts'], ...
[happeDir filesep 'scripts' filesep 'pipeline_scripts'], ...
eeglabDir, genpath([eeglabDir filesep 'functions'])) ;
rmpath(genpath([eeglabDir filesep 'functions' filesep 'octavefunc'])) ;
% ADD EEGLAB FOLDERS TO PATH
pluginDir = dir([eeglabDir filesep 'plugins']) ;
pluginDir = strcat(eeglabDir, filesep, 'plugins', filesep, {pluginDir.name}, ';') ;
addpath([pluginDir{:}]) ;
% ADD CLEANLINE FOLDERS TO PATH
if exist('cleanline', 'file')
cleanlineDir = which('eegplugin_cleanline.m') ;
cleanlineDir = cleanlineDir(1:strfind(cleanlineDir, 'eegplugin_cleanline.m')-1) ;
addpath(genpath(cleanlineDir)) ;
else; error('Please make sure cleanline is on your path') ;
end
%% DETERMINE AND SET PATH TO DATA
% Use input from the command line to set the path to the data. If an
% invalid path is entered, repeat until a valid path is entered.
while true
srcDir = input('Enter the path to the folder containing the dataset(s):\n> ','s') ;
if exist(srcDir, 'dir') == 7; break ;
else; fprintf(['Invalid input: please enter the complete path to the ' ...
'folder containing the dataset(s).\n']) ;
end
end
cd (srcDir) ;
%% DETERMINE IF REPROCESSING DATA
[reprocess, rerunExt] = isReprocessed() ;
%% DETERMINE IF USING PRESET PARAMETERS
ver = '2_3_0' ;
[preExist, params, changeParams] = isPreExist(reprocess, ver) ;
%% SET PARAMETERS
params = setParams(params, preExist, reprocess, changeParams, happeDir) ;
%% SAVE INPUT PARAMETERS
% If created new or changed parameter set, save as a new .mat file to a new
% folder (input_parameters) added to the source folder.
if ~preExist || changeParams
% CREATE "input_parameters" FOLDER AND ADD IT TO PATH, unless it
% already exists.
if ~isfolder([srcDir filesep 'input_parameters'])
mkdir([srcDir filesep 'input_parameters']) ;
end
addpath('input_parameters') ;
cd input_parameters ;
% DETERMINE PARAMETER FILE NAME: Prompt to use a default or custom name
% for parameter file. If file exists, ask to create new file with a
% different name or overwrite existing file.
fprintf(['Parameter file save name:\n default = Default name (input' ...
'Parameters_dd-mm-yyyy.mat).\n custom = Create a custom file name' ...
'\n']) ;
if choose2('custom', 'default')
paramFile = paramFile_validateExist(['inputParameters_' ...
datestr(now, 'dd-mm-yyyy') '.mat'], 'inputParameters_', 2) ;
else
fprintf('File name (Do not include .mat):\n') ;
paramFile = paramFile_validateExist([input('> ', 's') '.mat'], ...
'inputParameters_', 0) ;
end
% SAVE PARAMETERS: Save the params variable to a .mat file using the
% name created above.
params.HAPPEver = ver ;
fprintf('Saving parameters...') ;
save(paramFile, 'params') ;
fprintf('Parameters saved.') ;
end
%% CREATE OUTPUT FOLDERS
% Create the folders in which to store intermediate and final outputs.
% If you aren't running ERP analyses, don't add the ERP_filtered folder.
cd(srcDir) ;
fprintf('Creating output folders...\n') ;
allDirNames = {'intermediate_processing', 'wavelet_cleaned_continuous', ...
'muscIL', 'ERP_filtered', 'segmenting', 'processed', ...
'quality_assessment_outputs'} ;
if ~params.paradigm.ERP.on; allDirNames(ismember(allDirNames, 'ERP_filtered')) = []; end
if ~params.muscIL; allDirNames(ismember(allDirNames, 'muscIL')) = []; end
dirNames = cell(1,size(allDirNames,2)) ;
for i=1:length(allDirNames)
dirNames{i} = [num2str(i) ' - ' allDirNames{i}] ;
if ~isfolder([srcDir filesep num2str(i) ' - ' allDirNames{i}])
mkdir([srcDir filesep num2str(i) ' - ' allDirNames{i}]) ;
end
end
clear('allDirNames') ;
fprintf('Output folders created.\n') ;
%% COLLECT DATA TO RUN
fprintf('Gathering files...\n') ;
cd(srcDir) ;
% DETERMINE FILE EXTENSION USING USER LOAD INFO
if params.loadInfo.inputFormat == 1; inputExt = '.mat' ;
elseif params.loadInfo.inputFormat == 2; inputExt = '.raw' ;
elseif params.loadInfo.inputFormat == 3; inputExt = '.set' ;
elseif params.loadInfo.inputFormat == 4; inputExt = '.cdt' ;
elseif params.loadInfo.inputFormat == 5; inputExt = '.mff' ;
end
% COLLECT FILE NAMES: gather the names of all the files with the relevant
% file extension in the directory indicated by the user.
FileNames = {dir(['*' inputExt]).name} ;
% LOCATE STIM FILE AND NAMES
if params.loadInfo.inputFormat == 1 && params.paradigm.task
% ***
end
%% INITIALIZE QUALITY REPORT METRICS
fprintf('Initializing report metrics...\n') ;
% DATA QUALITY METRICS: create a variable holding all the names of each
% metric and a variable to hold the metrics for each file.
dataQCnames = {'File_Length_in_Seconds', 'Number_User-Selected_Chans', ...
'Number_Good_Chans_Selected', 'Percent_Good_Chans_Selected', 'Bad_Chan_IDs', ...
'Percent_Var_Retained_Post-Wav', 'Number_ICs_Rej', 'Percent_ICs_Rej', ...
'Chans_Interpolated_per_Seg', 'Number_Segs_Pre-Seg_Rej', ...
'Number_Segs_Post-Seg_Rej', 'Percent_Segs_Post-Seg_Rej'} ;
dataQC = cell(length(FileNames), length(dataQCnames)) ;
% If processing for tasks, create an additional variable to hold specific
% data metrics for each onset tag.
if params.paradigm.task
dataQC_task = cell(length(FileNames), length(params.paradigm.onsetTags)*3) ;
dataQCnames_task = cell(1, length(params.paradigm.onsetTags)*3) ;
for i=1:size(params.paradigm.onsetTags,2)
dataQCnames_task{i*3-2} = ['Number_' params.paradigm.onsetTags{i} ...
'_Segs_Pre-Seg_Rej'] ;
dataQCnames_task{i*3-1} = ['Number_' params.paradigm.onsetTags{i} ...
'_Segs_Post-Seg_Rej'] ;
dataQCnames_task{i*3} = ['Percent_' params.paradigm.onsetTags{i} ...
'_Segs_Post-Seg_Rej'] ;
end
% If grouping any tags by condition, create additional variable to hold
% specific data metrics for each condition.
if params.paradigm.conds.on
dataQC_conds = cell(length(FileNames), ...
size(params.paradigm.conds.groups,1)*3) ;
dataQCnames_conds = cell(1, size(params.paradigm.conds.groups,1)*3) ;
for i = 1:size(params.paradigm.conds.groups,1)
dataQCnames_conds{i*3-2} = ['Number_' params.paradigm.conds.groups{i, ...
1} '_Segs_Pre-Seg_Rej'] ;
dataQCnames_conds{i*3-1} = ['Number_' params.paradigm.conds.groups{i, ...
1} '_Segs_Post-Seg_Rej'] ;
dataQCnames_conds{i*3} = ['Percent_' params.paradigm.conds.groups{i, ...
1} '_Segs_Post-Seg_Rej'] ;
end
end
end
errorLog = {} ;
% PIPELINE QUALITY METRICS: initialize the variables to hold the pipeline
% quality metrics for assessing HAPPE performance on line noise reduction
% and waveleting.
if ~reprocess; lnMeans = [] ; wavMeans = [] ; end
if params.muscIL; icaMeans = [] ; end
%% RUN THE PROCESSING PIPELINE OVER EACH FILE
for currFile = 1:length(FileNames)
cd(srcDir) ;
try
if ~reprocess
%% LOAD AND VALIDATE RAW FILE
fprintf(['Loading ' FileNames{currFile} '...\n']) ;
try
% .MAT FILES
if params.loadInfo.inputFormat == 1
% Net Station:
if params.loadInfo.NSformat
load(FileNames{currFile}) ;
% Use known information about Net Station sampling
% rate variable names to determine the sampling
% rate of the file.
srate = intersect(who, {'samplingRate', 'EEGSamplingRate'}) ;
srate = double(eval(srate{1})) ;
% Use the potential variable names given by the user to
% determine the data variable.
eegVName = intersect(who, params.loadInfo.NSvarNames) ;
% Try to load the EEG using an EGI Net Station EEGLAB
% plugin. If unable to load the file, throws an error.
% Specifically looks for MATLAB:badsubscript.
try EEGraw = pop_importegimat(FileNames{currFile}, ...
srate, 0.00, eegVName{1}) ;
catch ME
if strcmp(ME.identifier, 'MATLAB:badsubscript')
fprintf(['Sorry, could not read the variable' ...
' name of the EEG data. Please check your' ...
' file.\n']) ;
end
rethrow(ME) ;
end
% Matrix:
else
% If all the files have the same sampling rate, use
% the single user specified value. Otherwise, load
% the table of sampling rates specified by the user
% and locate the corresponding value from the table
% using the file name.
if params.loadInfo.srate.same
srate = params.loadInfo.srate.val ;
else
srateTable = table2cell(readtable(params.loadInfo.srate.file)) ;
srate = srateTable{contains(list, ...
FileNames{currFile}), 2} ;
end
% Use the built-in EEGLAB function to load the
% data. Uses the sampling rate determined above and
% the user-specified channel locations. If the user
% did not specify any channel locations, is still
% able to load successfully.
EEGraw = pop_importdata('dataformat', 'matlab', ...
'data', FileNames{currFile}, 'srate', srate, ...
'chanlocs', params.loadInfo.chanlocs.file) ;
end
% If task data, import event information using the file
% indicated by the user.
if params.paradigm.task
cd(params.loadInfo.eventLoc) ;
EEGraw = pop_importevent(EEGloaded, 'append', 'no', ...
'event', StimNames{currFile}, 'fields', {'type', ...
'latency', 'status'}, 'skipline', 1, 'timeunit', ...
1E-3, 'align', NaN) ;
origEvents = EEGraw.urevent ;
events = EEGraw.event ;
end
else
% ---------------------------------------------------------
% .RAW FILES
if params.loadInfo.inputFormat == 2
EEGraw = pop_readegi(FileNames{currFile}, [], [], ...
'auto') ;
% ---------------------------------------------------------
% .SET FILES
elseif params.loadInfo.inputFormat == 3
EEGraw = load('-mat', FileNames{currFile}) ;
if isfield(EEGraw, 'EEG'); EEGraw = EEGraw.EEG; end
% ---------------------------------------------------------
% .CDT FILES
elseif params.loadInfo.inputFormat == 4
EEGraw = loadcurry([srcDir filesep ...
FileNames{currFile}], 'CurryLocation', 'false') ;
% ---------------------------------------------------------
% .MFF FILES
elseif params.loadInfo.inputFormat == 5
EEGraw = pop_mffimport(FileNames{currFile} , ...
params.loadInfo.typeFields, 0, 0) ;
end
% ---------------------------------------------------------
% Determine the sampling rate from the loaded EEG. If task
% related, save the events from the loaded file in its own
% variable.
srate = double(EEGraw.srate) ;
if params.paradigm.task; events = EEGraw.event ; end
end
% ---------------------------------------------------------
% Validate the loaded file using EEGLAB's checkset function
fprintf('Validating file...\n') ;
EEGraw.setname = 'rawEEG' ;
EEG = eeg_checkset(EEGraw) ;
dataQC{currFile, 1} = EEG.xmax ;
% If unable to load the file, indicate this to the user via the
% command line and throw an error.
catch ME
error('HAPPE:LoadFail', ...
['ERROR: Unable to load ' FileNames{currFile} '\n']) ;
end
%% CORRECT CHANNELS, IF NEEDED
% Channels need to be corrected for EGI files, regardless of
% input format. Please note that for .raw files, if exported
% incorrectly from netstation, this will not work!
if params.loadInfo.correctEGI
fprintf('Correcting Channels...\n') ;
% LOAD THE CHANNEL LOCATIONS FILE USING THE LAYOUT
% INFORMATION PARAMETERS
if params.loadInfo.layout(1) == 1
if params.loadInfo.layout(2) == 64
params.loadInfo.chanlocs.file = [] ;
params.loadInfo.chanlocs.file = [happeDir filesep ...
'acquisition_layout_information' filesep ...
'GSN65v2_0.sfp'] ;
else; error('Invalid sensor layout selection.') ;
end
elseif params.loadInfo.layout(1) == 2
if ismember(params.loadInfo.layout(2), [32, 64, 128, 256])
params.loadInfo.chanlocs.file = [happeDir filesep ...
'acquisition_layout_information' filesep ...
'GSN-HydroCel-' num2str(params.loadInfo.layout(2)+1) ...
'.sfp'] ;
else; error('Invalid sensor layout selection.') ;
end
end
% CORRECT THE CHANNELS USING EEGLAB FUNCTIONALITY
EEG = pop_chanedit(EEG, EEG.chanlocs, 'load', ...
{params.loadInfo.chanlocs.file 'filetype' 'autodetect'}) ;
% VALIDATE THAT THE EEG WAS CORRECTED CORRECTLY
EEG = eeg_checkset(EEG) ;
end
%% SET (AND CORRECT) THE REFERENCE CHANNEL
% If the reference channel is included in the data, remove it
% from the data channels. This will prevent it from being
% falsely flagged as a bad channel.
if params.loadInfo.chanlocs.inc && params.reref.on && ...
params.reref.flat
if ismember(params.reref.chan, {EEG.chanlocs.labels})
indx = find(strcmpi({EEG.chanlocs.labels}, ...
params.reref.chan)) ;
EEG.data(indx,:) = [] ;
EEG.nbchan = EEG.nbchan - 1;
EEG.chaninfo.refChan = EEG.chanlocs(indx) ;
EEG.chanlocs(indx) = [] ;
end
end
%% UPDATE CHANNEL IDS AND FILTER TO THE CHANNELS OF INTEREST
% IF CHANNEL LOCATIONS ARE INCLUDED, update the channel IDs in
% the EEG structure to match the channels of interest.
if params.loadInfo.chanlocs.inc
% Updates chanIDs (a variable holding the channels of
% interest for the current subject) using information
% provided by the user. If for some reason one or more
% channels are missing in a file, this will protect against
% the channel IDs being impacted in other files.
fprintf('Selecting channels of interest...\n') ;
if strcmpi(params.chans.subset, 'coi_exclude')
chanIDs = setdiff({EEG.chanlocs.labels}, params.chans.IDs) ;
elseif strcmpi(params.chans.subset, 'coi_include')
chanIDs = intersect({EEG.chanlocs.labels}, params.chans.IDs) ;
elseif strcmpi(params.chans.subset, 'all')
chanIDs = {EEG.chanlocs.labels} ;
end
% Filter to the channels of interest
fprintf('Filtering to channels of interest...\n') ;
EEG = pop_select(EEG, 'channel', chanIDs) ;
EEG.setname = [EEG.setname '_cs'] ;
% Report how many channels are present in the dataset as a
% data metric.
dataQC{currFile,2} = size(chanIDs,2) ;
% IF NO CHANNEL LOCATIONS ARE INCLUDED, report how many
% channels are present in the dataset as a data metric.
else
dataQC{currFile,2} = size(EEG.data,1) ;
chanIDs = 1:size(EEG.data,1) ;
end
%% REDUCE LINE NOISE
% Attempt to reduce line noise using the happe_reduceLN
% function (see function for further documentation). If HAPPE
% fails during this step, alert the user that processing failed
% during the line noise step in the command window and rethrow
% the error to proceed to the next file.
try
[EEG, lnMeans] = happe_reduceLN(EEG, params.lineNoise, ...
lnMeans) ;
EEG.setname = [EEG.setname '_ln'] ;
catch ME
fprintf(2, 'ERROR: Line noise reduction failed.\n') ;
rethrow(ME) ;
end
%% RESAMPLE DATA
% If the user enabled downsampling, resample the data to the
% user-specified frequency using EEGLAB functionality.
if params.downsample
fprintf(['Resampling the data to ' num2str(params.downsample) ...
' Hz...\n']) ;
EEG = pop_resample(EEG, params.downsample) ;
end
%% FILTER AND SAVE INTERMEDIATE FILE
% Enter the intermediate_processing folder
cd([srcDir filesep dirNames{contains(dirNames, ...
'intermediate_processing')}]) ;
% Filter the EEG data. For ERP paradigms, only use a 100 Hz
% bandpass. For all other paradigms, filter the data using both
% a 1 Hz highpass and a 100 Hz lowpass.
if params.paradigm.ERP.on; EEG = pop_eegfiltnew(EEG, [], 100, ...
[], 0, [], 0) ;
else; EEG = pop_eegfiltnew(EEG, 1, 100, [], 0, [], 0) ;
end
EEG.setname = [EEG.setname '_filt'] ;
% Save the filtered dataset to the intermediate_processing
% folder
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, '_filtered_lnreduced.set')) ;
%% DETECT BAD CHANNELS
% If bad channel detection is on, detect and reject bad
% channels using the happe_detectBadChans function (see
% function for further documentation). Regardless if detecting
% bad channels, update the number of channels included in
% processing (good channels), the percent of good channels
% selected (100% if bad channel detection is disabled), and the
% list of channels, if any, marked as bad.
origChans = EEG.chanlocs ;
if params.badChans.rej
[EEG, dataQC] = happe_detectBadChans(EEG, params, dataQC, ...
chanIDs, currFile) ;
EEG.setname = [EEG.setname 'badc'] ;
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, '_badchanrej.set')) ;
else
dataQC{currFile,3} = size(chanIDs,2) ;
dataQC{currFile,4} = dataQC{currFile,3}/dataQC{currFile,2}*100;
dataQC{currFile,5} = 'NA' ;
end
%% WAVELET THRESHOLDING
% Attempt to wavelet threshold using the happe_wavThresh
% function (see function for further documentation). If HAPPE
% fails during this step, alert the user that processing failed
% during wavelet thresholding in the command window and rethrow
% the error to proceed to the next file.
cd([srcDir filesep dirNames{contains(dirNames, ...
'wavelet_cleaned_continuous')}]) ;
try [EEG, wavMeans, dataQC] = happe_wavThresh(EEG, params, ...
wavMeans, dataQC, currFile) ;
catch ME
fprintf(2, 'ERROR: Wavelet thresholding failed.\n') ;
rethrow(ME) ;
end
% Save the wavelet-thresholded EEG as an intermediate output.
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, '_wavclean.set')) ;
else
%% LOAD AND VALIDATE PROCESSED FILE
try
fprintf(['Loading ' FileNames{currFile} '...\n']) ;
% LOAD FILTERED & LINE NOISE REDUCED FILE:
% Use the raw file name. Collect the original number of
% user-selected channels from this file, then clear the
% dataset.
EEGloaded = pop_loadset('filename', ...
strrep(FileNames{currFile}, inputExt, ...
'_filtered_lnreduced.set'), 'filepath', [srcDir filesep ...
dirNames{contains(dirNames, 'intermediate_processing')}]) ;
dataQC{currFile, 2} = size(EEGloaded.data,1) ;
origChans = EEGloaded.chanlocs ;
clear('EEGloaded') ;
% LOAD PRE-SEGMENTED PROCESSED DATA:
% Load the wavelet-thresholded data. Using the raw file to
% locate and load the correct file.
EEGloaded = pop_loadset('filename', strrep(FileNames{currFile}, ...
inputExt, '_wavclean.set'), 'filepath', [srcDir ...
filesep dirNames{contains(dirNames, ...
'wavelet_cleaned_continuous')}]) ;
% VALIDATE LOADED FILE
EEG = eeg_checkset(EEGloaded) ;
% FILL IN DATA QUALITY METRICS:
% Using the loaded data, fill in the following data quality
% metrics about the current file: length of the EEG
% recording, the number of good channels, and the percent
% of good channels.
dataQC{currFile,1} = EEG.xmax ;
dataQC{currFile,3} = size(EEG.data,1) ;
dataQC{currFile,4} = dataQC{currFile,3}/dataQC{currFile,2}*100 ;
dataQC{currFile,5} = 'NA' ;
dataQC{currFile,6} = 'NA' ;
% If unable to load the file, indicate this to the user via the
% command line and throw an error.
catch ME
error('HAPPE:LoadFail', ['ERROR: Unable to load ' ...
FileNames{currFile} '\n']) ;
end
end
%% MUSCIL
if ~params.paradigm.ERP.on && params.muscIL
cd([srcDir filesep dirNames{contains(dirNames, 'muscIL')}]) ;
try [EEG, dataQC] = happe_muscIL(EEG, dataQC, FileNames, ...
currFile, inputExt) ;
catch ME; rethrow(ME) ;
end
else
dataQC{currFile, 7} = 'NA' ;
dataQC{currFile, 8} = 'NA' ;
end
%% FILTER USING ERP CUTOFFS AND SAVE
% If performing preprocessing for ERPs, attempts to filter
% using the user-defined ERP cutoffs. If HAPPE fails during
% this step, alert the user that processing failed during
% filtering to ERP cutoffs in the command window and rethrow
% the error to proceed to the next file.
if params.paradigm.ERP.on
fprintf('Filtering using ERP cutoffs...\n') ;
cd([srcDir filesep dirNames{contains(dirNames, ...
'ERP_filtered')}]) ;
try
% If the user indicated to use ERPLAB's butterworth
% filter, filter using the butterFilt function (adapted
% from ERPLAB's functions - see butterFilt
% documentation for more information). Currently only
% available for ERP paradigms.
if params.filt.butter
EEG = butterFilt(EEG, params.paradigm.ERP.lowpass, ...
params.paradigm.ERP.highpass) ;
% If the user indicated to use EEGLAB's FIR filter,
% filter using the EEGLAB function.
else
EEG = pop_eegfiltnew(EEG, params.paradigm.ERP.highpass, ...
params.paradigm.ERP.lowpass, [], 0, [], 0) ;
end
EEG.setname = [EEG.setname '_forERP'] ;
catch ME
fprintf(2, 'ERROR: Filtering to ERP cutoffs failed.\n') ;
rethrow(ME) ;
end
% Save the filtered EEG as an intermediate output.
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, ['_wavcleaned_filteredERP' rerunExt '.set'])) ;
end
%% RE-ADD FLATLINE RE-REFERENCE
% If there was a flatline channel, return it to the data. If HAPPE
% is unable to find the designated channel, throw an error to
% proceed to the next file and alert the user via the command
% window.
if params.loadInfo.chanlocs.inc && params.reref.on && params.reref.flat
try
EEG.data(EEG.nbchan+1,:) = 0 ;
EEG.nbchan = EEG.nbchan + 1 ;
EEG.chanlocs(EEG.nbchan) = EEG.chaninfo.refChan ;
origChans(length(origChans)+1) = EEG.chaninfo.refChan ;
catch ME
fprintf(2, 'ERROR: No flatline channel detected.') ;
rethrow(ME) ;
end
end
%% SEGMENT DATA
% If the data is not already segmented, use the user-specified
% parameters to segment the data. If it is segmented, alert the
% user via the command window that the data could not be segmented.
cd([srcDir filesep dirNames{contains(dirNames, 'segmenting')}]) ;
if params.segment.on
try EEG = happe_segment(EEG, params) ;
% SAVE THE SEGMENTED DATASET AS AN INTERMEDIATE FILE
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, inputExt, ...
['_segmented' rerunExt '.set'])) ;
catch ME
fprintf(2, 'ERROR: Segmenting failed.\n') ;
rethrow(ME) ;
end
end
% Add the number of segments to the dataQC matrix. In the case of
% ERP paradigms, add the number of segments for each event tag and
% each condition to the dataQC_task and dataQC_conds matrices.
dataQC{currFile, 10} = EEG.trials ;
if params.paradigm.task
% Segments per Event Tag:
for i=1:length(params.paradigm.onsetTags)
try dataQC_task{currFile, i*3-2} = length(pop_selectevent(EEG, ...
'type', params.paradigm.onsetTags{i}).epoch) ;
catch
dataQC_task{currFile, i*3-2} = 0 ;
fprintf('No instances of "%s" in this file.\n', ...
params.paradigm.onsetTags{i}) ;
end
end
% Segments per Condition:
if params.paradigm.conds.on
for i=1:size(params.paradigm.conds.groups,1)
try
currGroup = params.paradigm.conds.groups(i,2:end) ;
currGroup = currGroup(~cellfun('isempty',currGroup)) ;
dataQC_conds{currFile, i*3-2} = length(pop_selectevent(EEG, ...
'type', currGroup).epoch) ;
catch
dataQC_conds{currFile, i*3-2} = 0 ;
end
end
end
end
%% BASELINE CORRECT AND SAVE (ERPs)
% If pre-processing ERP data and baseline correction is enabled,
% baseline correct the data using the baseline window specified by
% the user.
if params.paradigm.ERP.on && params.baseCorr.on
fprintf('Correcting baseline...\n') ;
try EEG = pop_rmbase(EEG, [params.baseCorr.start ...
params.baseCorr.end]) ;
% SAVE BASELINE CORRECTED DATA AS AN INTERMEDIATE FILE
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, inputExt, ...
['_segmented_BLcorr' rerunExt '.set'])) ;
catch ME
fprintf(2, 'ERROR: Baseline correction failed.\n') ;
rethrow(ME) ;
end
end
%% INTERPOLATE BAD DATA USING CHANNELS
% If enabled, interpolate data within segments by channels. Uses
% code adapted from the FASTER program (Nolan et al., 2010),
% interpolating channels scoring above/below z threshold of 3 for a
% segment. Not available if channels are not included or if bad
% channel detection was not performed. Adds the segments and
% interpolated channels, if any, to the dataquality metrics.
if params.segment.interp
try [EEG, dataQC] = happe_interpChanBySeg(EEG, dataQC, ...
currFile) ;
% SAVE INTERPOLATED DATA AS INTERMEDIATE FILE
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, inputExt, ...
['_segmented_interp' rerunExt '.set'])) ;
catch ME
fprintf(2, 'ERROR: Interpolation within segments failed.\n') ;
rethrow(ME) ;
end
else; dataQC{currFile, 9} = 'NA' ;
end
%% REJECT SEGMENTS
% If the data is segmented and segment rejection is enabled, reject
% segments containing artifact according to the criteria as set by
% the user. Can reject segments based on a region of interest or by
% the entire set of channels present in the data. If the data is
% continuous or is a single trial, cannot reject segments.
if params.segRej.on && EEG.trials > 1
try EEG = happe_segRej(EEG, params) ;
% SAVE FILE WITH SEGMENTS REJECTED
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, ['_segments_postrej' rerunExt '.set'])) ;
catch ME; rethrow(ME) ;
end
elseif params.segRej.on && EEG.trials == 1
error('HAPPE:rejOneSeg', ['ERROR: Cannot reject segments from' ...
' data with a single epoch.']) ;
end
% UPDATE DATA QUALITY METRICS: Add the number and percent of trials
% retained post-segment rejection.
dataQC{currFile, 11} = EEG.trials ;
dataQC{currFile, 12} = dataQC{currFile, 11}/dataQC{currFile, 10}*100 ;
%% INTERPOLATE BAD CHANNELS
% If the file has channel locations, interpolate the bad channels
% to enable proper re-referencing.
try
if params.loadInfo.chanlocs.inc
fprintf('Interpolating bad channels...\n') ;
EEG = pop_interp(EEG, origChans, 'spherical') ;
EEG.setname = [EEG.setname '_interp'] ;
end
catch ME
fprintf(2, 'ERROR: Bad channel interpolation failed.\n') ;
rethrow(ME) ;
end
%% RE-REFERENCE DATA
% If enabled, re-reference the data. Will perform an average
% rereference or re-reference to a subset based on the parameters
% set by the user. If a reference electrode is present in the data,
% will keep it.
if params.reref.on
fprintf('Re-referencing...\n') ;
if strcmpi(params.reref.method, 'average')
EEG = pop_reref(EEG, [], 'keepref', 'on') ;
elseif strcmpi(params.reref.method, 'subset')
EEG = pop_reref(EEG, intersect({EEG.chanlocs.labels}, ...
params.reref.subset, 'stable'), 'keepref', 'on') ;
EEG.setname = [EEG.setname '_reref'] ;
elseif strcmpi(params.reref.method, 'rest')
if strcmpi(params.reref.rest.ref, 'mastoid')
params.reref.rest.chanL = find(ismember({EEG.chanlocs.labels}, ...
params.reref.rest.chanL)) ;
params.reref.rest.chanR = find(ismember({EEG.chanlocs.labels}, ...
params.reref.rest.chanR)) ;
else
params.reref.rest.chanlist = 1:size({EEG.chanlocs.labels},2) ;
end
EEG = adaptedREST(EEG, params.reref.rest) ;
end
end
%% SPLIT EEG BY EVENT TAGS AND CONDITIONS
% In the case of multiple event tags in a dataset, create a unique
% EEG structure containing every instance of a single event tag
% for each event tag. If a tag does not appear in the original EEG
% dataset, alerts the user via the command window. Does the same
% per condition. Additionally updates the data quality metrics
% specific to tags and conditions with the number and percent
% trials retained.
if params.paradigm.task && size(params.paradigm.onsetTags,2) > 1
% Split by Event Tags:
fprintf('Creating EEGs by event tags...\n') ;
eegByTags = cell(1,size(params.paradigm.onsetTags,2)) ;
for i=1:size(params.paradigm.onsetTags,2)
try eegByTags{i} = pop_selectevent(EEG, 'type', ...
params.paradigm.onsetTags{i}) ;
dataQC_task{currFile, i*3-1} = length(eegByTags{i}.epoch) ;
dataQC_task{currFile, i*3} = dataQC_task{currFile, i*3-1} ...
/ dataQC_task{currFile, i*3-2}*100 ;
catch
fprintf('No instances of tag %s appear in this file.\n', ...
params.paradigm.onsetTags{i}) ;
dataQC_task{currFile, i*3-1} = 0 ;
dataQC_task{currFile, i*3} = 'NA' ;
end
end
% Split by Conditions:
if params.paradigm.conds.on
eegByConds = cell(1, size(params.paradigm.conds.groups,1)) ;
for i=1:size(params.paradigm.conds.groups,1)
eegByConds{i} = pop_selectevent(EEG, 'type', ...
params.paradigm.conds.groups(i,2:end)) ;
dataQC_conds{currFile, i*3-1} = length(eegByConds(i).epoch) ;
dataQC_conds{currFile, i*3} = dataQC_conds{currFile, i*3-1} ...
/ dataQC_conds{currFile, i*3-2}*100 ;
end
end
end
%% GENERATE VISUALIZATIONS
% If visualizations are enabled, generate topoplots for the current
% file. If preprocessing for ERPs, this will be the processed ERP
% spectrum.
cd([srcDir filesep dirNames{contains(dirNames, 'processed')}]) ;
if params.vis.enabled
if params.paradigm.ERP.on
figure; pop_timtopo(EEG, [params.vis.min params.vis.max], ...
params.vis.toPlot) ;
saveas(gcf, strrep(FileNames{currFile}, inputExt, ...
['_processedERPspectrum' rerunExt '.jpg'])) ;
else
figure; pop_spectopo(EEG, 1, [], 'EEG', 'freq', ...
[params.vis.toPlot], 'freqrange', [[params.vis.min] ...
[params.vis.max]], 'electrodes', 'off') ;
saveas(gcf, strrep(FileNames{currFile}, inputExt, ...
['_processedSpectrum' rerunExt '.jpg'])) ;
end
end
%% SAVE PRE-PROCESSED DATASET(S)
fprintf('Saving pre-processed dataset(s)...\n') ;
if params.outputFormat == 1 || params.outputFormat == 3
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, ['_processed' rerunExt '.set'])) ;
if params.paradigm.task && size(params.paradigm.onsetTags,2) > 1
for i=1:size(eegByTags, 2)
pop_saveset(EEG, 'filename', strrep(FileNames{currFile}, ...
inputExt, ['_processed_' params.paradigm.onsetTags{i} ...
'_' rerunExt '.set'])) ;
end
if params.paradigm.conds.on
for i=1:size(eegByConds,2)
pop_saveset(eegByConds{i}, 'filename', ...
strrep(FileNames{currFile}, inputExt, ...
['_' params.paradigm.conds.groups{i,1} ...
'_processed' rerunExt '.set'])) ;
end
end
end
end
if params.outputFormat == 2
save(strrep(FileNames{currFile}, inputExt, ['_processed' ...
rerunExt '.mat']), 'EEG') ;
if params.paradigm.task.on && size(params.paradigm.onsetTags, ...
2) > 1
for i=1:size(eegByTags, 2)
currEEG = eegByTags(i) ;
save(strrep(FileNames{currFile}, inputExt, ...
['_processed_' params.paradigm.onsetTags{i} ...
'_' rerunExt '.mat']), 'currEEG') ;
end
if params.paradigm.conds.on
for i=1:size(eegByConds,2)
save(strrep(FileNames{currFile}, inputExt, ...
['_' params.paradigm.conds.groups{i,1} ...
'_processed' rerunExt '.mat']), 'currEEG') ;
end
end
end
end
if params.outputFormat == 1
pop_export(EEG, strrep(FileNames{currFile}, inputExt, ...
['_processed_IndivTrial' rerunExt '.txt']), ...
'transpose', 'on', 'precision', 8) ;
pop_export(EEG, strrep(FileNames{currFile}, inputExt, ...
['_processed_AveOverTrials' rerunExt '.txt']), ...
'transpose', 'on', 'erp', 'on', 'precision', 8) ;
if size(params.paradigm.onsetTags, 2) > 1
for i=1:size(eegByTags, 2)
pop_export(eegByTags{i}, strrep(FileNames{currFile}, ...
inputExt, ['_processed_IndivTrial' ...
params.paradigm.onsetTags{i} rerunExt '.txt']), ...
'transpose', 'on', 'precision', 8) ;
pop_export(eegByTags{i}, strrep(FileNames{currFile}, ...
inputExt, ['_processed_AveOverTrials' ...
params.paradigm.onsetTags{i} rerunExt '.txt']), ...
'transpose', 'on', 'erp', 'on', 'precision', 8) ;
end
if params.paradigm.conds.on
for i=1:size(eegByConds,2)
pop_export(eegByConds{i}, strrep(FileNames{currFile}, ...
inputExt, ['_processed_IndivTrial' ...
params.paradigm.conds.groups{i,1} rerunExt, ...
'.txt']), 'transpose', 'on', 'precision', 8) ;
pop_export(eegByConds{i}, strrep(FileNames{currFile}, ...
inputExt, ['_processed_AveOverTrials' ...
params.paradigm.conds.groups{i,1} rerunExt, ...
'.txt']), 'transpose', 'on', 'erp', 'on', ...
'precision', 8) ;
end
end
end
end
%% ERRORS
% If HAPPE errors out, check first for common failures, and indicate
% the cause in the data and pipeline quality assessments. If an
% uncommon error, will simply fill with error. This allows HAPPE to
% continue to run over the whole dataset.
catch ME
% ADD ERROR TO THE ERROR LOG:
errorLog = [errorLog; {FileNames{currFile}, ME.message}] ; %#ok<AGROW>
% CHECK FOR COMMON ERRORS:
if strcmp(ME.identifier, 'HAPPE:LoadFail'); fill = 'LOAD_FAIL' ;
elseif strcmp(ME.identifier, 'HAPPE:allICsRej'); fill = 'ALL_IC_REJ' ;
elseif strcmp(ME.identifier, 'HAPPE:noTags'); fill = 'NO_TAGS' ;
elseif strcmp(ME.identifier, 'HAPPE:rejOneSeg'); fill = 'REJ_ONE_SEG' ;
elseif strcmp(ME.identifier, 'HAPPE:AllTrialsRej'); fill = 'ALL_SEG_REJ' ;
else; fill = 'ERROR' ;
end
% FILL DATA QC MATRICES WITH ERROR INDICATOR
for i=1:size(dataQC, 2)
if isempty(dataQC{currFile,i}); dataQC{currFile,i} = fill ; end
end
if params.paradigm.ERP.on
for i=1:size(dataQC_task, 2)
if isempty(dataQC_task{currFile, i})
dataQC_task{currFile, i} = fill ;
end
end
end
% FILL PIPELINE QC MATRICES WITH ERROR INDICATOR
if ~reprocess
lnMeans = qm_ERROR(lnMeans, 1+length(params.lineNoise.neighbors) + ...
length(params.lineNoise.harms.freqs), currFile) ;
wavMeans = qm_ERROR(wavMeans, 5+length(params.QCfreqs), ...
currFile) ;
end
end
end
%% GENERATE OUTPUT TABLES
fprintf('Generating quality assessment outputs...\n') ;
cd([srcDir filesep dirNames{contains(dirNames, ...
'quality_assessment_outputs')}]) ;
rmpath(genpath(cleanlineDir)) ;
try
% CREATE AND SAVE PIPELINE QUALITY ASSESSMENT
if ~reprocess
% Create line noise reduction names.
lnNames = {'r all freqs pre/post linenoise reduction'} ;
for i=2:size(params.lineNoise.neighbors, 2)+1
lnNames{i} = ['r ' num2str(params.lineNoise.neighbors(i-1)) ...
' hz pre/post linenoise reduction'] ;
end
for i=1:size(params.lineNoise.harms.freqs, 2)
lnNames{i+size(params.lineNoise.neighbors, 2)+1} = ['r ' num2str(params.lineNoise.harms.freqs) ...
' hz pre/post harmonic reduction'] ;
end
% Create wavelet thresholding names.
wavNames = {'RMSE post/pre waveleting', 'MAE post/pre waveleting', ...
'SNR post/pre waveleting', 'PeakSNR post/pre waveleting', ...
'r alldata post/pre waveleting'} ;
for i=1:size(params.QCfreqs,2)
wavNames{i+5} = ['r ' num2str(params.QCfreqs(i)) ' hz post/pre ' ...
'waveleting'] ;
end
% Concat the Names and QC matrices.
pipelineQC_names = [(lnNames) (wavNames)] ;
pipelineQC = [(lnMeans) (wavMeans)] ;
% Save the pipeline QC table.
pipelineQC_saveName = helpName(['HAPPE_pipelineQC' rerunExt '_' ...
datestr(now, 'dd-mm-yyyy') '.csv']) ;
writetable(array2table(pipelineQC, 'VariableNames', pipelineQC_names, ...
'RowNames', FileNames), pipelineQC_saveName, 'WriteRowNames', ...
true, 'QuoteStrings', true);
end
% CREATE AND SAVE DATA QUALITY ASSESSMENT
% Concat Names and QC matrices according to the presence or absence of
% multiple onset tags and conditions.