GO，kegg富集柱状图

转载至：https://blog.csdn.net/sinat_30623997/article/details/79250940

target_gene_id <- unique(read.delim("miRNA-gene interactions.txt")$EntrezID)
#富集的基因集

display_number = c(15, 15, 15)
## 使用clusterProfiler进行富集
library(clusterProfiler)
ego_MF <- enrichGO(OrgDb="org.Hs.eg.db",
             gene = target_gene_id,
             pvalueCutoff = 0.05,
             ont = "MF",
             readable=TRUE)
ego_result_MF <- as.data.frame(ego_MF)[1:display_number[1], ]
ego_result_MF <- ego_result_MF[order(ego_result_MF$Count),]
#获取MF的前15个并排序

ego_CC <- enrichGO(OrgDb="org.Hs.eg.db",
                   gene = target_gene_id,
                   pvalueCutoff = 0.05,
                   ont = "CC",
                   readable=TRUE)
ego_result_CC <- as.data.frame(ego_CC)[1:display_number[2], ]
ego_result_CC <- ego_result_CC[order(ego_result_CC$Count),]

ego_BP <- enrichGO(OrgDb="org.Hs.eg.db",
                   gene = target_gene_id,
                   pvalueCutoff = 0.05,
                   ont = "BP",
                   readable=TRUE)
ego_result_BP <- na.omit(as.data.frame(ego_BP)[1:display_number[3], ])
ego_result_BP <- ego_result_BP[order(ego_result_BP$Count),]

go_enrich_df <- data.frame(ID=c(ego_result_BP$ID, ego_result_CC$ID, ego_result_MF$ID),
                           Description=c(ego_result_BP$Description, ego_result_CC$Description, ego_result_MF$Description),
                           GeneNumber=c(ego_result_BP$Count, ego_result_CC$Count, ego_result_MF$Count),
                           type=factor(c(rep("biological process", display_number[1]), rep("cellular component", display_number[2]),
                                         rep("molecular function", display_number[3])), levels=c("molecular function", "cellular component", "biological process")))

#排序
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df)))


labels <- as.factor(rev(go_enrich_df$Description))


#设定颜色
CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")

library(ggplot2)
p <- ggplot(data=go_enrich_df, aes(x=number, y=GeneNumber, fill=type)) +
  geom_bar(stat="identity", width=0.8) + coord_flip() + 
  scale_fill_manual(values = CPCOLS) + theme_bw() + 
  scale_x_discrete(labels=labels) +
  xlab("GO term") + 
  theme(axis.text=element_text(face = "bold", color="gray50")) +
  labs(title = "The Most Enriched GO Terms")

p

pdf("go_enrichment_of_miRNA_targets.pdf")
p
dev.off()

svg("go_enrichment_of_miRNA_targets.svg")
p
dev.off()

###############################################################
##############################kegg#############################
###############################################################

kegg <- enrichKEGG(eg$ENTREZID, organism = 'hsa', keyType = 'kegg', pvalueCutoff = 0.05,pAdjustMethod = 'BH', 
                   minGSSize = 3,maxGSSize = 500,qvalueCutoff = 0.2,use_internal_data = FALSE)

kegg_all_diff <- kegg[order(kegg$Count),]

kegg_all_diff <- data.frame(ID=kegg_all_diff$ID,
                            Description=kegg_all_diff$Description,
                            GeneNumber=kegg_all_diff$Count
)


kegg_all_diff <- kegg_all_diff[(nrow(kegg_all_diff)-14):nrow(kegg_all_diff),]

kegg_all_diff$number <- factor(rev(1:nrow(kegg_all_diff)))

kegg_all_diff$Description <- str_split(kegg_all_diff$Description,'[,]',simplify = T)[,1]

labels <- as.factor(rev(kegg_all_diff$Description))



CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")

library(ggplot2)
p <- ggplot(data=kegg_all_diff, aes(x=number, y=GeneNumber, fill='purple')) +
  geom_bar(stat="identity", width=0.8) + coord_flip() + 
  scale_fill_manual(values = CPCOLS) + theme_bw() + 
  scale_x_discrete(labels=labels) +
  xlab("kegg term") + 
  theme(axis.text=element_text(face = "bold", color="gray50")) +
  labs(title = "The Most Enriched kegg Terms")

p

pdf("kegg_enrichment_all_diff.pdf")
p
dev.off()

svg("kegg_enrichment.svg")
p
dev.off()