考虑到很多做实验的小伙伴对很多生物信息学概念不是很了解，受实验小白的委托，我给大家写了一个非常简单的工具：KEGG富集分析

KEGG是干嘛的捏？

我这么跟你说吧：人类的七千多个基因组都是有已知功能的，KEGG把这七千多个基因分成了300个类，就是我们通常说的kegg通路；比如，我现在做了个实验，发现某细胞系里面的两万个基因里面有300个基因变化了，那这300个基因会涉及到KEGG数据库的哪几个通路？这时候就需要用到我的工具啦~将这300个基因加入工具里面，得出结果：有30个Cell
cycle通路。

就这么简单？可是富集分析怎么分析呐？

接下来还会用到我这个小工具：事实上，Cell Cycle KEGG 通路
hsa04110只有124个基因；而我处理了细胞系之后，在只有300个基因发生统计学显著变化的情况下，就有30个是Cell
Cycle通路？高达10%的概率，那到底这个Cell Cycle通路是不是被显著改变了呢？
首先，我会把用户的300个基因，都用KEGG数据库的300个通路注释，然后一个个通路循环做超几何分布检验，给出P值；比如刚才的Cell
Cycle 通路就很显著，因为124/7000就2%的概率，结果我
30/300有10%的概率~太可怕了，所以我的这个处理显著的改变了细胞系的Cell Cycle 通路。

使用方法其实没什么好说的，就是复制粘贴你感兴趣的基因到输入框即可。我这里主要讲解，这个网页如何写出来的。

一、代码编写

UI界面的代码：

1. 
   

   
   suppressPackageStartupMessages(library(shiny))
   

   
   


2. 
   

   
   suppressPackageStartupMessages(library(ggplot2))
   

   
   


3. 
   

   
   suppressPackageStartupMessages(library(DT))
   

   
   


4. 
   

   
   suppressPackageStartupMessages(library(stringr) )
   

   
   


5. 
   

   
   suppressPackageStartupMessages(library(clusterProfiler))
   

   
   


6. 
   

   
   suppressPackageStartupMessages(library(shinyjs))
   

   
   


7. 
   

   
   suppressPackageStartupMessages(library(shinyBS))
   

   
   


8. 
   

   
   ## 这个是侧边栏，就是你看到的网页坐标的输入框
   

   
   


9. 
   

   
   ## 分别是 下拉框，文字输入框，多选按钮，确定按钮。
   

   
   


10. 
    

    
    sp <- sidebarPanel(
    

    
    


11. 
    

    
    selectInput("IDtype",
    

    
    


12. 
    

    
    label = "choose the ID type",
    

    
    


13. 
    

    
    choices = c("HUGO Gene Symbol"= "symbols",
    

    
    


14. 
    

    
    "Entrez Gene ID"= "geneIds",
    

    
    


15. 
    

    
    "ENSEMBL Gene ID"= "geneEnsembl")),
    

    
    


16. 
    

    
    ## "symbols" "geneIds" "geneNames" "geneEnsembl" "geneMap" "geneAlias"
    

    
    


17. 
    

    
    h3("Type the Gene List:"),
    

    
    


18. 
    

    
    tags$style(type="text/css", "textarea {width:100%}"),
    

    
    


19. 
    

    
    tags$textarea(id = 'input_text',
    

    
    


20. 
    

    
    placeholder = "TP53nCBX6",
    

    
    


21. 
    

    
    rows = 8, ''),
    

    
    


22. 
    

    
    hr(),
    

    
    


23. 
    

    
    radioButtons("species", "Select a species:",
    

    
    


24. 
    

    
    c("human(Homo sapiens)"="human",
    

    
    


25. 
    

    
    "mouse(Mus Musculus)"="mouse"),
    

    
    


26. 
    

    
    selected = 'human'
    

    
    


27. 
    

    
    ),
    

    
    


28. 
    

    
    hr(),
    

    
    


29. 
    

    
    bsAlert("alert_ui_anchorId"),
    

    
    


30. 
    

    
    actionButton("do", "analysis", icon("paper-plane"),
    

    
    


31. 
    

    
    style="color: #fff; background-color: #337ab7; border-color: #2e6da4"
    

    
    


32. 
    

    
    )
    

    
    


33. 
    

    
    )
    

    
    


