$ 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) #rma関数によりバックグラウンド補正と正規化を行う。
source("http://bioconductor.org/biocLite.R")
biocLite("genefilter")
library(genefilter)
#8つのアレイのうち、発現量が7(適当に設定した。カットオフが緩いと後で計算処理が大変!)を越えるものが3つある遺伝子をえらぶ。(発現量によるフィルタリング)
f1 <- kOverA(3,7) #発現値に基づいたフィルタを作成
ffun <- filterfun(f1) #フィルタ関数を作る
which <- genefilter(exprs(eset), ffun)
sum(which)
[1] 18378
#ものとのesetからファルタにヒットしたものだけを集めて、exprSetオブジェクトを生成する
eset.filtering <- eset[which]
#t検定によってさらにフィルタリングしていく
eset.filtering$sample
[1] 1 2 3 4 5 6 7 8
eset.filtering$sample <- c(1,2,3,1,2,3,1,3) #$sampleのスロットを書き換える!!!
eset.filtering$sample #書き換わった!!
[1] 1 2 3 1 2 3 1 3
ttest.filter <- ttest(eset.filtering$sample, p = 0.0001)
ttest.filterfun <- filterfun(ttest.filter)
which.ttest <- genefilter(exprs(eset.filtering), ttest.filterfun)
以下にエラー t.test.formula(x ~ m) :
grouping factor must have exactly 2 levels
#2群比較になっていないから怒られた(T_T)
#GEOの情報によると、wt/wt, ko/wt, ko/ko, wt/wt, ko/wt, ko/ko, wt/wt, ko/ko
#ならば、ヘテロノックアウト以外のアレイ(カラム)を抽出して新にオブジェクトを作ってしまえ!!!!
eset.filtering.extract <- eset.filtering[, c(1,3,4,6,7,8)]
eset.filtering.extract$sample
[1] 1 3 1 3 1 3 #シャー!!うまくいった!!!
eset.filtering.extract$sample <- c(1,2,1,2,1,2)
eset.filtering.extract$sample
[1] 1 2 1 2 1 2
ttest.filter <- ttest(eset.filtering.extract$sample, p = 0.05) #t検定にに基づいたフィルタを作成
ttest.filterfun <- filterfun(ttest.filter) #フィルタ関数を作る
which.ttest <- genefilter(exprs(eset.filtering.extract), ttest.filterfun)
sum(which.ttest)
[1] 82
eset.ttest <- eset.filtering.extract[which.ttest] #有意差検定によりフィルタリングしたものをExprSetオブジェクトとして生成する
expr.ttest <- exprs(eset.ttest) #発現データのマトリックスを取り出す。
png("110203_heatmap.png")
heatmap(expr.ttest)
dev.off()