Rで遺伝子名をIDに変換する - アメリエフの技術ブログ

こんにちは、tsuyuhです。
最近は雨が降ったり止んだり、梅雨なのか梅雨じゃないのかよくわからない天気が続いています。個人的に雨は嫌いじゃないんですが、はっきりしない天気は嫌いです。

ところで、この間、とあるデータベースから遺伝子発現マトリックスを落としてきたら、それが "遺伝子名" と "発現量" の列で構成された表だったということがありました。
一見何も問題はないんですが、個人的には "遺伝子ID" は欲しいところかなと思います。

解析ツールによっては遺伝子IDを参照できないと不便なことがあります。
また、古い論文のデータだと、現在のデータベース上では遺伝子名が変更されていて、Synonymous扱いになっていることもよくあります。(""gene_A" not found" のように怒られます。)

そういうはっきりしない (?) 表は嫌いです。
という訳で、自分で好きな遺伝子IDを付与してみます。
あと、ついでに遺伝子の機能も追記してみようと思います。

RパッケージbiomaRtを使って遺伝子名に対応するIDや遺伝子機能を付与する

使用ツール
biomaRt v2.42.1
(Durinck et al., Bioinformatics 2005; Durinck et al., Nat. Protoc. 2009)

まず、こんな表を用意してみました。
とあるヒト遺伝子の発現量データです。
（※ 無作為に選んだ遺伝子と適当な数値です。）

gene_list.txt

gene_name	TPM
MIR421	5
SPATA21	10
CLEC4M	20
FUZ	100

（実際、これくらいの遺伝子数なら、Google先生に聞いた方が早いんですけどね。遺伝子が100個とか1000個とかの表でも同じことができますので...。）

ここでは、この表に "ensembl gene ID" と "entrez gene ID"、"遺伝子の機能詳細" を付与してみたいと思います。

まずは普通にRを起動して、

# ツールのインストール
> if (!requireNamespace("BiocManager", quietly = TRUE))  
       install.packages("BiocManager")
> BiocManager::install("biomaRt")

# ツールの読み込み
> library(biomaRt)

# データの読み込み
> data <- read.table("gene_list.txt", header = T)

次に、useMart() で参照するデータベースを選択します。今回はensemblでいい。
　※ というか、デフォルトでwww.ensembl.org を参照している模様。 (以前のバージョンだとensembl以外も参照できたんですが...。ensembl fungiとか。)

> db <- useMart("ensembl")

そして、listDatasets("データベース名") で、そのデータベースに登録されているデータを確認します。

> listDatasets(db)
       dataset                              description                 　　version
1     abrachyrhynchus_gene_ensembl          Pink-footed goose genes (ASM259213v1)         ASM259213v1
2         acalliptera_gene_ensembl          Eastern happy genes (fAstCal1.2)              fAstCal1.2
3       acarolinensis_gene_ensembl          Green anole genes (AnoCar2.0v2)               AnoCar2.0v2
4        acchrysaetos_gene_ensembl          Golden eagle genes (bAquChr1.2)               bAquChr1.2
(以下略)

ensemblでは、すべてのデータが "[生物種]_gene_ensembl" という名前になっていることがわかります。
ヒトなら "hsapiens_gene_ensembl" ですね。

ここまでわかったので、目的のデータを useDataset("データ名", mart = "データベース名") で読み込みます。

> hd <- useDataset("hsapiens_gene_ensembl", mart = db)

次に変換元と変換先の属性名を調べます。

今回は、変換元が "遺伝子名"、変換先が"ensembl gene ID" と "entrez gene ID"、"遺伝子の機能詳細" ですが、hsapiens_gene_ensemblではどんな名称で登録されているのかわかりません。

変換元の属性名を調べるのはlistFilters("データ名")

> listFilters(hd)
                                     name
1                         chromosome_name
2                                   start
3                                     end
4                              band_start
5                                band_end
6                            marker_start
7                              marker_end
(以下略)

ですが、種類がありすぎる上、名称に統一感がないので探すのは困難です。
そこで、grep() であたりをつけて検索します。
おそらく、遺伝子名は"gene_symbol" みたいな名称でしょう。

> filters <- listFilters(hd)

> filters[grep("symbol", filters[,1]),]
                 name                                description
83        hgnc_symbol                 HGNC symbol(s) [e.g. A1BG]
114 uniprot_gn_symbol UniProtKB Gene Name symbol(s) [e.g. 5HT1A]

検索した結果、"hgnc_symbol" という名前で登録されていることがわかりました。

次に、変換先の属性名はlistAttributes("データ名") でそれぞれ調べます。
検索ワードは適当に。

> attributes <- listAttributes(hd)

> attributes[grep("ensembl", attributes[,1]),]
                           name                  description         page
1               ensembl_gene_id               Gene stable ID feature_page
2       ensembl_gene_id_version       Gene stable ID version feature_page
3         ensembl_transcript_id         Transcript stable ID feature_page
4 ensembl_transcript_id_version Transcript stable ID version feature_page
5            ensembl_peptide_id            Protein stable ID feature_page
6    ensembl_peptide_id_version    Protein stable ID version feature_page
(以下略)

> attributes[grep("entrez", attributes[,1]),]
                     name                                 description          page
60  entrezgene_trans_name               EntrezGene transcript name ID      feature_page
80 entrezgene_description NCBI gene (formerly Entrezgene) description      feature_page
81   entrezgene_accession   NCBI gene (formerly Entrezgene) accession      feature_page
82          entrezgene_id          NCBI gene (formerly Entrezgene) ID      feature_page

