【R】遺伝子IDおよび遺伝子シンボルの相互変換【biomaRt】

pandas

遺伝子名のフォーマットはEnsemblやUniProt等のデータベース毎に異なります。例えばARF5は遺伝子シンボルと呼ばれるフォーマットですが、この遺伝子のEnsemblにおける名称はENSP00000000233、UniProt IDにおける名称はP84085となります。

Bioinformaticsの解析ではソースとなる実験データの遺伝子名の表記に異なるフォーマットが混在していることが多く、しばしば「Ensemble ID ↔︎ UniProt ID」のような遺伝子ID/シンボルの変換をして1つのフォーマットに統一するプロセスが必要です。

本記事では遺伝子ID/シンボルの変換に必要な「ID/シンボルの対応表」をRのbiomaRtで作成する方法を説明します。

1. biomaRtについて

biomaRtはRのライブラリで、

  1. 変換元の遺伝子IDの名称(ensembl_gene_idなど)
  2. 変換後の遺伝子IDの名称(uniprot_gn_symbolなど)
  3. 変換したい対象のIDの一覧(ENSP00000000233など)

を指定することで、①のフォーマットで記載された③のIDを②のフォーマットに変換するテーブルを出力できます。

2. biomaRtのインストール

まずは以下のようにRのコンソールを開き、BiocManagerを介してbiomaRtをインストールします。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("biomaRt")

3. データベースの選択

biomaRtでは最初に変換に利用するデータベースをuseEnsemblで選択します。

本記事の目的の範囲ならば、汎用性が高いgenesを選択しておけば問題ありません。

方法は以下の通りです。

library(biomaRt)
ensembl <- useEnsembl(biomart = "genes")

4. データセットの選択

次にデータセットの選択を行います。biomaRtにおける「データセット」とは「種」に相当するものです。つまり、「データセットの選択」とは「どの種を用いるか選択する」プロセスのことです。

前節で指定したデータベースで解析したい種がどのような名称で用いられているか分からないので、まずは解析したい種に対応するデータセット名をserachDatasetsで検索します。その際、目的の種の学名の一部をpatternの引数に指定すると検索が易しくなります。

例えば、マウスのデータセット名を検索する場合は以下のようになります。

searchDatasets(mart = ensembl, pattern = "mus")
                      dataset                                   description
18     bmusculus_gene_ensembl                Blue whale genes (mBalMus1.v2)
104 mmoschiferus_gene_ensembl Siberian musk deer genes (MosMos_v2_BIUU_UCD)
108    mmusculus_gene_ensembl                          Mouse genes (GRCm39)
158       psimus_gene_ensembl       Greater bamboo lemur genes (Prosim_1.0)
182     smaximus_gene_ensembl                   Turbot genes (ASM1334776v1)
209   umaritimus_gene_ensembl                 Polar bear genes (UrsMar_1.0)
               version
18         mBalMus1.v2
104 MosMos_v2_BIUU_UCD
108             GRCm39
158         Prosim_1.0
182       ASM1334776v1
209         UrsMar_1.0

ヒトのデータセット名を検索する場合は以下のようになります。

searchDatasets(mart = ensembl, pattern = "sapiens")
                 dataset              description    version
80 hsapiens_gene_ensembl Human genes (GRCh38.p14) GRCh38.p14

データセット名が判明したら、useDatasetで以下のようにデータセットを指定します。

ensembl <- useDataset(dataset = "mmusculus_gene_ensembl", mart = ensembl) # マウス
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl) # ヒト

5. IDの名称を検索

続いて、対応付けたいIDの名称を調べます。

前節で選択したデータセットにおいて、変換前後の遺伝子IDがどのような名称で登録されているかsearchFiltersで調べます。

例えば、

  1. Ensembl Gene ID
  2. Ensembl Protein ID
  3. NCBI refseq ID
  4. UniProt ID
  5. 遺伝子シンボル

