Microarrayの解析 - fukuitのブログ

とりあえず、NCBI GEOで公開されているMicroarrayのデータを使って、特徴的な遺伝子のリストを得るための、覚え書き。
例題は、GSE16515として公開されている、膵臓癌のデータを使うことにする。

library(GEOquery)
gse16515.eset <- getGEO("GSE16515")[[1]]
gse16515.eset.exp <- exprs(gse16515.eset)

で、このデータはTumorとNormalが整然と並んでいないので、面倒だが、何番目のデータがTumorなのかNormalなのかを保持する変数を持つ。なんか上手い方法があれば良いんだけれど、思いつかないので、地道に手で打つ。ちなみに、Tumorサンプルは、Tumor Positiveの略でTpos、逆にNormalサンプルはTnegということにする。

Tpos <- c(1,2,3, 4, 6, 8, 10, 12, 13, 14, 16, 18, 20, 21, 22, 23, 25, 26, 27, 28, 29, 31, 33, 35, 36, 37, 38, 39, 41, 42, 44, 45, 46, 48, 50, 51)
Tneg <- c(5, 7, 9, 11, 15, 17, 19, 24, 30, 32, 34, 40, 43, 47, 49, 52)

Tpos群とTneg群で全プローブセットについて、発現量の差を検定する。まあ、なんとなく、ここはWilcoxonの検定で。

pval.w <- apply(gse16515.eset.exp, 1, function(x)
 wilcox.test(x[Tpos], x[Tneg])$p.value)

Tpos群、Tneg群ごとに、各プローブの平均値を算出した上で、その倍差(Fold-change)を求める(fc)。logなので、比の算出は引き算で出る。

Fc <- 2^(rowMeans(gse16515.eset.exp[,Tpos]) - rowMeans(gse16515.eset.exp[,Tneg]))

発現差のあるプローブをピックアップする。ピックアップする基準は、p値が0.01以下で、かつFold-changeが5倍以上のもの。

gse16515.selected_probe.fc <- names(pval.w[pval.w < 0.01 & abs(log2(fc)) > log2(5)])

あるいは、siggensを使う。

gse16515.eset.cl <- c(rep(1,52))
for (i in Tpos) { gse16515.eset.cl[Tpos[i]] <- 0 }
sam_res <- sam(gse16515.eset.exp, gse16515.eset.cl, B=1000)
gse16515.selected_probe.siggenes <- list.siggenes(sam_res, 4.5)

それぞれのプローブが何なのか？の情報を得る。

getBioC("hgu133plus2.db")
library(hgu133plus2.db)
get(gse16515.selected_probe.fc[1], env=hgu133plus2SYMBOL) # プローブセットの遺伝子名
get(gse16515.selected_probe.fc[1], env=hgu133plus2GENENAME) # プローブセットの遺伝子名
get(gse16515.selected_probe.fc[1], env=hgu133plus2GO) # プローブセットの機能註釈
get(gse16515.selected_probe.fc[1], env=hgu133plus2ENTREZID) # プローブセットのENTREZID

heatmapも描いちゃう。

gse16515.scaled.sig <- t(scale(t(gse16515.eset.exp[gse16515.selected_probe.sig,])))
heatmap(gse16515.scaled.sig, scale="none", col=cm.colors(20))

以上のような作業のの結果、selected_probeの中にprobe「1554436_a_at」があり、これはREG4という遺伝子の発現量が特異的であることを示しているようだ。なるほど。

実際に、「ＲＥＧ４タンパク質を用いて膵癌を診断する方法」という公開特許もあるようなので、それなりに正しそうな結果が得られたようだ。