goseq是一个R包,用于寻找GO terms,即基因富集分析。
GO terms是标准化描述基因或基因产物的词汇,包括三方面,cellular component,molecular funciton,biological process。
每个GO term都有一个GO ID,比如 GO:006260,每个GO term背后都有一系列的相关基因。
GO分析的目的:在差异性基因分析后,我们可能得到很多差异基因,这些基因里的一部分可能跟某个生物过程相关,或几个生物过程相关。经过GO分析后,我们就能将差异性基因具体的生物功能展示出来,为下一步研究做准备。
GOseq需要输入的文件:
1.所有有count的genes。
2.差异性表达的genes。
3.genome信息,基因长度信息。#对于许多模式基因组来说,这些内容都被做成了独立的R包。
4.GO terms包。
>source("http://bioconductor.org/biocLite.R") >biocLite("goseq") >biocLite("geneLenDataBase")#genome,genes信息 >biocLite("org.Dm.eg.db")#果蝇的GO terms >library("goseq") >library("geneLenDataBase") >library("org.Dm.eg.db") >DEG<-read.table("DEG",header=FALSE) >ALL<-read.table("ALL",header=FALSE) #DEG:差异性基因表 ALL:所有基因表(数据框格式) >DEG.vector<-c(t(DEG)) >ALL.vector<-c(t(ALL)) #把数据格式转化为vector,便于下步操作 >gene.vector=as.integer(ALL.vector%in%DEG.vector) #生成二进制的gene vector(1代表差异性基因,0代表非差异性基因) >names(gene.vector)<-ALL.vector >pwf=nullp(gene.vector,"dm3","ensGene") #生成probability weighting function. >GO.wall=goseq(pwf,"dm3","ensGene") #生成GO terms ID 。这边的疑问:genes 没有mapping 到GO categories。 goseq函数有一个选项:gene2cat,如果gene2cat=NULL,则goseq会自动调用getgo函数实现mapping功能,并将输出值gene2cat。 >enriched.GO=GO.wall$category[GO.wall$over_represented_pvalue<.05] #生成差异性 GO terms ID >library(GO.db) >capture.output(for(go in enriched.GO[1:length(enriched.GO)]){ print(GOTERM[go]) cat("___________") } ,file="SigGo.txt") #生成具体的GO TERM详解
时间: 2024-11-08 19:51:15