$R
$ ls
GSM349024.CEL.gz GSM349027.CEL.gz GSM349030.CEL.gz
GSM349025.CEL.gz GSM349028.CEL.gz GSM349023.CEL.gz
GSM349026.CEL.gz GSM349029.CEL.gz
$R
library(affy)
GSE13869 <- ReadAffy() #カレントディレクトリにあるすべてのセルファイルを読み込んでAffyBatchオブジェクトを生成する
eset <- rma(GSE13869)
data <- eset[,c(1,3,4,6,7,8)]#GEOの情報によると、wt/wt, ko/wt, ko/ko, wt/wt, ko/wt, ko/ko, wt/wt, ko/koなので２群比較を視野にいれてwt/wtとko/koを抽出して新にオブジェクトを作った。
exprs <- exprs(data) #ExprsSetオブジェクトから発現データとアレイ名とスポット名だけを抜き取ってマトリックスを生成
mean_wt <- rowMeans(exprs[ ,c(1,3,5)])
mean_ko <- rowMeans(exprs[ ,c(2,4,6)])
difference <- mean_ko - mean_wt #RMAで正規化したときは、発現値がlog2値で変えされるので、Fold Changeは引き算で、log2(mean_ko/mean_wt)となる。
liobrary(genefilter)
<- mean_ko - mean_wt
a <- rowMeans(exprs)
# d <- rowMeans(exprs[ ,c(2,4,6)]) - rowMeans(exprs[ ,c(21,3,5)])ってこと。
library(genefilter)
label <- c(0,1,0,1,0,1)
#最後に、Valcanoプロットという手法で全遺伝子について、縦軸をp値のLOG10、横軸を平均との差を同時に描写することにした。
#俯瞰するために、
tt2 <- rowttests(exprs, factor(label)) #FC法を用いる前の段階のデータにt検定を行い、すべての遺伝子にp値を与えた
lod <- -log10(tt2$p.value)
png("110203_volcano.plot.1.png")
plot(d, lod,pch=".")
abline(v=1, col = "blue") #FCが2倍以上（log2で1以上）でカットオフラインを引く。vはvertical
abline(h=1, col = "red") #p値が0.1以下(-log10で1以上）でカットオフラインを引く。hはhorizontal
dev.off()
#カットオフを緩くして考えてみる。
png("110203_volcano.plot.2.png")
par(mfrow = c(2,2))
o1 <- order(abs(d), decreasing = TRUE)[1:10] #前回の解析で、FCが２倍以上のものが上位10個をマークしてみる。
o2 <- order(abs(tt2$statistic), decreasing = TRUE)[1:10] #p値についても上位10個をプロットしてみる。
o <- union(o1,o2)
plot(d[-o], lod[-o], xlim=c(-2,2), ylim = c(0, 4))
points(d[o1], lod[o1], pch = 18 , col = "blue")
points(d[o2], lod[o2], pch = 1 , col = "red")
o1 <- order(abs(d), decreasing = TRUE)[1:100] #前回の解析で、FCが２倍以上のものが上位100個をマークしてみる。
o2 <- order(abs(tt2$statistic), decreasing = TRUE)[1:100] #p値についても上位100個をプロットしてみる。
o <- union(o1,o2)
plot(d[-o], lod[-o], xlim=c(-2,2), ylim = c(0, 4))
points(d[o1], lod[o1], pch = 18 , col = "blue")
points(d[o2], lod[o2], pch = 1 , col = "red")
o1 <- order(abs(d), decreasing = TRUE)[1:1000] #前回の解析で、FCが２倍以上のものが上位1000個をマークしてみる。
o2 <- order(abs(tt2$statistic), decreasing = TRUE)[1:1000] #p値についても上位1000個をプロットしてみる。
o <- union(o1,o2)
plot(d[-o], lod[-o], xlim=c(-2,2), ylim = c(0, 4))
points(d[o1], lod[o1], pch = 18 , col = "blue")
points(d[o2], lod[o2], pch = 1 , col = "red")
o1 <- order(abs(d), decreasing = TRUE)[1:10000] #前回の解析で、FCが２倍以上のものが上位10000個をマークしてみる。
o2 <- order(abs(tt2$statistic), decreasing = TRUE)[1:10000] #p値についても上位10000個をプロットしてみる。
o <- union(o1,o2)
plot(d[-o], lod[-o], xlim=c(-2,2), ylim = c(0, 4))
points(d[o1], lod[o1], pch = 18 , col = "blue")
points(d[o2], lod[o2], pch = 1 , col = "red")
dev.off()
