我们做转录组分析,得到的结果可能是FPKM。但是FPKM确实存在不准确性,推荐使用TPM。
read count和FPKM结果都可以转成TPM,但是因为FPKM跟TPM的计算都考虑了基因长度,所以从FPKM转TPM最方便快捷。
假设原来的表达矩阵fpkm_expr.txt中行为基因,列为样本,中间数值是FPKM计算得到的值
先读取自己的表达矩阵
expMatrix<-read.table("fpkm_expr.txt",header = T,row.names = 1)
其实早已经有人帮我们整理好了
https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/
countToTpm <- function(counts, effLen)
{
rate <- log(counts) - log(effLen)
denom <- log(sum(exp(rate)))
exp(rate - denom + log(1e6))
}
countToFpkm <- function(counts, effLen)
{
N <- sum(counts)
exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
countToEffCounts <- function(counts, len, effLen)
{
counts * (len / effLen)
}
如果要计算TPM值,只需要用一下apply
函数
tpms <- apply(expMatrix,2,fpkmToTpm)
tpms[1:3,]
最后可以根据TPM的特征进行检查,看每列加和是否一致
colSums(tpms)
转自:刘小泽
网友评论