gafファイルを使って自身でGO enrichment解析をする [備忘録]

Gene Ontology(GO) とは、生物的プロセスなどの機能を整理し、それぞれの遺伝子に付与されるアノテーションのことをいう。GO enrichment解析はある生物がもつ全遺伝子に対して着目した遺伝子群において、どのGOタームが有意に濃縮されているかを調べる解析のことであり、着目した現象等がどのような意義をもつのかを見いだすのに有用なツールである。

線虫や酵母などでは解析ソフト中にデータベースが含まれているものがあるものの、データベースが用意されていない生物種も存在する。
本稿では、公開されている任意の生物種のGO annotation file (.gaf)をダウンロードし、GO enrichment解析を行う方法を備忘録としてまとめる。

## 解析に用いる生物種、データ、ソフト
本稿ではAspergillus oryzaeのRIB40を対象とした。
gafファイルはFungiDBの以下のURLより、FungiDB-64_AoryzaeRIB40_Curated_GO.gafを用いた。
また、GOの濃縮を見る遺伝子のリストは、後でsample関数でランダムに生成することにした。

https://fungidb.org/fungidb/app/downloads

解析はantiplastics(Koki Tsuyuzaki)さんのqiitaの記事に沿って行った(以下URL)。

https://qiita.com/antiplastics/items/add3f8438cbe6b2b594e

解析ソフトのversion等は以下の通り。パッケージは以下のものをダウンロードした前提で進める。

```
R 4.2.2 GUI 1.79 Big Sur ARM build (8160)

# 以下はパッケージ
library("GOstats")
library("GSEABase")

# パッケージをインストールする場合は以下
BiocManager::install("GOstats")

# versionは以下
GOstats: v2.64.0
GSEABase: v1.60.0
```
## 実際の解析
以下のスクリプトに従い、解析を進めた。

```
# gafファイルの読み込み, skipで先頭のいらない行をとばした
goframeData <- read.delim("FungiDB-64_AoryzaeRIB40_Curated_GO.gaf", skip=1, header=FALSE)

# GO、遺伝子ID、エビデンスコードの順番で並べる
goframeData <- goframeData[, c(5,7,2)]
colnames(goframeData) <- c("go_id", "Evidence", "gene_id")

# Go stats用のオブジェクト作成
goFrame <- GOFrame(goframeData, organism = "Aspergillus oryzae RIB40")
goAllFrame <- GOAllFrame(goFrame)
gsc <- GeneSetCollection(goAllFrame, setType = GOCollection())

# 全遺伝子を格納
all <- unique(goframeData$gene_id)

# GO解析したい遺伝子のリストの読み込み
# 読み込んだファイルがリストの場合には、unlist関数を使って変換する
gene<-as.vector(unlist(read.table(file="gene_interest.txt",header=F,sep=",")))
sample_gene<-intersect(gene,all) # ここで，gafファイルに記入されていない遺伝子名を除外する(ここを忘れると，この後動かないことがある)

# 今回は"sample_gene"の代わりに，ランダムに選んだ遺伝子を使って，GO解析する。
randam_gene <- sample(all, 200)

# p-valueは0.05、など各種のパラメータを設定。
# ここではMF (Molecular Function)を表示する。
p <- GSEAGOHyperGParams(
       name = "Paramaters",
       geneSetCollection = gsc,
       geneIds = randam_gene,
       universeGeneIds = all,
       ontology = "MF",
       pvalueCutoff = 0.05,
       conditional = FALSE,
       testDirection = "over"
     )

# Fisher's exact testのp値
result <- hyperGTest(p)

# 結果の表示
head(summary(result))
```

実際に上記のスクリプトで解析した結果が以下となる。

```
> head(summary(result))
      GOMFID      Pvalue OddsRatio   ExpCount Count Size                                                Term
1 GO:0043295 0.001940888  54.47260 0.07305936     2    4                                 glutathione binding
2 GO:1900750 0.001940888  54.47260 0.07305936     2    4                                oligopeptide binding
3 GO:1901681 0.003196058  36.31050 0.09132420     2    5                             sulfur compound binding
4 GO:0072341 0.004736734  27.22945 0.10958904     2    6                         modified amino acid binding
5 GO:0003962 0.018264840       Inf 0.01826484     1    1               cystathionine gamma-synthase activity
6 GO:0044016 0.018264840       Inf 0.01826484     1    1 histone acetyltransferase activity (H3-K4 specific)
```
パラメータのontologyの部分をCC (Cellular Component)やBP (Biological Process)にしても良い。

```
# CCの場合
      GOCCID       Pvalue  OddsRatio   ExpCount Count Size                                                             Term
1 GO:0043291 0.0009296942 111.836957 0.05326460     2    3                                                     RAVE complex
2 GO:0044613 0.0177548683        Inf 0.01775487     1    1                           nuclear pore central transport channel
3 GO:0017117 0.0177548683        Inf 0.01775487     1    1 single-stranded DNA-dependent ATP-dependent DNA helicase complex
4 GO:1902929 0.0177548683        Inf 0.01775487     1    1                              plasma membrane of growing cell tip
5 GO:0033202 0.0177548683        Inf 0.01775487     1    1                                             DNA helicase complex
6 GO:0022625 0.0237693322   5.255123 0.62142039     3   35                                cytosolic large ribosomal subunit

# BPの場合
      GOBPID      Pvalue OddsRatio   ExpCount Count Size                                              Term
1 GO:0070900 0.002102448  52.21341 0.07605085     2    4                   mitochondrial tRNA modification
2 GO:1900864 0.002102448  52.21341 0.07605085     2    4                    mitochondrial RNA modification
3 GO:0000959 0.002595646  13.11810 0.28519070     3   15               mitochondrial RNA metabolic process
4 GO:0090646 0.003460322  34.80488 0.09506357     2    5                     mitochondrial tRNA processing
5 GO:0070525 0.003460322  34.80488 0.09506357     2    5 tRNA threonylcarbamoyladenosine metabolic process
6 GO:0006418 0.003772471  11.24145 0.32321613     3   17       tRNA aminoacylation for protein translation

```