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")