2022/01/08追記, 1/13インストール追記
2022/12/25 説明追記
タイトルの通りの内容です。
質問があったのでそれに対応した記事になります。
１、アノテーションファイルの準備
TrinotateかeggNOG mapperを使ってGO termをアサインしていることを想定している。
A - Trinotateのアノテーション結果
Trinotate（Trinityやtransdecoderを開発したBrian Haas（参考）が作ったアノテーションツール）のアノテーション結果のSQLiteファイルを所有しているなら、orf_to_go.rbスクリプト（レポジトリ）を使って遺伝子（転写産物）とGO termのリストを抽出できる。
git clone https://bitbucket.org/njcaruana/trinotate-proteomics.git
cd trinotate-proteomics/
./orf_to_go.rb Trinotate.sqlite > orf2go.map
TrinotateはTrinityと連携しており色々連携可能だが、モデル生物から距離がある場合、GO termのアサイン率はeggNOG mapperと比べて劣る可能性がある。
B - eggNOG mapperのアノテーション結果
eggNOG mapperのアノテーション結果(.tsv)を所有しているなら、既にGO termやKO termがアサインされた結果を所有していることになる。コメント行を捨て、必要な列だけ保存する。
#先頭と末尾のコメント行(##)を捨てる
grep -v "##" eggnog_result.tsv > filtered.tsv
#GOは10列目（12列目はKO）
cut -f 1,10 filtered.tsv > orf2go.map
*下の補足１のような手順で1行に遺伝子名<tab>GO１つのみ、に変換してから使うことを考えるかもしれないが、topGOはGOの階層関係を読み込む。親子関係の情報を無くしてしまうと結果が変化するので気を付ける。また、GOがアサインされていない遺伝子は除かないとエラーになるので注意。
このようなファイルを用意する。ここではタブ区切り。１行目は#query<tab>GOs
orf2go.map
２、クエリの準備
トランスクリプトーム解析やプロテオーム解析などで得た関心のある遺伝子リスト（典型的には発現変動遺伝子のリスト）を準備する。１のアノテーション結果と同じ遺伝子名になっていないといけない。ここではファイル名”input”で保存した。
遺伝子名のみを1行１遺伝子ずつ記載したファイル。これで準備O.K
３、GOエンリッチメント解析
準備が出来たらGOエンリッチメント解析を行う。カスタムアノテーションに対応したいくつかの実装があり、webでもagriGO（紹介）のようなカスタムアノテーションが利用できるサービスがあるが、ここではカスタムアノテーションの使い方の詳しい説明があるBioconductorのtopGOパッケージを使う。
Bioconductor - topGO
https://bioconductor.riken.jp/packages/3.3/bioc/html/topGO.html
topGO manual
https://bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf
topGOをインストールする。これは自分の主観だが、Bioconductor のパッケージはBiocManagerで導入できるものの、ローカル環境にパッケージが増えてくると時間がかかる上に割と失敗したりと、非効率的な面が残っている。実はBioconductor のパッケージはbiocondaでも提供されているので、今回はmambaで仮想環境を作って素早く導入する（全てのパッケージ対応しているかどうかは不明）。base環境に入れないように注意する。
mamba create -n topGO -y
conda activate topGO
#conda (bioconda)
mamba install -c conda-forge -c bioconda bioconductor-topgo -y
mamba install -c conda-forge -c bioconda bioconductor-rgraphviz -y
準備が出来たらtopGOで解析する。流れはtopGOのマニュアルに則っている。
library(topGO)
library(Rgraphviz)
#readMappingsを使ってバックグラウンド（step1で準備）を読み込む。
geneID2GO <- readMappings(file = "orf2go.map", sep = "\t", IDsep = ",")
#check
str(head(geneID2GO))
#gene name
geneNames <- names(geneID2GO)
#DEGsなどのクエリのリスト（step2）を読み込む
list <- unlist(as.list(read.csv("input", header = TRUE)))
geneList <- factor(as.integer(geneNames %in% list ))
names(geneList) <- geneNames
#これで準備ができた。GODataオブジェクトを構築（ここではBPだがMFとCCも試す）
#nodeSize = 10 は、10未満の語彙からGO階層を刈り取るために使用される。
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO, nodeSize = 10)
#topGOdataクラスのオブジェクトを得たらエンリッチメント分析を始めることができる。
#ここでは遺伝子数に基づくフィッシャーの正確検定を行う
#runTest関数3つの引数（ topGOdata クラスのオブジェクト,GOグラフ構造を扱う方法,
検定統計量）を指定する
resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")
#（ほかにも遺伝子スコアに基づいて濃縮度を計算するKolmogorov-Smirnov ライクな検定が利用可能）classicかelimが利用可能。
resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")
#エンリッチメントテスト実施後、結果を分析し解釈するため、GenTableを使う。GenTable は最も有意な GO タームとそれに対応する p 値を分析する。
GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)
エンリッチメントテスト実施後、結果を分析し解釈するため、GenTableが用意されている。GenTable は最も有意な GO タームとそれに対応する p 値を分析する。
GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)
GO termの順位とp値を、classicとelim間で比較したりもできる。
GenTable(GOdata, classicFisher = resultFisher, classic = resultFis, elimKS = resultKS.elim, orderBy = "weight", ranksOf = "classic", topNodes = 10)
マニュアルでは各検定結果のp 値を比較するための流れも用意されている（マニュアル3.3）。
有意なGO termがどのようなものかGOグラフ上で調べることができる。
allRes <- GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)
#showSigOfNodesは誘導されたサブグラフを表示
showSigOfNodes(GOdata, score(resultFis), firstSigNodes = 5, useInfo = 'all')
#PDF保存（printGraphはshowSigOfNodesのワープ機能）
printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)
重要なノードは長方形で表される。プロットされたグラフは、これらの有意なノードによって生成される上位グラフである。色は相対的な重要度を表す。濃い赤（最も重要）から明るい黄色（最も重要でない）までの範囲。各ノードの中の文字は、最初の 2 行がGO 識別子とトリミングされた GO 名、3 行目には生の p-value 、4行目は有意な遺伝子数とその遺伝子にアノテーションされた総遺伝子数を示している。
topGOで利用できるどのような検定、メソッドを使うかについては、マニュアルのRunning the enrichment testを確認して下さい。
引用
Alexa A, Rahnenfuhrer J (2021). topGO: Enrichment Analysis for Gene Ontology. R package version 2.46.0.
Trinotate Proteomics
https://bitbucket.org/njcaruana/trinotate-proteomics/src/master/
参考
https://www.biostars.org/p/117294/
Enrichment analysis with topGO
https://gist.github.com/iracooke/73c58d11c921fae4fa7c
https://rpubs.com/aemoore62/TopGo_colMap_Func_Troubleshoot
補足１
#１遺伝子に複数GOがアサインされているが、これを1行に１遺伝子１語彙になるように仕分ける（変換後は同じ遺伝子が複数回登場する）。
eggnogmapper <- read.csv("filtered.tsv", header = TRUE, sep = "\t")
head(eggnogmapper ) #check
#GOは10列目、KOなら12列目
go <- eggnogmapper[c(1,10)]
head(go,3) #check
#変換
library(data.table)
go <- data.table(go)
go_sort <- go[, list(GOs = unlist(strsplit(GOs , ","))), by = X.query]
# 書き出す
write.table(go_sort, file = "orf2go.map", sep = "\t", quote=F)
補足２
上のリンクからダウンロードできるtopGOのソースコードにはカスタムアノテーションのデモデータとしてtopGO/inst/geneid2go.mapが含まれています。マニュアルの4.3でも使用されているファイルと思われます。
関連
