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

More than 1 year has passed since last update.


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ループで読んでいきます

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()で保存します

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

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

load(file = "GSE26714.Rdata")