Skip to content

Instantly share code, notes, and snippets.

@kokitsuyuzaki
Last active December 1, 2015 13:58
Show Gist options
  • Save kokitsuyuzaki/2a5694446e3e692fd5af to your computer and use it in GitHub Desktop.
Save kokitsuyuzaki/2a5694446e3e692fd5af to your computer and use it in GitHub Desktop.

MeSH ORA Frameworkの使い方

このgistの内容は、Bioconductorのmeshrパッケージのvignette(パッケージの使用方法が記されたドキュメント)を和訳したものです。

1. イントロ

このgistでは以下のMeSHに関連したパッケージの使い方を説明します。

  • MeSH.db : MeSHの情報を提供するパッケージ
  • MeSH.AOR.db : MeSHの祖先-子孫関係の情報を提供するパッケージ
  • MeSH.PCR.db : MeSHの親子関係の情報を提供するパッケージ
  • MeSH.XXX.eg.db型packages(旧org.MeSH.XXX.db) : 各生物種ごとのNCBI Gene ID - MeSH IDの対応関係を提供するパッケージ
  • MeSHDbi : MeSH関連のアノテーションパッケージのクラスの定義と、ユーザーが独自にMeSH関連のアノテーションパッケージを作成することを支援するパッケージ
  • meshr : MeSHを用いてORA(Over-representation analysis, enrichment analysis)を実行するパッケージ

MeSH(Medical Subject Headings)とはNLM(アメリカ国立医学区図書館)が管理している、MEDLINE/PubMed[1]に注釈情報を与える、網羅的なライフサイエンスの語彙です。 MeSHは25000もの臨床と生物学的な語彙を有しています。

図1

語彙数は、GO(Gene Ontology)[2]の二倍以上あり、カテゴリーもより幅広いです。 2015年現在で、MeSHは19ものカテゴリーを持っており、そのうちの16カテゴリーが実際にMeSH Termが割り当てられています。

各々のカテゴリーはNLMが定めた、1文字の大文字のアルファベットで略されます。

Abbreviation Category
A Anatomy
B Organisms
C Diseases
D Chemicals and Drugs
E Analytical, Diagnostic and Therapeutic Techniques and Equipment
F Psychiatry and Psychology
G Phenomena and Processes
H Disciplines and Occupations
I Anthropology, Education, Sociology and Social Phenomena
J Technology and Food and Beverages
K Humanities
L Information Science
M Persons
N Health Care
V Publication Type
Z Geographical Locations

MeSHはより詳細で網羅的な遺伝子アノテーションツールとして利用できると期待されます。 実際に、MeSHを使ったソフトウェアやデータベースもいくつか実装されています[3], [4], [5], [6]。

このgistではRでMeSHを扱うためのR/Bioconductorパッケージ群(MeSH ORA Framework)を紹介します。 オリジナルのMeSHのデータは、NLMのFTPサイトから取得できます。 データは、d2012.bin, q2015.binというプレーンテキストでダウンロードされます。 これらのファイルが図.2のデータ整形パイプラインで処理され、対応するデータがSQLite3ファイルとして保存されたのち、MeSH.db, MeSH.AOR.db, MeSH.PCR.dbパッケージとして提供されています。

1.1 MeSHについて

MeSH.dbパッケージは、MeSH ID、MeSHターム、MeSHカテゴリー、同義語、Qualifer ID、Qualiferタームの情報を提供します。 Qualiferタームは、MeSHよりもより大雑把な注釈情報を与えます。 MeSHはGOのように階層構造を持っています。 階層構造の情報は、MeSH.AOR.db(祖先-子孫関係)やMeSH.PCR.db(親子関係)が提供します。

図2

1.2 NLM MeSH IDとNCBI Gene IDの対応関係

MeSH.XXX.eg.db (XXXは何か生物の学名の略語、例. Hsa: Homo sapiens)型パッケージは、NCBI Gene IDとNLMのMeSH IDとの対応関係を提供します。 対応関係は、以下の3種類の方法で作成されています。

手法 Gene ID - MeSH ID対応関係作成方法
Gendoo テキストマイニング
gene2pubmed NCBIのキュレーターによる手作業でのキュレーション
RBBH (双方向ベストヒット) 双方向BLASTPによる、相同性検索

GendooはPubMedデータのテキストマイニングに基づいたWebアプリケーションです。 Gene ID - MeSH IDの共起関係は、PubMed文献は網羅的に取得され、情報科学的手法により、より関連性が高いものだけがフィルタリングされています。

