2295の生物種のアノテーションを扱うAnnotationHubパッケージ

バイオインフォのツールは、ヒト、マウスなど、メジャーな生物種を中心に開発される事が多い。

Bioconductorも例外ではなく、生物種ごとの遺伝子アノテーション情報を集めているOrgDb（生物種DBの一つ）は、19しか用意されておらず（http://bioconductor.org/packages/release/BiocViews.html#___OrgDb ）、ここに記載されていない生物種の解析をする研究者は苦労することになる。

例えば、GOstatsパッケージが使えなくなるなど（https://qiita.com/antiplastics/items/add3f8438cbe6b2b594e ）。

自分は、MeSH関連のパッケージで、できるだけ幅広い生物種に対応させようと頑張っていたので、74生物種はMeSHエンリッチメント解析ができるが（http://bioconductor.org/packages/release/BiocViews.html#___MeSHDb ）、それでもあらゆる生物種に対応しているわけではない。

ところで、最近、AnnotationHubというパッケージが相当幅広い生物種の情報を公開していることを知ったので紹介する。

準備

まずは、AnnotationHubパッケージをダウンロード&インストールする。

# パッケージインストール
BiocManager::install("AnnotationHub")
# ロード
library("AnnotationHub")

基本の使い方

まずAnnotationHubというそのままの関数があるので、これで設定を行う。

# まずAnnotationHubを構築
ah <- AnnotationHub()

ah
# AnnotationHub with 46429 records
# snapshotDate(): 2019-05-02
# $dataprovider: BroadInstitute, Ensembl, UCSC, ftp://ftp.ncbi.nlm.nih.gov/g...
# $species: Homo sapiens, Mus musculus, Drosophila melanogaster, Bos taurus,...
# $rdataclass: GRanges, BigWigFile, TwoBitFile, OrgDb, Rle, EnsDb, ChainFile...
# additional mcols(): taxonomyid, genome, description,
#   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#   rdatapath, sourceurl, sourcetype
# retrieve records with, e.g., 'object[["AH5012"]]'

#            title
#  AH5012  | Chromosome Band
#  AH5013  | STS Markers
#  AH5014  | FISH Clones
#  AH5015  | Recomb Rate
#  AH5016  | ENCODE Pilot
#  ...       ...
#  AH73982 | Ensembl 97 EnsDb for Xiphophorus couchianus
#  AH73983 | Ensembl 97 EnsDb for Xiphophorus maculatus
#  AH73984 | Ensembl 97 EnsDb for Xenopus tropicalis
#  AH73985 | Ensembl 97 EnsDb for Zonotrichia albicollis
#  AH73986 | Ensembl 79 EnsDb for Homo sapiens

46429ものデータベースが公開されているのがわかる。
個々のデータベースのメタ情報は、以下のように、ベクトルのようにアクセスできる。

ah[1]
# AnnotationHub with 1 record
# snapshotDate(): 2019-05-02
# names(): AH5012
# $dataprovider: UCSC
# $species: Homo sapiens
# $rdataclass: GRanges
# $rdatadateadded: 2013-03-26
# $title: Chromosome Band
# $description: GRanges object from UCSC track 'Chromosome Band'
# $taxonomyid: 9606
# $genome: hg19
# $sourcetype: UCSC track
# $sourceurl: rtracklayer://hgdownload.cse.ucsc.edu/goldenpath/hg19/database...
# $sourcesize: NA
# $tags: c("cytoBand", "UCSC", "track", "Gene", "Transcript",
#   "Annotation")
# retrieve record with 'object[["AH5012"]]'

