RPKM,FPKM,TPM等標(biāo)準(zhǔn)化方法還有那些問題？DESeq2的標(biāo)準(zhǔn)化方法的原理就是提高中等表達(dá)基因的地位一個(gè)例子

在上一節(jié)的StatQuest生物統(tǒng)計(jì)學(xué)專題中，我們簡單直白的討論了RPKM,FPKM,TPM的定義和生物學(xué)意義，明白了RPKM,FPKM,TPM標(biāo)準(zhǔn)化方法就是為了去除基因長度和測序深度對測序Read數(shù)的影響，見StatQuest生物統(tǒng)計(jì)學(xué)專題 - RPKM,FPKM,TPM，并且在一般意義上，更推薦使用TPM標(biāo)準(zhǔn)化方法。

然而不得不說明的是，諸如DESeq2、edgeR等差異表達(dá)分析軟件都不是使用的RPKM,FPKM,TPM方法，為什么呢？

RPKM,FPKM,TPM等標(biāo)準(zhǔn)化方法還有那些問題？

其實(shí)TPM之類的標(biāo)準(zhǔn)化方法雖然解決了基因長度和測序深度的影響的問題，但還是不能解決一個(gè)問題：那就是測序文庫組成不同造成的差異。

什么意思呢？

我們知道RPKM,FPKM,TPM是可以解決由于測序深度的差異而引起的Read數(shù)變異，如下圖所示，樣本2#的總Read數(shù)是樣本1#的2倍，每個(gè)基因的Read數(shù)也是1#的2倍。我們知道這種Read數(shù)差異并不是基因表達(dá)的不同，而是由于測序深度不同所致，只需要將樣本1#、2#除以各自的總Read數(shù)，那么這個(gè)Read數(shù)變異就會得到修正。

讓我們再考慮兩種新的情況，我們知道RNA-seq或其他高通量方法，往往是不同組織間的基因表達(dá)比較，比如肝細(xì)胞或脾細(xì)胞，他們會有一些各自特異性表達(dá)的基因，或者對于一些敲減基因的實(shí)驗(yàn)，有些基因在其中一個(gè)樣本中高表達(dá)，而在另一個(gè)樣本中不表達(dá)。

如下圖所示的兩個(gè)樣本，兩者的測序深度是一樣的，總Read數(shù)相同。但是由于A2M基因在樣本2#中不表達(dá)，導(dǎo)致563個(gè)Read會分配到樣本2#的其他5個(gè)基因中去：如樣本2#的基因A1BG的Read數(shù)是235，大于樣本1#中的30。然而實(shí)際上這種差異并不是生物學(xué)效應(yīng)所致，而是由于基因A2M被敲減所致，并不是A1BG等5個(gè)基因在樣本2#中表達(dá)增加了，而是基因A2M在樣本2#中表達(dá)減少了。

在這種情況下，這種測序文庫組成不同的差異是RPKM,FPKM,TPM等方法無法解決的。

DESeq2的標(biāo)準(zhǔn)化方法的原理就是提高中等表達(dá)基因的地位

那么如何解決這個(gè)問題呢？

本周先看一下DESeq2是如何進(jìn)行l(wèi)ibrary normalization的，DESeq2的標(biāo)準(zhǔn)化方法共有7步，看起來很繁瑣，但是原理很簡單，它有一個(gè)貫穿始終的基本思想——提高中等表達(dá)基因的地位。

而且這7步只是為了得到一個(gè)標(biāo)準(zhǔn)化因子，并進(jìn)行變換。

首先以下述數(shù)據(jù)集為例，共有3個(gè)樣本，每個(gè)樣本有3個(gè)基因：

第一步 對Read矩陣取對數(shù)變換

DESeq2默認(rèn)是使用自然對數(shù)，也可以使用log2或log10。

第二歩 取各基因的平均數(shù)

要說明的是，由于樣本1的基因1的Read數(shù)是0，所以在取對數(shù)時(shí)，它的值是-Inf（負(fù)無窮大），因此對基因1取平均數(shù)時(shí)就直接得出是-Inf即可。

為何要取對數(shù)？

其實(shí)是為了減少高表達(dá)異常值對標(biāo)準(zhǔn)化的影響，需要注意的是異常值不代表是錯(cuò)誤值，只是說它相比數(shù)據(jù)趨勢比較異常。

以原始表達(dá)矩陣為例，基因3的3個(gè)Read數(shù)分別是33、55、200，尤其是200，相比較整個(gè)表達(dá)矩陣來說是一個(gè)高表達(dá)異常值。

如果對原始表達(dá)矩陣求基因均值，那么基因3的均值是(33+55+200)/3 = 96，而對于對數(shù)變換后的表達(dá)矩陣來說，基因3的均值是(3.5+4.0+5.3)/3 = 4.3； e^4.3 =73.7，73.7<96，也就是說對數(shù)變換后基因3的均數(shù)受到200的影響更?。ㄟ@種取對數(shù)求得的平均數(shù)是幾何平均數(shù)）。

第三步 過濾掉-Inf基因

將存在-Inf值的基因過濾掉，過濾掉的基因不再參與標(biāo)準(zhǔn)化因子的計(jì)算。

實(shí)際上，這一步是把在一個(gè)或多個(gè)樣本中存在零表達(dá)的基因剔除。假定本實(shí)驗(yàn)是在比較不同組織細(xì)胞如肝細(xì)胞和脾細(xì)胞的表達(dá)量差異，那么這一步會剔除掉組織特異性表達(dá)的基因，而只保留管家基因——在不同細(xì)胞中都或多或少會表達(dá)的基因。