gene2pubmedはNCBI Gene IDとNLM PubMed IDとの対応関係をまとめたファイルです。 これらの対応関係は、NCBIのキュレーターチームによって、手作業で割り当てられています。 私達は、独自にライセンス版のPubMedのデータを取得しており、このデータとgene2pubmedのデータをマージすることで、Gene IDとMeSH IDの対応関係を作成しています。

幾つかのマイナーな生物種に対しては、Gendooでもgene2pubmedでも十分なアノテーションが得られなかったため、私達はあらかじめ、よくアノテーションがついている15生物種と(メジャー生物種)、アノテーションが不十分な100生物種(マイナー生物種)を定義し、これらの間の全組み合わせに対してBLASTP (E-value < 10^−50)を実行し、ヒットしたメジャー生物種に割り当てられたMeSHタームを、マイナー生物種のMeSHタームとして二次利用しました。

図3

1.3 MeSH ORA Frameworkのための、データベースインターフェース

私達は、MeSHDbiというデータベースインターフェース(DBI)を実装しました。 このパッケージは二つの重要な働きをします。 一つ目に、MeSH.db、MeSH.AOR.db、MeSH.PCR.db、and MeSH.XXX.eg.db型パッケージが継承するためのクラスを定義し、全てのパッケージの挙動を統一しています。 二つ目に、ユーザーがオリジナルのMeSH.XXX.eg.db型パッケージを作成することを支援します。 DNAシーケンサーの急速な普及により、様々な生物種のゲノム配列が次々と解読されています。 また、Gene IDとMeSH IDの対応関係は、Gendoo以外にも数多くのデータベースが作成・公開しています。 そのため、私達は、ユーザーが、なんらかの方法で、Gene ID - MeSH IDの対応関係を取得した状況を考え、MeSH.XXX.eg.db型パッケージの作成支援をする関数を実装しました。

1.4 MeSH ORA

実際のオミックス解析でMeSHパッケージを利用する状況を想定し、私達はmeshrという、MeSHデータを使ってORAを実行するパッケージを開発しました。 このパッケージは、内部で、MeSH.db、MeSH.AOR.db、MeSH.PCR.db、MeSH.XXX.eg.dbをインポートし、興味がある遺伝子リストに有意に含まれるMeSH IDを検出します。

図4

2. 演習

あらかじめ、このgistで必要なBioconductorパッケージをインストールしておきます。

source("http://bioconductor.org/biocLite.R")
biocLite("MeSH.db")
biocLite("MeSH.AOR.db")
biocLite("MeSH.PCR.db")
biocLite("MeSH.Hsa.eg.db")
biocLite("MeSH.Aca.eg.db")
biocLite("MeSH.Bsu.168.eg.db")
biocLite("MeSH.Syn.eg.db")
biocLite("MeSHDbi")
biocLite("meshr")
biocLite("fdrtool")
biocLite("cummeRbund")

2.1 MeSH.db、MeSH.AOR.db, MeSH.PCR.dbパッケージ

2.1.1 columns, keytypes, keys, select関数

MeSH ORA Frameworkでは全てのデータは、AnnotationDbiパッケージが規定した、columns, keytypes, keys, selectの4種類の関数を使ってアクセスします。MeSH.dbパッケージを使って、これらの関数の使い方を練習しましょう。まずMeSH.dbパッケージをロードします。

library("MeSH.db")

ls関数で、このパッケージ内で定義されたオブジェクトを見ることができます。

ls("package:MeSH.db")

MeSH.dbという名前のオブジェクトがあることがわかります。MeSH.dbと打ってみます。

MeSH.db

これは、MeSH.dbパッケージの内部にあるsqliteデータベースとのコネクションです。 このコネクションに対して、columns, keytypes, keys, select関数を実行することができます。 まずはcolumnsを使ってみます。この関数は、データベースがどのような項目を持っているかを返します。

columns(MeSH.db)

"CATEGORY", "MESHID", "MESHTERM", "QUALIFIER", "QUALIFIERID", "SYNONYM"という6種類のデータを持っていることがわかります。

次はkeytypesを使います。

keytypes(MeSH.db)

この関数は、データベースの行を一意に定めることができるデータ(主キー)を返します。

次は、keysを使います。

k <- keys(MeSH.db, keytype="MeSHID")
length(k)
head(k)

この関数は、キーの具体的な値をMeSH.dbから取得します。 最後に、selectです。 この関数では、columns,keys, keytypeという3種類のパラメーターを定める必要があります。columnsは取得したいカラム、keytypeは主キーとして指定するカラム、keysは主キーの具体的な値を意味します。さきほどのkeysで取得した、kを使ってデータを取得してましょう。

select(MeSH.db, columns = c("MESHID", "MESHTERM"), keytype = "MESHID", keys = k[1:10])

これにより、MESHIDがD000001〜D000010のMESHIDを持つ、MESHTERMを取得することができました。

