やりたいこと

遺伝子名を羅列したリストに、ensembl gene IDの情報も追加したい。

方法

人のゲノムデータベース（Ensembl）を使って遺伝子名（inputデータ）を、ensembl gene ID（outputデータ）に変換し、リストに追加する。

使用するツール

R
RパッケージbiomaRt

使用するデータをRに取り込み

以下のようにヒトの遺伝子名を羅列したリストをエクセルで作ってtxt形式で保存しました(gene_list.txt)。
目標は、gene_name列の隣にensembl gene IDを追加する。

リストのあるディレクトリに移動してデータを確認。

Terminal

% 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を起動する。

Terminal

% 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)

出来上がりはこんな感じ。

こちらの記事参考にしました