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)
Leave a Reply