GOstats

GOstats は GO enrichment analysis test を行う R / Bioconductor パッケージである。

GOstats パッケージは AnnotationForge パッケージに保存されていアノテーションるデータを利用して解析を行う。次のコマンドを R 上で実行すると、GO 解析がサポートされている全生物種を確認できる。

library(AnnotationForge)
available.dbschemas()
##  [1] "ARABIDOPSISCHIP_DB" "BOVINE_DB"          "BOVINECHIP_DB"     
##  [4] "CANINE_DB"          "CANINECHIP_DB"      "CHICKEN_DB"        
##  [7] "CHICKENCHIP_DB"     "ECOLI_DB"           "ECOLICHIP_DB"      
## [10] "FLY_DB"             "FLYCHIP_DB"         "GO_DB"             
## [13] "HUMAN_DB"           "HUMANCHIP_DB"       "HUMANCROSSCHIP_DB" 
## [16] "INPARANOIDDROME_DB" "INPARANOIDHOMSA_DB" "INPARANOIDMUSMU_DB"
## [19] "INPARANOIDRATNO_DB" "INPARANOIDSACCE_DB" "KEGG_DB"           
## [22] "MALARIA_DB"         "MOUSE_DB"           "MOUSECHIP_DB"      
## [25] "PFAM_DB"            "PIG_DB"             "PIGCHIP_DB"        
## [28] "RAT_DB"             "RATCHIP_DB"         "WORM_DB"           
## [31] "WORMCHIP_DB"        "YEAST_DB"           "YEASTCHIP_DB"      
## [34] "ZEBRAFISH_DB"       "ZEBRAFISHCHIP_DB"  

サポートされていない生物種は、blast2GO などを利用して独自のアノテーションファイルを作成する。アノテーションファイルを自作する場合は、3 列からなるタブ区切りのテキスト形式で作成する。1 列目は GO ID、2 列目はエビデンスコード、3 列目には Ensembl Gene ID を書き入れる。R に読み込むときは read.table で読み込みます。

アノテーションファイル(ann.txt)の例。

GO_id	Evidence	Gene_id
GO:0008150	ND	1
GO:0001869	IDA	2
GO:0007597	TAS	2

自作のアノテーションファイルを R に読み込む例。

goframeData <- read.table("ann.txt", header = TRUE)

GOstats の使い方

GOstats による GO 解析例。まずアノテーションデータを用意する。

library(org.Hs.eg.db)

frame <- toTable(org.Hs.egGO)
goframeData <- data.frame(frame$go_id, frame$Evidence, frame$gene_id)

# アノテーションファイルから読み込む場合
# goframeData <- read.table("annotation.txt")

次に、アノテーションデータを GOstats で利用できるように変換する。

library(GSEABase)

goFrame <- GOFrame(goframeData, organism = "Homo sapiens")
goAllFrame <- GOAllFrame(goFrame)

gsc <- GeneSetCollection(goAllFrame, setType = GOCollection())

準備が整えた段階で、GOstats を利用して GO 解析を行う。

library(GOstats)

all <- Lkeys(org.Hs.egGO)
degs <- all[1:500]
p <- GSEAGOHyperGParams(
       name = "Paramaters",
       geneSetCollection = gsc,
       geneIds = degs,
       universeGeneIds = all,
       ontology = "MF",
       pvalueCutoff = 0.05,
       conditional = FALSE,
       testDirection = "over"
     )
result <- hyperGTest(p)

head(summary(result))
##       GOMFID       Pvalue OddsRatio    ExpCount Count Size                                   Term
## 1 GO:0016805 0.0001767507 127.57500 0.020074301     2   14                   dipeptidase activity
## 2 GO:0019239 0.0007253218  58.82692 0.040148602     2   28                     deaminase activity
## 3 GO:0004180 0.0012680042  43.67429 0.053053510     2   37              carboxypeptidase activity
## 4 GO:0033882 0.0014338786       Inf 0.001433879     1    1         choloyl-CoA hydrolase activity
## 5 GO:0050129 0.0014338786       Inf 0.001433879     1    1 N-formylglutamate deformylase activity
## 6 GO:0004060 0.0028657946 729.52381 0.002867757     1    2 arylamine N-acetyltransferase activity


# GO解析結果を保存
write.table(summary(result), file = "results.txt")