の対応表を作成をする場合、以下のように進めます。

5.1. Ensembl ID

1と2のEnsembl ID (Ensembl Gene ID (ENSMUSG00000000001など) / Protein ID (ENSP00000000233など))の名称は以下のように調べます。

# マウス
ensembl <- useDataset(dataset = "mmusculus_gene_ensembl", mart = ensembl)
searchFilters(mart = ensembl, pattern = "ensembl.*id")
                            name
43               ensembl_gene_id
44       ensembl_gene_id_version
45         ensembl_transcript_id
46 ensembl_transcript_id_version
47            ensembl_peptide_id
48    ensembl_peptide_id_version
49               ensembl_exon_id
                                                        description
43                      Gene stable ID(s) [e.g. ENSMUSG00000000001]
44       Gene stable ID(s) with version [e.g. ENSMUSG00000000001.5]
45                Transcript stable ID(s) [e.g. ENSMUST00000000001]
46 Transcript stable ID(s) with version [e.g. ENSMUST00000000001.5]
47                   Protein stable ID(s) [e.g. ENSMUSP00000000001]
48    Protein stable ID(s) with version [e.g. ENSMUSP00000000001.5]
49                             Exon ID(s) [e.g. ENSMUSE00000097910]
# ヒト
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
searchFilters(mart = ensembl, pattern = "ensembl.*id")
                            name
53               ensembl_gene_id
54       ensembl_gene_id_version
55         ensembl_transcript_id
56 ensembl_transcript_id_version
57            ensembl_peptide_id
58    ensembl_peptide_id_version
59               ensembl_exon_id
                                                      description
53                       Gene stable ID(s) [e.g. ENSG00000000003]
54       Gene stable ID(s) with version [e.g. ENSG00000000003.16]
55                 Transcript stable ID(s) [e.g. ENST00000000233]
56 Transcript stable ID(s) with version [e.g. ENST00000000233.10]
57                    Protein stable ID(s) [e.g. ENSP00000000233]
58     Protein stable ID(s) with version [e.g. ENSP00000000233.5]
59                              Exon ID(s) [e.g. ENSE00000000001]

descriptionにIDの例が掲載されているので、変換したい対象のIDとパターンがマッチするか目視で確認して名称を同定します。

マウスを例にすると、ENSMUSG00000000001のようなIDと最も近いフォーマットは、descriptionの43行目のensembl_gene_idであるため、Ensembl Gene ID の名称はensembl_gene_idであると分かります。また、ENSP00000000233のようなIDの名称も同様に比較することでEnsembl Protein IDの名称はensembl_peptide_idとなることが分かります。

5.2. NCBI refseq ID

NCBI refseq IDについては以下のように検索できます。

# マウス
ensembl <- useDataset(dataset = "mmusculus_gene_ensembl", mart = ensembl)
searchFilters(mart = ensembl, pattern = "refseq")
                            name
27              with_refseq_mrna
28    with_refseq_mrna_predicted
29             with_refseq_ncrna
30   with_refseq_ncrna_predicted
31           with_refseq_peptide
32 with_refseq_peptide_predicted
76                   refseq_mrna
77         refseq_mrna_predicted
78                  refseq_ncrna
79        refseq_ncrna_predicted
80                refseq_peptide
81      refseq_peptide_predicted
                                          description
27                             With RefSeq mRNA ID(s)
28                   With RefSeq mRNA predicted ID(s)
29                            With RefSeq ncRNA ID(s)
30                  With RefSeq ncRNA predicted ID(s)
31                          With RefSeq peptide ID(s)
32                With RefSeq peptide predicted ID(s)
76              RefSeq mRNA ID(s) [e.g. NM_001001130]
77    RefSeq mRNA predicted ID(s) [e.g. XM_001002281]
78                RefSeq ncRNA ID(s) [e.g. NR_000002]
79      RefSeq ncRNA predicted ID(s) [e.g. XR_001538]
80           RefSeq peptide ID(s) [e.g. NP_001001130]
81 RefSeq peptide predicted ID(s) [e.g. XP_001004117]
# ヒト
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
searchFilters(mart = ensembl, pattern = "refseq")
                             name