ah[2]
# AnnotationHub with 1 record
# snapshotDate(): 2019-05-02
# names(): AH5013
# $dataprovider: UCSC
# $species: Homo sapiens
# $rdataclass: GRanges
# $rdatadateadded: 2013-03-26
# $title: STS Markers
# $description: GRanges object from UCSC track 'STS Markers'
# $taxonomyid: 9606
# $genome: hg19
# $sourcetype: UCSC track
# $sourceurl: rtracklayer://hgdownload.cse.ucsc.edu/goldenpath/hg19/database...
# $sourcesize: NA
# $tags: c("stsMap", "UCSC", "track", "Gene", "Transcript",
#   "Annotation")
# retrieve record with 'object[["AH5013"]]'

このAH*というのは、個々のデータベースのIDらしい。

mcolsという関数を使うと、このメタ情報を、データフレーム状に閲覧することができる。

# mcolsでデータフレーム状に一覧
mcols(ah)
# DataFrame with 46429 rows and 15 columns
#                                               title dataprovider
#                                         <character>  <character>
# AH5012                              Chromosome Band         UCSC
# AH5013                                  STS Markers         UCSC
# AH5014                                  FISH Clones         UCSC
# ...省略...

# AH*は行名になっているらしい
head(rownames(mcols(ah)))

# 15種類のメタ情報が用意されている
dim(mcols(ah))
# 46429    15

# 列名
colnames(mcols(ah))
# [1] "title"              "dataprovider"       "species"
# [4] "taxonomyid"         "genome"             "description"
# [7] "coordinate_1_based" "maintainer"         "rdatadateadded"
# [10] "preparerclass"      "tags"               "rdataclass"
# [13] "rdatapath"          "sourceurl"          "sourcetype"

そして、対応生物種が異常に多く、2295も用意してある！
菌は株ごとに用意しているレベル。

# 対応生物種の多さ（2295種）
sort(unique(ah$species))

データベースも様々なものが用意されているようである。

# 公開データのタイプの多さ（24種類）
unique(ah$rdataclass)

キャッシュについて

先ほどまでの操作は、あくまでデータベースのメタ情報を見ていただけで、実際のデータベースを手元にダウンロードしているわけではない。

以下のように指定すると、実際にダウンロードを始める。

ah[["AH5012"]]

そもそも46429ものデータベースを勝手にダウンロードし始めたら大変なことになる。

AnnotationHubの良くできている点は、キャッシュの概念をうまく取り入れた点である。

つまり、一度ダウンロードしたデータベースはキャッシュとして、保存する方針にしている（この仕組み自体は、BiocFileCacheパッケージを利用）。

これにより、ユーザーが必要になった時に初めて、ダウンロードを開始し、一度ダウンロードした同じファイルは、再度ダウンロードしなくて済む。

データベースファイル自体（sqlite3ファイルなど）は、~/Library/Caches/AnnotationHub以下に保存されるらしい。

# 個々のデータ自体は、ダウンロードされるたびに
# ~/Library/Caches/AnnotationHub以下にキャッシュされる
hubCache(ah)
system("ls ~/Library/Caches/AnnotationHub")

# mcolでは、~/.AnnotationHub以下に置かれたキャッシュを見ている
#（まだデータ本体はダウンロードしていない）
system("ls ~/.AnnotationHub")

# 2019-05-02にできたAnnotationHubらしい（古いVerを利用することも可能）
snapshotDate(ah)

非モデルの生物種データベースを利用

生物種データベースの中でも、学名がそのままの名前になっているタイプのパッケージ（OrganismDb）は、select関数でシンプルにデータにアクセスできて、使いやすいため重宝していたのだが、この類いのパッケージは、陽には、Homo.sapiens, Mus.musculus, Rattus.norvegicusしか作っていないと自分は思っていた。

しかし、AnnotationHubベースで本体のsqlite3は別のサーバーに置いて、必要な時だけ、取りに行くという方針で、実は相当な量を公開していることを最近知った。
例えば、以下のコードは、biomaRtでトラブっていたEnsemblが返すNCBI Gene IDが一意じゃなかった件についてのコードを、Homo.sapiensパッケージで実行したものだが、

library("Homo.sapiens")
select(Homo.sapiens, columns=c("GENEID", "SYMBOL"),
    keytype="GENEID", keys=c("6013", "6019"))
