どの生物種までGO解析ができるか

GO解析について

ゲノム解析では、扱う遺伝子数が膨大なので、自分の興味がある遺伝子リストに関連した遺伝子の機能タームを抽出する、エンリッチメント解析がよく利用される

機能タームの情報源としては、Gene Ontology（GO）、BioCyc、Reactome、KEGG、BioCarta、MeSH、...など多岐にわたる

とにかく何らかの観点で、機能タームと遺伝子IDの対応が付いていれば、この解析は適用できるわけで、最近ではenrichr、MetaScape、GeneTrail2、GeneSetDBなど、一つのインターフェースで、複数の機能タームに対してのエンリッチメント解析を横断的に行うツールも登場している

特にenrichrは、132種類の機能タームDBに対して計算するので、ヒト・マウス・ラットなどの生物種を対象とした解析であれば、もう決定版に近い気がするが、家畜動物であったり、バクテリア、植物であったり、ちょっとマイナーな生物種にはまだ適用できない

というか、そもそも、検索するためのデータベースがこういった生物種ではそれほど整備されてない

GOはそういった中でも、比較的幅広い生物種を対象として遺伝子の情報をまとめているので、ここでは、GOのエンリッチメント解析（GO解析）が現在までのところ、どの程度の生物種で適用できるのかを調べた

ちなみに、自分が作ったMeSH Termでエンリッチメント解析するmeshrパッケージは74生物種までなら使えるので（2018/8/4現在）、マイナー生物種ではこちらがおすすめ（ https://www.bioconductor.org/packages/release/BiocViews.html#___MeSHDb ）

org.XX.eg.db型パッケージを利用した場合

BioconductorでGO解析するツールとしては、有名どころとして、GOstats、clusterProfiler、topGOなどが挙げられる

GO解析には、生物種ごとの遺伝子IDとGO IDの対応表が必要で、それらの情報はBioconductorでは、org.XX.eg.db（例 : org.Hs.eg.db）という名前のパッケージの、org.XX.egGOというオブジェクト（例 : org.Hs.egGO）として、整備されている

このタイプのパッケージは、以下のように、現在のところ19件Bioconductorに登録されている
https://bioconductor.org/packages/release/BiocViews.html#___OrgDb

なので、ここに書かれた19生物種で、GOstatsの動作確認をしてみた

なお、シロイヌナズナ（At）の時は拡張子が.eg.dbではなく、.tair.dbなので注意

（egはEntrez Gene IDの略だが、これだけTAIRデータベースのIDだからだと思われる）

同様にマラリア原虫（Pf）も、拡張子がplasmo.db、出芽酵母（Sc）も拡張子がsgd.dbになっている

# パッケージのダウンロード & インストール
install.packages("BiocManager")
BiocManager::install("org.Ag.eg.db")
BiocManager::install("org.At.tair.db") # これはhair.db
BiocManager::install("org.Bt.eg.db")
BiocManager::install("org.Ce.eg.db")
BiocManager::install("org.Cf.eg.db")
BiocManager::install("org.Dm.eg.db")
BiocManager::install("org.Dr.eg.db")
BiocManager::install("org.EcK12.eg.db")
BiocManager::install("org.EcSakai.eg.db")
BiocManager::install("org.Gg.eg.db")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("org.Mm.eg.db")
BiocManager::install("org.Mmu.eg.db")
BiocManager::install("org.Pf.plasmo.db") # これはplasmo.db
BiocManager::install("org.Pt.eg.db")
BiocManager::install("org.Rn.eg.db")
BiocManager::install("org.Sc.sgd.db") # これはsgd.db
BiocManager::install("org.Ss.eg.db")
BiocManager::install("org.Xl.eg.db")
BiocManager::install("GOstats")

# パッケージのロード
library("org.Ag.eg.db")
library("org.At.tair.db")
library("org.Bt.eg.db")
library("org.Ce.eg.db")
library("org.Cf.eg.db")
library("org.Dm.eg.db")
library("org.Dr.eg.db")
library("org.EcK12.eg.db")
library("org.EcSakai.eg.db")
library("org.Gg.eg.db")
library("org.Hs.eg.db")
library("org.Mm.eg.db")
library("org.Mmu.eg.db")
library("org.Pf.plasmo.db")
library("org.Pt.eg.db")
library("org.Rn.eg.db")
library("org.Sc.sgd.db")
library("org.Ss.eg.db")
library("org.Xl.eg.db")
library("GOstats")