36               with_refseq_mrna
37     with_refseq_mrna_predicted
38              with_refseq_ncrna
39    with_refseq_ncrna_predicted
40            with_refseq_peptide
41  with_refseq_peptide_predicted
95                    refseq_mrna
96          refseq_mrna_predicted
97                   refseq_ncrna
98         refseq_ncrna_predicted
99                 refseq_peptide
100      refseq_peptide_predicted
                                           description
36                              With RefSeq mRNA ID(s)
37                    With RefSeq mRNA predicted ID(s)
38                             With RefSeq ncRNA ID(s)
39                   With RefSeq ncRNA predicted ID(s)
40                           With RefSeq peptide ID(s)
41                 With RefSeq peptide predicted ID(s)
95                  RefSeq mRNA ID(s) [e.g. NM_000014]
96     RefSeq mRNA predicted ID(s) [e.g. XM_003403597]
97                 RefSeq ncRNA ID(s) [e.g. NR_000005]
98    RefSeq ncRNA predicted ID(s) [e.g. XR_001736914]
99               RefSeq peptide ID(s) [e.g. NP_000005]
100 RefSeq peptide predicted ID(s) [e.g. XP_003403645]

NCBI refseq IDの名称はrefseq_peptideであることが分かります。

5.3. UniProt ID

UniProt IDは以下のようにして検索します。

# ヒト
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
searchFilters(mart = ensembl, pattern = "uniprot")
                     name                                  description
48        with_uniprot_gn               With UniProtKB Gene Name ID(s)
49   with_uniprot_isoform                 With UniProtKB isoform ID(s)
50  with_uniprotswissprot              With UniProtKB/Swiss-Prot ID(s)
51   with_uniprotsptrembl                  With UniProtKB/TrEMBL ID(s)
107         uniprot_gn_id  UniProtKB Gene Name ID(s) [e.g. A0A023T6R1]
108     uniprot_gn_symbol   UniProtKB Gene Name symbol(s) [e.g. 5HT1A]
109       uniprot_isoform  UniProtKB isoform ID(s) [e.g. A0A096LP49-1]
110      uniprotswissprot UniProtKB/Swiss-Prot ID(s) [e.g. A0A024R1R8]
111       uniprotsptrembl     UniProtKB/TrEMBL ID(s) [e.g. A0A023T6R1]

UniProt IDの名称はuniprot_gn_idであることが分かります。

5.4. 遺伝子シンボル

遺伝子シンボル表記については以下のように検索します。

# マウス
# マウス・ヒト共通
searchFilters(mart = ensembl, pattern = "symbol")
                name                                        description
62       hgnc_symbol                      HGNC symbol(s) [e.g. OR5BS1P]
65        mgi_symbol                 MGI symbol(s) [e.g. 0610005C13Rik]
87 uniprot_gn_symbol UniProtKB Gene Name symbol(s) [e.g. 0610009B22Rik]

uniprot_gn_symbolが名称です。

6. 変換のテスト

一度に全ての対応表を作成する前に、まずは1件について対応表を作成してみます。

対応表の作成はgetBMで実施します。以下のようにattributesに前節で調べたID名称をベクトルで指定します。filtersには変換元のID名称を、valuesには変換対象のIDを指定します。

# ヒト
getBM(
  attributes=c("ensembl_peptide_id", "ensembl_gene_id", "refseq_peptide", "uniprot_gn_id", "uniprot_gn_symbol"),
  filters = "ensembl_peptide_id",
  values = c("ENSP00000000233"),
  mart= ensembl
)
 ensembl_peptide_id ensembl_gene_id refseq_peptide uniprot_gn_id
