いろいろ倉庫

KNIME、EXCEL、R、Pythonなどの備忘録

【R】遺伝子のIDを変換したり検索したい

・お題:データセットによって遺伝子のIDの記載方法が異なることがある。くっつけて解析する際に困るので変換・統一する方法を調べてみたい。

 

・biomaRtというパッケージが使えるらしい。公式サイトは以下。

bioconductor.org

・今回参考にするチュートリアルは以下。インストールは割愛。

http://127.0.0.1:15967/library/biomaRt/doc/accessing_ensembl.html

 

・ライブラリを読み込む。

> library(biomaRt)

・まずは必要なデータベースを選ぶために、利用可能なBioMartのサービスを確認する。

> listEnsembl()

・今回は、genesを使うことにする。useEnsembl関数を使って読み込んでみる。

> ensembl <- useEnsembl(biomart = "genes")
> ensembl 
Object of class 'Mart':
  Using the ENSEMBL_MART_ENSEMBL BioMart database
  No dataset selected.

・「データセットを指定していないよ」といわれてしまった。データセットの一覧は、以下のコマンドで確認することができる。

> datasets <- listDatasets(ensembl)
> datasets

中略

後略

・ヒトのデータセットhsapiens_gene_ensemblが存在することが分かった。

・データセット名が分かったので、データセットを読み込むことができる。

> ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
Ensembl site unresponsive, trying useast mirror

・うまくデータにアクセスできなかったので、ミラーを変更したらしい。これ以上のエラーが出なかったので、多分うまくデータセットをダウンロードできたのだろう。

> ensembl
Object of class 'Mart':
  Using the ENSEMBL_MART_ENSEMBL BioMart database
  Using the hsapiens_gene_ensembl dataset

・次に、クエリを作成してgetBM()関数でEnsembl BioMartサーバーに投げる。getBM()関数には、 最低3つの引数(values、filters、attribute)を入れ込む必要がある。

・valuesは検索する元のIDを入力する。

・filterは検索対象のデータセットを制限するための引数らしい。filterには以下のような選択肢がある。今回は、AFFY HG U133A 2 probeのIDを意味が分かる遺伝子名などに変換したいので、filterはaffy_hg_u133a_2になる。

> listFilters(ensembl)

中略

後略
・attributesは変換後の遺伝子名などの属性を指定するための引数らしい。attributeには以下のような選択肢がある。今回はentrezgene_idとwikigene_nameにする。

> listAttributes(ensembl)

後略
・これらで得られた引数をgetBMに渡して、IDを変換してみる。

> affyids <- c("202763_at","209310_s_at","207500_at")

> getBM(attributes = c('affy_hg_u133_plus_2', 'entrezgene_id', 'wikigene_name'),
+       filters = 'affy_hg_u133_plus_2',
+       values = affyids, 
+       mart = ensembl)

・なんとかうまくいったっぽい。

・ここまでで、何度か表から入力する引数を頑張って探す作業があった。serch●●関数で検索できるらしい。

> searchDatasets(mart = ensembl, pattern = "hsapiens")

> searchAttributes(mart = ensembl, pattern = "hgnc")

> searchFilters(mart = ensembl, pattern = "ensembl.*id")

 

・せっかくなので、何パターンか試してみる。

・EntrezGene IDにGOアノテーションをくっつける。

> entrez=c("673","837")
> goids = getBM(attributes = c('entrezgene_id', 'go_id'), 
+               filters = 'entrezgene_id', 
+               values = entrez, 
+               mart = ensembl)

> goids

・“MAP kinase activity”というGOタームに関連する遺伝子のentrezgene_idとhgnc_symbolをとって来る。

> getBM(attributes = c('entrezgene_id','hgnc_symbol'), 
+       filters = 'go', 
+       values = 'GO:0004707', 
+       mart = ensembl)

・EntrezGene IDから、100bp上流のプロモーター配列を釣ってくる。
> entrez=c("673","7157","837")
> getSequence(id = entrez, 
+             type="entrezgene_id",
+             seqType="coding_gene_flank",
+             upstream=100, 
+             mart=ensembl


・いろいろな使い方ができそう。

 

 

おわり