R(とBioConductor)とGEO
R(とBioConductor)で、GEOにて公開されているマイクロアレイのデータを読み込んでみる企画。
準備
Rを起動して、GEOqueryを使えるようにしておく。
> Sys.setenv("http_proxy"="http://proxy.example.com/") > Sys.setenv("ftp_proxy"="http://proxy.example.com/") > source("http://bioconductor.org/getBioC.R") > biocLite("GEOquery") > library(GEOquery)
GEOqueryによるマイクロアレイデータの取得
ひとまず、サンプルとしてGSE7505のデータを取得してみることにする。
> gse <- getGEO( 'GSE7505' ) 以下にエラー curlPerform(curl = curl, .opts = opts, .encoding = .encoding) : FTP response reading failed
というワケで実行不可能。proxyを超えられないようだ。
> getGEO function (GEO = NULL, filename = NULL, destdir = tempdir(), GSElimits = NULL, GSEMatrix = TRUE, AnnotGPL = FALSE) { filename <- filename con <- NULL if (!is.null(GSElimits)) { if (length(GSElimits) != 2) { stop("GSElimits should be an integer vector of length 2, like (1,10) to include GSMs 1 through 10") } } if (is.null(GEO) & is.null(filename)) { stop("You must supply either a filename of a GEO file or a GEO accession") } if (is.null(filename)) { GEO <- toupper(GEO) geotype <- toupper(substr(GEO, 1, 3)) if (GSEMatrix & geotype == "GSE") { return(getAndParseGSEMatrices(GEO)) } filename <- getGEOfile(GEO, destdir = destdir, AnnotGPL) } if (length(grep("\\.gz$", filename, perl = TRUE)) > 0) { con <- gzfile(filename, open = "rt") } else { con <- file(filename, "r") } ret <- parseGEO(con, GSElimits) close(con) return(ret) } <environment: namespace:GEOquery>
どうやら、filename= というオプションがあるので、ファイルをダウンロードしておけば良さそうだ。
> gse <- getGEO( filename = "GSE7505_family.soft.gz" ) Parsing.... ^PLATFORM = GPL5080 ^SAMPLE = GSM181868 (中略) ^SAMPLE = GSM181929 50 件以上の警告がありました (警告を見るには warnings() を使って下さい) > warnings() 警告メッセージ： 1: In grep("\\.gz$", filename, perl = TRUE) : perl=TRUE は UTF-8 ロケールに対してのみ完全実装されています (略)
ということで、どういうわけだか、今度はちゃんと読み込むことができたようだ。
> names(GSMList(gse7505)) [1] "GSM181867" "GSM181868" "GSM181869" "GSM181870" "GSM181871" "GSM181872" [7] "GSM181873" "GSM181874" "GSM181875" "GSM181876" "GSM181877" "GSM181878" [13] "GSM181879" "GSM181880" "GSM181881" "GSM181882" "GSM181883" "GSM181884" [19] "GSM181885" "GSM181886" "GSM181887" "GSM181888" "GSM181889" "GSM181890" [25] "GSM181891" "GSM181892" "GSM181893" "GSM181894" "GSM181895" "GSM181896" [31] "GSM181897" "GSM181898" "GSM181899" "GSM181900" "GSM181901" "GSM181902" [37] "GSM181903" "GSM181904" "GSM181905" "GSM181906" "GSM181907" "GSM181908" [43] "GSM181909" "GSM181910" "GSM181911" "GSM181912" "GSM181913" "GSM181914" [49] "GSM181915" "GSM181916" "GSM181917" "GSM181918" "GSM181919" "GSM181920" [55] "GSM181921" "GSM181922" "GSM181923" "GSM181924" "GSM181925" "GSM181926" [61] "GSM181927" "GSM181928" "GSM181929" >
データの変換
gse7505をRで発現データを扱う時に使われるexprSetに変換する。
もうGEOqueryのドキュメント(pdf)に載ってるそのまんま、ワケも分からずに入力する。
> gsmplatforms <- lapply(GSMList(gse7505), function(x) { + Meta(x)$platform + }) > Table(GSMList(gse7505)[[1]])[1:5, ] > Columns(GSMList(gse7505)[[1]])[1:5, ] > probesets <- Table(GPLList(gse7505)[[1]])$ID > data.matrix <- log2(do.call("cbind", lapply(GSMList(gse7505), function(x) { + tab <- Table(x) + mymatch <- match(probesets, tab$ID_REF) + return(tab$VALUE[mymatch]) + }))) > require(Biobase) > rownames(data.matrix) <- probesets > colnames(data.matrix) <- names(GSMList(gse7505)) > pdata <- data.frame(samples = names(GSMList(gse7505))) > rownames(pdata) <- names(GSMList(gse7505)) > pheno <- new("phenoData", pData = pdata, varLabels = as.list("samples")) エラー： The phenoData class is defunct, use AnnotatedDataFrame (with ExpressionSet) instead
いやなところでエラーが出る。phenoDataクラスは廃止になったからAnnotatedDataFrameを使え、と。
> pheno <- new("AnnotatedDataFrame", pData = pdata, varLabels = as.list("samples")) 以下にエラー .nextMethod(.Object, ...) : クラス "AnnotatedDataFrame" のスロットに対しては不正な名前です: pData, varLabels >
ちょっと、考えないといけないな。
あれ？[Bioc-devel] Refactored Biobase in BioC 2.2 - exprSet and phenoData classes defunct, use ExpressionSet and AnnotatedDataFrameというのがあって、BioConductor2.2以降では、exprSetとphenoDataを使うこと自体が良くないようだ。ここで、今使っているBioConductor2.4を2.1以前に戻すというのは、後向きすぎる解決手段だ。
ExpressionSetとAnnotatedDataFrameを使うにはどうするか？を、ちょっと調べてみることにしよう。
とりあえず、検索してみると、[BioC] biobase phenodataに次のように書いてある。
This will get you an ExpressionSet:
https://stat.ethz.ch/pipermail/bioconductor/2008-September/024099.html
gse6901.eset <- getGEO("GSE6901")[[1]]
Now, gse6901.eset is an ExpressionSet.
なんか、簡単じゃん。やりなおしてみることに。
> gse7505 <- getGEO( filename="GSE7505_family.soft.gz" )[[1]]; Parsing.... ^PLATFORM = GPL5080 ^SAMPLE = GSM181868 (中略) ^SAMPLE = GSM181929 以下にエラー getGEO(filename = "GSE7505_family.soft.gz")[[1]] : この S4 クラスは部分代入可能ではありません
どうせいっちゅうんじゃ。
自宅で実行
自宅なら、proxyを使わないので、もっと簡単にできるはず。
library(GEOquery) gse7505.eset <- getGEO('GSE7505')[[1]] Found 1 file(s) GSE7505_series_matrix.txt.gz URL 'ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE7505/GSE7505_series_matrix.txt.gz' を試しています ftp data connection made, file length 2036446 bytes 開かれた URL ================================================== downloaded 1.9 Mb URL 'http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?targ=self&acc=GPL5080&form=text&view=full' を試しています Content type 'geo/text' length unknown 開かれた URL .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... .......... ... downloaded 603 Kb File stored at: /var/folders/XX/XXaRkr7wFNma1xCWGoA3fU+++TI/-Tmp-//RtmpbgXdK8/GPL5080.soft > gse7505.eset ExpressionSet (storageMode: lockedEnvironment) assayData: 6728 features, 63 samples element names: exprs phenoData sampleNames: GSM181867, GSM181868, ..., GSM181929 (63 total) varLabels and varMetadata description: title: NA geo_accession: NA ...: ... data_row_count: NA (41 total) featureData featureNames: 1, 2, ..., 6728 (6728 total) fvarLabels and fvarMetadata description: ID: NA CLONE.ID: NA ...: ... SPOT_ID: NA (8 total) additional fvarMetadata: Column, Description experimentData: use 'experimentData(object)' Annotation: GPL5080 >
あっさり完了。会社で数時間悩んだのはなんだったんだ。どうやら、getGEO()に読み込ませるファイルを間違えていたようだ。
BioConductorサイコー！
(多分、続く)
