やりたいこと
遺伝子名を羅列したリストに、ensembl gene IDの情報も追加したい。
方法
人のゲノムデータベース（Ensembl）を使って遺伝子名（inputデータ）を、ensembl gene ID（outputデータ）に変換し、リストに追加する。
使用するツール
R
RパッケージbiomaRt
使用するデータをRに取り込み
以下のようにヒトの遺伝子名を羅列したリストをエクセルで作ってtxt形式で保存しました(gene_list.txt)。
目標は、gene_name列の隣にensembl gene IDを追加する。
リストのあるディレクトリに移動してデータを確認。
% cd Desktop
% ls
gene_list.txt
% cat gene_list.txt
pmid result gene_name
33794282 down-regurated PDYN
33794282 ND OPRK1
33794282 no-changes TAC3
33794282 ND TACR3
33794282 ND KISS1R
33794282 down-regurated KISS1
33794282 ND EMR5P
33794282 ND EMR4P%
ターミナルでRを起動する。
% R
R version 4.1.1 (2021-08-10) -- "Kick Things"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
R は、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。
R は多くの貢献者による共同プロジェクトです。
詳しくは 'contributors()' と入力してください。
また、R や R のパッケージを出版物で引用する際の形式については
'citation()' と入力してください。
'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。
>
使用する遺伝子リスト（gene_list.txt ）を読み込み、変数dataに代入。
遺伝子リストにヘッダー(gene_name)がついているのでheader = T(TRUE)とする。
> data <- read.table("gene_list.txt", header = T)
使用するRパッケージbiomaRtをインストールし読み込み。
## ツールのインストール
> if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
> BiocManager::install("biomaRt")
# ツールの読み込み
> library(biomaRt)
biomaRtツールのuseMart() 関数を使って、参照するデータベース（今回はensembl）に接続する。
biomaRtツールのlistDatasets()関数でデータベースに登録されているデータセットを確認する。
> listDatasets(useMart("ensembl"))
dataset
1 abrachyrhynchus_gene_ensembl
2 acalliptera_gene_ensembl
3 acarolinensis_gene_ensembl
4 acchrysaetos_gene_ensembl
：
：(省略)
：
79 hhucho_gene_ensembl
80 hsapiens_gene_ensembl
81 ipunctatus_gene_ensembl
(以下省略)
今回使うヒトのdataset名は "hsapiens_gene_ensembl" （80行目にあった）ということがわかった。
useMart() については以下参照。
https://www.rdocumentation.org/packages/biomaRt/versions/2.28.0/topics/useMart
listDatasets() については以下参照。
https://www.rdocumentation.org/packages/biomaRt/versions/2.28.0/topics/listDatasets
目的のdatasetを読み込む。
> hd <- useDataset("hsapiens_gene_ensembl", mart = useMart("ensembl"))
useDataset() については以下参照。
https://www.rdocumentation.org/packages/biomaRt/versions/2.28.0/topics/useDataset
今回は、遺伝子名（inputデータ）を、ensembl gene ID（outputデータ）に変換したい。
biomaRtでの変換可能なinputデータ属性はlistFilters()関数で、outputデータ属性はlistAttributes()関数で確認できる。
まずは、inputデータとして使用できるデータの属性名をlistFilters()関数で表示させてみる。
目標：今回inputデータとして使う遺伝子名が入っているデータの属性名が知りたい。
> listFilters(hd)
name
1 chromosome_name
2 start
3 end
4 band_start
5 band_end
6 marker_start
7 marker_end
8 encode_region
9 strand
10 chromosomal_region
11 with_biogrid
12 with_ccds
(以下省略)
今回inputデータとして使用する遺伝子名の属性はhgnc_symbol。
83 hgnc_symbol
次に、outputデータとして使用できるデータの属性名をlistAttributes()関数で表示させてみる。
目標：今回outputデータとして使うensembl gene IDの属性名が知りたい。
> listAttributes(hd)
name
1 ensembl_gene_id
2 ensembl_gene_id_version
3 ensembl_transcript_id
4 ensembl_transcript_id_version
5 ensembl_peptide_id
(以下省略)
今回目的としているoutputデータのensembl gene IDの属性名はensembl_gene_id。
ついでに、遺伝子名の属性名はhgnc_symbolであることも確認しました。
1 ensembl_gene_id
65 hgnc_symbol
最後に、ここまでで調べてきた属性名を使用して遺伝子目からensembl IDを取得。
遺伝子名からensembl IDを取得する方法はgetBM（）関数を使う。使い方は以下の通り。
getBM(attributes = "outputの属性名", filters = "inputの属性名", values = "inputデータ", mart = "参照データベース", useCache = FALSE)
※ キャッシュの読み込みと書き込みを無効 (useCache = FALSE) にしないとエラーが出るそうです。
> res <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id"), filters = "hgnc_symbol", values = data$gene_name, mart = hd, useCache = FALSE)
> res
hgnc_symbol ensembl_gene_id
1 KISS1 ENSG00000170498
2 KISS1R ENSG00000116014
3 OPRK1 ENSG00000082556
4 PDYN ENSG00000101327
5 TAC3 ENSG00000166863
6 TACR3 ENSG00000169836
元のinputデータ（data）と比較してみるとoutputデータ(res)の方が数が少ない。
(今回inputの遺伝子名の中でデータベースにid登録されていなかった遺伝子があるためことが原因です)
inputデータ（data）とoutputデータ（res）を遺伝子名で並び替えてから精査する。
> data_sort <- data[order(data$gene_name),]
> data_sort
pmid result gene_name
8 33794282 ND EMR4P
7 33794282 ND EMR5P
6 33794282 down-regurated KISS1
5 33794282 ND KISS1R
2 33794282 ND OPRK1
1 33794282 down-regurated PDYN
3 33794282 no-changes TAC3
4 33794282 ND TACR3
> res_sort <- res[order(res$hgnc_symbol),]
> table <- res_sort[match(data_sort$gene_name, res_sort$hgnc_symbol),]
> table
hgnc_symbol ensembl_gene_id
NA <NA> <NA>
NA.1 <NA> <NA>
1 KISS1 ENSG00000170498
2 KISS1R ENSG00000116014
3 OPRK1 ENSG00000082556
4 PDYN ENSG00000101327
5 TAC3 ENSG00000166863
6 TACR3 ENSG00000169836
検索結果とinputデータを合わせて完成品を眺めてみる。
> summary <- cbind(data_sort,table$ensembl_gene_id)
> summary
pmid result gene_name table$ensembl_gene_id
8 33794282 ND EMR4P <NA>
7 33794282 ND EMR5P <NA>
6 33794282 down-regurated KISS1 ENSG00000170498
5 33794282 ND KISS1R ENSG00000116014
2 33794282 ND OPRK1 ENSG00000082556
1 33794282 down-regurated PDYN ENSG00000101327
3 33794282 no-changes TAC3 ENSG00000166863
4 33794282 ND TACR3 ENSG00000169836
完璧ですね！
gene_nameの隣の列に見事ensembl_gene_idが追加されました。
データベースに載ってないものはNAとなっています。
最後にcsvがファイルとして出力。
> write.csv(summary,"res.txt",quote=F)
出来上がりはこんな感じ。
こちらの記事参考にしました