蝴蝶图代码

library(clusterProfiler)
library(org.Hs.eg.db)
#####  筛选条件设置 #######
log2FC_cutoff = 0.5
#pvalue_cutoff = 0.005
padj_cutoff = 0.05

####获取DEG结果的上下调基因 ####
#DEG_DESeq2[,c(2,5,6)]  DEG_limma_voom[,c(1,4,5)] DEG_edgeR[,c(1,4,5)] 
need_DEG <- nrDEG[,c(1,4,5)]  
head(need_DEG)
colnames(need_DEG) <- c('log2FoldChange','pvalue','padj')

gene_up=rownames(need_DEG[with(need_DEG,log2FoldChange>log2FC_cutoff  & padj<padj_cutoff),])
gene_down=rownames(need_DEG[with(need_DEG,log2FoldChange < -log2FC_cutoff & padj<padj_cutoff),])

#### 转化基因名为entrez ID ####
#org.Hs.eg.db\org.Mm.eg.db包含着各大主流数据库的数据，如entrez ID和ensembl,
#keytypes(org.Hs.eg.db) #查看所有支持及可转化类型 常用 "ENTREZID","ENSEMBL","SYMBOL"
gene_up_entrez <- as.character(na.omit(bitr(gene_up, #数据集
                                            fromType="SYMBOL", #输入格式
                                            toType="ENTREZID", # 转为ENTERZID格式
                                            OrgDb="org.Hs.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"
gene_down_entrez <- as.character(na.omit(bitr(gene_down, #数据集
                                              fromType="SYMBOL", #输入格式
                                              toType="ENTREZID", # 转为ENTERZID格式
                                              OrgDb="org.Hs.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"
gene_diff_entrez <- unique(c(gene_up_entrez ,gene_down_entrez ))
##或者使用mapIds
# gene_up_ENTREZID <- as.character(na.omit(mapIds(x = org.Mm.eg.db,
#                                                  keys =  gene_up,
#                                                  keytype = "SYMBOL",
#                                                  column = "ENTREZID")))


#### KEGG、GO富集 ####
kegg_enrich_results_up <- enrichKEGG(gene  = gene_up_entrez,
                                  organism  = "hsa", #物种人类 hsa 小鼠mmu
                                  pvalueCutoff = 0.05,
                                  qvalueCutoff = 0.2)
kegg_enrich_results_up <- DOSE::setReadable(kegg_enrich_results_up, 
                                         OrgDb="org.Hs.eg.db", 
                                         keyType='ENTREZID')#ENTREZID to gene Symbol
write.csv(kegg_enrich_results_up@result,'KEGG_gene_up_enrichresults.csv') 
save(kegg_enrich_results, file ='KEGG_gene_up_enrichresults.Rdata')

kegg_enrich_results_down <- enrichKEGG(gene  = gene_down_entrez,
                                     organism  = "hsa", #物种人类 hsa 小鼠mmu
                                     pvalueCutoff = 0.05,
                                     qvalueCutoff = 0.2)
kegg_enrich_results_down <- DOSE::setReadable(kegg_enrich_results_down, 
                                         OrgDb="org.Hs.eg.db", 
                                         keyType='ENTREZID')#ENTREZID to gene Symbol
write.csv(kegg_enrich_results_down@result,'KEGG_gene_down_enrichresults.csv')


go_enrich_results_up <- enrichGO(gene = gene_up_entrez,
                              OrgDb = "org.Hs.eg.db",
                              ont   = "ALL"  ,     #One of "BP", "MF"  "CC"  "ALL" 
                              pvalueCutoff  = 0.05,
                              qvalueCutoff  = 0.2,
                              readable      = TRUE)
write.csv(go_enrich_results_up@result, 'GO_gene_up_BP_enrichresults.csv') 
go_enrich_results_down <- enrichGO(gene = gene_down_entrez,
                                 OrgDb = "org.Hs.eg.db",
                                 ont   = "ALL"  ,     #One of "BP", "MF"  "CC"  "ALL" 
                                 pvalueCutoff  = 0.05,
                                 qvalueCutoff  = 0.2,
                                 readable      = TRUE)
write.csv(go_enrich_results_down@result, 'GO_gene_down_BP_enrichresults.csv') 
save(go_enrich_results, file ='GO_gene_up_enrichresults.Rdata')




cnetplot(kegg_enrich_results, categorySize="pvalue", 
         circular = TRUE, colorEdge = TRUE, 
         node_label = "category")




barplot(kegg_enrich_results, showCategory = 10,color = "pvalue")
dev.off()

up=kegg_enrich_results_up@result
down=kegg_enrich_results_down@result


up=up[order(up$pvalue,decreasing = F),]# 
up$Description=factor(up$Description,levels = up$Description)
id_up=levels(up$Description)
up=up[1:10,]
up=up[order(up$pvalue,decreasing = T),]
up$Description=factor(up$Description,levels = up$Description)
id_up=levels(up$Description)



down=down[order(down$pvalue,decreasing = F),]
down$Description=factor(down$Description,levels = down$Description)
id_down=levels(down$Description)
down=down[1:10,]
down=down[order(down$pvalue,decreasing = T),]
down$Description=factor(down$Description,levels = down$Description)
id_down=levels(down$Description)


##ggplot2 作图
library(ggplot2)

#上调 GO
p_up <- ggplot(up, aes( Description, log(pvalue, 10))) +
  geom_col(fill = 'red2', color = 'black', width = 0.6) +
  theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent')) +
  theme(axis.line.x = element_line(colour = 'black'), axis.line.y = element_line(colour = 'transparent'), axis.ticks.y = element_line(colour = 'transparent')) +
  theme(plot.title = element_text(hjust = 0.5, face = 'plain')) +
  coord_flip() +
  geom_hline(yintercept = 0) +
  labs(x = '', y = '', title = 'UP') +
  scale_y_continuous(expand = c(0, 0), breaks = c(-12, -8, -4, 0), labels = as.character(abs(c(-12, -8, -4, 0)))) +     #这儿更改间距设置
  scale_x_discrete(labels = id_up)

#下调 GO
p_down <- ggplot(down, aes(Description, -log(pvalue, 10))) +
  geom_col(fill = 'green4', color = 'black', width = 0.6) +
  theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent')) +
  theme(axis.line.x = element_line(colour = 'black'), axis.line.y = element_line(colour = 'transparent'), axis.ticks.y = element_line(colour = 'transparent')) +
  theme(plot.title = element_text(hjust = 0.5, face = 'plain')) +
  geom_hline(yintercept = 0) +
  coord_flip() +
  labs(x = '', y = '', title = 'DOWN') +
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 2, 4, 6), labels = as.character(c(0, 2, 4, 6))) +       #这儿更改间距设置
  scale_x_discrete(labels = id_down, position = 'top')


library(cowplot)

pdf('butterfly_kegg.pdf', width = 14, height = 5)
plot_grid(p_up, p_down, nrow = 2, ncol = 2, rel_heights = c(9, 1), labels = 'KEGG_Enrichment Score (-log10(p-value))', label_x = 0.5, label_y = 0, label_fontface = 'plain')
dev.off()


up_go= go_enrich_results_up@result
down_go=go_enrich_results_down@result

up_go=up_go[order(up_go$pvalue,decreasing = F),]#这里做的是一个p值的降序排列，因为p值肯定越小越好
up_go$Description=factor(up_go$Description,levels = up_go$Description)#将ID转化为因子，因为画图用到的是因子
id_up_go=levels(up_go$Description)
up_go=up_go[1:10,]#选取前10个ID
#up_go=up_go[order(up_go$pvalue,decreasing = T),]#这个不用，那时候想错了
#up_go$Description=factor(up_go$Description,levels = up_go$Description)
#id_up_go=levels(up_go$Description)



down_go=down_go[order(down_go$pvalue,decreasing = F),]
down_go$Description=factor(down_go$Description,levels = down_go$Description)
id_down_go=levels(down_go$Description)
down_go=down_go[1:10,]
#down_go=down_go[order(down_go$pvalue,decreasing = T),]
#down_go$Description=factor(down_go$Description,levels = down_go$Description)
#id_down_go=levels(down_go$Description)

##ggplot2 作图
library(ggplot2)

#上调 GO
p_up <- ggplot(up_go, aes( Description, log(pvalue, 10))) +
  geom_col(fill = 'red2', color = 'black', width = 0.6)# fill后面可以改颜色，其他的就不用动了
  +
  theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent')) +
  theme(axis.line.x = element_line(colour = 'black'), axis.line.y = element_line(colour = 'transparent'), axis.ticks.y = element_line(colour = 'transparent')) +
  theme(plot.title = element_text(hjust = 0.5, face = 'plain')) +
  coord_flip() +
  geom_hline(yintercept = 0) +
  labs(x = '', y = '', title = 'UP') +
  scale_y_continuous(expand = c(0, 0), breaks = c(-12, -8, -4, 0), labels = as.character(abs(c(-12, -8, -4, 0)))) +     #这儿更改间距设置
  scale_x_discrete(labels = id_up_go)

#下调 GO
p_down <- ggplot(down_go, aes(Description, -log(pvalue, 10))) +
  geom_col(fill = 'green4', color = 'black', width = 0.6) +
  theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent')) +
  theme(axis.line.x = element_line(colour = 'black'), axis.line.y = element_line(colour = 'transparent'), axis.ticks.y = element_line(colour = 'transparent')) +
  theme(plot.title = element_text(hjust = 0.5, face = 'plain')) +
  geom_hline(yintercept = 0) +
  coord_flip() +
  labs(x = '', y = '', title = 'DOWN') +
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 2, 4, 6), labels = as.character(c(0, 2, 4, 6))) +       #这儿更改间距设置
  scale_x_discrete(labels = id_down_go, position = 'top')


library(cowplot)

pdf('butterfly_go.pdf', width = 14, height = 5)
plot_grid(p_up, p_down, nrow = 2, ncol = 2, rel_heights = c(9, 1), labels = 'GO_Enrichment Score (-log10(p-value))', label_x = 0.5, label_y = 0, label_fontface = 'plain')
dev.off()