-
Notifications
You must be signed in to change notification settings - Fork 0
/
ex04.R
56 lines (36 loc) · 1.12 KB
/
ex04.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
library("Biobase")
library("ALL")
library("genefilter")
data("ALL")
ALL$BT
bcell = grep("^B", as.character(ALL$BT))
types = c("NEG", "BCR/ABL")
moltyp = which(as.character(ALL$mol.biol) %in% types)
ALL_bcrneg = ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)
library("genefilter")
sds = rowSds(exprs(ALL_bcrneg))
sh = shorth(sds)
sh
hist(sds, breaks=50, col="mistyrose",
xlab="desviación estándar",main="")
abline(v=sh, col="blue", lwd=3, lty=2)
ALLsfilt = ALL_bcrneg[sds>=sh, ]
dim(exprs(ALLsfilt))
dim(exprs(ALL_bcrneg))
table(ALLsfilt$mol.biol)
tt = rowttests(ALLsfilt, "mol.biol")
names(tt)
hist(tt$p.value, breaks=50, col="mistyrose", xlab="valor p",main="")
library("multtest")
mt = mt.rawp2adjp(tt$p.value, proc="BH")
g = featureNames(ALLsfilt)[mt$index[1:10]]
library("hgu95av2.db")
links(hgu95av2SYMBOL[g])
mb = ALLsfilt$mol.biol
y = exprs(ALLsfilt)[g[1],]
ord = order(mb)
plot(y[ord], pch=c(1,16)[mb[ord]],
col=c("black", "red")[mb[ord]],
main=hgu95av2SYMBOL[g[1]], ylab="Expresión",
xlab="Muestras")