34. 
    

    
    ## 这个是主栏，就显示两个表格即可。
    

    
    


35. 
    

    
    mp <- mainPanel(
    

    
    


36. 
    

    
    h3('KEGG enrichment:'),
    

    
    


37. 
    

    
    DT::dataTableOutput('KEGG_df'),
    

    
    


38. 
    

    
    hr() ,
    

    
    


39. 
    

    
    h3('Gene Annotation result:'),
    

    
    


40. 
    

    
    DT::dataTableOutput('gene_df'),
    

    
    


41. 
    

    
    hr()
    

    
    


42. 
    

    
    )
    

    
    


43. 
    

    
    ## 主程序入口，是一个侧边栏格式的网页，所以需要定义侧栏和主栏！
    

    
    


44. 
    

    
    shinyUI(
    

    
    


45. 
    

    
    fluidPage(
    

    
    


46. 
    

    
    titlePanel('KEGG enrichment'),
    

    
    


47. 
    

    
    sidebarLayout(
    

    
    


48. 
    

    
    sidebarPanel = sp,
    

    
    


49. 
    

    
    mainPanel = mp
    

    
    


50. 
    

    
    )
    

    
    


51. 
    

    
    )
    

    
    


52. 
    

    
    )
    

    
    

server服务端代码

1. 
   

   
   suppressPackageStartupMessages(library(shiny))
   

   
   


2. 
   

   
   suppressPackageStartupMessages(library(ggplot2))
   

   
   


3. 
   

   
   suppressPackageStartupMessages(library(DT))
   

   
   


4. 
   

   
   suppressPackageStartupMessages(library(stringr) )
   

   
   


5. 
   

   
   suppressPackageStartupMessages(library(clusterProfiler))
   

   
   


6. 
   

   
   suppressPackageStartupMessages(library(shinyjs))
   

   
   


7. 
   

   
   suppressPackageStartupMessages(library(shinyBS))
   

   
   


8. 
   

   
   suppressMessages(library(RSQLite))
   

   
   


9. 
   

   
   sqlite <- dbDriver("SQLite")
   

   
   


10. 
    

    
    createLink <- function(base,val) {
    

    
    


11. 
    

    
    sprintf('<a href="%s" class="btn btn-link" target="_blank" >%s</a>',base,val) ##target="_blank"
    

    
    


12. 
    

    
    }
    

    
    


13. 
    

    
    log_cat <- function(info='hello world~',file='log.txt'){
    

    
    


14. 
    

    
    cat(as.character(Sys.time()),info ,"n",file=file,append=TRUE)
    

    
    


15. 
    

    
    }
    

    
    