# GO解析では、このorg.XX.egGOというオブジェクトが使われる
head(toTable(org.Ag.egGO))
head(toTable(org.At.tairGO))
head(toTable(org.Bt.egGO))
head(toTable(org.Ce.egGO))
head(toTable(org.Cf.egGO))
head(toTable(org.Dm.egGO))
head(toTable(org.Dr.egGO))
head(toTable(org.EcK12.egGO))
head(toTable(org.EcSakai.egGO))
head(toTable(org.Gg.egGO))
head(toTable(org.Hs.egGO))
head(toTable(org.Mm.egGO))
head(toTable(org.Mmu.egGO))
head(toTable(org.Pf.plasmoGO))
head(toTable(org.Pt.egGO))
head(toTable(org.Rn.egGO))
head(toTable(org.Sc.sgdGO))
head(toTable(org.Ss.egGO))
head(toTable(org.Xl.egGO))

# 全Gene IDを抽出
GeneID_Ag <- unique(toTable(org.Ag.egGO)$gene_id)
GeneID_At <- unique(toTable(org.At.tairGO)$gene_id)
GeneID_Bt <- unique(toTable(org.Bt.egGO)$gene_id)
GeneID_Ce <- unique(toTable(org.Ce.egGO)$gene_id)
GeneID_Cf <- unique(toTable(org.Cf.egGO)$gene_id)
GeneID_Dm <- unique(toTable(org.Dm.egGO)$gene_id)
GeneID_Dr <- unique(toTable(org.Dr.egGO)$gene_id)
GeneID_EcK12 <- unique(toTable(org.EcK12.egGO)$gene_id)
GeneID_EcSakai <- unique(toTable(org.EcSakai.egGO)$gene_id)
GeneID_Gg <- unique(toTable(org.Gg.egGO)$gene_id)
GeneID_Hs <- unique(toTable(org.Hs.egGO)$gene_id)
GeneID_Mm <- unique(toTable(org.Mm.egGO)$gene_id)
GeneID_Mmu <- unique(toTable(org.Mmu.egGO)$gene_id)
GeneID_Pf <- unique(toTable(org.Pf.plasmoGO)$gene_id)
GeneID_Pt <- unique(toTable(org.Pt.egGO)$gene_id)
GeneID_Rn <- unique(toTable(org.Rn.egGO)$gene_id)
GeneID_Sc <- unique(toTable(org.Sc.sgdGO)$gene_id)
GeneID_Ss <- unique(toTable(org.Ss.egGO)$gene_id)
GeneID_Xl <- unique(toTable(org.Xl.egGO)$gene_id)

# 何件がGO Termと関連付けられているか
length(GeneID_Ag) # 13件
length(GeneID_At) # 26862件
length(GeneID_Bt) # 5722件
length(GeneID_Ce) # 13963件
length(GeneID_Cf) # 9616件
length(GeneID_Dm) # 13730件
length(GeneID_Dr) # 20720件
length(GeneID_EcK12) # 3499件
length(GeneID_EcSakai) # 0件
length(GeneID_Gg) # 2200件
length(GeneID_Hs) # 19456件
length(GeneID_Mm) # 24423件
length(GeneID_Mmu) # 12255件
length(GeneID_Pf) # 4584件
length(GeneID_Pt) # 12797件
length(GeneID_Rn) # 19711件
length(GeneID_Sc) # 0件
length(GeneID_Ss) # 7812件
length(GeneID_Xl) # 18820件

# 疑似的なDEGを作成（全遺伝子の1/100をランダムに選択）
sig.GeneID_Ag <- sample(GeneID_Ag, length(GeneID_Ag)/100)
sig.GeneID_At <- sample(GeneID_At, length(GeneID_At)/100)
sig.GeneID_Bt <- sample(GeneID_Bt, length(GeneID_Bt)/100)
sig.GeneID_Ce <- sample(GeneID_Ce, length(GeneID_Ce)/100)
sig.GeneID_Cf <- sample(GeneID_Cf, length(GeneID_Cf)/100)
sig.GeneID_Dm <- sample(GeneID_Dm, length(GeneID_Dm)/100)
sig.GeneID_Dr <- sample(GeneID_Dr, length(GeneID_Dr)/100)
sig.GeneID_EcK12 <- sample(GeneID_EcK12, length(GeneID_EcK12)/100)
sig.GeneID_EcSakai <- sample(GeneID_EcSakai, length(GeneID_EcSakai)/100)
sig.GeneID_Gg <- sample(GeneID_Gg, length(GeneID_Gg)/100)
sig.GeneID_Hs <- sample(GeneID_Hs, length(GeneID_Hs)/100)
sig.GeneID_Mm <- sample(GeneID_Mm, length(GeneID_Mm)/100)
sig.GeneID_Mmu <- sample(GeneID_Mmu, length(GeneID_Mmu)/100)
sig.GeneID_Pf <- sample(GeneID_Pf, length(GeneID_Pf)/100)
sig.GeneID_Pt <- sample(GeneID_Pt, length(GeneID_Pt)/100)
sig.GeneID_Rn <- sample(GeneID_Rn, length(GeneID_Rn)/100)
sig.GeneID_Sc <- sample(GeneID_Sc, length(GeneID_Sc)/100)
sig.GeneID_Ss <- sample(GeneID_Ss, length(GeneID_Ss)/100)
sig.GeneID_Xl <- sample(GeneID_Xl, length(GeneID_Xl)/100)

