关于Count，FPKM，TPM，RPKM等表达量的计算

原文链接：关于Count，FPKM，TPM，RPKM等表达量的计算及转换 | 干货

写在前面

今天使用count值转化TPM，或是使用FPKM转换成TPM。这样的教程，我们在前面已经出国一起相对比较详细的教程了，一文了解Count、FPKM、RPKM、TPM | 相互间的转化，在这个教程中，我们也归纳了各个数值的含义。但是，也许你看到后，会添加到自己的收藏夹中，但是，后面就没有看了。自己也是这样的，一个人的时间和精力是有限的，我们不可能有那么多的精力。因此，做学习笔记就有很大的帮助，当自己使用的时候有地方找寻。

本教程涉及的数据、代码和文件等在社群中可获得！！

回顾一下知识点

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按长度进行标准化之后数值的和。

如何准换呢？？

方法一

使用GenomicFeature包导入对应的gtf文件或gff文件，获得基因外显子长度。
**优点：**方便，快捷。
**缺点：**对与gtf或gff文件的要求比较高，我自己使用stringtie组装的注释文件，是无法获得导入的。此外，这个方法只能获得是gene的表达量，若你想获得transcript的表达量，自己未成功。

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)

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))

方法二

使用featureCount等计算出Count值，获得结果中就有对应的外显子长度，因此，你可以直接使用其进行转化。

1. 导入数据

##
count_df <- read.csv("count.csv",header = T, row.names = 1)
dim(count_df)
names(count_df)

##'@提取Count值
expr_df <- count_df[,2:ncol(count_df)]
head(expr_df)
###'@基因表达量之和大于0
expr_df <- expr_df[rowSums(expr_df) > 0,]
dim(expr_df)
##'@保存count值矩阵
write.csv(expr_df,"过滤后_count.csv")

计算TPM

##'@计算的TPM
##'@提取基因长度，基因长度需要转化成kb
gene_length_kb <- count_df$Length / 1000
head(gene_length_kb)

### 每千碱基reads（per million scaling factor)长度标准化
data_rpk <- expr_df /gene_length_kb
##'@每百万
TPM <- t(t(data_rpk) / colSums(data_rpk) * 1000000)
head(TPM)
## 求均值
avg_tmp <- data.frame(avg_tmp = rowMeans(TPM))
head(avg_tmp)
##'@保存数据
write.csv(TPM,"Tomato_TPM.csv")

计算FPKM

FPKM <- t(t(data_rpk) / colSums(expr_df) * 10^6)
##'@保存数据
write.csv(TPM,"Tomato_FPKM.csv")

参考：

1. https://mp.weixin.qq.com/s/JOVYGh7pO8SiW-9iIeMO6Q
2. https://www.jianshu.com/p/aec488f358d2
3. https://www.jianshu.com/p/6b6bb306b76e
4. https://www.bioinfo-scrounger.com/archives/407/

方法三

我们直接使用的perl脚本进行转化，很是方便。

perl CountToFPKM.pl Count.txt > FPKM.txt 

原文链接：关于Count，FPKM，TPM，RPKM等表达量的计算及转换 | 干货

代码和数据链接：

本教程涉及的数据、代码和文件等在社群中可获得！！

若我们的分享对你有用，希望您可以点赞+收藏+转发，这是对小杜最大的支持。

往期文章：

1. 复现SCI文章系列专栏

2. 《生信知识库订阅须知》,同步更新，易于搜索与管理。

3. 最全WGCNA教程（替换数据即可出全部结果与图形）

• WGCNA分析 | 全流程分析代码 | 代码一

• WGCNA分析 | 全流程分析代码 | 代码二

• WGCNA分析 | 全流程代码分享 | 代码三

• WGCNA分析 | 全流程分析代码 | 代码四

• WGCNA分析 | 全流程分析代码 | 代码五(最新版本)

---

4. 精美图形绘制教程

• 精美图形绘制教程

5. 转录组分析教程

转录组上游分析教程[零基础]

一个转录组上游分析流程 | Hisat2-Stringtie

小杜的生信筆記 ，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料!!