# 'select()' returned 1:1 mapping between keys and columns
#   GENEID SYMBOL
# 1   6013   RLN1
# 2   6019   RLN2

全く同じ操作を、実はAnnotationHubでも実行できる。

query関数というものを使うと、実際のデータベースをダウンロードして、キャッシュとして蓄えてくれる。

# query関数
hs <- query(ah, c("OrgDb", "Homo sapiens"))[[1]]

hs
# OrgDb object:
# | DBSCHEMAVERSION: 2.1
# | Db type: OrgDb
# | Supporting package: AnnotationDbi
# | DBSCHEMA: HUMAN_DB
# | ORGANISM: Homo sapiens
# | SPECIES: Human
# | EGSOURCEDATE: 2019-Apr26
# | EGSOURCENAME: Entrez Gene
# | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
# | CENTRALID: EG
# | TAXID: 9606
# | GOSOURCENAME: Gene Ontology
# | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
# | GOSOURCEDATE: 2019-Apr24
# | GOEGSOURCEDATE: 2019-Apr26
# | 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: 2018-Dec3
# | ENSOURCEDATE: 2019-Apr08
# | ENSOURCENAME: Ensembl
# | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
# | UPSOURCENAME: Uniprot
# | UPSOURCEURL: http://www.UniProt.org/
# | UPSOURCEDATE: Fri Apr 26 20:12:58 2019

select(hs, columns=c("ENTREZID", "SYMBOL"),
    keytype="ENTREZID", keys=c("6013", "6019"))
# 'select()' returned 1:1 mapping between keys and columns
#   ENTREZID SYMBOL
# 1     6013   RLN1
# 2     6019   RLN2

GENEIDという名前ではなく、ENTREZIDという名前でNCBI Gene IDを指定する必要はあるが、操作感はほとんど同じである。

なお、query関数の最後の[[1]]の部分が不思議な感じがするが、query関数の第二引数が、キーワード検索になっている？らしく、必ずしも、1つが検索でヒットするわけではない場合に、[[2]]みたいに、検索でヒットした別のデータベースをダウンロードしたりもできるらしい。

[[1]]の部分を書かない場合は、データをとりにはいかず、メタ情報だけを返す。

AHのIDで指定した方が、一意で安全かもしれない。
RのVersionが変わると、IDが変わってしまうらしい。
例えば、ヒトのIDはR-3.5.0では、AH66156だが、R-3.6.0では、AH70572になっていた。

以下とやっていることは多分同じ。

hs <- ah[["AH70572"]]

Homo sapiensのところを、別の生物種の名前に入れ替えれば、これは実質、あらゆる生物種でのOrganismDbを作ったのと同じである。

mm <- query(ah, c("OrgDb", "Mus musculus"))[[1]]
rn <- query(ah, c("OrgDb", "Rattus norvegicus"))[[1]]
at <- query(ah, c("OrgDb", "Arabidopsis thaliana"))[[1]]
bt <- query(ah, c("OrgDb", "Bos taurus"))[[1]]
ce <- query(ah, c("OrgDb", "Caenorhabditis elegans"))[[1]]
dm <- query(ah, c("OrgDb", "Drosophila melanogaster"))[[1]]
dr <- query(ah, c("OrgDb", "Danio rerio"))[[1]]
gg <- query(ah, c("OrgDb", "Gallus gallus"))[[1]]
pa <- query(ah, c("OrgDb", "Pongo abelii"))[[1]]
xs <- query(ah, c("OrgDb", "Xenopus", "Silurana"))[[1]]
ss <- query(ah, c("OrgDb", "Sus scrofa"))[[1]]

これは、非モデル生物種の研究者の解析を、相当捗らせるのではないだろうか。

自分のscTensorパッケージも、biomaRtパッケージ関係のトラブルが最近絶えないので、AnnotationHubベースに切り替えようか検討している。