• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      老俊俊的生信笔记

      • 首页
      • 博客
      • 老俊俊的生信笔记
      • shiny CountToTPM/FPKM

      shiny CountToTPM/FPKM

      • 发布者 weinfoadmin
      • 分类 老俊俊的生信笔记
      • 日期 2021年9月10日
      • 评论 0评论

      感谢老俊俊的大力支持。我们会每日跟新,欢迎您关注老俊俊的生信笔记。

      定量

      转录组数据分析上游分析的最后步骤就是对比对的bam文件进行定量,用的比较多的是HTseq和Subread软件里面的featureCounts函数也可以定量。

      下面是它们文章的引用量:


      featureCounts使用方法可参考:转录组定量-featureCounts[1]

      • 基本用法
      featureCounts -T 15  # 使用线程数
                    -t exon  # 对什么feature进行定量(gene、transcript、exon、cds等)
                    -p  # 测序是双端就加上,单端不加,paired的意思
                    -g gene_id  # 根据gtf文件提取meta信息
                    --extraAttributes gene_name  # 把基因名也一起输出,推荐
                    -a gencode.gtf  # 注释文件
                    -o test.count.txt  #输出文件
                    A.bam B.bam C.bam # bam文件
      • 输出文件内容

      counts to TPM/FPKM

      定量count结果得到后一边直接走差异分析,一边将count值转为TPM或者FPKM值进行下游可视化分析,比如像热图heatmap,每次想得到tpm/fpkm值都得写代码,还得把gene_id转为基因名,于是为了解放一下双手,就干脆写个shiny小程序就行了。

      基本思路:

      • 1.直接输入featureCounts输出的原始count文件
      • 2.点击submit一键获得tpm/fpkm标准化值
      • 3.点击download一键下载

      ui部分:

      library(shiny)
      library(DT)
      library(ggsci)
      library(colourpicker)
      library(stringr)
      library(data.table)
      library(shinythemes)
      library(shinycssloaders)
      # 设置上传数据大小,最大100M
      options(shiny.maxRequestSize=100*1024^2)

      ui <- navbarPage(
          theme = shinytheme("simplex"),
          
          title = tags$strong('ZhouLab Normalization'),
          tabPanel(
              'Normalization',icon = icon('balance-scale'),
              
              
              tabsetPanel(
                  tabPanel(
                      'upload',
                      fileInput('data',
                                h4('Upload Count data'),
                                accept = c("text/csv",
                                           "text/comma-separated-values,text/plain",
                                           ".csv"))
                  ),
                  tabPanel('example data',
                           hr(),
                           # 这里可以在www文件夹下放一张上传数据的示例文件图片
                           tags$img(src='count.PNG'))
                  
              ),
              wellPanel(style = "background: #faf0af",
              h4('1、The raw data showed here:',h5('the Chr:Start:End:Strand is too long,not showed here)',style='color:#01a9b4')),
              hr(style="border-color: #29bb89"),
              DT::dataTableOutput('raw_table')),
                       
              hr(),
              h4('2、The normalized data showed here:'),
              hr(),
              actionButton('nor_submit','submit',icon = icon('android')),
              hr(style="border-color: #29bb89"),
              
              fluidRow(column(6,
                              wellPanel(style = "background: #fcdada",
                              h4(tags$strong('TPM table',style='color:#f39189')),
                              withSpinner(type = 7,size = 0.5,DT::dataTableOutput('tpm_table')),
                              hr(),
                              downloadButton('download_tpm_table','Download'))),
                       column(6,
                              wellPanel(style = "background: #99f3bd",
                              h4(tags$strong('FPKM table',style='color:#0278ae')),
                              withSpinner(type = 7,size = 0.5,DT::dataTableOutput('fpkm_table')),
                              hr(),
                              downloadButton('download_fpkm_table','Download')))),
          )
      )

      server部分:

      server <- function(input, output) {

          raw_tab <- reactive({
              infile <- input$data$datapath
              if (is.null(infile))
                  return(NULL)
              
              d <- infile
              type <- str_sub(d,-3)
              
              if(type=='csv'){
                  count <- fread(infile,header = T,sep = ',',skip = 1)
                  }else{
                      count <- fread(infile,header = T,sep = 't',skip = 1)
                  }
          })
          #' @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
          output$raw_table <- renderDataTable(options = list(scrollX = TRUE,
                                                             aLengthMenu = c(5, 10, 15), 
                                                             iDisplayLength = 5),{
              filter <- raw_tab()[,c(-2:-5)]
              filter
          })
          
          # ---------------------------------------------------------------
          nor_tpm_dat <- eventReactive(input$nor_submit,{
              infile <- input$data$datapath
              if (is.null(infile))
                  return(NULL)
              
              d <- infile
              type <- str_sub(d,-3)
              
              if(type=='csv'){
                  counts <- read.table(infile,skip = '#',header = T,row.names = 1,sep = ',')
              }else{
                  counts <- read.table(infile,skip = '#',header = T,row.names = 1,sep = 't')
              }
              
              my_couts <- counts[,-1:-6]
              
              gene_anno=data.frame(gene_id=rownames(counts),gene_name=counts[,6])
              ##############################################################
              #计算tpm
              kb <- counts$Length/1000
              rpk <- my_couts/kb
              tpm <- t(t(rpk)/colSums(rpk)*1000000)
              
              tpm=as.data.frame(tpm)
              tpm$gene_id <- rownames(tpm)
              
              tpm_anno=merge(tpm,gene_anno,by='gene_id')
              # write.csv(tpm_anno,file = 'TPM-anno.csv',row.names = F)
              
          })
          
          # @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
          nor_fpkm_dat <- eventReactive(input$nor_submit,{
              infile <- input$data$datapath
              if (is.null(infile))
                  return(NULL)
              
              d <- infile
              type <- str_sub(d,-3)
              
              if(type=='csv'){
                  counts <- read.table(infile,skip = '#',header = T,row.names = 1,sep = ',')
              }else{
                  counts <- read.table(infile,skip = '#',header = T,row.names = 1,sep = 't')
              }
              
              my_couts <- counts[,-1:-6]
              
              gene_anno=data.frame(gene_id=rownames(counts),gene_name=counts[,6])
              ##############################################################
              #计算rpk
              kb <- counts$Length/1000
              rpk <- my_couts/kb
              
              ##############################################################
              #计算fpkm
              fpkm <- t(t(rpk)/colSums(my_couts)*1000000)
              fpkm=as.data.frame(fpkm)
              fpkm$gene_id <- rownames(fpkm)
              
              fpkm_anno=merge(fpkm,gene_anno,by='gene_id')
              # write.csv(fpkm_anno,file = 'FPKM-anno.csv',row.names = F)
          })
          
          # -------------------------------------fpkm_table tpm_table
          output$tpm_table <- renderDT(options = list(scrollX = TRUE),{
              nor_tpm_dat()
          })
          
          output$fpkm_table <- renderDT(options = list(scrollX = TRUE),{
              nor_fpkm_dat()
          })
          
          #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
          # download_tpm_table download_fpkm_table
          output$download_tpm_table <- downloadHandler(
              filename = function() {
                  paste('Normalized-TPM', '.csv',sep = "")
              },
              
              content = function(file) {
                  write.csv(nor_tpm_dat(),file,row.names = FALSE)
              })
          
          output$download_fpkm_table <- downloadHandler(
              filename = function() {
                  paste('Normalized-FPKM', '.csv',sep = "")
              },
              
              content = function(file) {
                  write.csv(nor_fpkm_dat(),file,row.names = FALSE)
              })
      }

      启动app:

      # Run the application 
      shinyApp(ui = ui, server = server)

      运行App

      示例数据:

      拿个测试数据试一试:


      完成!大功告成。

      欢迎小伙伴留言评论!

      点击我留言!

      今天的分享就到这里了,敬请期待下一篇!

      最后欢迎大家分享转发,您的点赞是对我的鼓励和肯定!

      如果觉得对您帮助很大,打赏一下吧!

      参考资料

      [1]

      转录组定量-featureCounts: https://www.bioinfo-scrounger.com/archives/407/


      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      shiny VennDiagram
      2021年9月10日

      下一篇文章

      自己模仿画个-- m6A distribution on transcript
      2021年9月10日

      你可能也喜欢

      8-1651542331
      跟着Nature学绘图(2) 箱线图-累积分布曲线图
      2 5月, 2022
      9-1651542322
      Julia 笔记之字符串
      2 5月, 2022
      0-1651542343
      Julia 笔记之数学运算和初等函数
      1 5月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?