定量
转录组数据分析上游分析的最后步骤就是对比对的bam文件进行定量，用的比较多的是HTseq和Subread软件里面的featureCounts函数也可以定量。
下面是它们文章的引用量：
捕获1.PNG
捕获2.PNG
featureCounts使用方法可参考：转录组定量-featureCounts
- 基本用法
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文件
-
输出文件内容
count.PNG
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
捕获3.PNG
示例数据：
捕获4.PNG
拿个测试数据试一试：
捕获5.PNG
捕获6.PNG
完成！
欢迎交流留言评论，关注我的公众号：老俊俊的生信笔记
网友评论