【生信基础】深入理解FPKM/TPM，只有GTF文件如何将cou

2022-05-26 本文已影响0人 xizzy

相信学习生信的小伙伴不难发现，在定量的时候大多的文章或者公司的结果都有这几种定量的方式：FPKM（RPKM），TPM，CPM(RPM)还有count等等，这些究竟该如何使用，如何相互转换，该用哪个？今天就和大家谈谈这个问题。

一、基础知识

1. count

什么是count呢？实际上count就是原始序列比对到参考基因组上后，对应的基因有多少条reads命中到这个基因。是后续一切其它归一化方法的基础
适用范围：通常count可以用于后续的DESeq2，edgeR等软件进行差异分析，因为他们会对count进行另一种归一化的方法——TMM后，默认使用负二项分布检验进行差异分析。

2.FPKM/RPKM

FPKM和RPKM分别对应Fragments Per Kilobase of exon model per Million mapped fragments(每千个碱基的转录每百万映射读取的fragments)和Reads Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的reads)，两者的计算方法是一致的，只是应用的场景有所不同，通常前者用于双端测序，后者用于单端测序，其余的内涵是一致的。
FPKM的计算方法如下图，C为比对到基因的fragments数（count），N为比对到参考基因的总fragments数，L为该基因的有效长度。FPKM的初衷是为了能消除基因长度和测序量差异对计算基因表达的影响。

fpkm计算公式

3.TPM

TPM代表Transcripts Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts)是FPKM的一种改进算法，如果数学敏感的读者应该会发现在FPKM的公式中，当比较同一个基因时，除了他们的C可能不同，测序总量带来的N同样的是不同的，两个变量都不同的情况进行比较是可笑的，所以TPM的作者认为FPKM不适合在不同组样本之间进行比较于是提出了TPM。
公式如下：

TPM计算公式

这里看着有点复杂，其实说白了就是先消除长度的影响，先把每个基因除去他们的长度，在求和然后用对一个基因走正常的FPKM的运算后除去刚才的求和的值后乘百万。这样的好处就使得刚才N不等的问题消除掉了，理论上就更适合不同组的样本之间的比较了。
适用范围：同FPKM，同时也可以粗略的比较不同组的基因的表达量（不推荐！）

4.CPM/RPM

这里这里的CPM或者RPM(read counts per million) ，其实就是不考虑长度不考虑长度直接把count除总count的N后乘百万就完事了，很多公司很坑爹的因为miRNA的长度基本相同不考虑后直接把CPM写成TPM，带来了严重的误导！
适用范围：很难评估有效基因长度或基因长度基本相当的组学，如circRNA-seq、ChIP-seq、CUT&Tag、ATAC、MeRIP-seq等。

二、几种归一化方法的比较

看了上述的几种归一化方法，是否会让你觉得TPM是一种完美的归一化方法，其实非也，2020年的一个研究结果如下：

归一化方法排名

图中y轴代表的归一化方法的排名，可以看到TMM法基本吊打全局，和TMM相比其它方法都谈不上优秀，看似很牛的TPM中位数还不如FPKM！
关于TMM归一化算法和优势感兴趣的我们留一期单独谈谈。

三、count转FPKM、TPM

这里首先引入一个概念，上面谈到的基因长度都是指有效基因长度，通常认为有效基因长度等于所有非冗余的外显子的长度总和。明白了这一点我们就可以计算FPKM/TPM了，以R为例代码如下：

首先，得到用htseq等工具或者TCGA下载到的count文件，以及对应物种的gtf文件（Ensembl下载），读到R中，这里以hg38.gtf和count.tsv为例子

library(tidyverse)
#读gtf文件，计算所有外显子的长度
gtf <- read_tsv("hg38.gtf", comment="#", col_names=c('chr','source','type','start','end','score','strand','phase','attributes')) %>% filter(type=='exon') %>% mutate(len = end - start + 1) %>% select(start, end, attributes,len)
#计算基因的非冗余外显子的长度，即获得有效基因长度
gtf$attributes %>% str_extract(., "gene_id \"[\\w|\\.]+") %>% str_remove(., "gene_id \"") -> gtf$gene_id
gtf %>% select(start, end, gene_id, len) %>% distinct(start,end,gene_id, .keep_all = T) %>% select(gene_id,len) %>% group_by(gene_id) %>% summarise(est_len=sum(len)) -> gtf
#读取count的表达量矩阵，列名为gene_id和count，其中gene_id是和gtf文件一致的，然后和刚才计算得到的有效基因长度合并
expmat <- read_tsv(count.tsv) %>% inner_join(gtf , by = 'gene_id' ) %>% drop_na() 

#各种转换的方法
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))
}
 
countToCPM <- function( counts)
{
    N <- sum(counts)
    exp( log(counts) + log(1e6) - log(N) )
}


expmat %>% 
    mutate( FPKM = countToFpkm(.$count, .$est_len) ) %>% #转FPKM
    mutate( TPM = countToTpm(.$count, .$est_len) )  %>% #转TPM
    mutate( CPM = countToCPM(.$count) ) %>% #转CPM
   select(-est_len) %>% write_tsv("out.xls") #输出结果

输出效果