🤖

[R] 遺伝子ID/遺伝子名の変換方法

2023/03/13に公開・約16,000字

biomart

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にまとめてくれる。

multiVals = "list"

mapIds(x = org.Hs.eg.db, keys = keys, keytype = "ENSEMBL", column = c("REFSEQ"),
	multiVals = "list")

CharacterListというIRangesパッケージの専用クラスで返すこともできる。

multiVals = "CharacterList"

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

ログインするとコメントできます

目次

AnnotationDataパッケージ

Genome wide annotationパッケージ

項目確認

データ取り出し

遺伝子変換

select

mapIds

clusterProfilerのbitr機能

Ensdbパッケージ

項目確認

遺伝子変換

bioMart

使用するMartの選択

データベースの選択

項目確認

情報取得

ホモログ変換

ヒト遺伝子名をマウス遺伝子名に変換

Discussion