Laboratorium EEG/EEGLAB

Laboratorium_EEG/EEGLAB

Zajęcia bazują na podręczniku praktycznym użytkowania pakietu EEGLAB dostępnym na stronach:

• rozdziały 1-8 z części I: http://sccn.ucsd.edu/wiki/EEGLAB#The_EEGLAB_Tutorial_Outline

• rozdziały z części: http://sccn.ucsd.edu/wiki/Chapter_02:_Writing_EEGLAB_Scripts

Analiza potencjałów wywołanych

Wprowadzenie

W tej części zajęć naszym celem jest przećwiczenie pracy z danymi z eksperymentu w którym badamy potencjały wywołane na grupie osób. Dane pochodzą z eksperymentu opisanego w artykule: Effects of Valence and Origin of Emotions in Word Processing Evidenced by Event Related Potential Correlates in a Lexical Decision Task

Skupimy się na badaniu różnic pomiędzy przetwarzaniem słów i pseudo-słów.

Przebieg eksperymentu

Uczestnicy zostali poinformowani o celu eksperymentu i charakterze pomiaru EEG. Zachęciliśmy ich, aby utrzymywali wygodną postawę i kontrolowali mruganie. Protokół zapewniał 3-s przerwy na normalne mruganie co 10 prób, jak również dwie dłuższe przerwy, których czas trwania jest kontrolowany przez uczestnika, dla odpoczynku i poprawy postawy. Długie przerwy występowały co 270 prób.

Zadaniem było odczytanie bodźców, które pojawiły się na środku ekranu, i zaklasyfikowanie ich jako słów lub pseudo-słów przez naciśnięcie oznaczonych klawiszy na klawiaturze. Rejestrowano typ i opóźnienie odpowiedzi. Pojedynczy blok eksperymentalny zawierał 135 słów i 135 pseudo-słów; ten blok powtórzono trzy razy. Słowa i pseudo słowa były wyświetlane w losowej kolejności we wszystkich blokach. Próby przebiegały w następujący sposób: (1) punkt fiksacji wyświetlany przez 500 ms; (2) bodziec wyświetlany, dopóki uczestnik nie odpowie; (3) pusty ekran wyświetlany przez losowo zmienny odstęp czasu między 1000 a 1100 ms.

Przygotujmy katalogi do analizy. Przydadzą się następujące:

• DANE
• SKRYPTY
• RYSUNKI

Dane dostępne są tutaj.

Proszę je ściągnąć i wypakować do katalogu DANE. Powinno pojawić się 66 plików, po 2 na każdą osobę. Pliki *.set to pliki matlabowe zawierające struktury EEG, pliki *fdt zawierają właściwe sygnały.

SKRYPTY

Etapy analizy, które musimy wykonać to:

• FILTROWANIE i CIECIE sygnałów
• DOBÓR PROGÓW DLA ARTEFAKTOW na podstawie oglądania sygnału i prób z różnymi progami do autoamtu (abnormal values / abnorlmal trends)
• Automatyczne usuwanie artefaktów
• Przygotowanie danych do analiz wg. przynależności do czynnika grupującego
• Wybór obszarów zainteresowania (grup elektrod do analizy)
• Analiza wariancji
• Zilustrowanie wyników analizy

Aby nasze analizy były łatwo powtarzalne przygotujemy główny skrypt np: run_me.m, do którego będziemy dodawać stopniowo skrypty i funkcje kolejnych etapów analizy. Dla czytelności kodu i możliwości uruchamiania osobno poszczególnych etapów posłużymy się technologią komórek w Matlabie (komórkę rozpoczynają dwa procenty '%%');

Zacznijmy zatem pisać run_me.m. Ustawiamy zmienne ze ścieżkami, wybieramy, które dane wchodzą do analizy (część została odrzucona ze względu na nadmiar artefaktów), częstości odcięcia filtrów dolno- i górno-przepustowego, czasy wokół momentu wyświetlenia bodźca, informacja, które kanały zawierają EEG, zakres czasu do korekty poziomu "0":