> attributes[grep("description", attributes[,1]),]
                     name                                 description         page
8             description                            Gene description      feature_page
38  phenotype_description                       Phenotype description      feature_page
50 goslim_goa_description                      GOSlim GOA Description      feature_page
73   mim_gene_description                        MIM gene description      feature_page
75 mim_morbid_description                      MIM morbid description      feature_page
80 entrezgene_description NCBI gene (formerly Entrezgene) description      feature_page
(以下略)

それぞれ、 "ensembl_gene_id"、"entrezgene_id" 、"description" という名前を使えばよさそうです。

ここまでわかったので、実際に遺伝子名をそれぞれに変換します。

変換するコマンドは getBM(attributes = "変換先の属性名", filters = "変換元の属性名", values = "変換元", mart = "参照データ") です。
ここでは、"遺伝子名" と "ensembl gene ID"、"entrez gene ID"、"遺伝子の機能詳細" を横並びにしたいので、attributes = c(hgnc_symbol", "ensembl_gene_id", "entrezgene_id", "description") としています。
また、変換元のデータは用意した表に記載されている遺伝子名なので、data変数内のgene_name列を使います。
　　※ なお、キャッシュの読み込みと書き込みを無効 (useCache = FALSE) にしないとエラーが出ます。

> res <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id", 
         "entrezgene_id", "description"), 
         filters = "hgnc_symbol", values = data$gene_name, 
         mart = hd, useCache = FALSE)

# 結果の確認
> res
  hgnc_symbol ensembl_gene_id entrezgene_id
1      CLEC4M ENSG00000104938         10332
2         FUZ ENSG00000010361         80199
3      MIR421 ENSG00000202566        693122
4     SPATA21 ENSG00000187144        374955
                                                                 description
1 C-type lectin domain family 4 member M [Source:HGNC Symbol;Acc:HGNC:13523]
2     fuzzy planar cell polarity protein [Source:HGNC Symbol;Acc:HGNC:26219]
3                           microRNA 421 [Source:HGNC Symbol;Acc:HGNC:32793]
4          spermatogenesis associated 21 [Source:HGNC Symbol;Acc:HGNC:28026]

できました！
あとは、この結果と発現マトリックスを合体させるだけ。

# 遺伝子名でソート
> data <- data[order(data$gene_name),]
> res <- res[order(res$hgnc_symbol),]

# dataとresで遺伝子名が一致しているか確認
> table_1 <- data[match(data$gene_name, res$hgnc_symbol),]
> table_2 <- res[match(data$gene_name, res$hgnc_symbol),]
> table_1
  gene_name TPM
3    CLEC4M  20
4       FUZ 100
1    MIR421   5
2   SPATA21  10

> table_2
  hgnc_symbol ensembl_gene_id entrezgene_id
1      CLEC4M ENSG00000104938         10332
2         FUZ ENSG00000010361         80199
3      MIR421 ENSG00000202566        693122
4     SPATA21 ENSG00000187144        374955
                                                                 description
1 C-type lectin domain family 4 member M [Source:HGNC Symbol;Acc:HGNC:13523]
2     fuzzy planar cell polarity protein [Source:HGNC Symbol;Acc:HGNC:26219]
3                           microRNA 421 [Source:HGNC Symbol;Acc:HGNC:32793]
4          spermatogenesis associated 21 [Source:HGNC Symbol;Acc:HGNC:28026]

# 遺伝子名が一致していたので、合体
> cbind(table_2, table_1$TPM)
  hgnc_symbol ensembl_gene_id entrezgene_id
1      CLEC4M ENSG00000104938         10332
2         FUZ ENSG00000010361         80199
3      MIR421 ENSG00000202566        693122
4     SPATA21 ENSG00000187144        374955
                                                                 description
1 C-type lectin domain family 4 member M [Source:HGNC Symbol;Acc:HGNC:13523]
2     fuzzy planar cell polarity protein [Source:HGNC Symbol;Acc:HGNC:26219]
3                           microRNA 421 [Source:HGNC Symbol;Acc:HGNC:32793]
4          spermatogenesis associated 21 [Source:HGNC Symbol;Acc:HGNC:28026]
  table_1$TPM
1          20
2         100
3           5
4          10

完成！

ちなみに、dataとresの遺伝子名が一致しないという場合もよくあります。
特に、non-coding RNA遺伝子などアノテーションがはっきり付与されていない遺伝子はbiomaRtでデータベースを検索しても見つかりません。
そんなときは、以下のようにします。

# 例えば、FUZ遺伝子がデータベースになかった場合、
> table_3 <- data[match(data$gene_name, res$hgnc_symbol),]
   gene_name TPM
3     CLEC4M  20
NA      <NA>  NA
1     MIR421   5
2    SPATA21  10

# このように、その行がNAと表示されるので、na.omit()関数でその行を削除します。
> table_3 <- na.omit(table_3)
> table_3
  gene_name TPM
3    CLEC4M  20
1    MIR421   5
2   SPATA21  10

これで、もう片方の表でも同様のことをして合体すればOK。
ただ、あまりにも検索で抜け落ちる場合には、属性名が間違っている可能性があるので、別の属性名で試してみるのがいいかもですね。

ちなみに、同様の方法でGene Ontologyも付与することができます。
当ブログの過去記事↓ staffblog.amelieff.jp