16. 
    

    
    shinyServer(function(input, output,session) {
    

    
    


17. 
    

    
    ## 定义几个全局变量
    

    
    


18. 
    

    
    glob_values <- reactiveValues(
    

    
    


19. 
    

    
    gene_df=NULL,
    

    
    


20. 
    

    
    kegg_df=NULL,
    

    
    


21. 
    

    
    species=NULL
    

    
    


22. 
    

    
    )
    

    
    


23. 
    

    
    reactiveValues.reset <-function(){
    

    
    


24. 
    

    
    glob_values$gene_df=NULL
    

    
    


25. 
    

    
    glob_values$species=NULL
    

    
    


26. 
    

    
    glob_values$kegg_df=NULL
    

    
    


27. 
    

    
    }
    

    
    


28. 
    

    
    ## 主程序入口，用户点击了确定运行程序的按钮，就开始判断文字输入框是否有基因列表。
    

    
    


29. 
    

    
    observeEvent(input$do,{
    

    
    


30. 
    

    
    reactiveValues.reset()
    

    
    


31. 
    

    
    db=ifelse(input$species == 'human','human_gene_info','mouse_gene_info')
    

    
    


32. 
    

    
    gene_list=NULL
    

    
    


33. 
    

    
    ## first check the upload file:
    

    
    


34. 
    

    
    inFile <- input$file1
    

    
    


35. 
    

    
    if( ! is.null(inFile) ){
    

    
    


36. 
    

    
    gene_list=read.table(inFile$datapath, header=input$header)[,1]
    

    
    


37. 
    

    
    }
    

    
    


38. 
    

    
    ## Then check the text input area:
    

    
    


39. 
    

    
    if( nchar(input$input_text) >5){
    

    
    


40. 
    

    
    gene_list = strsplit(input$input_text,'n')[[1]]
    

    
    


41. 
    

    
    }
    

    
    


42. 
    

    
    if( !is.null(gene_list)){
    

    
    


43. 
    

    
    if(length(gene_list) < 20){
    

    
    


44. 
    

    
    createAlert(session, "alert_ui_anchorId", "exampleAlert", title = "Oops",
    

    
    


45. 
    

    
    content = " 你给的基因数量少于20，没啥意思，不给你富集了 ！", append = FALSE)
    

    
    


46. 
    

    
    }elseif(length(gene_list) > 2000){
    

    
    


47. 
    

    
    createAlert(session, "alert_ui_anchorId", "exampleAlert", title = "Oops",
    

    
    


48. 
    

    
    content = " 给我多于2000个基因，我的服务器受不了呀，要不你先捐赠一下呗 ", append = FALSE)
    

    
    


49. 
    

    
    }else{
    

    
    


50. 
    

    
    closeAlert(session, "exampleAlert")
    

    
    


51. 
    

    
    gene_list=unique(gene_list)
    

    
    


52. 
    

    
    con <- dbConnect(sqlite,"hg19_bioconductor.sqlite")
    

    
    


53. 
    

    
    sql <- paste0('select * from ',db," where ",input$IDtype,
    

    
    


54. 
    

    
    " in ('",paste0(gene_list,collapse = "','"),"')")
    

    
    


55. 
    

    
    glob_values$gene_df=dbGetQuery(con, sql)
    

    
    


56. 
    

    
    dbDisconnect(con)
    

    
    


57. 
    

    
    if(T ){ ## for kegg
    

    
    


58. 
    

    
    suppressPackageStartupMessages(library(org.Mm.eg.db))
    

    
    


59. 
    

    
    suppressPackageStartupMessages(library(org.Hs.eg.db))
    

    
    


60. 
    

    
    gene_df = glob_values$gene_df
    

    
    


61. 
    

    
    ############################################################
    

    
    


62. 
    

    
    ############ gene ID transfer #############################
    

    
    


63. 
    

    
    ############################################################
    

    
    


64. 
    

    
    gene_list <- gene_df$symbol
    

    
    


65. 
    

    
    gene.df=''
    

    
    


66. 
    

    
    if(input$species == 'human'){
    

    
    


67. 
    

    
    gene.df <- bitr(gene_list, fromType = "SYMBOL",
    

    
    


68. 
    

    
    toType = c("ENSEMBL", "ENTREZID"),
    

    
    


69. 
    

    
    OrgDb= org.Hs.eg.db )
    

    
    


70. 
    

    
    ego <- enrichKEGG(gene = gene.df$ENTREZID,
    

    
    


71. 
    

    
    organism = 'hsa',
    

    
    


72. 
    

    
    pvalueCutoff = 0.99,
    

    
    


73. 
    

    
    qvalueCutoff=0.99
    

    
    


74. 
    

    
    )
    

    
    


75. 
    

    
    ego <- setReadable(ego, OrgDb= org.Hs.eg.db, keytype="ENTREZID")
    

    
    


76. 
    

    
    }else{
    

    
    


77. 
    

    
    gene.df <- bitr(gene_list, fromType = "SYMBOL",
    

    
    


78. 
    

    
    toType = c("ENSEMBL", "ENTREZID"),
    

    
    


79. 
    

    
    OrgDb= org.Mm.eg.db )
    

    
    


80. 
    

    
    ego <- enrichKEGG(gene = gene.df$ENTREZID,
    

    
    


81. 
    

    
    organism = 'mmu',
    

    
    


82. 
    

    
    pvalueCutoff = 0.99,
    

    
    


83. 
    

    
    qvalueCutoff=0.99
    

    
    


84. 
    

    
    )
    

    
    


85. 
    

    
    ego <- setReadable(ego, OrgDb= org.Mm.eg.db, keytype="ENTREZID")
    

    
    


86. 
    

    
    }
    

    
    


87. 
    

    
    glob_values$kegg_df <- as.data.frame(ego)
    

    
    


88. 
    

    
    } ## for kegg
    

    
    


89. 
    

    
    }} ## if( !is.null(gene_list)){
    

    
    


90. 
    

    
    })
    

    
    


91. 
    

    
    ## 显示第一个表格，基因注释表格
    

    
    


92. 
    

    
    output$gene_df <- DT::renderDataTable({
    

    
    


93. 
    

    
    if(! is.null(glob_values$gene_df)){
    

    
    


94. 
    

    
    dat=glob_values$gene_df
    

    
    


95. 
    

    
    dat$geneIds=createLink(paste0("http://www.ncbi.nlm.nih.gov/gene/",dat$geneIds),dat$geneIds)
    

    
    


96. 
    

    
    dat
    

    
    


97. 
    

    
    }
    

    
    


98. 
    

    
    } ,rownames= FALSE,escape = FALSE,options = list(
    

    
    


99. 
    

    
    pageLength = 10,
    

    
    


100. 
     

     
     lengthMenu = list(c(10, 50, 100,-1), c('10', '50','100', 'All'))
     

     
     


101. 
     

     
     )## endforoptions
     

     
     


102. 
     

     
     )
     

     
     


103. 
     

     
     ## 显示第二个表格，kegg富集结果表格
     

     
     


104. 
     

     
     output$KEGG_df <- DT::renderDataTable({
     

     
     


105. 
     

     
     if(! is.null(glob_values$kegg_df))
     

     
     


106. 
     

     
     glob_values$kegg_df
     

     
     


107. 
     

     
     } ,rownames= FALSE,options = list(
     

     
     


108. 
     

     
     pageLength = 10,
     

     
     


109. 
     

     
     lengthMenu = list(c(10, 50, 100,-1), c('10', '50','100', 'All'))
     

     
     


110. 
     

     
     ,
     

     
     


111. 
     

     
     scrollX = TRUE,
     

     
     


112. 
     

     
     fixedHeader = TRUE,
     

     
     


113. 
     

     
     fixedColumns = TRUE ,
     

     
     


114. 
     

     
     deferRender = TRUE
     

     
     


115. 
     

     
     ),
     

     
     


116. 
     

     
     #filter = 'top',
     

     
     


117. 
     

     
     escape = FALSE
     

     
     


118. 
     

     
     )
     

     
     


119. 
     

     
     })
     

     
     

组合两个程序，在Rstudio里面运行即可

二、安装及运行

1、首先安装R，再安装Rstudio，然后安装shiny等R包。

2、把上面UI端代码拷贝成ui.R 的文件，再把服务端代码拷贝成server.R文件，放在同一个文件夹，命名为yourAPP里面。

3、进入R里面，该文件夹上层目录，用runAPP('yourAPP')来执行即可

直通车寄语

当然，最简单的就是去https://github.com/jmzeng1314/myShiny 我的GitHub里面把代码下载，或者查看其它工具咯

UI界面如下：

点击阅读原文，试用小工具，赶快输入一串基因去调戏一下这个工具吧！

比如：

LGALS1

DUSP10

LOC402644

SOX11

LOC100131940

LSM5

HEATR2

UCP2

ASNS

RNU6-15

RNU6-1

DUSP10

CDCA7L

DBNL

AXL

TXNRD2

HNRNPA2B1

DHRS2