goseq

統計解析用フリーソフトウェア R では、GO 解析関連のパッケージがいくつか提供されている。goseq (Bioconductor)はその一つである。RNA-seq よる発現変動遺伝子の同定では、トランスクリプトームの長さが発現変動遺伝子の同定結果に影響を与える。goseq はこのような影響を取り除くように作られている。

以下はサンプルデータに対して goseq による GO 解析例である。サンプルデータは 2 列からなる。1 列目は遺伝子の Ensembl ID であり、2 列目は edgeR などのパッケージにより発現変動遺伝子を同定した結果として得られる p-value である。

goseq パッケージが GO 解析を行う際に、その遺伝子が発現変動遺伝子かどうかの情報とその塩基長の 2 種類の情報を与える必要がある。発現変動遺伝子かどうかの情報はベクトルとして与える。ベクトルの要素は 0 と 1 からなり、0 は非発現変動遺伝子を表し、1 は発現変動遺伝子を表す。また、塩基長に関してはヒトゲノムのアノテーション情報(hg19)を与えれば、プログラムが自ら塩基長を調べてくる。

library(goseq)

# データを読み込む(1 列目は Ensembl ID、2 列目は p-value)
x <- read.table("https://bi.biopapyrus.net/data/go.sample.txt")


# p-value から FDR を計算する
fdr <- p.adjust(x[, 2], method = "BH")

# 発現変動遺伝子判定(FDR < 0.01 ならば発現変動遺伝子とする)
is.DEG <- as.numeric(fdr < 0.01)
names(is.DEG) <- x[,1]

# GO 解析
# ヒトゲノムアノテーション(hg19)を利用して遺伝子の塩基長を調べる
# その際に Ensembl ID を利用するため ensGene を指定する 
pwf <- nullp(is.DEG, "hg19", "ensGene")
result <- goseq(pwf, "hg19", "ensGene", test.cats = "GO:MF")  # MF のみ
head(result)
##        category over_represented_pvalue under_represented_pvalue numDEInCat numInCat
## 1403 GO:0015291            3.641180e-08                        1         18      101
## 1405 GO:0015293            6.489392e-08                        1         15       73
## 1183 GO:0008509            1.552203e-07                        1         20      122
## 1408 GO:0015296            2.732143e-07                        1          9       27
## 1406 GO:0015294            2.970684e-07                        1         11       43
## 1326 GO:0015103            4.022130e-07                        1         12       49

複数の GO 分類についても解析を行うことが可能。

# MF と BP
pwf <- nullp(is.DEG, "hg19", "ensGene")
result <- goseq(pwf, "hg19", "ensGene", test.cats = c("GO:MF", "GO:BP")) 

goseq を利用した GO 解析では、解析結果として GO の accession number と p-value は表記されるが、その定義などは表記されない。統計量とともにその定義も合わせて、表形式で出力するには、GO.db パッケージを呼び出して、調整する必要がある。

library(GO.db)

# GO 解析
result <- goseq(pwf, "hg19", "ensGene", test.cats = "GO:MF")
go.pvalue <- result$over_represented_pvalue

# GO 解析結果の p-value から FDR を計算
go.fdr <- p.adjust(go.pvalue, method = "BH")

# FDR < 0.01 の GO だけを選ぶ
gotable <- result[go.fdr < 0.01, ]

# Term と Definition を保存するための変数を作成
godef <- data.frame(
           Term = rep("", length = nrow(gotable)),
           Definition = rep("", length = nrow(gotable)),
           stringsAsFactors = FALSE
        )

# GO accession から Term と Definition を探す
for (i in 1:nrow(gotable)) {
  acc <- gotable[i, "category"]
  godef[i, "Term"] <- GOTERM[[acc]]@Term[1]
  godef[i, "Definition"] <- GOTERM[[acc]]@Definition[1]
}

# GO 解析結果と GO 定義を結合する
table <- cbind(gotable, godef)
#write.table(table, file = "GO.result.txt", row.names = FALSE, sep = "\t")