%% USTAWIENIA:
directoryIn  = '../DANE/';                       % lokalizacja folderu z surowymi danymi
directoryOut = '../DANE/';                       % lokalizacja do zapisu przygotowanych danych
addName      = '_epoch';                        % dodatek do nazwy pliku przy zapisywaniu przygotowanych danych
subjects     = [1:6,8:11 ,15:24,26:37];         % lista badanych do wczytania (z poprawnymi danymi)
LowFreqCut   = 0.1;                             % granica odciecia filtru gornoprzepustowego
HighFreqCut  = 30;                              % granica odciecia filtru dolnoprzepustowego
startT       = -0.3;                            % start do wyciecia epok [s]
stopT        = 1.0;                             % stop do wyciecia epok [s]
eegChan      = 1:19;                            % numery kanalow EEG
baseline     = [-200, 0];                       % baseline [od,do] [ms] do odjecia


%% FILTROWANIE i CIECIE


eeglab
filtruj_i_tnij(EEG , ALLEEG , CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)

filtruj_i_tnij.m

Musimy teraz napisać funkcję filtruj_i_tnij. Musi ona kolejno:

• wczytać nasze zbiory danych,
• stworzyć montaż kanałów do średniej ze wszystkich kanałów, tzw. common average reference (CAR),
• przefiltrować je,
• kanał po kanale,
• wyciąć interesujące nas epoki,
• wykonać dla nich korektę poziomu "0",
• zapisać tak przetworzone dane do nowych "*.set"-ów.

Szkielet tej funkcji wygląda następująco, proszę uzupełnić brakujące fragmenty:

