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関数でランダムに生成することにした。
解析はantiplastics(Koki Tsuyuzaki)さんのqiitaの記事に沿って行った(以下URL)。
解析ソフトの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