# パラメーター設定
para_Ag <- new("GOHyperGParams",geneIds=sig.GeneID_Ag,universeGeneIds=GeneID_Ag,annotation="org.Ag.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_At <- new("GOHyperGParams",geneIds=sig.GeneID_At,universeGeneIds=GeneID_At,annotation="org.At.tair.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Bt <- new("GOHyperGParams",geneIds=sig.GeneID_Bt,universeGeneIds=GeneID_Bt,annotation="org.Bt.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Ce <- new("GOHyperGParams",geneIds=sig.GeneID_Ce,universeGeneIds=GeneID_Ce,annotation="org.Ce.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Cf <- new("GOHyperGParams",geneIds=sig.GeneID_Cf,universeGeneIds=GeneID_Cf,annotation="org.Cf.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Dm <- new("GOHyperGParams",geneIds=sig.GeneID_Dm,universeGeneIds=GeneID_Dm,annotation="org.Dm.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Dr <- new("GOHyperGParams",geneIds=sig.GeneID_Dr,universeGeneIds=GeneID_Dr,annotation="org.Dr.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_EcK12 <- new("GOHyperGParams",geneIds=sig.GeneID_EcK12,universeGeneIds=GeneID_EcK12,annotation="org.EcK12.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_EcSakai <- new("GOHyperGParams",geneIds=sig.GeneID_EcSakai,universeGeneIds=GeneID_EcSakai,annotation="org.EcSakai.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Gg <- new("GOHyperGParams",geneIds=sig.GeneID_Gg,universeGeneIds=GeneID_Gg,annotation="org.Gg.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Hs <- new("GOHyperGParams",geneIds=sig.GeneID_Hs,universeGeneIds=GeneID_Hs,annotation="org.Hs.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Mm <- new("GOHyperGParams",geneIds=sig.GeneID_Mm,universeGeneIds=GeneID_Mm,annotation="org.Mm.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Mmu <- new("GOHyperGParams",geneIds=sig.GeneID_Mmu,universeGeneIds=GeneID_Mmu,annotation="org.Mmu.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Pf <- new("GOHyperGParams",geneIds=sig.GeneID_Pf,universeGeneIds=GeneID_Pf,annotation="org.Pf.plasmo.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Pt <- new("GOHyperGParams",geneIds=sig.GeneID_Pt,universeGeneIds=GeneID_Pt,annotation="org.Pt.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Rn <- new("GOHyperGParams",geneIds=sig.GeneID_Rn,universeGeneIds=GeneID_Rn,annotation="org.Rn.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Sc <- new("GOHyperGParams",geneIds=sig.GeneID_Sc,universeGeneIds=GeneID_Sc,annotation="org.Sc.sgd.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Ss <- new("GOHyperGParams",geneIds=sig.GeneID_Ss,universeGeneIds=GeneID_Ss,annotation="org.Ss.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")
para_Xl <- new("GOHyperGParams",geneIds=sig.GeneID_Xl,universeGeneIds=GeneID_Xl,annotation="org.Xl.eg.db",ontology="BP",pvalueCutoff=0.05,conditional=FALSE,testDirection="over")

# エンリッチメント解析を実行
(resAg <- hyperGTest(para_Ag)) # 対応するGOなし
(resAt <- hyperGTest(para_At))
(resBt <- hyperGTest(para_Bt))
(resCe <- hyperGTest(para_Ce))
(resCf <- hyperGTest(para_Cf))
(resDm <- hyperGTest(para_Dm))
(resDr <- hyperGTest(para_Dr))
(resEcK12 <- hyperGTest(para_EcK12))
(resEcSakai <- hyperGTest(para_EcSakai)) # 対応するGOなし
(resGg <- hyperGTest(para_Gg))
(resHs <- hyperGTest(para_Hs))
(resMm <- hyperGTest(para_Mm))
(resMmu <- hyperGTest(para_Mmu))
(resPf <- hyperGTest(para_Pf))
(resPt <- hyperGTest(para_Pt))
(resRn <- hyperGTest(para_Rn))
(resSc <- hyperGTest(para_Sc))
(resSs <- hyperGTest(para_Ss)) # 対応するGOなし
(resXl <- hyperGTest(para_Xl))

以上から、パッケージとしてはまとめられているものの、遺伝子IDとGO IDの対応が極端に少ないものが3件あり、まともに使えそうなのは、それらを除いた以下の16生物種だと思われる。

• Arabidopsis thaliana（シロイヌナズナ）
• Bos taurus（ウシ）
• Caenorhabditis elegans（線虫）
• Canis familiaris（イヌ）
• Drosophila melanogaster（ショウジョウバエ）
• Danio rerio（ゼブラフィッシュ）
• Escherichia coli K12（大腸菌）
• Gallus gallus（ニワトリ）
• Homo sapiens（ヒト）
• Mus musculus（マウス）
• Macaca mulatta（アカゲザル）
• Plasmodium falciparum（マラリア原虫）
• Pan troglodytes（チンパンジー）
• Rattus norvegicus（ラット）
• Saccharomyces cerevisiae（出芽酵母）
• Xenopus laevis（アフリカツメガエル）

自力でGene ID-GO ID対応表から行う場合

ところで、Gene Ontologyのウェブサイトを見ると、生物種でみても27件ほどあったので、GO IDがGene IDに割り当てているものの、orgパッケージになっていない生物種が幾つも存在している（ http://www.geneontology.org/page/download-go-annotations ）

例えば、緑膿菌、分裂酵母、イネ、真菌系のorgパッケージは、Bioconductorに登録されていそうな気もするが、実際はまだ存在しない（MeSHパッケージでは作ったが）

おそらく、今回の3件のように、GO IDが十分な数割り当てられていない場合や、Bioconductor Core Teamがまだ対応していないものだと思われる

上記のように楽はできないものの、biopapyrusの例（ https://bi.biopapyrus.jp/pathway/go/r/gostats.html ）に習うように、GSEABaseパッケージと組み合わせれば、orgパッケージがない場合でもGOstatsを実行することはできるらしい

ここでは試しに緑膿菌で、GOstatsを利用してみる

まず緑膿菌はorgパッケージにはなってないので、Gene Ontologyのウェブサイトから直接、遺伝子IDとGO IDの対応表をダウンロードする

download.file(url="http://geneontology.org/gene-associations/gene_association.pseudocap.gz", destfile="pseudocap_valid.gaf")

このファイルをデータフレームとして読み込み、GOID、Evidence Code、Gene IDの順に並び替える。

# 2カラム目は遺伝子ID、5カラム目がGO、7カラム目がエビデンスコード
goframeData <- read.delim("pseudocap_valid.gaf", skip=5, header=FALSE)
goframeData <- goframeData[, c(5,7,2)]
colnames(goframeData) <- c("go_id", "Evidence", "gene_id")

次に、GSEABaseパッケージを使って、設定

# GSEABaseで、GOstats用のオブジェクトを作成する
library("GSEABase")
goFrame <- GOFrame(goframeData, organism = "Pseudomonas aeruginosa PAO1")
goAllFrame <- GOAllFrame(goFrame)
gsc <- GeneSetCollection(goAllFrame, setType = GOCollection())

最後に、GO解析のパラメーターの設定を、GSEAGOHyperGParamsで行い、hyperGTestでFisher's exact testのp値を算出する

# 全遺伝子
all <- unique(goframeData$gene_id)
# 疑似的なDEG
degs <- sample(all, 100)
p <- GSEAGOHyperGParams(
       name = "Paramaters",
       geneSetCollection = gsc,
       geneIds = degs,
       universeGeneIds = all,
       ontology = "MF",
       pvalueCutoff = 0.05,
       conditional = FALSE,
       testDirection = "over"
     )
result <- hyperGTest(p)
head(summary(result))

このように、GO ID、Evidence、Gene IDの対応表さえあれば、自力でなんとかGO解析ができることがわかった

なお、MeSH関連パッケージの一つであるMeSHDbiパッケージでは、makeGeneMeSHPackageという関数を用意しており、ユーザーがCSVでMeSH IDとGene IDの対応表を持っていれば、この関数でユーザーオリジナルのMeSHパッケージを作成でき、このパッケージとmeshrパッケージを組み合わせて使うことで、Bioconductorに登録されていない生物種でも、MeSHエンリッチメント解析ができるようにしている

興味があれば是非こちらも使ってみていただきたい（ https://www.bioconductor.org/packages/release/bioc/vignettes/meshr/inst/doc/MeSH.pdf ）