10、差異基因topGO富集

參考：http://www.biotrainee.com/thread-558-1-1.html

http://bioconductor.org/packages/3.7/bioc/

http://www.bioconductor.org/packages/release/bioc/html/topGO.html

https://www.jianshu.com/p/9e21f2196178

https://rpubs.com/aemoore62/TopGo_colMap_Func_Troubleshoot

構建topGOdata對象的3個數據

1. 基因某種ID的列表（可以有另一個對應的分數值，如p值或t統計量，或者是差異表達值）
2. 基因的這種ID與GO的映射表，在ID為芯片的探針ID時，可以直接使用bioconductor的芯片注釋包如hgu95av2.db包
3. GO的層次關系數據，這個結果可以從GO.db包獲得，topGO也只支持GO.db包定義的層次結

topGO使用：

首先，我們制作准備文件，CC、BP、MF三個注釋文件，格式為：基因ID\t GO:xxx,GO:yyy（topGO.map）

然后准備我們的差異基因列表，兩列差異基因ID\t FDR。（topGO.list）

 

library("topGO")
geneID2GO<-readMappings(choose.files())　　##讀取所有基因注釋信息
geneNames<- names(geneID2GO)

data<-read.table(choose.files(), row.names = 1, header=TRUE,check.names =F)　　##讀取差異基因的ID
geneList<-data[,1]
names(geneList) <- rownames(data)
topDiffGenes<-function(allScore){return(allScore<0.05)}
1.1、###BP
sampleGOdata <- new("topGOdata",nodeSize = 6,ontology="BP", allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO,geneSel=topDiffGenes)
resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
allRes <- GenTable(sampleGOdata,KS = resultKS.elim,ranksOf = "classic", topNodes = attributes(resultKS.elim)$geneData[4])
write.table(allRes, file="T01_vs_T02.topGO_BP.xls", sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)

pdf("T01_vs_T02.topGO_BP.pdf")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")　　##作圖
dev.off()

png("T01_vs_T02.topGO_BP.png")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()
1.2、##MF
sampleGOdata <- new("topGOdata",nodeSize = 6,ontology="MF", allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO,geneSel=topDiffGenes)
resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
allRes <- GenTable(sampleGOdata,KS = resultKS.elim,ranksOf = "classic", topNodes = attributes(resultKS.elim)$geneData[4])
write.table(allRes, file="T01_vs_T02.topGO_MF.xls", sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
pdf("T01_vs_T02.topGO_MF.pdf")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()

png("T01_vs_T02.topGO_MF.png")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()
1.3、##CC
sampleGOdata <- new("topGOdata",nodeSize = 6,ontology="CC", allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO,geneSel=topDiffGenes)
resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")
allRes <- GenTable(sampleGOdata,KS = resultKS.elim,ranksOf = "classic", topNodes = attributes(resultKS.elim)$geneData[4])
write.table(allRes, file="T01_vs_T02.topGO_CC.xls", sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
pdf("T01_vs_T02.topGO_CC.pdf")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()

png("T01_vs_T02.topGO_CC.png")
showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()

1.4、##輸出每個go的基因以及注釋到這個go的差異基因[可以不做這個]

allGO =usedGO(object = sampleGOdata)

for (gos in allGO){

goID <-gos;

gene.universe <- genes(sampleGOdata);

go.genes <- genesInTerm(sampleGOdata,goID)[[1]];

sig.genes <- sigGenes(sampleGOdata);

file1=paste("GO-TMP_BP_sig_",gos,sep="");

write.table(sig.genes,file=file1);

file2=paste("GO-TMP_BP_go_",gos,sep="");

write.table(go.genes,file=file2);

}

 

 

 

 

2、來自生信技能樹

###BP
sampleGOdata <- new("topGOdata",nodeSize = 6,ontology="BP", allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO,geneSel=topDiffGenes)

allGO =usedGO(object = sampleGOdata)

resultCFis<-runTest(sampleGOdata,algorithm="classic",statistic="fisher")

gtFis<-GenTable(sampleGOdata,classicFisher=resultCFis,orderBy="classic",ranksOf="classicFisher",topNodes=length(allGO))

fdr<-p.adjust(p=gtFis[,"classicFisher"],method="fdr")

r <-cbind(gtFis,fdr)

write.table(r,file="topGO_BP.xls",sep="\t")

showSigOfNodes(GOdata,score(resultCFis),firstSigNodes= 5, useInfo = "all")