上下调基因各自独立进行GO数据库的3分类富集(求美图代码)

## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
if(T){
  ###   over-representation test
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      universe     = gene_all,
                      pvalueCutoff = 0.9,#根据实际情况调整
                      qvalueCutoff =0.9)#根据实际情况调整
  head(kk.up)[,1:6]
  dotplot(kk.up );ggsave('kk.up.dotplot.png')
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        universe     = gene_all,
                        pvalueCutoff = 0.9,#根据实际情况调整
                        qvalueCutoff =0.9)#根据实际情况调整
  head(kk.down)[,1:6]
  dotplot(kk.down );ggsave('kk.down.dotplot.png')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.05)#根据实际情况调整
  head(kk.diff)[,1:6]
  dotplot(kk.diff );ggsave('kk.diff.dotplot.png')
  
  kegg_diff_dt <- as.data.frame(kk.diff)
  kegg_down_dt <- as.data.frame(kk.down)
  kegg_up_dt <- as.data.frame(kk.up)
  down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group=-1
  up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group=1
  ## 原来的画图函数
  source('functions.R') 
  g_kegg=kegg_plot(up_kegg,down_kegg)
  print(g_kegg)
  
  ggsave(g_kegg,filename = 'kegg_up_down.png')
  
  ###  GSEA 
  kk_gse <- gseKEGG(geneList     = geneList,
                    organism     = 'hsa',
                    nPerm        = 1000,#根据实际情况调整
                    minGSSize    = 10,#根据实际情况调整
                    pvalueCutoff = 0.9,##根据实际情况调整
                    verbose      = FALSE)
  head(kk_gse)[,1:6]
  gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
  
  down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
  up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
  
  g_kegg=kegg_plot(up_kegg,down_kegg)
  print(g_kegg)
  ggsave(g_kegg,filename = 'kegg_up_down_gsea.png')
  
  
}

### Update: Xiaojie Liang
### Date: 2021-07-12 23:38:07
### Email: liangxiaojiecom@163.com
### ---------------
library(ggthemes)
library(ggplot2)
go.kegg_plot <- function(up.data,down.data){
  # up.data <- up.data
  # down.data <- down.data
  dat=rbind(up.data,down.data)
  colnames(dat)
  dat$p.adjust = -log10(dat$p.adjust)
  dat$p.adjust=dat$p.adjust*dat$group 
  
  dat=dat[order(dat$p.adjust,decreasing = F),]
  
  gk_plot <- ggplot(dat,aes(reorder(Description, p.adjust), y=p.adjust)) +
    geom_bar(aes(fill=factor((p.adjust>0)+1)),stat="identity", width=0.7, position=position_dodge(0.7)) +
    coord_flip() +
    scale_fill_manual(values=c("#0072B2", "#B20072"), guide=FALSE) +
    labs(x="", y="" ) +
    theme_pander()  +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          #axis.ticks.x = element_blank(),
          axis.line.x = element_line(size = 0.3, colour = "black"),#x轴连线
          axis.ticks.length.x = unit(-0.20, "cm"),#修改x轴刻度的高度，负号表示向上
          axis.text.x = element_text(margin = margin(t = 0.3, unit = "cm")),##线与数字不要重叠
          axis.ticks.x = element_line(colour = "black",size = 0.3) ,#修改x轴刻度的线                         
          axis.ticks.y = element_blank(),
          axis.text.y  = element_text(hjust=0),
          panel.background = element_rect(fill=NULL, colour = 'white')
    )
}

library(org.Hs.eg.db)
#把SYMBOL转换ENTREZID，这里会损失一部分无法匹配到的
gene_up=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
                                    keys = gene_up,
                                    columns = 'ENTREZID',
                                    keytype = 'SYMBOL')[,2]))
gene_down=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
                                      keys = gene_down,
                                      columns = 'ENTREZID',
                                      keytype = 'SYMBOL')[,2]))

library(clusterProfiler) 
#对正相关基因进行富集分析
go <- enrichGO(gene_up, OrgDb = "org.Hs.eg.db", ont="all") 
library(ggplot2)
library(stringr)
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")+ 
  ggsave('gene_up_GO_all_barplot.png') 
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_up_GO_all_barplot.csv')

#对负相关基因进行富集分析
go <- enrichGO(gene_down, OrgDb = "org.Hs.eg.db", ont="all") 
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")+ 
  ggsave('gene_down_GO_all_barplot.png') 
go=DOSE::setReadable(go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv( go@result,file = 'gene_down_GO_all_barplot.csv')

GSM2560200 plasmaRNA_health H13
GSM2560201 plasmaRNA_health H5
GSM2560202 plasmaRNA_health H8
GSM2560203 plasmaRNA_OSCC X474
GSM2560204 plasmaRNA_OSCC X574
GSM2560205 plasmaRNA_OSCC X602