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

Wracamy na chwilę do run_me.m: zebranie danych>

Dalsze analizy będzie wygodniej prowadzić gdy dane od wszystkich osób zbierzemy w jednej macierzy o rozmiarze [kategoria x osoby x kanaly x probki]. Będziemy w niej przechowywać dane dla każdej z osób uśrednione po realizacjach należących do tej samej kategorii, tzn. 2 kategorie: (1.Words|2.Pseudo).

W run_me.m dopisujemy:

%% przygotowanie danych

data = preparujDane(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects);

Funkcja: preparujDane.m

Tym razem wczytując dane musimy iterować się po kolejnych zbiorach danych 'subIdx' abstrachując od numeru osoby w nazwie pliku. Chcemy aby w powstającej macierzy 'data' nie było zerowych wpisów z powodu odrzucenia (braku pliku) jakieś osoby. Proponuję po uzupełnieniu kodu zatrzymać się debugerm wewnątrz pętli sumującej triale i upewnić, że dokładnie rozumiemy jak są agregowane dane.

function [data] = preparujDane(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects)
 % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
       ... 
    data    = zeros(...);       % macierz z danymi o rozmiarze [2(1.Words|2.Pseudo) x subjects x kanaly x probki]
    trialNb = zeros(...);                                         % macierz z informacja o liczbie triali na warunek [2(1.Words|2.Pseudo) x subjects]
                    

    for subIdx = 1:length(subjects)
        
        %    składamy nazwę pliku
        if subjects(subIdx) < 10
            nameOfSubj = ['B0' num2str(subjects(subIdx))];
        else
            nameOfSubj =...
        end

        % wczytujemy plik:
        EEG = ...
        [ALLEEG, EEG, CURRENTSET] = ...
        EEG = eeg_checkset( EEG);
        
        % dla danej osoby iterujemy się po trialach
    
        for trial = 1:EEG.trials
            if ~(EEG.reject.rejconst(1,trial) || EEG.reject.rejthresh(1,trial))   % wybieramy tylko triale bez artefaktow 
                type = strcmp(EEG.epoch(1,trial).eventtype,'Pseudo')+1; % określamy kategorię trialu: 
                                                         % dla triali  Pseudo, strcmp da wartość True, czyli 1 więc w wyniku dostaniemy 2,
                                                         % cała reszta to słowa i dla nich wynik porównania strcmp da fałsz , czyli 0 
                data(type,subIdx,:,:) = squeeze(data(type,subIdx,:,:)) + squeeze(EEG.data(:,:,trial)); % i sumujemy ich 
                                                         % przebiegi dla każdej z dwoch kategorii: Words/Pseudo, żeby wymiary
                                                         % dodawanych macierzy się zgadzały musimy wpierw usunąć wymiary singletowe
                trialNb(type,subIdx) = trialNb(type,subIdx) + 1;
            end
        end
        
        % dzielimy przebiegi przez zliczenia triali w kategorii aby uzyskac sredni przebieg
        data(1,subIdx,:,:) = data(1,subIdx,:,:)/squeeze(trialNb(1,subIdx));
        data(2,subIdx,:,:) = data(2,subIdx,:,:)/squeeze(trialNb(2,subIdx));
        disp(['preparujDane ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end
    % przywracamy stan wejściowy zmiennych globalnych eeglaba
       ... 
    
end

Wybór okien czasowych dla komponentów ERP: metoda GFP

Wybór okienek do analizy ERP można oprzeć o krzywą globalnej mocy pola (ang. global field power GFP). GFP jest obliczany jako przestrzenne odchylenie standardowe i określa sumę aktywności elektrycznej wszystkich elektrod w danym punkcie czasowym. Opóźnienia maksymalnych wartości GFP wskazują na opóźnienia potencjalnych składników potencjału wywołanego (Lehmann i Skrandies, 1980; Skrandies, 1990). Intuicja, która za tym stoi mówi, że dla pewnego stanu aktywności mamy zróżnicowaną topografię amplitudy (rozkład dipolowy , kwadrupolowy itp. ) i przejście do innego stanu powinno wieść przez brak takiej różnorodności.

Do pliku

%% GFP1

timeMarks = [...]; %wpiszmy tu czasy z artykułu
plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs) 

%% GFP2 - widać że powyższe czasy można poprawić:
timeMarks = [...];
plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs)

