topGO是一個半自動的GO富集包旗吁，該包的主要優(yōu)勢是集中了好幾種統(tǒng)計檢驗的方法理肺，目前支持的統(tǒng)計方法如下：

一、安裝

BiocManager::install('topGO')
需要R的版本為>=2.10，但biocmanager安裝需要的R版本更高礁遵，現(xiàn)在應該是3.6脱惰。

二搏嗡、數(shù)據(jù)準備

富集工作主要包括3個步驟：
1、準備相關數(shù)據(jù)拉一；
2采盒、進行富集統(tǒng)計檢驗；
3蔚润、分析結(jié)果磅氨。
所以最重要的工作就是數(shù)據(jù)的準備。需要的數(shù)據(jù)包括包含全部geneID（背景基因名抽碌，一般是研究物種的全部基因）的文件悍赢，需要進行富集分析的geneID（差異表達基因或感興趣的基因）文件，還有g(shù)ene-to-GO的注釋文件货徙。

物種全部的geneID和差異基因ID比較容易獲得左权，比較費勁的是gene-to-GO文件。
topGO提供了一些函數(shù)來幫助我們自動獲取注釋信息：
annFUN.db：用于Bioconductor上有注釋包的物種的芯片數(shù)據(jù)痴颊；
annFUN.org：用于Bioconductor上有“org.XX.XX”注釋包的數(shù)據(jù)赏迟；
annFUN.gene2GO：用戶自己提供gene-to-GO文件；
annFUN.GO2gene：用戶提供的GO-to-gene文件也可以蠢棱；
annFUN.file：讀取有g(shù)ene2GO或GO2gene的txt文件锌杀。
一般Bioconductor提供的注釋物種并不多，我的方法主要是用AnnotationHub的select函數(shù)或biomaRt的getBM函數(shù)來獲取泻仙，具體操作見:https://github.com/xianyu426/gene_annotation

自己提供gene2GO文件時糕再，格式應該為：
gene_ID<TAB>GO_ID1, GO_ID2, GO_ID3, ....

三、數(shù)據(jù)導入

library(topGO)
# 讀取gene-to-GO mapping文件
gene2go <- readMapping(file = "gene-to-GO文件") # 這里我用的是物種全部的基因?qū)狦O文件
# 讀取差異基因文件
DEGs <- read.table("差異基因文件", header = TRUE)

# 定義背景基因和感興趣基因
genenames <- names(gene2go)
genelist <- factor(as.integer(genenames %in% DEGs$geneid)) 
# 這里會生成一個factor玉转，有兩個levels：0和1突想，其中1表示感興趣的基因。
names(genelist) <- genenames
GOdata <- new("topGOdata", ontology="MF", allGenes = genelist, 
              annot = annFUN.gene2GO, gene2GO = gene2go)

這樣就定義了一個topGOdata對象。

四猾担、統(tǒng)計檢驗

test.stat <- new("classicCount", testStatistic = GOFisherTest, name = "Fisher test")
resultFisher <- getSigGroups(GOdat, test.stat)
test.stat <- new("elimScore", testStatistic = GOKSTest, name = "Fisher test", cutOff = 0.01)
resultElim <- getSigGroups(GOdata, test.stat)
test.stat <- new("weightCount", testStatistic = GOFisherTest, name="Fisher test", sigRatio = "ratio")
resultWeight <- getSigGroups(GOdata, test.stat)
test.stat <- new("classicScore", testStatistic = GOKSTest, name = "KS tests")
resultKS <- getSigGroups(GOdata, test.stat)
elim.ks <- runTest(GOdata, algorithm = "elim", statistic = "ks")
allRes <- GenTable(GOdat, classic=elim.ks, KS=resultKS, weight = resultWeight,
                   orderBy = "weight", ranksOf = "classic", topNodes =10)
write.table(allRes, file = "sig_GO_result.txt",
            row.name = FALSE, col.names=TRUE)

結(jié)果可以作氣泡富集圖袭灯。

五、顯示結(jié)果

showSigOfNodes(GOdata, score(resultWeight), firstSigNodes = 10, useInfo = "all")