RNA 数据下机后,如果处理成read counts matrix的话,是一定要进行基于基因长度的标准化的(TMP,RPKM,TPKM等)。目前最常用的是TPM,网上已经有很多关于这三个标准的计算方法了,在此不赘述,主要说一下这几个数据的计算公式和相互转换。
前提知识点
RPKM, FPKM, TPM区别
www.plob.org
图标
计算公式
FPKM、RPKM
Reads Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的reads)
FPKM: Fragments Per Kilobase of exon model per Million mapped fragments(每千个碱基的转录每百万映射读取的fragments)。与RPKM计算过程类似。只有一点差异:RPKM计算的是reads,FPKM计算的是fragments。single-end/paired-end测序数据均可计算reads count,fragments count只能通过paired-end测序数据计算。paired-end测序数据时,两端的reads比对到相同区域,且方向相反,即计数1个fragments;如果只有单端reads比对到该区域,则一个reads即计数1个fragments。所以fragments count接近且小于2 * reads count
参考:https://www.jianshu.com/p/c25e84383ae3
TPM
Transcripts Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts)
i为比对到第i个exon的reads数; Li为第i个exon的长度;sum()为所有 (n个)exon按长度进行标准化之后数值的和
CPM
RPM/CPM: Reads/Counts of exon model per Million mapped reads (每百万映射读取的reads),多进行样本间比较,无法进行样本内差异表达分析
相互转换代码
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))
}
################################################################################
# An example
################################################################################
# count convert
cnts <- c(4250, 3300, 200, 1750, 50, 0)
lens <- c(900, 1020, 2000, 770, 3000, 1777)
countDf <- data.frame(count = cnts, length = lens)
## assume a mean(FLD) = 203.7
#countDf$effLength <- countDf$length - 203.7 + 1
countDf$effLength=countDf$length
countDf$tpm <- with(countDf, countToTpm(count, effLength))
countDf$fpkm <- with(countDf, countToFpkm(count, effLength))
countDf (INPUT FORMAT)
本来还有一个effect length 的计算,即校正实验误差后的序列长度,同时由effect length 产生effect counts,为了方便理解,此处把原始数据当成effect并进行后续计算,详细见下方英文文章说明
结果输出
fpkmToTpm
expMatrix<-read.table("fpkm_expr.txt",header = T,row.names = 1)
tpms <- apply(expMatrix,2,fpkmToTpm)
tpms[1:3,]
最后可以根据TPM的特征进行检查,看每列加和是否一致
colSums(tpms)
fpkm_expr.txt:行为基因,列为样本,中间数值是FPKM计算得到的值
转换后的TPM
超全英文版参考资料:
https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/
haroldpimentel.wordpress.com
本篇多为资料整合,方便自行查看的,感谢大神们分享指教,侵删,谢谢~