【R】遺伝子のIDを変換したり検索したい
・お題:データセットによって遺伝子のIDの記載方法が異なることがある。くっつけて解析する際に困るので変換・統一する方法を調べてみたい。
・biomaRtというパッケージが使えるらしい。公式サイトは以下。
・今回参考にするチュートリアルは以下。インストールは割愛。
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)
・いろいろな使い方ができそう。
おわり