GEOから遺伝子発現データを取ってくる【R】

#GEOとは
https://www.ncbi.nlm.nih.gov/geo/
遺伝子発現データ（マイクロアレイ,次世代シークエンス)のレポジトリです
transcriptomeだけではなく, メチル化アレイのデータもあります

#データの取得
Rを使ってダウンロードします

##データの概要
GSE26174
ラットの網膜に電気刺激(transcorneal electrical stimulation, TES)を与えた時の効果を調べたもの

Willmann G, Schäferhoff K, Fischer MD, Arango-Gonzalez B et al. Gene expression profiling of the retina after transcorneal electrical stimulation in wild-type Brown Norway rats. Invest Ophthalmol Vis Sci 2011 Sep 29;52(10):7529-37. PMID: 21873684

処理群4サンプル, 対照群4サンプル　合計8サンプル

##コード
ライブラリの読み込み

library(Biobase)
library(GEOquery)
library(data.table)

GEOのウェブサイトからデータを取得します
getGEO()の引数 GSEMatrix = F とするとサンプルごとの数値データ, GSEMatrix = T とするとデータセットの概要(GEOのaccessionのページではseries matrix file)を取得します

gsem <- getGEO("GSE26174",GSEMatrix = T)
gse <- getGEO("GSE26174",GSEMatrix = F) 

###サンプル名とか実験条件の取得
gsemはlistオブジェクトなので, listの1番目を取り出します

x <- pData(phenoData(gsem[[1]]))

得られた表はこんな感じ

> names(x)
 [1] "title"                   "geo_accession"           "status"                  "submission_date"        
 [5] "last_update_date"        "type"                    "channel_count"           "source_name_ch1"        
 [9] "organism_ch1"            "characteristics_ch1"     "characteristics_ch1.1"   "characteristics_ch1.2"  
[13] "molecule_ch1"            "extract_protocol_ch1"    "label_ch1"               "label_protocol_ch1"     
[17] "taxid_ch1"               "hyb_protocol"            "scan_protocol"           "data_processing"        
[21] "platform_id"             "contact_name"            "contact_email"           "contact_institute"      
[25] "contact_address"         "contact_city"            "contact_zip/postal_code" "contact_country"        
[29] "supplementary_file"      "supplementary_file.1"    "data_row_count"          "age:ch1"                
[33] "gender:ch1"              "strain:ch1"             

> x[1:3]
                     title geo_accession                status
GSM642729 TES Control 1 4h     GSM642729 Public on Nov 21 2011
GSM642730 TES Control 2 4h     GSM642730 Public on Nov 21 2011
GSM642801 TES Control 3 4h     GSM642801 Public on Nov 21 2011
GSM642802 TES Control 4 4h     GSM642802 Public on Nov 21 2011
GSM642803         TES 1 4h     GSM642803 Public on Nov 21 2011
GSM642804         TES 2 4h     GSM642804 Public on Nov 21 2011
GSM642805         TES 3 4h     GSM642805 Public on Nov 21 2011
GSM642806         TES 4 4h     GSM642806 Public on Nov 21 2011

```
ここで必要になるのは各サンプルのGEOaccessionと実験条件です
このデータセットでは"title"列に実験条件が書かれていますが, それぞれのデータセットによって表記がいろいろなので, 必ずseries matrixを確認します

###数値データの取得
gseは各サンプルのデータが入ったGSEインスタンスです
gse@gsemにマイクロアレイの数値データ（プロセシング済み)が入っています
Table()でプローブIDと対応する数値データを取り出し, サンプルごとにforループで読んでいきます

```R
y <- gse@gsms

dat <- Table(y[[1]])
dat <- dat[,1:2]

for (i in 2:length(y)){
  sub_dat <- Table(y[[i]])
  dat <- cbind(dat,sub_dat[,2])
}
dat <- data.frame(dat,row.names = 1)
colnames(dat) <- x$geo_accession
```
サンプルごとの数値データをdata.frameとして得ることができました

```
> dat[1:3,]
         GSM642729 GSM642730 GSM642801 GSM642802 GSM642803 GSM642804 GSM642805 GSM642806
10700001  10.34770  10.22240  10.23050  10.37600  10.55910  10.29450  10.30840  10.16820
10700002   2.77305   2.85680   2.80159   2.64945   2.72091   2.77305   2.72123   2.77305
10700003   8.86587   8.91388   8.85737   8.87960   9.09890   8.91388   8.85572   8.82685
```

###保存と次回の読み込み
Rdata形式にしておくと読み込みがはやいのでsave()で保存します

```R
save(dat,file = "GSE26714.Rdata")
```

次回からはload()で読み込むとgetGEO()から戻らなくてよくなります（しかも速い）

```R
load(file = "GSE26714.Rdata")
```