UNIVERZA V LJUBLJANI ˇ BIOTEHNISKA FAKULTETA ODDELEK ZA ZOOTEHNIKO
Gregor GORJANC1
Mikromreˇ ze in hierarhiˇ cno razvrˇ sˇ canje
Domaˇca naloga pri predmetu “genomika” Podpiplomski ˇstudij genetike
Ljubljana, 2006 1
[email protected]
Gorjanc G. Mikromreˇ ze in hierarhiˇ cno razvrˇsˇ canje. Domaˇ ca. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
1
2
Priprava podatkov
Za delo sem izbral odprtokodni programski jezik in okolje za analizo podatkov R (Ihaka and Gentleman, 1996; R Development Core Team, 2005). Prikazujem tako potek (programsko kodo), kakor tudi rezultate. Najprej sem uvozil podatke in jih ustrezno uredil za nadaljno analizo. R> podatki <- read.csv(file = "podatki.csv") R> rownames(podatki) <- podatki$gene R> podatki <- t(as.matrix(podatki[, -1])) R> podatki C D E
F G
H
I J K L
M
N
X0 X2
1 1 1 1.00 1 1.00 1.0 1 1 1 1.00 1.00 8 3 4 1.00 2 0.50 4.0 2 1 2 0.33 0.13
X4 X6
12 4 8 1.00 3 0.33 8.0 1 1 3 0.25 0.08 16 4 8 0.25 4 0.25 4.0 2 1 4 0.25 0.06
X8 12 3 8 0.25 3 0.33 1.0 1 3 3 0.33 0.08 X10 8 2 8 0.10 2 0.50 0.5 2 3 2 0.50 0.13 Pripravil sem tudi set podatkov, kjer sem vrednosti transformiral z funkcijo log2 . R> podatkiLog <- podatki R> podatkiLog <- log(podatki, base = 2) R> round(podatkiLog, digits = 2)
X0
C D E 0.00 0.00 0
F G 0.00 0.00
H 0.0
I J K L 0 0 0.00 0.00
M 0.0
N 0.00
X2
3.00 1.58 2
0.00 1.00 -1.0
2 1 0.00 1.00 -1.6 -2.94
X4 X6
3.58 2.00 3 0.00 1.58 -1.6 4.00 2.00 3 -2.00 2.00 -2.0
3 0 0.00 1.58 -2.0 -3.64 2 1 0.00 2.00 -2.0 -4.06
X8 3.58 1.58 3 -2.00 1.58 -1.6 0 0 1.58 1.58 -1.6 -3.64 X10 3.00 1.00 3 -3.32 1.00 -1.0 -1 1 1.58 1.00 -1.0 -2.94 Sledil je izraˇcun korelacijskih koeficientov po Pearsonu za oba seta podatkov: R> korelacije <- cor(podatki, method = "pearson") R> korelacijeLog <- cor(podatkiLog, method = "pearson")
2
Hierarhiˇ cno razvrˇ sˇ canje/zdruˇ zevanje
2.1 Transformirani podatki Izraˇcunane korelacijske koeficiente za transformirane podatke sem uporabil za hierarhiˇcno razvrˇsˇcanje/zdruˇzevanje. Ker je matrika korelacijskih koef. simetriˇcna, sem za laˇzje delo izbral le vrednosti nad diagonalo. Koeficienti so na zaˇcetku imeli tako sledeˇce vrednosti.
Gorjanc G. Mikromreˇ ze in hierarhiˇ cno razvrˇsˇ canje. Domaˇ ca. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
3
C D E F C 0 0.94 0.96 -0.40
G H 0.95 -0.95
I 0.41
J 0.36
K 0.23
L M N 0.95 -0.94 -1.00
D 0 0.00 0.84 -0.10 E 0 0.00 0.00 -0.57
0.94 -0.94 0.89 -0.89
0.68 0.21
0.24 -0.07 0.30 0.43
0.94 -1.00 -0.95 0.89 -0.84 -0.96
F 0 0.00 0.00 G 0 0.00 0.00
0.00 -0.35 0.35 0.00 0.00 -1.00
0.60 -0.43 -0.79 -0.35 0.10 0.40 0.48 0.22 0.11 1.00 -0.94 -0.96
H 0 0.00 0.00 I 0 0.00 0.00
0.00 0.00
0.00 0.00
0.00 -0.48 -0.21 -0.11 -1.00 0.94 0.96 0.00 0.00 0.00 -0.75 0.48 -0.68 -0.42
J 0 0.00 0.00 K 0 0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.22 -0.24 -0.33 0.11 0.07 -0.22
L 0 0.00 0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00 -0.94 -0.96
M 0 0.00 0.00 N 0 0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.00 0.00
0.95 0.00
Najviˇsji korelacijski koeficient na zaˇcetku zdruˇzevanja je bil med genoma G in L, kjer je znaˇsal koeficient kar 1. Za laˇzje delo sem napisal funkcijo, ki poiˇsˇce lokacijo najveˇcjega korelacijski koeficient v matriki, s tem pa tudi katera dva gena/skupine moramo zdruˇziti v posameznem koraku. R> poisciMaxKoef <- function(x) { + k <- nrow(x) + +
tmp <- which.max(x) tmp1 <- floor(tmp/k)
+ +
i <- tmp - (tmp1 * k) j <- tmp1 + 1
+ +
list(coef = x[i, j], dim = c(i, j), group = c(rownames(x)[i], colnames(x)[j]))
+
}
Za prvi korak je kot smo ˇze povedali najviˇsja korelacija med genoma G in L. R> poisciMaxKoef(korelacijeLog) $coef [1] 1 $dim [1]
5 10
$group [1] "G" "L" Nadalje sem napisal funkcijo, ki sprejme osnovne podatke in rezultat funkcije poisciMaxKoef ter izraˇcuna povpreˇcje novo zdruˇzene skupine. R> zdruzi <- function(x, y) { + if (!ncol(x) > 3)
Gorjanc G. Mikromreˇ ze in hierarhiˇ cno razvrˇsˇ canje. Domaˇ ca. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
+ +
cols <- colnames(x) x <- cbind(rowMeans(x[, y$dim, drop = TRUE]), x[, -y$dim, drop = TRUE])
+ +
if (ncol(x) > 2) { colnames(x)[1] <- paste("[", paste(y$group, collapse = ""),
+ +
}
+ +
else { colnames(x) <- c(paste("[", paste(y$group, collapse = ""),
4
"]", collapse = "", sep = "")
+ +
"]", collapse = "", sep = ""), cols[-y$dim]) }
+
x
+ } R> round(zdruzi(x = podatkiLog, y = poisciMaxKoef(korelacijeLog)), +
digits = 2) [GL]
C
D E
F
H
I J
K
M
N
X0 X2
0.00 0.00 0.00 0 1.00 3.00 1.58 2
0.00 0.0 0.00 -1.0
0 0 0.00 0.0 0.00 2 1 0.00 -1.6 -2.94
X4 X6
1.58 3.58 2.00 3 0.00 -1.6 2.00 4.00 2.00 3 -2.00 -2.0
3 0 0.00 -2.0 -3.64 2 1 0.00 -2.0 -4.06
X8 1.58 3.58 1.58 3 -2.00 -1.6 0 0 1.58 -1.6 -3.64 X10 1.00 3.00 1.00 3 -3.32 -1.0 -1 1 1.58 -1.0 -2.94 Sedaj sem samo ˇse napisal funkcijo, ki izvede zdruˇzevanje in preraˇcunavanje za vse gene in sproti izpisuje rezultate vse dokler ni v podatkih samo ˇse en stolpec - kar nakazuje, da smo konˇcali z zdruˇzevanjem. R> hierarhicnoRazvrscanje <- function(x) { +
while (ncol(x) > 1) {
+ +
kor <- cor(x, method = "pearson") kor <- urediMatrikoKorelacij(kor)
+ +
maxKoef <- poisciMaxKoef(kor) cat(maxKoef$group, "\n")
+ +
x <- zdruzi(x = x, y = maxKoef) }
+ } R> hierarhicnoRazvrscanje(podatkiLog) G L H N C E [HN] M [GL] D [[GL]D] [CE] F I [[[GL]D][CE]] J
Gorjanc G. Mikromreˇ ze in hierarhiˇ cno razvrˇsˇ canje. Domaˇ ca. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
5
[[[[GL]D][CE]]J] K [FI] [[HN]M] [[FI][[HN]M]] [[[[[GL]D][CE]]J]K]
2.2 Ne-transformirani podatki Enako lahko storimo tudi za ne-transformirane podatke: R> hierarhicnoRazvrscanje(podatki) G L [GL] C H M [HM] N [[GL]C] D [[[GL]C]D] E F I [[[[GL]C]D]E] K [[[[[GL]C]D]E]K] J [[[[[[GL]C]D]E]K]J] [FI] [[[[[[[GL]C]D]E]K]J][FI]] [[HM]N] Opazimo, da se rezultati hierarhiˇcnega razvrˇsˇcanja na ne-transformiranih podatkih ne razlikujejo bistveno od rezultatov pri transformiranih podatkih. Omembe vredna razlika je le v zdruˇzevanju skupine [FL], ki se pri transformiranih podatkih zdruˇzuje s skupino [[HN]M], pri ne-transformiranih podatkih pa se le ta skupina zdruˇzuje z [[[[[[GL]C]D]E]K]J] skupino.
3
Dendrogrami
Program R ima tudi funkcije za hierarhiˇcno razvrˇsˇcanje in risanje dendrogramov. R> tmpLog <- hclust(d = as.dist((1 - cor(podatkiLog))/2), method = "average") R> tmp <- hclust(d = as.dist((1 - cor(podatki))/2), method = "average") Iz neznanih razlogov funkcija hclust vrne ravno obratno zdruˇzevanje skupine [FI] glede na transformacijo podatkov kot sem izraˇcunal sam.
Gorjanc G. Mikromreˇ ze in hierarhiˇ cno razvrˇsˇ canje. Domaˇ ca. naloga. Ljubljana, Univ. v Ljubljani, Biotehniˇska fakulteta, 2006
6
R> par(mfrow = c(1, 2)) R> plot(tmpLog, axes = FALSE, ylab = "", xlab = "", main = "transformacija") R> plot(tmp, axes = FALSE, ylab = "", xlab = "", main = "brez transformacije")
brez transformacije
C G L N H M
D
C E D G L
M H N
E
F I
F I
J
K
J
K
transformacija
ˇ Slika 1: Stevilo zadetkov v bazi PubMed glede na statistiˇcni pristop in leto
4
VIRI
Ihaka, R. and Gentleman, R. (1996). R: A language for data analysis and graphics. J. Comput. Graph. Stat., 5(3):299–314. R Development Core Team (2005). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0.