將有興趣的基因列表使用org.Hs.eg.db查找最新的對應GO terms,再使用topGO做GO analysis

topGO package在bioconductor裡已經超過10年了,算是最早R裡面可以用來做GO term enrichment analysis的工具,但也因此裡面處理輸入資料的方式相對的來說頗不方便,語法比較接近R base,所以對於用習慣tidyverse體系的package的人來說會有點頭痛,這邊則來分享一下,如何使用topGO來做GO analysis

使用topGO做GO analysis主要有五個步驟:
1. 將基因列表轉換為Entrez ID
2. 使用org.Hs.eg.db資料庫分別找出每個Enrez ID所對應到的GO terms
3. 將資料結構處理成topGO中所使用的topGOdata object
4. 使用runTest進行GO analysis
5. 使用GenTable來獲取結果成data.frame

最麻煩的地方是理解topGOdata object建立所需要的常數,和每個輸入的結構!

這邊補充一下,各個下面代碼的重要資料結構長得樣子

screenshot.png

library(annotationDbi)
library(org.Hs.eg.db)
library(purrr)
library(dplyr)
library(topGO)

all.up.GO.Sphere_CLS.entrezID  <-  up.tukey.aov.Sphere_CLS %>% with(EntrezID)

#使用annotationDbi中提供的select方法,query的結果為data.frame的形式
all.up.GO.Sphere_CLS <- AnnotationDbi::select(org.Hs.eg.db,                                            keys = all.up.GO.Sphere_CLS.entrezID [!is.na(all.up.GO.Sphere_CLS.entrezID )] %>% as.character,
                                           columns = "GO",
                                           keytype = "ENTREZID")

#將entrezID對應GO term的data.frame格式轉換成list-like結構

all.up.GO.Sphere_CLS.list <- all.up.GO.Sphere_CLS %>%
                        filter(EVIDENCE != "IEA") %>% split(.$ONTOLOGY, .$ENTREZID)
all.up.GO.BP.Sphere_CLS   <- all.up.GO.Sphere_CLS.list$BP
all.up.GO.BP.Sphere_CLS   <- all.up.GO.BP.Sphere_CLS %>% split(.$ENTREZID) %>% map(with, GO)
all.up.GO.BP.Sphere_CLS.GO2geneID <- inverseList(all.up.GO.BP.Sphere_CLS)

#建立topGO所需要輸入參數allGenes, gene2GO兩個輸入變數geneList,all.up.GO.BP.Sphere_CLS
geneNames  <- names(all.up.GO.BP.Sphere_CLS)
geneList <- factor(as.integer(geneNames %in% (up.Sphere_CLS$EntrezID[!is.na(up.Sphere_CLS$EntrezID)] %>% as.character)))
names(geneList) <- geneNames

# construct topGOdata
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,
              annot = annFUN.gene2GO, gene2GO = all.up.GO.BP.Sphere_CLS)

# calculate the result
resultFisher <- runTest(GOdata , algorithm = "classic", statistic = "fisher")
resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")
resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks")

# query the result
allRes  <- GenTable(GOdata, classicFisher= resultFisher,
                    classicKS = resultKS, elimKS = resultKS.elim,
                    orderBy = "classicFisher", ranksOf = "classicFisher", topNodes = 200)
allRes  <- allRes %>% transmute(GO.ID=GO.ID,
                                Term=Term,
                                Annotated=Annotated,
                                Significant=Significant,
                                Expected=Expected,
                                classicFisher=as.numeric(classicFisher),
                                classicKS=as.numeric(classicKS),
                                elimKS=as.numeric(elimKS))

allRes.up.Sphere_CLS <- allRes %>% filter(classicFisher < 0.05)

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 /  變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 /  變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 /  變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 /  變更 )

w

連結到 %s