PS：如果你需要本教程的练习代码和文档，可以在公众号回复“20220122”即可获得。
前言：
今早看到一篇博文，提到了FPKM与TPM间转化。我自己也系统的再次进行整理一下（PS：自己前期的基础不是很牢固，基本只是使用Count和FPKM，其它的表达丰度基本没涉及）。因此，本教程博文也是自己的一个学习笔记。也是结合很多篇博文组装出来的，记录记录！！！
---- Du
•转录组使用hisat2比对后，我们会使用featureCounts、HTseq-count等软件计算每个基因Count值（每个基因比对上的reads数），count值是最原始的，也是最接近真实的基因表达情况，是没被标准化的数值，因此，很多的差异表达分析，输入文件（input data)使用Count值。以及，后面所有的FPK、RPKM、TPM等都是依据Count值转换出来的。
•计算FPKM值，可以根据Count值进行计算，此步需要我们后期自己计算，但也是使用Stringtie软件进行计算。该软件也可以使用其脚本prepDE.py进行转化，由FPKM To Count，使用也是相对比较方便。详情到网址：StringTie (jhu.edu) - http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual
Count
定义：高通量测序中比对到exon上的reads数。可以使用featureCounts、HTseq-count等软件进行计算。
优点：可以有效说明该区域是否真的有表达及真实的表达丰度。能够近似呈现真实的表达情况。
缺点：由于exon长度不同，难以进行不同exon丰度比较；由于测序总数不同，难以对不同测序样本间比较。
FPKM
FPKM: FPKM的全称为Fragments Per Kilobase Million，Fragments Per Kilobase of exon model per Million mapped fragments(每千个碱基的转录每百万映射读取的fragments)。通俗讲，把比对到的某个基因的Fragment数目，除以基因的长度，其比值再除以所有基因的总长度。注意，这里的基因长度是指基因外显子的总长度。
RPKM
RPKM: Reads Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的reads)；
FPKM与RPKM的区别
RPKM通常用于单端测序，FPKM常用于双端测序
如果是单端测序，那么一个fragmetns就对应了一条read，如下所示：
如果是双端测序，那么一条fragments就对应两条reads，当然，有时候双端测序也有可能出现一条fragment对应一条read（另外一条read有可能会因为质量低而被剔除），FPKM就保证了，一条fragment的两条reads不会被统计2次，如下所示：
FPKM是以fragment为准，而不是以reads数为准，它们的计算方式是一样的
RPM
定义：RPM/CPM: Reads/Counts of exon model per Million mapped reads (每百万映射读取的reads)
公式：RPM = ExonMappedReads * 10^6 /TotalMappedReads
优点：利于进行样本间比较。根据比对到基因组上的总reads count，进行标准化。即：不论比对到基因组上的总reads count是多少，都将总reads count标准化为10^6。sRNA_seq等测序长度较短的高通量测序经常采用RPM进行标准化，因为sRNA长度差异较小，18-35 nt较多，所以长度对不同的small RNAs相互比较影响较小 (优点：计算简单、方便。)。
缺点：未消除exon长度造成的表达差异，难以进行样本内exon差异表达的比较。
TPM
定义：TPM的全称为Transcripts per million，Transcripts Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts)
解释：Ni为比对到第i个exon的reads数；Li为第i个exon的长度；sum(N1/L1+N2/L2 + ... + Nn/Ln)为所有 (n个)exon按长度进行标准化之后数值的和。
获得gene外显子长度!
library(GenomicFeatures)
## 导入gff3文件
txdb <- makeTxDbFromGFF("ITAG4.1_gene_models.gff", format = "gff")
## 获取外显子位置
exons_gene <- exonsBy(txdb, by = "gene")
## 去除外显子重叠部分，计算外显子长度
exons_gene_len <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
exons_gene_len <- as.matrix(t(exons_gene_len))
write.csv(exons_gene_len,"tomato_gene_length_4.1.csv", row.names = T)
我使用番茄基因组4.1的注释文件，但是提取只得到100多个基因的外显子长度，不知道是哪一步出现问题。你可以使用你的物种注释文件操作一下。如果你会写编程，也是使用脚本进行提取，也是相对容易的循环就可以提取得到。
Code | 各表达量间的转化
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)
}
# An example
################################################################################
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$tpm <- with(countDf, countToTpm(count, effLength))
countDf$fpkm <- with(countDf, countToFpkm(count, effLength))
with(countDf, all.equal(tpm, fpkmToTpm(fpkm)))
countDf$effCounts <- with(countDf, countToEffCounts(count, length, effLength))
参考：
关于readsCount、RPKM/FPKM、RPM(CPM)、TPM的理解
RNA-Seq的count、RPKM/FPKM、CPM、TPM的关系
RPKM, FPKM and TPM, clearly explained | RNA-Seq Blog
StatQuest学习笔记24——RPKM FPKM TPM
Htseq Count To Fpkm | KeepNotes blog
“小杜的生信筆記”公众号、知乎、简书平台，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！
网友评论