第四步 將對數(shù)矩陣減去對數(shù)均值，得到對數(shù)比值矩陣

將對數(shù)矩陣的每個(gè)Read分別減去此基因的對數(shù)矩陣。

對數(shù)相減的意義何在？

對數(shù)相減其實(shí)是真數(shù)相除，也就是：

log(reads for gene X)-log(average for gene X) = log(reads for gene X/average for gene X)

注意：此時(shí)的average for gene X是幾何平均數(shù)。

第五步 計(jì)算每個(gè)樣本的對數(shù)比值矩陣的中位數(shù)

取中位數(shù)，而不是均值，也是為了進(jìn)一步降低異常值的影響，具有較大表達(dá)差異的基因?qū)χ形粩?shù)的影響甚微?？紤]到絕大部分情況下，表達(dá)差異大的基因都是很少的，所以這個(gè)“中位數(shù)”更能代表的是中等表達(dá)基因或管家基因的情況。

第六步 將對數(shù)中位數(shù)轉(zhuǎn)換為其相應(yīng)的真數(shù)，得到各個(gè)樣本的標(biāo)準(zhǔn)化因子

第七步 將原始表達(dá)矩陣除以這個(gè)標(biāo)準(zhǔn)化因子

原始矩陣的每個(gè)樣本的全部Read數(shù)均除以各自的標(biāo)準(zhǔn)化因子（包括較大表達(dá)差異的基因）。

我們可以看到標(biāo)準(zhǔn)化之后，對于基因3來說，樣本1#的Read數(shù)提升，而樣本3#的Read數(shù)下降，基因3在3個(gè)樣本中的表達(dá)其實(shí)是接近的。

做一下總結(jié)：

1. 對數(shù)變換可以減少只表達(dá)在某些樣本中的基因的影響，同時(shí)還可以減少異常值的影響（幾何平均數(shù)）；

2. 將在某些樣本中表達(dá)量為0的基因剔除，不參與中位數(shù)計(jì)算，可以剔除特異性表達(dá)基因的影響；

3. 中位數(shù)算法可以進(jìn)一步降低高差異表達(dá)基因的影響，而提高中等表達(dá)的基因的地位；

標(biāo)準(zhǔn)化因子的生物學(xué)意義

其實(shí)這個(gè)標(biāo)準(zhǔn)化因子算法就是選出一個(gè)有代表性的gene X（其實(shí)是每個(gè)樣本一個(gè)代表性gene X），而這個(gè)gene X的reads for gene X/average for gene X比值就是標(biāo)準(zhǔn)化因子。

只不過選取gene X的時(shí)候，通過對數(shù)變換和中位數(shù)的方法，更多的參考了中等表達(dá)基因和管家基因的數(shù)據(jù)趨勢，而剔除了特異性表達(dá)基因和高差異表達(dá)基因的影響。

相比較RPKM,FPKM,TPM標(biāo)準(zhǔn)化方法是除以總Read數(shù)，DESeq2標(biāo)準(zhǔn)化方法是除以一個(gè)有代表性基因的Read數(shù)，只不過這個(gè)Read數(shù)進(jìn)行了變換（它除以了幾何平均Read數(shù)， reads for gene X/average for gene X）。因?yàn)楦芴幚泶嬖谔禺愋员磉_(dá)基因和高差異表達(dá)基因的數(shù)據(jù)。

一個(gè)例子

我們按照這個(gè)DESeq2的標(biāo)準(zhǔn)化方法的思想，對圖2中的數(shù)據(jù)進(jìn)行一個(gè)簡單的標(biāo)準(zhǔn)化，沒有完全按照上述的7步法，只是體會這種標(biāo)準(zhǔn)化的意思（結(jié)果沒有大的差異，但是算法并不完全正確）。

#Gene        Sample#1    Saple#2

A1BG        30          235
A1BG-AS1    24          188
A1CF        0           0
A2M            563         0
A2M-AS1        5           39
A2ML1        13          102

• 首先找到參照基因Gene X

  對于Sample 1#來說，A1BG-AS1和A2ML1是中位數(shù)基因（不計(jì)算A1CF、A2M），而樣本 2#也是同樣。

  為了簡便，就用A2ML1基因了，其實(shí)結(jié)果是類似的。

原本應(yīng)該根據(jù)ln(reads for gene X/average for gene X)的值計(jì)算中位數(shù)，這里直接根據(jù)原始Read值計(jì)算，會有一定的差異。

• 求解reads for gene X/average for gene X比值

  由于，

  average for gene A2ML1= (ln13+ln102)/2 = 3.12

  于是，

  reads for gene A2ML1/average for gene A2ML1= (13,102)/3.12 = (4.16,32.66)

  也就是說兩個(gè)樣本的標(biāo)準(zhǔn)化因子分別是4.16和32.66。

• 進(jìn)行標(biāo)準(zhǔn)化變換

  將原始矩陣按照樣本的不同除以各自的標(biāo)準(zhǔn)化因子，得下表：

  可以發(fā)現(xiàn)，除了A2M基因有表達(dá)差異外，其他基因表達(dá)無明顯差異。

#Gene        Sample#1        Saple#2

A1BG        7.205869671     7.194095374
A1BG-AS1    5.764695737     5.755276299
A1CF        0               0
A2M            135.2301542     0
A2M-AS1        1.200978278     1.1939137
A2ML1        3.122543524     3.122543524

參考資料

StatQuest課程：https:///video-index/