[R] 遺伝子ID/遺伝子名の変換方法
目次
AnnotationDataパッケージ
各生物種の遺伝子アノテーション情報をまとめたデータベースがRのパッケージとして提供されている。
ここから遺伝子アノテーションに加え、パスウェイデータベースなどのパッケージも探すことができる。
Genome wide annotationパッケージ
上記リンクからorg.生物種.eg.db
というパターンのパッケージを見つけてインストールしておく。ヒトの場合はorg.Hs.eg.db、マウスの場合はorg.Mm.eg.dbを使用する。
(egはEntrez geneの略??)
library(org.Hs.eg.db)
一度パッケージを読み込むと、パッケージ名のオブジェクトが参照できるようになる。
試しにオブジェクト名を打ってみると、データベース内の各情報ソースの名前やversionなどが確認できる。
org.Hs.eg.db
データベースオブジェクトの中身
OrgDb object:
| DBSCHEMAVERSION: 2.1
| Db type: OrgDb
| Supporting package: AnnotationDbi
| DBSCHEMA: HUMAN_DB
| ORGANISM: Homo sapiens
| SPECIES: Human
| EGSOURCEDATE: 2022-Mar17
| EGSOURCENAME: Entrez Gene
| EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
| CENTRALID: EG
| TAXID: 9606
| GOSOURCENAME: Gene Ontology
| GOSOURCEURL: http://current.geneontology.org/ontology/go-basic.obo
| GOSOURCEDATE: 2022-03-10
| GOEGSOURCEDATE: 2022-Mar17
| GOEGSOURCENAME: Entrez Gene
| GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
| KEGGSOURCENAME: KEGG GENOME
| KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
| KEGGSOURCEDATE: 2011-Mar15
| GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
| GPSOURCEURL:
| GPSOURCEDATE: 2022-Nov23
| ENSOURCEDATE: 2021-Dec21
| ENSOURCENAME: Ensembl
| ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
| UPSOURCENAME: Uniprot
| UPSOURCEURL: http://www.UniProt.org/
| UPSOURCEDATE: Fri Apr 1 14:42:16 2022
アノテーションパッケージを読み込むと、合わせてAnnotationDbiパッケージが読み込まれる。
データベースからのデータ取り出しにはAnnotationDbiに用意されてた関数を使用することになる。
項目確認
データベースに存在する項目名を確認するにはkeytypes機能かcolumns機能を使用する。
keytypes(org.Hs.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" "ENZYME"
[8] "EVIDENCE" "EVIDENCEALL" "GENENAME" "GENETYPE" "GO" "GOALL" "IPI"
[15] "MAP" "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
[22] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" "UNIPROT"
ENSEMBLデータベースでの遺伝子ID, EntreZ ID, 遺伝子名 (SYMBOL)だけでなく、Gene ontologyなどの情報もある。
各項目の説明はここにあった。
項目の説明
Keytype | Description |
---|---|
ACCNUM | GenBank accession numbers |
ALIAS | Commonly used gene symbols |
ARACYC | KEGG Identifiers for arabidopsis as indicated by aracyc |
ARACYCENZYME | Aracyc enzyme names as indicated by aracyc |
CHR | Chromosome (deprecated for Bioc > 3.1) For this information you should look at a TxDb or OrganismDb object and search for an appropriate field like TXCHROM, EXONCHROM or CDSCHROM. This information can also be retrieved from these objects using an appropriate range based accesor like transcripts, transcriptsBy etc. |
CHRLOC | Chromosome and starting base of associated gene (deprecated for Bioc > 3.1) For this information you should look at a TxDb or OrganismDb object and search for an appropriate field like TXSTART etc. or even better use the associated range based accessors like transcripts or transcriptsBy to get back GRanges objects. |
CHRLOCEND | Chromosome and ending base of associated gene (deprecated for Bioc > 3.1) For this information you should look at a TxDb or OrganismDb object and search for an appropriate field like TXEND etc. or even better use the associated range based accessors like transcripts or transcriptsBy to get back GRanges objects. |
COMMON | Common name |
DESCRIPTION | The description of the associated gene |
ENSEMBL | The ensembl ID as indicated by ensembl |
ENSEMBLPROT | The ensembl protein ID as indicated by ensembl |
ENSEMBLTRANS | The ensembl transcript ID as indicated by ensembl |
ENTREZID | Entrez gene Identifiers |
ENZYME | Enzyme Commission numbers |
EVIDENCE | Evidence codes for GO associations with a gene of interest |
EVIDENCEALL | Evidence codes for GO (includes less specific terms) |
GENENAME | The full gene name |
GO | GO Identifiers associated with a gene of interest |
GOALL | GO Identifiers (includes less specific terms) |
INTERPRO | InterPro identifiers |
IPI | IPI accession numbers |
MAP | cytoband locations |
OMIM | Online Mendelian Inheritance in Man identifiers |
ONTOLOGY | For GO Identifiers, which Gene Ontology (BP, CC, or MF) |
ONTOLOGYALL | Which Gene Ontology (BP, CC, or MF), (includes less specific terms) |
ORF | Yeast ORF Identifiers |
PATH | KEGG Pathway Identifiers |
PFAM | PFAM Identifiers |
PMID | Pubmed Identifiers |
PROBEID | Probe or manufacturer Identifiers for a chip package |
PROSITE | Prosite Identifiers |
REFSEQ | Refseq Identifiers |
SGD | Saccharomyces Genome Database Identifiers |
SMART | Smart Identifiers |
SYMBOL | The official gene symbol |
TAIR | TAIR Identifiers |
UNIGENE | Unigene Identifiers |
UNIPROT | Uniprot Identifiers |
データ取り出し
各項目の中身を確認するにはkeys機能を使用する。
keys(データベースオブジェクト, keytype = "項目名")
試しにorg.Hs.eg.dbの全項目をチラ見してみる。
for(i in columns(org.Hs.eg.db)){
print(i)
str(keys(org.Hs.eg.db,i))
}
各項目の中身
[1] "ACCNUM"
chr [1:806658] "AA484435" "AAH35719" "AAL07469" "AB073611" "ACJ13639" "AF414429" "AI022193" "AK055885" ...
[1] "ALIAS"
chr [1:131663] "A1B" "ABG" "GAB" "HYST2477" "A1BG" "A2MD" "CPAMD5" "FWP007" "S863-7" "A2M" "A2MP" "A2MP1" ...
[1] "ENSEMBL"
chr [1:38636] "ENSG00000121410" "ENSG00000175899" "ENSG00000256069" "ENSG00000171428" "ENSG00000156006" ...
[1] "ENSEMBLPROT"
chr [1:21427] "ENSP00000286479" "ENSP00000428416" "ENSP00000483300" "ENSP00000482269" "ENSP00000478547" ...
[1] "ENSEMBLTRANS"
chr [1:38893] "ENST00000543404" "ENST00000566278" "ENST00000545343" "ENST00000544183" "ENST00000286479" ...
[1] "ENTREZID"
chr [1:66102] "1" "2" "3" "9" "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" ...
[1] "ENZYME"
chr [1:975] "2.3.1.5" "2.3.1.87" "6.1.1.7" "2.6.1.22" "2.6.1.19" "2.7.10.2" "1.4.3.22" "2.4.1.37" "2.4.1.40" ...
[1] "EVIDENCE"
chr [1:20] "ND" "HDA" "IDA" "TAS" "IEA" "IBA" "IPI" "IMP" "NAS" "ISS" "IEP" "IGI" "IC" "EXP" "HEP" "HMP" "RCA" ...
[1] "EVIDENCEALL"
chr [1:20] "ND" "HDA" "IDA" "TAS" "IEA" "IBA" "IPI" "IMP" "NAS" "ISS" "IEP" "IC" "IGI" "EXP" "HMP" "HEP" "RCA" ...
[1] "GENENAME"
chr [1:65232] "alpha-1-B glycoprotein" "alpha-2-macroglobulin" "alpha-2-macroglobulin pseudogene 1" ...
[1] "GENETYPE"
chr [1:11] "protein-coding" "pseudo" "other" "unknown" "ncRNA" "tRNA" "rRNA" "scRNA" "snoRNA" "snRNA" ...
[1] "GO"
chr [1:18644] "GO:0003674" "GO:0005576" "GO:0005615" "GO:0008150" "GO:0031093" "GO:0034774" "GO:0062023" ...
[1] "GOALL"
chr [1:22834] "GO:0003674" "GO:0005575" "GO:0005576" "GO:0005615" "GO:0005622" "GO:0005737" "GO:0008150" ...
[1] "IPI"
chr [1:37945] "IPI00644018" "IPI00022895" "IPI00478003" "IPI00644361" "IPI00930394" "IPI00004421" ...
[1] "MAP"
chr [1:2003] "19q13.43" "12p13.31" "8p22" "14q32.13" "3q25.1" "2q35" "17q25.1" "16q22.1" "19q13" "19q13-qter" ...
[1] "OMIM"
chr [1:22739] "138670" "103950" "108345" "243400" "612182" "107280" "600338" "603488" "600950" "601065" ...
[1] "ONTOLOGY"
chr [1:3] "MF" "CC" "BP"
[1] "ONTOLOGYALL"
chr [1:3] "MF" "CC" "BP"
[1] "PATH"
chr [1:229] "04610" "00232" "00983" "01100" "00380" "00970" "00250" "00280" "00410" "00640" "00650" "02010" ...
[1] "PFAM"
chr [1:6280] "PF13895" "PF07678" "PF07677" "PF17791" "PF00207" "PF01835" "PF07703" "PF17789" "PF00797" ...
[1] "PMID"
chr [1:737868] "2591067" "3458201" "3610142" "8889549" "12477932" "14702039" "14760718" "15221005" "15461460" ...
[1] "PROSITE"
chr [1:1787] "PS50835" "PS00477" "PS00284" "PS01174" "PS50082" "PS00678" "PS50294" "PS51186" "PS50860" ...
[1] "REFSEQ"
chr [1:285470] "NM_130786" "NP_570602" "NM_000014" "NM_001347423" "NM_001347424" "NM_001347425" "NP_000005" ...
[1] "SYMBOL"
chr [1:66091] "A1BG" "A2M" "A2MP1" "NAT1" "NAT2" "NATP" "SERPINA3" "AADAC" "AAMP" "AANAT" "AARS1" "AAVS1" ...
[1] "UCSCKG"
chr [1:222852] "ENST00000596924.1" "ENST00000263100.8" "ENST00000600123.5" "ENST00000595014.1" ...
[1] "UNIPROT"
chr [1:32129] "P04217" "V9HWD8" "P01023" "P18440" "Q400J6" "F5H5R8" "A4Z6T7" "P11245" "A0A024R6P0" "P01011" ...
項目内の登録数が揃っていないのが分かるだろう。
ALIAS
もSYMBOL
も遺伝子名の情報であるが、実際にkeysで取り出して違いを見てみる。
alias <- keys(org.Hs.eg.db, "ALIAS")
symbol <- keys(org.Hs.eg.db, "SYMBOL")
登録されている遺伝子数が全然違う。
length(alias)
length(symbol)
SYMBOLの遺伝子名はALIASにも登録されている。
length(intersect(alias, symbol))
setdiff(symbol, alias)
遺伝子変換
複数のデータベースが個別にIDをつけているため1つの遺伝子に複数の表記が存在する。
上記のうち、ENSEMBL
, ENTREZ
, SYMBOL
辺りをよく使うのではないだろうか。
遺伝子名 (SYMBOL)で解析を進めたいところであるが、解析パッケージによって求められる遺伝子表記が異なり、その度に遺伝子名 <-> 遺伝子ID 変換が求められる。
select
データベースのある項目から他の項目の情報を取得するにはselect機能を使用する。selectだけだと他のパッケージのselectと機能名が重複していることがあるので、AnnotationDbi::select
と書くと間違い無い。
具体的には次のように記載する。
select(x = データベースオブジェクト
, keys = 変換したい遺伝子のベクトル
, keytype = keysに指定したベクトルの表記法
, colums = 変換先の項目
)
試しにENSEMBL IDをSYMBOLとENTREZIDに変換してみる。
keys <- head(keys(org.Hs.eg.db, keytype = "ENSEMBL"))
select(x = org.Hs.eg.db,
keys = keys,
keytype = "ENSEMBL",
columns = c("SYMBOL","ENTREZID"))
GENENAMEやGENETYPEの情報も遺伝子を理解するのに有用である。
select(x = org.Hs.eg.db,
keys = keys,
keytype = "ENSEMBL",
columns = c("SYMBOL","GENENAME","GENETYPE"))
select機能の注意点として、一つのkeyに対して複数の情報がある場合である。
以下の例では1つのENSEMBL IDにmRNAのRefseq IDとproteinのRefseq IDがあり、指定したkeysの数よりも大きなdata.frameが返されてしまう。
select(x = org.Hs.eg.db,
keys = keys,
keytype = "ENSEMBL",
columns = c("REFSEQ"))
mapIds
似た機能にmapIdsがあるが、こちらは項目を一つしか指定できない。
デフォルトの返り値はkeysが名前属性に入ったベクトルである。
mapIds(x = org.Hs.eg.db, keys = keys, keytype = "ENSEMBL", column = "SYMBOL")
しかし、mapIdsにはmultiVals
オプションがあり、1つのkeyに対して複数のhitがある場合の処理が選べる。
デフォルトではmultiVals="first"
で、最初のhitが返されている。
multiVals="list"
にすると1つのkeyに対する複数のhitをlistにまとめてくれる。
mapIds(x = org.Hs.eg.db, keys = keys, keytype = "ENSEMBL", column = c("REFSEQ"),
multiVals = "list")
CharacterListというIRangesパッケージの専用クラスで返すこともできる。
mapIds(x = org.Hs.eg.db, keys = keys, keytype = "ENSEMBL", column = c("REFSEQ"),
multiVals = "CharacterList")
clusterProfilerのbitr機能
clusterProfilerパッケージのbitr機能でもID変換が可能であるが、これも内部でAnnotationDbiを使用している。
引数名が違うぐらいでselect機能と殆ど使い方が変わらない。
library(clusterProfiler)
genes <- head(keys(x = org.Hs.eg.db, keytype = "SYMBOL"))
bitr(geneID = genes,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
ただデフォルトで変換が旨くいかなかった場合はその遺伝子を除くようになっている。
指定遺伝子を減らしたくない場合や、変換がうまくいかなかった遺伝子を知りたい場合は、drop=FALSE
オプションをつけておくよい。
bitr(geneID = genes,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db,
drop=FALSE)
Ensdbパッケージ
Ensemblデータベースを基にした遺伝子データベース。
ヒトの物でも複数のversionのパッケージが提供されているが、genome versionやensembl versionの違いにより登録されている遺伝子数、転写産物数が違うようである。
# BiocManager::install("EnsDb.Hsapiens.v75")
library(EnsDb.Hsapiens.v75)
EnsDb.Hsapiens.v75
# BiocManager::install("EnsDb.Hsapiens.v86")
library(EnsDb.Hsapiens.v86)
EnsDb.Hsapiens.v86
項目確認
AnnotationDbiのkeytypes機能で項目名が確認できる。
keytypes(EnsDb.Hsapiens.v86)
[1] "ENTREZID" "EXONID" "GENEBIOTYPE" "GENEID" "GENENAME"
[6] "PROTDOMID" "PROTEINDOMAINID" "PROTEINDOMAINSOURCE" "PROTEINID" "SEQNAME"
[11] "SEQSTRAND" "SYMBOL" "TXBIOTYPE" "TXID" "TXNAME"
[16] "UNIPROTID"
注意点は、ENSEMBLという項目名が無く、このデータベースではGENEIDがEnsembl IDである。
遺伝子変換
AnnotationDbiのselectが上述と同様に使用可能である。
AnnotationDbi::select(x = EnsDb.Hsapiens.v86,
keys = keys,
keytype = "GENEID",
columns = "ENTREZID")
bioMart
おそらくこれが一番、多量の情報を取得できる。その分、欲しいデータまでたどり着くのが大変。ここではあくまで遺伝子名 <-> 遺伝子ID変換 に必要な情報のみ扱う。
bioMartはウェブ上からも同様のデータを取得することができる。上部タブのBioMart
から移行できる。
これからfiltersやattributesといった概念で情報の選別をしていくのだが、コマンドラインだと情報が膨大でわかりづらいのでブラウザーで感覚をつかんでおくと良いだろう。
togotvでもブラウザ操作が紹介されている。
# BiocManager::install("biomaRt")
library(biomaRt)
使用するMartの選択
listMarts()
db <- useMart(biomart = "ENSEMBL_MART_ENSEMBL") # useMart("ensembl")でも良いみたい。
db
Object of class 'Mart':
Using the ENSEMBL_MART_ENSEMBL BioMart database
No dataset selected.
このようにしてミラーサイトを指定することも可能。server errorが出る際にはこちらの方がいいかも。
db <- useEnsembl(biomart = "ENSEMBL_MART_ENSEMBL", mirror = "asia", dataset = "hsapiens_gene_ensembl")
データベースの選択
使用できるデータベース一覧を表示
head(listDatasets(db))
213種の生物種のデータベースが登録されている。
ヒトのデータベースを取得するには以下を実行する。
hg <- useDataset("hsapiens_gene_ensembl", mart = db)
項目確認
bioMartでは、filtersで欲しい情報の制限 (特定の染色体番号や、queryとする遺伝子一覧、、、など) し、attributesで欲しい情報 (遺伝子名や遺伝子ID、、、など)を選択する。
listFilters機能、listAttributes機能で指定できるtermを確認できるが、膨大な条件が選択可能であることを見てみると良いだろう。(Rのコマンドの出力から探すよりブラウザで見る方が現実的である。)
ただブラウザでの表記名とRで実際に指定するときの指定名は異なっている。listFiltersやlistAttributesの出力のdescription列がブラウザでの表記名なので、これを基に探すとよいだろう。
RでのlistFiltersの出力
ブラウザでのFilters画面
Attributesの方は、ホモログの情報も含んでおり、他の生物種の情報も入れると膨大な項目が得られる。
page = "feature_page"にするとまだ見やすいだろう。
listAttributes(mart = hg, page = "feature_page")
情報取得
ENSEMBL IDをkeyにして、遺伝子名を取得してみる。
keys <- head(keys(org.Hs.eg.db, keytype = "ENSEMBL"))
keys
[1] "ENSG00000121410" "ENSG00000175899" "ENSG00000256069" "ENSG00000171428" "ENSG00000156006" "ENSG00000196136"
入力する遺伝子の表記はEnsembl ID。上記で確認した通りRでの指定名は"ensembl_gene_id"になる。
attributesには遺伝子名っぽい項目が幾つかあるので、ひとまず全て使用してみる。
実際の情報取得にはgetBM機能を使用する。attributes
に欲しい情報、 mart
に使用するデータベース、filters
にvaluesに指定する遺伝子表記法, values
にqueryとなる遺伝子ベクトルを指定した。
getBM(attributes = c("ensembl_gene_id","hgnc_symbol","external_gene_name","wikigene_name"),
mart = hg,
filters = "ensembl_gene_id",
values = keys)
data.frameの出力が得られるので、適宜成型して使用できる。
ちなみに、attributesとmart引数だけを使用すると、全遺伝子の指定attributesが取得可能である。
ホモログ変換
biomaRtパッケージのgetLDS機能を使って2つのdataset間でhomology mappingを行うことができる。
ヒト遺伝子名をマウス遺伝子名に変換
ヒト、マウスそれぞれのMartオブジェクトを作成。
human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = c("hgnc_symbol"),
filters = "hgnc_symbol",
values = keys,
mart = human,
attributesL = c("mgi_symbol"),
martL = mouse,
uniqueRows=T)
(記事作成時にはserver errorが出たため出力例は無し。)
Discussion