macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

カスタムアノテーションを使った GO enrichment解析の例

2022/01/08追記, 1/13インストール追記

 

タイトルの通りの内容です。

質問があったのでそれに対応した記事になります。

 

1、アノテーションファイルの準備

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のような手順で1行に遺伝子名<tab>GO1つのみ、に変換してから使うことを考えるかもしれないが、topGOはGOの階層関係を読み込む。親子関係の情報を無くしてしまうとおそらく結果が変化するので気を付ける。また、GOがアサインされていない遺伝子は除かないとエラーになるので注意。

 

 

2、クエリの準備

トランスクリプトーム解析やプロテオーム解析などで得た関心のある遺伝子リスト(典型的には発現変動遺伝子のリスト)を準備する。1のアノテーション結果と同じ遺伝子名になっていないといけない。ここではファイル名”input”で保存した。

f:id:kazumaxneo:20220107235623p:plain

遺伝子名のみを1行1遺伝子ずつ記載したファイル。これで準備O.K

 

 

3、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 bioconda bioconductor-topgo -y
mamba install -c bioconda bioconductor-rgraphviz -y

 

準備が出来たらtopGOで解析する。流れはtopGOのマニュアルに則っている。

library(topGO)

#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)

f:id:kazumaxneo:20220108140251p:plain

 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行目は有意な遺伝子数とその遺伝子にアノテーションされた総遺伝子数を示している。

f:id:kazumaxneo:20220108003029p:plain

 

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

 

補足1

#1遺伝子に複数GOがアサインされているが、これを1行に1遺伝子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)

 

補足2

上のリンクからダウンロードできるtopGOのソースコードにはカスタムアノテーションのデモデータとしてtopGO/inst/geneid2go.mapが含まれています。マニュアルの4.3でも使用されているファイルと思われます。

 

関連