読者です 読者をやめる 読者になる 読者になる

アメリエフのブログ

Amelieff Staff Blog

SRAをRで検索

毎日寒いですね。
家ではこたつから出られなくなって、こたつの中で生活している方もいらっしゃるのではないでしょうか。

今日は、RでNCBI SRAのデータを検索するRパッケージ「SRAdb」をご紹介します。

最近出たばかりのSRAdbの論文はこちらです。
Zhu Y, Stephens RM, Meltzer PS, Davis SR.
SRAdb: query and use public next-generation sequencing data from within R.
BMC Bioinformatics. 2013 Jan 17;14(1):19.


・インストール
source("http://bioconductor.org/biocLite.R")
biocLite("SRAdb")

・ライブラリ読み込み
library(SRAdb)

・検索
sqlfile = getSRAdbFile()  # 10分くらいかかります
sra_con = dbConnect(SQLite(), sqlfile)
rs <- getSRA(search_terms = "Human RNA-Seq", out_types = c("run", "study"), sra_con = sra_con)
# ↑"Human RNA-Seq"で検索する場合の例
head(rs[, c(1:5, 7)])

・結果
run_alias run run_date updated_date spots run_center
1 poly_1 DRR000008 2012-09-06 19153889
2 poly_3 DRR000010 2012-09-06 18079621
3 poly_2 DRR000009 2012-09-06 17668118
4 cyto_2 DRR000012 2012-09-06 15075547
5 cyto_1 DRR000011 2009-06-12 2012-09-06 15788579 UT-MGS
6 cyto_3 DRR000013 2009-06-15 2012-09-06 15490013 UT-MGS

結果が返ってきました!


面白いと思ったのは、SRAdbには、起動中のIGVの表示領域を切り替える機能もあることです。

※予めIGVを起動しておき、「View」→「Preferences」→「Advanced」の「Enable port」にチェックが入っていて、「60151」になっていることを確認してください。

・IGVの表示領域切り替え
myBam = file.path("path/to/bam")  
# bamと同じ場所にbam.baiも置いておきます
sock <- IGVsocket()
IGVgenome(sock, "hg19")
IGVload(sock, myBam)
IGVgoto(sock, "chr1:1-1000")

今はDBCLS SRAなど便利なサイトがたくさんあるので、SRAデータの検索をあえてRで行いたい場面はあまり思いつかないのですが、面白いパッケージだと思いました。

「こたつ」ならぬ「Rに入ったままなんでもやりたい」方にはよさそうです。