clusterProfiler富集分析

BiocManager::install('clusterProfiler')
BiocManager::install('org.Hs.eg.db')
BiocManager::install('DOSE')
library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)

entrezID <- read.table("11.xls",header=F,sep="\t")
entrezID <- entrezID$V1
BP <- enrichGO(entrezID,"org.Hs.eg.db",ont="BP",keyType = "ENSEMBL",pAdjustMethod = "BH",pvalueCutoff = 0.05,qvalueCutoff = 0.1,readable = T)
dotplot(BP, x = "GeneRatio", color = "p.adjust", showCategory = 20, size = NULL, split = NULL, font.size = 12, title="Dotplot for Gene Ontology Analysis")
write.table(BP, 'go_tmp.txt', sep = '\t', row.names = FALSE, quote = FALSE)
id_with_fc <- read.table("22.xls",header=T,sep="\t")
id_with_fc2 <- id_with_fc[,2]
names(id_with_fc2) <- id_with_fc[,1]
geneList=id_with_fc2
cnetplot(BP,showCategory = 10,foldChange = geneList)
emapplot(BP,showCategory = 20)
cnetplot(BP,showCategory = 50,circular=TRUE,colorEdge=TRUE)
barplot(BP, x = "GeneRatio", color = "p.adjust", showCategory = 20, size = NULL, split = NULL, font.size = 12, title="Dotplot for Gene Ontology Analysis")
#KEGG
gene.df <- bitr(entrezID,fromType="ENSEMBL",toType=c("ENTREZID","ENSEMBL"),OrgDb = org.Dr.eg.db)
BP <- enrichKEGG(gene.df$ENTREZID,organism = 'dre',keyType = "kegg",pAdjustMethod = "BH",pvalueCutoff = 1,qvalueCutoff = 0.9)
#y <- setReadable(BP,OrgDb = org.Mm.eg.db, keyType = "ENTREZID")#将KEGG结果ENTREZID转换成gene symble
barplot(BP, showCategory = 10)
write.table(BP, 'kegg_tmp.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.csv(as.data.frame(kegg_enrich),"KEGG_enrichment.csv",row.names = F)

非模式物种关键步骤是mp_term2gene、mp_term2name两个库文件的获得,也就是term discription和term id的对应关系,enricher函数不局限与做KEGG或GO分析,类似的功能富集都行,比如Mammalian phenotypes (MP terms)

gene <- read.table("ctide_PSG_id.txt")
mp_term2gene <- data.frame(MP$V2,MP$V1)
mp_term2name <- data.frame(MP$V2,MP$V3)
names(mp_term2gene) <- c("MP_term","Gene")
names(mp_term2name) <- c("MP_term","Name")
gene <- as.vector(gene$V1)
mp_enrich <- enricher(gene=gene,pvalueCutoff = 1,qvalueCutoff = 1,pAdjustMethod = "BH",TERM2GENE = mp_term2gene,TERM2NAME = mp_term2name)