Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
Celem laboratorium jest zapoznanie Państwa z metodami analizy danych pochodzących z
eksperymentów mikromacierzowych. W trakcie laboratorium zostanie przeprowadzona
analiza danych opisanych przez Goluba – są to dane chorych na dwa róŜne typy białaczki
– dane z mikromacierzy oznaczonych symbolem ALL to dane chorych na ostrą białaczkę
limfoblastyczną, natomiast mikromacierzy oznaczone symbolem AML pochodzą z badań
chorych na ostrą białaczkę szpikową. Wynikiem analiz będzie opisana (zanotowana) lista
50 genów, które najbardziej róŜnicują obydwa zbiory danych.
Opis danych :
http://www.broad.mit.edu/
cgi-bin/cancer/publications/pub_paper.cgi?mode=view&paper_id=43
Odnośnik do artykułu Goluba, w którym przedstawiono wyniki obliczeń dla danych
wykorzystywanych w trakcie laboratorium.
http://www.sciencemag.org/cgi/reprint/286/5439/531.pdf
Utwórz katalog roboczy i skopuj do niego skrypty Normalize.R i getGOIDs.R
Uruchom program R.
Wykorzystując polecenie setwd lub z menu (File->Change Dir) zmień domyślny katalog
roboczy na katalog utworzony przez Ciebie
Pobierz oraz skopiuj do przestrzeni danych zbior golubEsets.
library(golubEsets)
data(Golub_Merge)
Jakiego typu jest to obiekt? Sprawdź jak wygląda struktura wczytanych danych –
polecenie str(Golub_Merge). Jakie dane moŜna przechowywać w obiektach klasy
ExpressionSet?
Wyświetl informacje na temat wczytanej zmiennej poprzez podanie jej nazwy:
Golub_Merge
W tym wypadku zostanie wyświetlona tylko ogólna informacja o zmiennej
zadeklarowanej metodą show.
Wykonaj na danych funckcję featureNames: featureNames(Golub_Merge)
Co otrzymałeś wyniku?
Preprocesing danych.
Aby znormalizować dane wykorzystaj funkcję normalize.R dostarczoną w ramach
laboratorium. Postaraj się zrozumieć działanie tej funkcji. Skomentuj w raporcie kaŜdą
jej linijkę.
Za pomocą polecenia source wczytaj skrypt do przestrzeni roboczej R:
source("normalize.R")
Wykorzystaj skrypt do znormalizowania danych
Strona 1 z 6
Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
expressedData =normalize(Golub_Merge)
expr.Golub=expressedData$expressedData
expr.ProbeSetsNames=expressedData$probeSetsNames
Aby w prosty sposób moc odwoływać się do grup danych dla typu białaczki AML oraz
ALL, wygodnie jest utworzyć zmienne zawierające numery indeksów dla obydwu grup
danych.
ALL_indexes =which(Golub_Merge$ALL.AML=="ALL")
AML_indexes =which(Golub_Merge$ALL.AML=="AML")
Ile było mikromacierzy dla białaczki typu AML w analizowanych danych, a ile
mikromacierzy dla białaczki typu ALL? Wykorzystaj funkcję length.
Wyznaczanie zbioru genów róŜnicujących pomiędzy dwoma typami białaczki:
W celu wyznaczenia genów róŜnicujących grupy białaczki zostanie wykorzystany prosty
test t. Rezultaty tego testu mogą zostać wykorzystane do wyznaczenia, które geny uległy
ekspresji w sposób istotny odbiegający od średniej. Takie geny mogą później zostać
wykorzystane jako geny róŜnicujące pomiędzy grupą AML oraz ALL
Przykładowy test t dla pierwszej sondy:
t.test(expr.Golub[1,ALL_indexes], expr.Golub[1,AML_indexes])
Sprawdź jak wygląda struktura zwracanej zmiennej - poleceniem str
Dla wykonywanych przez nas analiz interesująca jest wartość poziomu istotności testu
(p-value). W jaki sposób moŜemy otrzymać tą wartość?
Test dla kaŜdej sondy macierzy moŜe zostać przeprowadzony na dwa sposoby: z
wykorzystaniem pętli for lub teŜ poprzez uŜycie funkcji sapply, która w poniŜszym zapisie jest równoznaczna z pętlą for, lecz jest prostsza w zapisie:
expr.Golub.pVal<-sapply(1:nrow(expr.Golub), function(i)
t.test(expr.Golub[i,AML_indexes],expr.Golub[i,ALL_indexes])$p.value)
Mając tak wyliczony wektor p-value moŜna sprawdzić jaka liczba sond rozróŜnia
wskazane grupy na poziome istotności mniejszym od 5%. Podaj tą liczbę.
table(expr.Golub.pVal<0.05)
PoniewaŜ wykonana powyŜej analiza jest wielokrotnym testowaniem konieczne jest
wprowadzenie odpowiedniej korekcji p-value uzyskanych w teście t. Wykorzystana
zostanie tu metoda oszacowania odsetka wyników fałszywie dodatnich (FDR)
Benjaminiego i Hochberga. Do tego celu zostanie wykorzystana funkcja p.adjust:
expr.Golub.pVal.bh=p.adjust(expr.Golub.pVal,method="BH")
Strona 2 z 6
Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
Jaka liczba sond ma wartość p-value poniŜej 5% dla wyników skorygowanych metodą
BH? Podaj otrzymaną wartość.
Przygotuj zmienną zawierającą raport z otrzymanych wyników:
W raporcie (zmiennej0 powinny znaleźć się następujące informacje:
- nazwy sond Affymetrixa - dla genów które uległy ekspresji (results$AffyID)
- wartości poziomu istotności p (p-value) (results$pVal)
- skorygowane wartości p-value (results$pVal.bh)
- średnie wartości w całym przedziale danych (results$A)
- krotność ekspresji (fold change) (results$M)
- średnie wartości ekspresji dla danych ALL (results$M.ALL)
- średnie wartości ekspresji dla danych AML (results$M.AML)
Wyniki wpisz do zmiennej results – listy :
results=list()
results$AffyID<
results$pVal<-
results$pVal.bh<-
results$A<-
results$M<-
results$M.ALL<-
results$M.AML<-
Do obliczenia średnich wartości wykorzystaj funkcję apply. Funkcja apply wykonuje
na zadanej macierzy funkcję względem zadanego wymiaru. Obliczenie średniej dla
kaŜdego z wierszy tablicy moŜna wykonać następująco:
tmp.A<-apply(expr.Golub,1,mean)
Sprawdź strukturę otrzymanej zmiennej str(tmp.A)
PoniewaŜ otrzymana zmienna posiada atrybuty nazwy, który w naszym przypadku jest
zbędny, naleŜy go usunąć wykorzystując polecenie unname:
results$A<-unname(tmp.A)
Aby obliczyć wartość krotności M ekspresji (dla danych zlogarytmowanych) oblicz
róŜnicę średnich ekspresji w grupie ALL i AML. Usuń zbędne nazwy.
apply(expr.Golub[,AML_indexes],1,mean)-
apply(expr.Golub[,ALL_indexes],1,mean)
Na podstawie powyŜszych przykładów oblicz średnie wartości ekspresji w grupach AML
i ALL.
Ogranicz raport tylko do genów (sond) mających znamienną wartość p-value
skorygowaną (przy załoŜeniu Ŝe za istotne uznamy geny dla których FDR<0.01,
uzyskana lista zawierać będzie wtedy 1% genów fałszywie dodatnich).
Aby osiągnąć ten cel posortuj wyniki.
Strona 3 z 6
Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
PoniewaŜ najprościej wykonać sortowanie dla obiektu typu data.frame, przekonwertuj
zmienną results typiu lista na zmienną typu data frame (ramka danych):
results=as.data.frame(results)
Posortuj wszystkie kolumny ramki danych według kolumny results$pVal.bh
res.sorted=results[order(results$pVal.bh),]
Zostaw tylko te wartości, dla których p-value po skorygowaniu jest mniejsza od 0.01
res.sorted<-res.sorted[res.sorted$pVal.bh<0.01,]
Adnotacja otrzymanych wyników
Dla tak otrzymanego zbioru danych dokonaj adnotacji.
Adnotacja będzie składać się z dwóch kroków. W pierwszym kroku przetłumacz nazwy
sond Affymetryxa na zrozumiałe dla biologów symbole genów. Najprościej zrobić to
wykorzystując gotowe biblioteki adnotacji przygotowanych dla kaŜdego typu macierzy.
W jaki sposób moŜna sprawdzić jakiego typu macierz została wykorzystana w
eksperymencie przeprowadzonym przez Goluba?
Biblioteki niezbędne do przeprowadzenia adnotacji sond:
library(annotate)
library(hu6800.db)
- biblioteka zawierająca adnotacje dla macierzy HG8600
Dopisanie symbolu genu do listy wynikowej obciętej do p.value<0.01 moŜna uzyskać
stosując funkcję getSYMBOL. Funkcja ta wymaga podania parametrów: numeru sondy
(w tym przypadku jest to wektor numerów sond) oraz typu macierzy.
res.sorted $AffyID =as.character(res.sorted $AffyID)
res.sorted $Symbol=getSYMBOL(res.sorted $AffyID,”hu6800”)
W bazie danych Gene Ontology znajdują się terminy (GO Terms), które mogą być
wykorzystywane do adnotacji genów – czyli opisywania tych genów pod względem ich
funkcji iologicznej. Terminy te dzielą się na trzy podstawowe ontologie: Proces
Biologiczny, Komponent Komórkowy oraz Funkcja Molekularna.
Drugi krok adnotacji będzie polegał na pobraniu dla kaŜdego genu jego zbioru
identyfikatorów GO dla Procesu Biologicznego, przetłumaczenie ich na terminy GO i opisanie tymi terminami kaŜdego genu. Do tego celu konieczne będzie wykorzystanie biblioteki „GO”
Pobranie listy terminów GO dla wszystkich genów w zbiorze:
gg=getGO(res.sorted$AffyID,"hu6800")
Aby pobrać listy GO dla wszystkich sond znajdujących się w wekorze res.sorted$AffyID
wykorzystaj funkcję (getGOIDs):
source("getGOIDs.R",echo=TRUE)
GOIDsList= getGOIDs(res.sorted$AffyID)
Strona 4 z 6
Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
Napisz skrypt który dla podanego wektora sond zwróci listy terminów GO
MoŜe się przydać:
Lista wszystkich identyfikatorów GO dla pierwszej sondy:
oneGeneGOIDs= GOIDsList[[1]]
Pobranie opisu GO dla pierwszego terminu GO z listy oneGeneGOIDs -wszystkie
ontologie. Funkcja getGOdesc przyjmuje dwa argumenty: identyfikator GO oraz symbol
typu ontolgi który chcemy uzyskać (ANY – wszystkie ontologie, BP – proces
biologiczny)
GOdesc=getGOdesc(oneGeneGOIDs[1],"ANY")
Pobranie terminu GO dla pierwszego terminu GO z listy oneGeneGOIDs
goTerm=GOdesc[[1]]@Term
Wynikiem końcowym laboratorium (poza skryptem) powinien być plik tekstowy
zawierający listę dla 50 genów – kryteria doboru genów zostaną podane na zajęciach.
KaŜdy gen powinien znajdować się w osobnym wierszu, a wartości dla kaŜdego genu powinny być oddzielone znakiem tabulacji.
Dla kaŜdego genu w pliku powinny znajdować się następujące informacje:
-nazwa sondy Affymetrixa
-symbol genu
-wartość p-value
-skorygowana wartość p-value
-krotność ekspresji
-lista terminów GO dla wszystkich genów – typ ontologii zostanie podany na zajęciach Materiały źródłowe:
1.Gala G.: „Podstawy analizy danych uzyskanych techniką mikromacierzy
oligonukleotydowych
wysokiej
gęstości
Affymetrix”
w
środowisku
R/Bioconductor. Warsztaty Bioinformatyczne Projeku PBZ-MNiI-2/1/2005
Analiza Danych w Genomie Funkcjonalnej
2. http://www.r-project.org/
3. http://www.bioconductor.org/
Strona 5 z 6
Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia
Zadanie dodatkowe:
Dla otrzymanych wyników stworzyć wykres (volcano plot) rozrzutu pomiędzy wynikami
krotności ekspresji, a wartością poziomu istotności.
Geny dla których –log10 z poziomu istotności jest większy od 10 i równocześnie
absolutna wartość krotności ekspresji jest większa od 2 zaznaczyć oraz podpisać tak jak
na poniŜszym wykresie.
Informacje pomocnicze:
plot(M,-log10(expr.Golub.pVal)) – funkcja która narysuje wykres rozrzutu
wartości M oraz –log10 z poziomu istotności dla danych z laboratorium.
Do opisania wykresu wykorzystać funkcje: lines, points, text.
Aby znaleźć indeksy genów spełniające obydwa warunki (-log10 z poziomu istotności >
10 oraz wartość absolutna krotności ekspresji >2) wykorzystaj funkcję which oraz
funkcję intersect.
Strona 6 z 6