function filtruj_i_tnij(EEG, ALLEEG, CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)
    
  % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
    EEG0 = EEG; 
    ALLEEG0 = ALLEEG; 
    CURRENTSET0 = CURRENTSET;
    events = {'Pseudo','ANeg','ANeu','APos','ONeg','ONeu','OPos','RNeg','RNeu','RPos'};   % typ eventow do wyciecia
    
    for sub = subjects
        
        % tworzymy nazwę setu dla osoby z numerem 'sub'
        if sub < 10
            nameOfSubj = ['B0' num2str(sub)];
        else
            nameOfSubj = ['B' num2str(sub)];
        end

        % wczytujemy set z danymi
        EEG = pop_loadset('filename',[nameOfSubj '.set'],'filepath',directoryIn);

        % uaktualniamy zestaw danych globalnych dla eeglaba:
       [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );
        EEG = eeg_checkset( EEG);
        
        % zostawiamy tylko potrzebne kanaly EEG , te które są wymienione w zmiennej eegChan:
        EEG.data = EEG.data(....); % wycianmy dane
        EEG.nbchan = .... ; % uaktualniamy liczbę kanałów
        EEG.chanlocs = EEG.chanlocs(1,...);  % obcianamy lstę położeń kanałów
        
        % referencja - odjecie sredniej ze wszystkich kanalow
        EEG = pop_reref(.... );
        EEG = eeg_checkset(EEG);

        % filtracja danych - zastosujemy filtr Butterwartha  2-go rzędu
        [b1,a1] = butter(2, .... , 'high');
        [b2,a2] = butter(2, ... , 'low');
        [b3,a3] = butter(2, 2*[49.5 50.5]/EEG.srate, 'stop'); % filtr pasmowozaporowy  "notch" na sieć
        for channel = 1:size(EEG.data,1)
            EEG.data(channel,:) = ... % górnoprzepustowy
            EEG.data(channel,:) = ... % dolno
            EEG.data(channel,:) = ... % notch
        end
        EEG = eeg_checkset( EEG );
        
        % ekstrakcja epok: 
        % zauważmy, że podajemy w tej funkcji listę eventów wokół których mają być robione wycinki i okresy wycinków
        % przy okazji nadajemy setowi nazwę z dodanym stringiem addName: 
        EEG = pop_epoch( EEG, events, [startT, stopT], 'newname', [nameOfSubj addName], 'epochinfo', 'yes');

        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'gui','off'); 
        EEG = eeg_checkset( EEG );
        
        % usuniecie baseline
        EEG = pop_rmbase(...);
        EEG = eeg_checkset( EEG );
        
        % zapisanie setu
        EEG = pop_saveset( EEG, 'filename',[nameOfSubj addName],'filepath',directoryOut);
        disp(['filtrowanie i cięcie ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end

    % przywracamy  stan wejściowy zmiennych globalnych eeglaba
    EEG0 = EEG; 
    EEG = EEG0; 
    ALLEEG = ALLEEG0; 
    CURRENTSET = CURRENTSET0;
end

DOBÓR PROGOW DLA ARTEFAKTÓW

Usuwanie artefaktów stanowi jeden z bardziej żmudnych ale niezbędnych etapów przygotowania danych do dalszej analizy. Jeśli mamy dostatecznie dużo danych to najbardziej bezpiecznym podejściem jest usunięcie tych epok (spośród wcześniej wyciętych), które mają artefakty. Do najczęstszych i najbardziej szkodliwych w analizie potencjałów wywołanych należą:

• mrugnięcia
• płynięcie sygnału
• mięśnie

Zwykle usuwanie artefaktów wykonujemy dwuetapowo:

• Dla ustalonych progów na wykrycie artefaktu oznaczamy automatycznie złe fragmenty
• Dokonujemy wizualnego przeglądu tak oznaczonych danych i wykonujemy ewentualne korekty, w miejscach nietypowych zaburzeń.

Tu dla oszczędzenia czasu zrobimy tylko procedurę uproszczoną:

• wczytamy jeden set, np.: B01_epoch.set
• popracujemy chwilę z narzędziem: Tools -> Reject data epochs -> Reject data (all methods)
• link do tutoriala EEGlab na ten temat: https://sccn.ucsd.edu/wiki/Chapter_01:_Rejecting_Artifacts#Rejecting_artifacts_in_epoched_data

Naszym celem jest takie dobranie progów

• abnormalValue: wartosc [uV] powyżej/poniżej(na minusie) której sygnał uznawany jest za artefakt
• abnormalTrend: maksymalne nachylenie trendu [uV/epokę], R squared limit (0-1)]

Jak już zdecydujemy się na konkretne wartości to ...

Wracamy na chwilę do run_me.m: Automatyczne usuwanie artefaktów

Zaaplikujmy ustalone przed chwilą parametry do procedury automatycznego ich usuwania. W pliku run_me.m dopisujemy następujące linijki:

%% Automatyczne usowanie artefaktow

abnormalValue = ...                            % wartosc [uV] powyzej/ponizej(na minusie) ktorej oznaczane sa artefakty
abnormalTrend = ...                     % [max slope [uV/epoch], R squared limit (0-1)]
oznaczArtefakty(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects, abnormalValue, abnormalTrend)

Funkcja: oznaczArtefakty.m

W pętli po zbiorach danych aplikujemy ustalone progi:

function [] = oznaczArtefakty(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects, abnormalValue, abnormalTrend)

      % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
       ... 
    for sub = ...
        % składamy nazwę pliku
          ....

        % wczytujemy plik:
        EEG = ...
        [ALLEEG, EEG, CURRENTSET] = ...
        EEG = eeg_checkset( EEG);
        
        % automatyczne oznaczanie artefaktow
               
        EEG = pop_eegthresh(...); % abnormal value
        EEG = eeg_checkset(EEG);
        EEG = pop_rejtrend(... ); % abnormalTrend
        EEG = eeg_checkset(EEG);
        close(gcf)
        
        % zapisanie setu
        EEG = pop_saveset(EEG, 'filename', [nameOfSubj addName], 'filepath', directoryOut);
        disp(['oznaczArtefakt ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end
    % przywracamy stan wejściowy zmiennych globalnych eeglaba
       ... 
end

Funkcja: preparujDane.m

%% przygotowanie danych

data = preparujDane(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects);     
% data - macierz z danymi o rozmiarze [2(1.Words|2.Pseudo) x subjects x kanaly x probki]

Funkcja: preparujDane.m

%% GFP

timeMarks = []; plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs) timeMarks = [70,120,225,290,375,670]; plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs)

%% ANOVA

ROI.channels = [[18,13];[15,10];[19,17];[9,4];[11,6]]; ROI.labels = {'LF','CF','RF','LP','RP'}; dataANOVA = prepareDataForANOVA(data,timeMarks,ROI,startT,stopT,EEG.srate); % dataANOVA - [5 okien czasowych x obserwacje x |1.srednia amplituda 2. subject 3.slowo(0)/pseudo(1) 4.ROI|]

for window = 1:(length(timeMarks)-1)

   dANOVA = squeeze(dataANOVA(window,:,:));
   stats  = anova(dANOVA(:,1),dANOVA(:,2),dANOVA(:,3),dANOVA(:,4),{'type','ROI'});
   stats2 = PostHocStats(stats,dANOVA,ROI);
   
   disp(['------------- Time window: ' num2str(timeMarks(window)) '-' num2str(timeMarks(window+1)) '-------------'])
   disp(stats)
   disp(stats2)

end

%% rysowanie ERP w ROI

plotERP_ROI(data,startT,stopT,EEG.srate,timeMarks,ROI)

</source>

Laboratorium_EEG/EEGLAB