2.1.2 白血病のMeSH Termを取得

次は、MeSH.dbを使って、白血病(Leukemia)関連の情報をひきだしてみましょう。 select関数を使って"Leukemia"というMeSH Termがあるか検索をかけます。

LEU <- select(MeSH.db, keys = "Leukemia", columns = c("MESHID", "MESHTERM", "CATEGORY", "SYNONYM"), keytype = "MESHTERM")
LEU

これにより、LeukemiaはCカテゴリー(Diseases)に分類され、MeSH IDはD007938だとわかりました。 またLeukemiaはLeucocythaemias、Leucocythaemia、Leucocythemias、Leukemiasといった同義語があることもわかりました。

上で説明したように、MeSHは階層構造を持っています。 MeSH.AOR.db, MeSH.PCR.dbがこれら階層構造の情報を提供します。 MeSH.AOR.dbにより、以下のようにLeukemiaより上の階層のタームを調べることができます。

library("MeSH.AOR.db")
ANC <- select(MeSH.AOR.db, keys = "D007938", columns = c("ANCESTOR", "OFFSPRING"), keytype = "OFFSPRING")
ANC

これにより、D009370というMeSH IDがLeukemiaより上くらいに位置していることがわかります。 このMeSH IDを以下のようにMeSH Termに変換します。

select(MeSH.db, keys = ANC[1, 1], columns = c("MESHTERM"), keytype = "MESHID")

これにより、LeukemiaはNeoplasms by Histologic Typeの一つとして分類されていることがわかりました。 MeSH.AOR.dbに対して、keytypeパラメーターを逆に(OFFSPRINGからANCESTORに)書き換えることで、より下位の階層のMeSH IDを取得することもできます。

OFF <- select(MeSH.AOR.db, keys = "D007938", columns = c("ANCESTOR", "OFFSPRING"), keytype = "ANCESTOR")
OFF
select(MeSH.db, keys = OFF[, 2], columns = c("MESHTERM"), keytype = "MESHID")

MeSH.PCR.dbは直下(または直上)のタームだけを返します。

library("MeSH.PCR.db")
 CHI <- select(MeSH.PCR.db, keys = LEU[1, 1], columns = c("PARENT", "CHILD"), keytype = "PARENT")
head(CHI)
head(select(MeSH.db, keys = CHI[, 2], columns = c("MESHTERM"), keytype = "MESHID"))

これにより、Leukemiaは、AvianLeukosis、BlastCrisis、Leukemia、Erythroblastic、Acuteといった、幾つものサブタイプを持っていることがわかりました。

2.1.3 その他の関数

その他の関数はより複雑なデータ取得をする際に利用できます。 この節では、SQLの知識が要求されます。 dbInfo関数は、このパッケージの情報を返します。 dbFileはコンピューター内でのsqliteファイルの位置を返します。 dbschemaはデータベースのスキーマ(設計図)を返します。 dbconnは、RSQLiteで作られたデータベースのコネクションを返します。

dbInfo(MeSH.db)
dbfile(MeSH.db)
dbschema(MeSH.db)
dbconn(MeSH.db)

dbschemaを見ると、データはDATAというテーブルでsqliteデータベース内に保存されており、MESHID, MESHTERM, CATEGORY, SYNONYM, QUALIFIERID, QUALIFIERという6つのカラムを持っていることがわかります。 そこで以下のように、より複雑なSQLクエリを投げることも可能です。

library("RSQLite")
SQL1 <- paste("SELECT MESHTERM, QUALIFIERID, QUALIFIER FROM DATA", "WHERE MESHID = 'D000001'", "AND QUALIFIERID = 'Q000494'")
dbGetQuery(dbconn(MeSH.db), SQL1)
SQL2 <- paste("SELECT ANCESTOR, OFFSPRING FROM DATA", "WHERE OFFSPRING = 'D000002'", "OR OFFSPRING = 'D000003'", "OR OFFSPRING = 'D000004'", "OR ANCESTOR = 'D009275'") > dbGetQuery(dbconn(MeSH.AOR.db), SQL2)
SQL3 <- paste("SELECT PARENT, CHILD FROM DATA", "WHERE PARENT = 'D000005'", "AND NOT CHILD = 'D004312'")
dbGetQuery(dbconn(MeSH.PCR.db), SQL3)

2.2 MeSH.XXX.eg.db型パッケージ

2.2.1 89生物種のアノテーション

MeSH.db, MeSH.AOR.db, MeSH.PCR.dbと同様に、MeSH.XXX.eg.db型パッケージも、データを取得するために、keytypes, columns, keys, selectを実行することができます。

