首页 > 科研教程 > KEGG富集分析从未如此简单
2022
06-29

KEGG富集分析从未如此简单

考虑到很多做实验的小伙伴对很多生物信息学概念不是很了解,受实验小白的委托,我给大家写了一个非常简单的工具: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



最后编辑:
作者:萌小白
一个热爱网络的青年!

发布评论

表情