劉小澤寫于19.2.21
或許這還是個比較常用的需求
我們做轉(zhuǎn)錄組分析撬码，得到的結(jié)果可能是FPKM。但是FPKM確實存在不準確性版保，推薦使用TPM呜笑。
read count和FPKM結(jié)果都可以轉(zhuǎn)成TPM，但是因為FPKM跟TPM的計算都考慮了基因長度彻犁，所以從FPKM轉(zhuǎn)TPM最方便快捷叫胁。
假設(shè)原來的表達矩陣fpkm_expr.txt中行為基因，列為樣本汞幢，中間數(shù)值是FPKM計算得到的值
先讀取自己的表達矩陣
expMatrix<-read.table("fpkm_expr.txt",header = T,row.names = 1)
其實早已經(jīng)有人幫我們整理好了
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函數(shù)
tpms <- apply(expMatrix,2,fpkmToTpm)
tpms[1:3,]
最后可以根據(jù)TPM的特征進行檢查，看每列加和是否一致
colSums(tpms)
