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

2024-03-22 14:12
文章标签 计算 count 表达 tpm fpkm rpkm

本文主要是介绍关于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脚本,获得FPKM

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的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!

这篇关于关于Count,FPKM,TPM,RPKM等表达量的计算的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/835456

相关文章

使用C#代码计算数学表达式实例

《使用C#代码计算数学表达式实例》这段文字主要讲述了如何使用C#语言来计算数学表达式,该程序通过使用Dictionary保存变量,定义了运算符优先级,并实现了EvaluateExpression方法来... 目录C#代码计算数学表达式该方法很长,因此我将分段描述下面的代码片段显示了下一步以下代码显示该方法如

如何用Java结合经纬度位置计算目标点的日出日落时间详解

《如何用Java结合经纬度位置计算目标点的日出日落时间详解》这篇文章主详细讲解了如何基于目标点的经纬度计算日出日落时间,提供了在线API和Java库两种计算方法,并通过实际案例展示了其应用,需要的朋友... 目录前言一、应用示例1、天安门升旗时间2、湖南省日出日落信息二、Java日出日落计算1、在线API2

poj 1113 凸包+简单几何计算

题意: 给N个平面上的点,现在要在离点外L米处建城墙,使得城墙把所有点都包含进去且城墙的长度最短。 解析: 韬哥出的某次训练赛上A出的第一道计算几何,算是大水题吧。 用convexhull算法把凸包求出来,然后加加减减就A了。 计算见下图: 好久没玩画图了啊好开心。 代码: #include <iostream>#include <cstdio>#inclu

uva 1342 欧拉定理(计算几何模板)

题意: 给几个点,把这几个点用直线连起来,求这些直线把平面分成了几个。 解析: 欧拉定理: 顶点数 + 面数 - 边数= 2。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#inc

uva 11178 计算集合模板题

题意: 求三角形行三个角三等分点射线交出的内三角形坐标。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <

XTU 1237 计算几何

题面: Magic Triangle Problem Description: Huangriq is a respectful acmer in ACM team of XTU because he brought the best place in regional contest in history of XTU. Huangriq works in a big compa

音视频入门基础:WAV专题(10)——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

一、引言 从文章《音视频入门基础:WAV专题(6)——通过FFprobe显示WAV音频文件每个数据包的信息》中我们可以知道,通过FFprobe命令可以打印WAV音频文件每个packet(也称为数据包或多媒体包)的信息,这些信息包含该packet的pts、dts: 打印出来的“pts”实际是AVPacket结构体中的成员变量pts,是以AVStream->time_base为单位的显

计算数组的斜率,偏移,R2

模拟Excel中的R2的计算。         public bool fnCheckRear_R2(List<double[]> lRear, int iMinRear, int iMaxRear, ref double dR2)         {             bool bResult = true;             int n = 0;             dou

GPU 计算 CMPS224 2021 学习笔记 02

并行类型 (1)任务并行 (2)数据并行 CPU & GPU CPU和GPU拥有相互独立的内存空间,需要在两者之间相互传输数据。 (1)分配GPU内存 (2)将CPU上的数据复制到GPU上 (3)在GPU上对数据进行计算操作 (4)将计算结果从GPU复制到CPU上 (5)释放GPU内存 CUDA内存管理API (1)分配内存 cudaErro

Java - BigDecimal 计算分位(百分位)

日常开发中,如果使用数据库来直接查询一组数据的分位数,就比较简单,直接使用对应的函数就可以了,例如:         PERCENT_RANK() OVER(PARTITION BY 分组列名 ORDER BY 目标列名) AS 目标列名_分位数         如果是需要在代码逻辑部分进行分位数的计算,就需要我们自己写一个工具类来支持计算了 import static ja