このgistの内容は、Bioconductorのmeshrパッケージのvignette（パッケージの使用方法が記されたドキュメント）を和訳したものです。
この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もの臨床と生物学的な語彙を有しています。
語彙数は、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パッケージとして提供されています。
MeSH.dbパッケージは、MeSH ID、MeSHターム、MeSHカテゴリー、同義語、Qualifer ID、Qualiferタームの情報を提供します。 Qualiferタームは、MeSHよりもより大雑把な注釈情報を与えます。 MeSHはGOのように階層構造を持っています。 階層構造の情報は、MeSH.AOR.db（祖先-子孫関係）やMeSH.PCR.db（親子関係）が提供します。
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タームとして二次利用しました。
私達は、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型パッケージの作成支援をする関数を実装しました。
実際のオミックス解析でMeSHパッケージを利用する状況を想定し、私達はmeshrという、MeSHデータを使ってORAを実行するパッケージを開発しました。 このパッケージは、内部で、MeSH.db、MeSH.AOR.db、MeSH.PCR.db、MeSH.XXX.eg.dbをインポートし、興味がある遺伝子リストに有意に含まれるMeSH IDを検出します。
あらかじめ、この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")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を取得することができました。
次は、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といった、幾つものサブタイプを持っていることがわかりました。
その他の関数はより複雑なデータ取得をする際に利用できます。 この節では、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)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)ほとんどの人は意識せずに使っていますが、MeSHDbiはMeSH ORA Frameworkで重要な働きをしています。 このパッケージは、MeSH関連のアノテーションパッケージのクラス定義をしています。 また、makeGeneMeSHPackage関数により、ユーザーが独自にMeSH.XXX.eg.db型パッケージを作成することを支援します。詳しくは、example("makeGeneMeSHPackage")で出力されるコードを参考にしてください。
library("MeSHDbi")
example("makeGeneMeSHPackage")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