相信學(xué)習(xí)生信的小伙伴不難發(fā)現(xiàn)，在定量的時(shí)候大多的文章或者公司的結(jié)果都有這幾種定量的方式：FPKM（RPKM）纠炮，TPM月趟，CPM(RPM)還有count等等，這些究竟該如何使用恢口，如何相互轉(zhuǎn)換孝宗，該用哪個(gè)？今天就和大家談?wù)勥@個(gè)問題耕肩。
一因妇、基礎(chǔ)知識(shí)
1. count
什么是count呢？實(shí)際上count就是原始序列比對(duì)到參考基因組上后猿诸，對(duì)應(yīng)的基因有多少條reads命中到這個(gè)基因婚被。是后續(xù)一切其它歸一化方法的基礎(chǔ)
適用范圍：通常count可以用于后續(xù)的DESeq2，edgeR等軟件進(jìn)行差異分析梳虽，因?yàn)樗麄儠?huì)對(duì)count進(jìn)行另一種歸一化的方法——TMM后址芯，默認(rèn)使用負(fù)二項(xiàng)分布檢驗(yàn)進(jìn)行差異分析。
2.FPKM/RPKM
FPKM和RPKM分別對(duì)應(yīng)Fragments Per Kilobase of exon model per Million mapped fragments(每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的fragments)和Reads Per Kilobase of exon model per Million mapped reads (每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的reads)，兩者的計(jì)算方法是一致的谷炸，只是應(yīng)用的場(chǎng)景有所不同北专，通常前者用于雙端測(cè)序，后者用于單端測(cè)序淑廊，其余的內(nèi)涵是一致的逗余。
FPKM的計(jì)算方法如下圖，C為比對(duì)到基因的fragments數(shù)（count）季惩，N為比對(duì)到參考基因的總fragments數(shù)录粱，L為該基因的有效長(zhǎng)度。FPKM的初衷是為了能消除基因長(zhǎng)度和測(cè)序量差異對(duì)計(jì)算基因表達(dá)的影響画拾。
3.TPM
TPM代表Transcripts Per Kilobase of exon model per Million mapped reads (每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的Transcripts)是FPKM的一種改進(jìn)算法啥繁，如果數(shù)學(xué)敏感的讀者應(yīng)該會(huì)發(fā)現(xiàn)在FPKM的公式中，當(dāng)比較同一個(gè)基因時(shí)青抛，除了他們的C可能不同旗闽，測(cè)序總量帶來的N同樣的是不同的，兩個(gè)變量都不同的情況進(jìn)行比較是可笑的蜜另，所以TPM的作者認(rèn)為FPKM不適合在不同組樣本之間進(jìn)行比較于是提出了TPM适室。
公式如下：
這里看著有點(diǎn)復(fù)雜，其實(shí)說白了就是先消除長(zhǎng)度的影響举瑰，先把每個(gè)基因除去他們的長(zhǎng)度捣辆，在求和然后用對(duì)一個(gè)基因走正常的FPKM的運(yùn)算后除去剛才的求和的值后乘百萬。這樣的好處就使得剛才N不等的問題消除掉了此迅，理論上就更適合不同組的樣本之間的比較了汽畴。
適用范圍：同F(xiàn)PKM，同時(shí)也可以粗略的比較不同組的基因的表達(dá)量（不推薦Ｋ市颉）
4.CPM/RPM
這里這里的CPM或者RPM(read counts per million) 忍些，其實(shí)就是不考慮長(zhǎng)度不考慮長(zhǎng)度直接把count除總count的N后乘百萬就完事了，很多公司很坑爹的因?yàn)閙iRNA的長(zhǎng)度基本相同不考慮后直接把CPM寫成TPM坎怪，帶來了嚴(yán)重的誤導(dǎo)罢坝！
適用范圍：很難評(píng)估有效基因長(zhǎng)度或基因長(zhǎng)度基本相當(dāng)?shù)慕M學(xué)，如circRNA-seq芋忿、ChIP-seq炸客、CUT&Tag、ATAC戈钢、MeRIP-seq等。
二是尔、幾種歸一化方法的比較
看了上述的幾種歸一化方法殉了，是否會(huì)讓你覺得TPM是一種完美的歸一化方法，其實(shí)非也拟枚，2020年的一個(gè)研究結(jié)果如下：
圖中y軸代表的歸一化方法的排名薪铜，可以看到TMM法基本吊打全局众弓，和TMM相比其它方法都談不上優(yōu)秀，看似很牛的TPM中位數(shù)還不如FPKM隔箍！
關(guān)于TMM歸一化算法和優(yōu)勢(shì)感興趣的我們留一期單獨(dú)談?wù)劇?/p>
三谓娃、count轉(zhuǎn)FPKM、TPM
這里首先引入一個(gè)概念蜒滩，上面談到的基因長(zhǎng)度都是指有效基因長(zhǎng)度滨达，通常認(rèn)為有效基因長(zhǎng)度等于所有非冗余的外顯子的長(zhǎng)度總和。明白了這一點(diǎn)我們就可以計(jì)算FPKM/TPM了俯艰，以R為例代碼如下：
首先捡遍，得到用htseq等工具或者TCGA下載到的count文件，以及對(duì)應(yīng)物種的gtf文件（Ensembl下載）竹握，讀到R中画株，這里以hg38.gtf和count.tsv為例子
library(tidyverse)
#讀gtf文件，計(jì)算所有外顯子的長(zhǎng)度
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)
#計(jì)算基因的非冗余外顯子的長(zhǎng)度啦辐，即獲得有效基因長(zhǎng)度
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的表達(dá)量矩陣谓传，列名為gene_id和count，其中g(shù)ene_id是和gtf文件一致的芹关，然后和剛才計(jì)算得到的有效基因長(zhǎng)度合并
expmat <- read_tsv(count.tsv) %>% inner_join(gtf , by = 'gene_id' ) %>% drop_na()
#各種轉(zhuǎn)換的方法
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) ) %>% #轉(zhuǎn)FPKM
mutate( TPM = countToTpm(.$count, .$est_len) ) %>% #轉(zhuǎn)TPM
mutate( CPM = countToCPM(.$count) ) %>% #轉(zhuǎn)CPM
select(-est_len) %>% write_tsv("out.xls") #輸出結(jié)果
