Zadania

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 1 z 6

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

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

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 2 z 6

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

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 3 z 6

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.

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 4 z 6

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)

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 5 z 6

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

sto

Affymetrix”

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/

Laboratorium Bioinformatyka autor: Aleksandra Gruca
Zadnia

Strona 6 z 6

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.