ExprSetからmatixに変換したのち、FC法やらt検定やらでフィルタリングする方法

#ExpressionSetオブジェクトではなく、マトリックスに対してFold Change法によりフィルタリングするとともに、t検定により発現変動遺伝子を同定する。

$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なので2群比較を視野にいれて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)  #2倍以上高発現しているプローブ数をしらべる。ちなみに、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

#よし!!発現が2倍以上上昇している遺伝子の発現データだけを取り出すことができたぞ!!!



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と0を交換してもよい。しなくてもよい。