码迷,mamicode.com
首页 > 其他好文 > 详细

10、差异基因topGO富集

时间:2018-08-19 20:57:16      阅读:1158      评论:0      收藏:0      [点我收藏+]

标签:stat   pack   mil   family   info   run   not   release   rar   

参考: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

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)}
###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()
##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()
##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()

10、差异基因topGO富集

标签:stat   pack   mil   family   info   run   not   release   rar   

原文地址:https://www.cnblogs.com/renping/p/9502514.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!