library("MeSH.Hsa.eg.db")
columns(MeSH.Hsa.eg.db)
keytypes(MeSH.Hsa.eg.db)
key_HSA <- keys(MeSH.Hsa.eg.db, keytype = "MESHID")
select(MeSH.db, keys = key_HSA[1:10], columns = c("MESHID", "MESHTERM"), keytype = "MESHID")

また、これらのパッケージは、その他にも、species, nomenclature, listDatabasesという関数を持っています。 speciesは生物種の一般名を、nomenclatureは学名を返します。

library("MeSH.Aca.eg.db")
library("MeSH.Bsu.168.eg.db")
library("MeSH.Syn.eg.db")
species(MeSH.Hsa.eg.db)
species(MeSH.Aca.eg.db)
species(MeSH.Bsu.168.eg.db)
species(MeSH.Syn.eg.db)
nomenclature(MeSH.Hsa.eg.db)
nomenclature(MeSH.Aca.eg.db)
nomenclature(MeSH.Bsu.168.eg.db)
nomenclature(MeSH.Syn.eg.db)

listDatabasesは、どのような情報源でGene ID - MeSH IDの対応関係を作成したのかを返します。 RBBHに関しては、BLASTPを投げた相手先の生物種名を返します。 これらの値は2.4.1 MeSH ORAで必要になります。

listDatabases(MeSH.Hsa.eg.db)
listDatabases(MeSH.Aca.eg.db)
listDatabases(MeSH.Bsu.168.eg.db)
listDatabases(MeSH.Syn.eg.db)

2.3 MeSHDbiパッケージ

2.3.1 ユーザー独自のMeSH.XXX.eg.dbパッケージの作成

ほとんどの人は意識せずに使っていますが、MeSHDbiはMeSH ORA Frameworkで重要な働きをしています。 このパッケージは、MeSH関連のアノテーションパッケージのクラス定義をしています。 また、makeGeneMeSHPackage関数により、ユーザーが独自にMeSH.XXX.eg.db型パッケージを作成することを支援します。詳しくは、example("makeGeneMeSHPackage")で出力されるコードを参考にしてください。

library("MeSHDbi")
example("makeGeneMeSHPackage")

2.4 meshrパッケージ

2.4.1 MeSH ORA

meshrパッケージは、MeSHを利用してORAを実行するように開発されています。 実装は、GOstatsと類似するように書かれているため、GOstatsユーザーは容易にMeSH ORA Frameworkを利用することができます。

以下のように実行すると、cummeRbundパッケージで利用されていた、ヒトのiPS細胞とES細胞データを取得することができます。

 library("meshr")
data(geneid.cummeRbund)
data(sig.geneid.cummeRbund)

これにより303の遺伝子(ここでは全遺伝子とみたてる)と104のそこから選ばれた遺伝子(ここでは発現変動遺伝子とみたてる)が出てくるので、これを使ってMeSH ORAを実行してみましょう。

dim(geneid.cummeRbund)[1]
dim(sig.geneid.cummeRbund)[1]

ヒトのデータなので、MeSH.Hsa.eg.dbパッケージを呼び出します。 またfdrtoolパッケージも必要です。

library("fdrtool")
library("MeSH.Hsa.eg.db")

最初に以下のように、"MeSHHyperGParams"という名前のオブジェクト内でパラメーターを設定します。 geneIdsには、先ほどの104の発現変動遺伝子を、universeGeneIdsには、303の全遺伝子を、annotationには、使用するMeSH.XXX.eg.db型パッケージの名前を、categoryには、MeSH ORAを実行したいMeSHのカテゴリーを、databaseには使用したいデータベースを、pvalueCutoffには、有意とするp値(またはFDR制御をする場合q値)の閾値を設定します。 pAdjustではFDR制御法を設定することができ、

meshParams <- new("MeSHHyperGParams", geneIds = sig.geneid.cummeRbund[,2], universeGeneIds = geneid.cummeRbund[, 2], annotation = "MeSH.Hsa.eg.db",
category = "C", database = "gendoo", pvalueCutoff = 0.05,
pAdjust = "none")

meshHyperGTest関数により、超幾何分布のp値が計算されます。

meshR <- meshHyperGTest(meshParams)
meshR

結果の詳細は、summary関数で取得することができます。これにより、有意にエンリッチしたMeSH ID、MeSHターム、p値を取得することができます。

head(summary(meshR))

他のMeSHカテゴリーやデータベースに切り替えることも簡単です。 例えば、カテゴリーはG(Phenomena and Processes)、データベースはgene2pubmedにしたい場合、以下のように実行するだけです。

category(meshParams) <- "G"
database(meshParams) <- "gene2pubmed"
meshR <- meshHyperGTest(meshParams)
meshR
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment