根據(jù)b站基因課視頻修改舷暮，按照輸入文件名輸出四個富集圖像。

knitr::opts_chunk$set(echo = TRUE)

library(readr)
library(dplyr)
library(clusterProfiler)
library(tidyr)
library(ggplot2)
library(org.hsa.eg.db)
#-------------------------------------------------------------------------------
#文件讀取(只改變輸入文件即可)
#-------------------------------------------------------------------------------
doc_name <- "DESeq2.DE_results" //兩兩比對的DE_results結(jié)果

#準(zhǔn)備注釋文件黎茎，根據(jù)教程生成 //http://www.reibang.com/p/f2e4dbaae719
pathway2gene <- AnnotationDbi::select(org.hsa.db,
                                      keys = keys(org.hsa.db),
                                      columns = c("Pathway","KO")) %>%
  na.omit() %>%
  dplyr::select(Pathway, GID)

#準(zhǔn)備基因列表
DE_results <- read_delim(doc_name, 
                        delim = "\t", escape_double = FALSE, 
                        trim_ws = TRUE)
#篩選前200個差異基因作為kegg富集對象
#-------------------------------------------------------------------------------
DE_results_kegg <- arrange(DE_results,desc(abs(log2FoldChange)))%>%
    dplyr::slice(1:2000)%>%
    dplyr::select(gene,log2FoldChange)

ekp <- enricher(DE_results_kegg$gene[1:200],
                TERM2GENE = pathway2gene,
                TERM2NAME = pathway2name,
                pvalueCutoff = 0.05,  # 表示全部保留，可以設(shè)為0.05作為閾值
                qvalueCutoff = 0.05, # 表示全部保留当悔，可以設(shè)為0.05作為閾值
                pAdjustMethod = "BH",
                minGSSize = 1)
png(filename = paste(doc_name, "_kegg.png", sep = ""), width = 700, height = 700, res = 100)
dotplot(ekp, showCategory = 10) # 只顯示前10個類別
dev.off()
# 基因差異列表
geneList <- DE_results_kegg$log2FoldChange
names(geneList) <- DE_results_kegg$gene
geneList <- sort(geneList, decreasing = T)

png(filename = paste(doc_name, "_kegg_cent.png", sep = ""), width = 2000, height = 2000, res = 200)  
cnetplot(ekp, 
         color.params = list(foldChange = geneList, edge = TRUE),
         showCategory = 3,
         node_label = "all", # category | gene | all | none
         circular = TRUE)
dev.off()

#gos富集分析
#-------------------------------------------------------------------------------
ego <- enrichGO(gene=na.omit(DE_results_kegg$gene[1:200]),
                OrgDb=org.hsa.eg.db,
                keyType="GID",
                ont="ALL",   #CC/BP/MF可選
                qvalueCutoff = 0.05,
                pvalueCutoff =0.05)
png(filename = paste(doc_name, "_go.png", sep = ""), width = 700, height = 700, res = 100)
dotplot(ego, showCategory = 10) # 只顯示前10個類別
dev.off()
#gos富集分析ggplot作圖
go.res <- data.frame(ego)
goBP <- subset(go.res,subset = (ONTOLOGY == "BP"))[1:10,]
goCC <- subset(go.res,subset = (ONTOLOGY == "CC"))[1:10,]
goMF <- subset(go.res,subset = (ONTOLOGY == "MF"))[1:10,]
go.df <- rbind(goBP,goCC,goMF)
go.df$Description <- factor(go.df$Description,levels = rev(go.df$Description))
go_bar <- ggplot(data = go.df, # 繪圖使用的數(shù)據(jù)
  aes(x = Description, y = Count,fill = ONTOLOGY))+ # 橫軸坐標(biāo)及顏色分類填充
      geom_bar(stat = "identity",width = 0.9)+ # 繪制條形圖及寬度設(shè)置
      coord_flip()+theme_bw()+ # 橫縱坐標(biāo)反轉(zhuǎn)及去除背景色
      scale_x_discrete(labels = function(x) str_wrap(x,width = 120))+ # 設(shè)置term名稱過長時換行
      labs(x = "GO terms",y = "GeneNumber",title = "Barplot of sx_jx Enriched GO Terms")+ # 設(shè)置坐標(biāo)軸標(biāo)題及標(biāo)題
      theme(axis.title = element_text(size = 13), # 坐標(biāo)軸標(biāo)題大小
          axis.text = element_text(size = 11), # 坐標(biāo)軸標(biāo)簽大小
          plot.title = element_text(size = 14,hjust = 0.5,face = "bold"), # 標(biāo)題設(shè)置
          legend.title = element_text(size = 13), # 圖例標(biāo)題大小
          legend.text = element_text(size = 11), # 圖例標(biāo)簽大小
          plot.margin = unit(c(0.5,0.5,0.5,0.5),"cm")) # 圖邊距
    ggsave(go_bar,filename = paste(doc_name, "_go_barplot.png", sep = ""),width = 15,height = 7)
#-------------------------------------------------------------------------------