Funkcja: plotGFP.m

function plotGFP(data,startT,stopT,Fs,timeMarks,chanlocs)
    
    fontSize  = 14;
    colors    = [[1 0 0];[0 0 1]];
    colorMark = [0.5 0.5 0.5];
    marginesX = 0.1;
    marginesY = 0.1;
    labels    = {'Words','Pseudo'};
    
    figure()
    time = (startT:1/Fs:stopT);
    scaleTopo = [min(min(min(mean(data,2)))) max(max(max(mean(data,2))))];
    positionPlot = [marginesX, 1-marginesY-(1-2*marginesY)*1/2, 1-2*marginesX, (1-2*marginesY)*1/2];
    positionTopo = zeros(length(timeMarks)-1,4); 
    positionTopo(:,2:4) = repmat([1-marginesY-(1-2*marginesY)*3/4, (1-2*marginesX)/length(timeMarks), (1-2*marginesY)*1/6],(length(timeMarks)-1),1);
    for tim = 1:(length(timeMarks)-1)
        positionTopo(tim,1) = marginesX+(tim-1)*(1-2*marginesX)/length(timeMarks)+(tim-1)*(1-2*marginesX)/length(timeMarks)/(length(timeMarks)-2);
    end
    positionTopo2 = positionTopo;
    positionTopo2(:,2) = 1-marginesY-(1-2*marginesY);
    positionCbar = [1-1.1*marginesX, 1-marginesY-(1-2*marginesY)*7/8, marginesX/4, (1-2*marginesY)*1/6];
    
    % rysowanie krzywej GFP i pionowych linii przedzialow czasowych
    axes('position',positionPlot)
    GFP = std(squeeze(mean(data(:,:,:,:),2)),[],2);
    YLIM = [floor(min(min(GFP))), ceil(max(max(GFP)))];
    for type=1:size(data,1)
        hold on
        plot(time,GFP(type,:)','Color',colors(type,:))
    end 
    if ~isempty(timeMarks)
        for l=1:length(timeMarks)
            line([timeMarks(l)/1000,timeMarks(l)/1000], [YLIM(1) YLIM(2)*0.9], 'Color', colorMark)
            legend off
            hold on
        end
    end
    xlim([startT,stopT])
    ylim(YLIM)
    set(gca,'box','off');
    set(gca, 'FontSize', fontSize)
    drawaxis(gca,'x',0,'y',0);
    axis off
    text(stopT*0.9,-0.1*YLIM(1)-0.3,'[ms]','FontSize',fontSize,'Interpreter','Latex')
    text(startT*0.5,0.8*YLIM(2),'[\textit{$\mu$}V]','FontSize',fontSize,'Interpreter','Latex')
    legend(labels,'Box','off')
    
    if ~isempty(timeMarks)
        % rysowanie glowek topo dla przedzialow czasowych
        for l=1:(length(timeMarks)-1)
            axes('position',positionTopo(l,:))
            idx = (time>timeMarks(l)/1000)&(time<timeMarks(l+1)/1000);
            topoplot(mean(data(1,:,:,idx),2), chanlocs, 'plotrad' ,0.6 ,'maplimits',scaleTopo);
            if l==1
                text(-1.5,0,labels{1},'FontSize',fontSize);
            end
        end
        
        for l=1:(length(timeMarks)-1)
            axes('position',positionTopo2(l,:))
            idx = (time>timeMarks(l)/1000)&(time<timeMarks(l+1)/1000);
            topoplot(mean(data(2,:,:,idx),2), chanlocs, 'plotrad' ,0.6 ,'maplimits',scaleTopo);
            if l==1
                text(-1.5,0,labels{2},'FontSize',fontSize);
            end
        end
        
        % rysowanie skali kolorow
        axes('position',positionCbar)
        c = colorbar;
        set(gca,'box','off');
        axis off
        caxis(scaleTopo);
        set(c,'FontSize',fontSize)
        text(9,0.5,'[\textit{$\mu$}V]','FontSize',fontSize,'Interpreter','Latex')
    end
    
    orient landscape
    set(findall(gcf,'-property','FontName'),'FontName','Liberation Serif')
    set(gcf,'DefaultFigurePaperType','A4');
end

ANOVA

%% 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)

Laboratorium_EEG/EEGLAB