$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)])
#mean_wtのデータ構造を覗いてみる
class(mean_wt)
[1] "numeric"
summary(mean_wt)
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.653 4.394 6.177 6.401 8.124 14.530
length(mean_wt)
[1] 45101
difference <- mean_ko - mean_wt #RMAで正規化したときは、発現値がlog2値で変えされるので、Fold Changeは引き算で、log2(mean_ko/mean_wt)となる。
sum(difference > 1) #２倍以上高発現しているプローブ数をしらべる。ちなみに、sum(differnce)とすると、differenceの成分の総和が計算される。
liobrary(genefilter)
difference_list <- difference[difference>1]
class(difference_list)
[1] "matrix"
difference_list
1417065_at 1417210_at 1419011_at 1420453_at 1420536_at 1424903_at
1.368660 2.924009 1.734660 1.066935 1.147841 2.096864
1426438_at 1426439_at 1426598_at 1440509_at 1443505_at 1449161_at
4.188947 2.455281 1.716046 1.069932 1.331611 1.931633
1450671_at 1452077_at 1452486_a_at 1457582_at 1457945_at 1459565_at
1.103439 2.344211 1.744559 1.380012 2.022879 1.830987
names(difference_list)
[1] "1417065_at" "1417210_at" "1419011_at" "1420453_at" "1420536_at"
[6] "1424903_at" "1426438_at" "1426439_at" "1426598_at" "1440509_at"
[11] "1443505_at" "1449161_at" "1450671_at" "1452077_at" "1452486_a_at"
[16] "1457582_at" "1457945_at" "1459565_at"
difference_exprs <- as.matrix(exprs[names(difference_list),])
difference_exprs
GSM349023.CEL.gz GSM349025.CEL.gz GSM349026.CEL.gz
1417065_at 8.068466 11.100941 10.818559
1417210_at 3.464524 7.376194 3.316209
1419011_at 7.912762 10.150693 6.194791
1420453_at 9.438416 10.468867 6.213329
1420536_at 9.049456 10.713780 6.812691
1424903_at 3.054923 6.006283 3.089801
1426438_at 2.961528 8.951067 2.902040
1426439_at 3.649746 7.373851 4.163506
1426598_at 2.981470 5.657400 3.295904
1440509_at 6.155688 7.915514 6.490324
1443505_at 4.283419 6.262506 4.264163
1449161_at 4.382851 7.634153 5.062025
1450671_at 9.424871 10.975197 6.938640
1452077_at 3.351106 6.911309 3.467514
1452486_a_at 9.342876 11.715102 7.086912
1457582_at 2.838758 4.631070 2.846881
1457945_at 3.962957 6.409204 3.932521
1459565_at 4.332023 6.348979 4.300042
GSM349028.CEL.gz GSM349029.CEL.gz GSM349030.CEL.gz
1417065_at 12.346885 13.181451 12.726629
1417210_at 9.140852 8.731609 7.767322
1419011_at 5.671617 5.056715 8.545940
1420453_at 4.412419 4.533255 8.504519
1420536_at 6.077600 6.804226 9.318515
1424903_at 6.743316 6.141642 5.827358
1426438_at 9.874542 9.239822 8.844623
1426439_at 8.089183 7.646497 7.362559
1426598_at 6.389696 6.154807 5.533222
1440509_at 6.881806 6.787127 7.845616
1443505_at 4.491399 4.102569 5.891079
1449161_at 4.955703 6.033426 8.683346
1450671_at 4.323317 4.933257 9.308570
1452077_at 8.174362 8.357335 7.122919
1452486_a_at 5.447898 5.774073 10.274539
1457582_at 5.535733 5.168933 4.827805
1457945_at 8.010864 7.414966 6.959013
1459565_at 7.661426 6.794672 6.909294
#よし！！発現が２倍以上上昇している遺伝子の発現データだけを取り出すことができたぞ！！！
d <- mean_ko - mean_wt
a <- rowMeans(exprs)
# d <- rowMeans(exprs[ ,c(2,4,6)]) - rowMeans(exprs[ ,c(21,3,5)])ってこと。
png(filename="110203_MA-plot.png", width = 600, height = 600, )
par(mfrow = c(3,2))
plot(a,d, ylim = c(-1.5,1.5), main = "MA-plot")
plot(a,d, ylim = c(-1.5,1.5), main = "MA-plot",pch = ".")
plot(a,d, ylim = c(-1.5,1.5), main = "MA-plot",pch="☆")
plot(a,d, ylim = c(-1.5,1.5), main = "MA-plot",pch="生")
plot(a,d, ylim = c(-3.0,3.0), xlim =c(0,20),main = "MA-plot")
plot(a,d, ylim = c(-3.0,3.0), xlim =c(0,20),main = "MA-plot",pch=".")
dev.off()
#スポットをいろんなもので打ってみて遊んでみました。
library(genefilter)
label <- c(0,1,0,1,0,1)
tt <- rowttests(difference_exprs,factor(label)) #t検定を行う
tt
statistic dm p.value
1417065_at -0.8791187 -1.368660 0.42897866
1417210_at -1.5724074 -2.924009 0.19095943
1419011_at -1.1183894 -1.734660 0.32602918
1420453_at -0.4654789 -1.066935 0.66580241
1420536_at -0.7342557 -1.147841 0.50351255
1424903_at -1.9766008 -2.096864 0.11926466
1426438_at -1.9684666 -4.188947 0.12038079
1426439_at -1.9208674 -2.455281 0.12714448
1426598_at -1.6433605 -1.716046 0.17565329
1440509_at -2.8145537 -1.069932 0.04809363
1443505_at -2.4556116 -1.331611 0.07002112
1449161_at -1.5980379 -1.931633 0.18527690
1450671_at -0.4629578 -1.103439 0.66746078
1452077_at -1.3827829 -2.344211 0.23891398
1452486_a_at -0.8066202 -1.744559 0.46511093
1457582_at -1.6776282 -1.380012 0.16872379
1457945_at -1.6213698 -2.022879 0.18025610
1459565_at -2.0130724 -1.830987 0.11439743
#符号が逆になっているのでlabelの1と０を交換してもよい。しなくてもよい。
