专栏首页生物信息云RNA-Seq的Counts和FPKM数据如何转换成TPM?

RNA-Seq的Counts和FPKM数据如何转换成TPM?

我们做转录组分析,得到的数据通常是raw counts 的数据,raw counts 的数据有很多R包进行归一化。在TCGA数据库中下载的RNA-Seq的数据就有2种形式,raw counts 和FPKM,尽管有很多文章是直接利用FPKM进行分析的,但是FPKM存在不准确性,通常我们会使用TPM。关于什么是FPKM?什么是TPM?我在前面的文章中就有介绍:RNA-seq的counts,RPM, RPKM, FPK值到底有什么区别?。如果从原始的下机数据开始分析,那就根据自己需要进行转换,但通常我们大多数拿到的是raw counts数据,一般送测序,也会要求返回raw counts的数据,从数据库下载的数据我们通常也是选择raw counts数据或者FPKM的数据。那么我们如何将这些数据进行转换成TPM的数据呢?read count和FPKM结果都可以转成TPM,但是因为FPKM跟TPM的计算都考虑了基因长度,所以从FPKM转TPM最方便快捷。只需要按照下面公式就可以计算:

具体可参考前面的文章:RNA-seq的counts,RPM, RPKM, FPK值到底有什么区别?,这里提供的是R代码。

首先我们得有FPKM的数据,这里我以之前TCGA数据库的数据为例。数据可在文章【TCGA数据库33个Project的RNA-Seq转录组数据为你整理打包好了】中下载。

load("F:/TCGA/HTSeq-FPKM/Rdata/data/TCGA-COAD-Exp.Rdata")
exp <- transomeData[["proteinCodingExpData"]][["Exp"]]

之前上传的数据是字符串,需要转换成数值。

library(dplyr)
exp <- exp %>% data.matrix() %>% as.data.frame()
head(exp)[,1:2]

首先,我们定义个函数,也就是上面的公式。

FPKM2TPM <- function(fpkm){
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

然后我们利用apply函数进行遍历,就可以转换啦。

TPMs <- apply(exp,2,FPKM2TPM)

除了FPKM转换成TPM外,其他的数据也可以进行转换。

  • Counts转TPM
Counts2TPM <- function(counts, effLen){
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}
  • Counts转FPKM
Counts2FPKM <- function(counts, effLen){
  N <- sum(counts)
  exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
  • Counts转Effective counts
Counts2EffCounts <- function(counts, len, effLen){
  counts * (len / effLen)
}

参考:https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-07-25

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言基础教程——第3章:数据结构——数组

    数组(array)与矩阵类似,但是维度可以大于2。数组可通过array函数创建,形式如下:

    DoubleHelix
  • 转录组分析 | 使用RSeQC软件对生成的BAM文件进行质控

    RSeQC是发表于2012年的一个RNA-Seq质控工具,属于python包。它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据,比如一些基本...

    DoubleHelix
  • PanGPCR | 预测多个潜在的GPCR靶标及其在组织中的表达位置,副作用以及GPCR药物的可能用途

    靶向G蛋白偶联受体(GPCR)(已知的最大治疗靶标)的药物发现具有挑战性。为了促进GPCR药物的快速发现和开发,Yufeng J Tseng等人构建了PanGP...

    DoubleHelix
  • 4 个有效提升 Jupyter Notebooks 效果的非凡技巧

    链接 | https://towardsdatascience.com/4-awesome-tips-for-enhancing-jupyter-noteboo...

    AI研习社
  • 聊聊flink的ParameterTool

    flink-core-1.7.1-sources.jar!/org/apache/flink/api/common/ExecutionConfig.java

    codecraft
  • 聊聊flink的ParameterTool

    flink-core-1.7.1-sources.jar!/org/apache/flink/api/common/ExecutionConfig.java

    codecraft
  • MySQL 每秒 570000 的写入,如何实现?

    一个朋友接到一个需求,从大数据平台收到一个数据写入在20亿+,需要快速地加载到MySQL中,供第二天业务展示使用。

    数据和云
  • SDK、API和OPEN API有啥区别,这是最为形象的比喻

    SDK 就是 Software Development Kit 的缩写,中文意思就是“软件开发工具包”。

    IT大咖说
  • MySQL 每秒 570000 的写入,如何实现?

    一个朋友接到一个需求,从大数据平台收到一个数据写入在20亿+,需要快速地加载到MySQL中,供第二天业务展示使用。

    芋道源码
  • MySQL 每秒 570000 的写入,如何实现?

    一个朋友接到一个需求,从大数据平台收到一个数据写入在20亿+,需要快速地加载到MySQL中,供第二天业务展示使用。

    Java团长

扫码关注云+社区

领取腾讯云代金券