1    ENSP00000000233 ENSG00000004059      NP_001653        A4D0Z3
2    ENSP00000000233 ENSG00000004059      NP_001653        P84085
3    ENSP00000000233 ENSG00000004059      NP_001653        C9J1Z8
  uniprot_gn_symbol
1              ARF5
2              ARF5
3              ARF5

指定した遺伝子ID/シンボルの対応が得られました。

7. IDの変換

予め変化したい遺伝子リストを作成しておき、例えばnamelist.txtのようにして保存しておきます。

namelist.txtの内容は、以下のような1行目がヘッダで2行目以降に遺伝子IDリストが掲載されている形式にしておきます。

ensembl_id
ENSP00000000233
ENSP00000000234
ENSP00000000235
...

次に以下のようにしてnamelist.txtを読み込んで、遺伝子IDの部分だけ選択します。

d <- read.csv("namelist.txt")
v <- d[, 1]

続いて、以下のようにして変換を実行します。

res <- getBM(
 	attributes=c("ensembl_peptide_id", "ensembl_gene_id", "refseq_peptide", "uniprot_gn_id", "uniprot_gn_symbol"),
 	filters = "ensembl_peptide_id",
  	values = v,
  	mart= ensembl
)

最後に以下のようにして保存をします。

write.csv(res, "id_conversion/ensemblP_to_symbols.csv", quote=FALSE, row.names=FALSE)

8. ここまでのまとめ

上記を1つにまとめたスクリプトは以下です。

8.1. マウス

# install
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("biomaRt")

# select database and datasets
ensembl <- useEnsembl(biomart = "genes")
ensembl <- useDataset(dataset = "mmusculus_gene_ensembl", mart = ensembl)

# id conversion
d <- read.csv("namelist.txt")
v <- d[, 1]
res <- getBM(
 	attributes=c("ensembl_peptide_id", "ensembl_gene_id", "refseq_peptide", "uniprot_gn_id", "uniprot_gn_symbol"),
 	filters = "ensembl_peptide_id",
  	values = v,
  	mart= ensembl
)

# save
write.csv(res, "id_conversion/ensemblP_to_symbols.csv", quote=FALSE, row.names=FALSE)

8.2. ヒト

# install
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("biomaRt")

# select database and datasets
ensembl <- useEnsembl(biomart = "genes")
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)

# id conversion
d <- read.csv("namelist.txt")
v <- d[, 1]
res <- getBM(
 	attributes=c("ensembl_peptide_id", "ensembl_gene_id", "refseq_peptide", "uniprot_gn_id", "uniprot_gn_symbol"),
 	filters = "ensembl_peptide_id",
  	values = v,
  	mart= ensembl
)

# save
write.csv(res, "id_conversion/ensemblP_to_symbols.csv", quote=FALSE, row.names=FALSE)

9. Pythonから結果を読み込む

以下のようにしてRのbiomaRtで作成した遺伝子ID/シンボルの対応表をPythonから読み込んでみます。

import pandas as pd
df_eid_sym = pd.read_csv('id_conversion/ensemblP_to_symbols.csv', index_col=0)
df_eid_sym.head()
ensembl_peptide_iduniprot_gn_iduniprot_gn_symbol
ENSP00000354687U5Z754ND1
ENSP00000354687P03886MT-ND1
ENSP00000354499U5YWV7COX1
ENSP00000354499P00395MT-CO1
ENSP00000354665U5Z977
以下のようにして辞書にすれば、Ensembl Protein IDから遺伝子シンボルを取得することができるようになります。
eid_to_sym = df_eid_sym['uniprot_gn_symbol'].to_dict()

10. まとめ

本記事では、RのbiomaRtを用いて遺伝子ID/シンボルの変換テーブルを作成する手順を示しました。

以上になります。

スポンサーリンク

コメント

ホーム
トップ