ENTREZID SYMBOL ENSEMBLE 的相互转换
suppressMessages(library(org.Hs.eg.db))# 载入包
keytypes(org.Hs.eg.db) #查看支持对选项
rt=read.table("gene.txt",sep="\t",check.names=F,header=T)
rt<-as.data.frame(rt[,1])
class(rt)
#keys是自己的基因,columns是输出的类型,keytype是输入的类型
gene_list<-select(org.Hs.eg.db, keys=as.character(rt$`rt[, 1]`), columns=c("SYMBOL","ENTREZID"), keytype="ENSEMBL")
gene_list[1:4,1:3]
write.table(gene_list,file="a.txt",sep="\t",quote=F,row.names=F)
1 ID的转换
library("org.Hs.eg.db")
rt=read.table("symbol.txt",sep="\t",check.names=F,header=T)
genes=as.vector(rt[,1])
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
dim(out)
out<-out[(out$entrezID!='NA'),] #删除NA值
dim(out)
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F)
2 GO分析
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]
gene=rt$entrezID
#GO富集分析
kk <- enrichGO(gene = gene,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)
#柱状图
tiff(file="barplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
barplot(kk, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#气泡图
tiff(file="dotplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
dotplot(kk,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#热图
tiff(file="heatplot.tiff",width = 40,height = 20,units ="cm",compression="lzw",bg="white",res=600)
heatplot(kk,showCategory =20, foldChange=cor)
dev.off()
3 KEGG分析
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]
gene=rt$entrezID
#kegg富集分析
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05)
write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F)
#柱状图
tiff(file="barplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
barplot(kk, drop = TRUE, showCategory = 20)
dev.off()
#气泡图
tiff(file="dotplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
dotplot(kk, showCategory = 20)
dev.off()
#热图
tiff(file="heatplot.tiff",width = 25,height = 15,units ="cm",compression="lzw",bg="white",res=600)
heatplot(kk,showCategory =20, foldChange=cor)
dev.off()
#通路图
library("pathview")
keggxls=read.table("KEGG.txt",sep="\t",header=T)
for(i in keggxls$ID){
pv.out <- pathview(gene.data = cor, pathway.id = i, species = "hsa", out.suffix = "pathview")
}
网友评论