RPKM,FPKM,TPM等標準化方法還有那些問題？DESeq2的標準化方法的原理就是提高中等表達基因的地位一個例子
在上一節(jié)的StatQuest生物統(tǒng)計學(xué)專題中，我們簡單直白的討論了RPKM,FPKM,TPM的定義和生物學(xué)意義，明白了RPKM,FPKM,TPM標準化方法就是為了去除基因長度和測序深度對測序Read數(shù)的影響，見StatQuest生物統(tǒng)計學(xué)專題 - RPKM,FPKM,TPM，并且在一般意義上，更推薦使用TPM標準化方法。
然而不得不說明的是，諸如DESeq2、edgeR等差異表達分析軟件都不是使用的RPKM,FPKM,TPM方法，為什么呢？
其實TPM之類的標準化方法雖然解決了基因長度和測序深度的影響的問題，但還是不能解決一個問題：那就是測序文庫組成不同造成的差異。
什么意思呢？
我們知道RPKM,FPKM,TPM是可以解決由于測序深度的差異而引起的Read數(shù)變異，如下圖所示，樣本2#的總Read數(shù)是樣本1#的2倍，每個基因的Read數(shù)也是1#的2倍。我們知道這種Read數(shù)差異并不是基因表達的不同，而是由于測序深度不同所致，只需要將樣本1#、2#除以各自的總Read數(shù)，那么這個Read數(shù)變異就會得到修正。
讓我們再考慮兩種新的情況，我們知道RNA-seq或其他高通量方法，往往是不同組織間的基因表達比較，比如肝細胞或脾細胞，他們會有一些各自特異性表達的基因，或者對于一些敲減基因的實驗，有些基因在其中一個樣本中高表達，而在另一個樣本中不表達。
如下圖所示的兩個樣本，兩者的測序深度是一樣的，總Read數(shù)相同。但是由于A2M基因在樣本2#中不表達，導(dǎo)致563個Read會分配到樣本2#的其他5個基因中去：如樣本2#的基因A1BG的Read數(shù)是235，大于樣本1#中的30。然而實際上這種差異并不是生物學(xué)效應(yīng)所致，而是由于基因A2M被敲減所致，并不是A1BG等5個基因在樣本2#中表達增加了，而是基因A2M在樣本2#中表達減少了。
在這種情況下，這種測序文庫組成不同的差異是RPKM,FPKM,TPM等方法無法解決的。
那么如何解決這個問題呢？
本周先看一下DESeq2是如何進行l(wèi)ibrary normalization的，DESeq2的標準化方法共有7步，看起來很繁瑣，但是原理很簡單，它有一個貫穿始終的基本思想——提高中等表達基因的地位。
而且這7步只是為了得到一個標準化因子，并進行變換。
首先以下述數(shù)據(jù)集為例，共有3個樣本，每個樣本有3個基因：
第一步 對Read矩陣取對數(shù)變換
DESeq2默認是使用自然對數(shù)，也可以使用log2或log10。
第二歩 取各基因的平均數(shù)
要說明的是，由于樣本1的基因1的Read數(shù)是0，所以在取對數(shù)時，它的值是-Inf（負無窮大），因此對基因1取平均數(shù)時就直接得出是-Inf即可。
為何要取對數(shù)？
其實是為了減少高表達異常值對標準化的影響，需要注意的是異常值不代表是錯誤值，只是說它相比數(shù)據(jù)趨勢比較異常。
以原始表達矩陣為例，基因3的3個Read數(shù)分別是33、55、200，尤其是200，相比較整個表達矩陣來說是一個高表達異常值。
如果對原始表達矩陣求基因均值，那么基因3的均值是
(33+55+200)/3 = 96，而對于對數(shù)變換后的表達矩陣來說，基因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值的基因過濾掉，過濾掉的基因不再參與標準化因子的計算。
實際上，這一步是把在一個或多個樣本中存在零表達的基因剔除。假定本實驗是在比較不同組織細胞如肝細胞和脾細胞的表達量差異，那么這一步會剔除掉組織特異性表達的基因，而只保留管家基因——在不同細胞中都或多或少會表達的基因。
第四步 將對數(shù)矩陣減去對數(shù)均值，得到對數(shù)比值矩陣
將對數(shù)矩陣的每個Read分別減去此基因的對數(shù)矩陣。
對數(shù)相減的意義何在？
對數(shù)相減其實是真數(shù)相除，也就是：
log(reads for gene X)-log(average for gene X) = log(reads for gene X/average for gene X)注意：此時的average for gene X是幾何平均數(shù)。
第五步 計算每個樣本的對數(shù)比值矩陣的中位數(shù)
取中位數(shù)，而不是均值，也是為了進一步降低異常值的影響，具有較大表達差異的基因?qū)χ形粩?shù)的影響甚微?？紤]到絕大部分情況下，表達差異大的基因都是很少的，所以這個“中位數(shù)”更能代表的是中等表達基因或管家基因的情況。
第六步 將對數(shù)中位數(shù)轉(zhuǎn)換為其相應(yīng)的真數(shù)，得到各個樣本的標準化因子
第七步 將原始表達矩陣除以這個標準化因子
原始矩陣的每個樣本的全部Read數(shù)均除以各自的標準化因子（包括較大表達差異的基因）。
我們可以看到標準化之后，對于基因3來說，樣本1#的Read數(shù)提升，而樣本3#的Read數(shù)下降，基因3在3個樣本中的表達其實是接近的。
做一下總結(jié)：
對數(shù)變換可以減少只表達在某些樣本中的基因的影響，同時還可以減少異常值的影響（幾何平均數(shù)）；
將在某些樣本中表達量為0的基因剔除，不參與中位數(shù)計算，可以剔除特異性表達基因的影響；
中位數(shù)算法可以進一步降低高差異表達基因的影響，而提高中等表達的基因的地位；
標準化因子的生物學(xué)意義
其實這個標準化因子算法就是選出一個有代表性的gene X（其實是每個樣本一個代表性gene X），而這個gene X的reads for gene X/average for gene X比值就是標準化因子。
只不過選取gene X的時候，通過對數(shù)變換和中位數(shù)的方法，更多的參考了中等表達基因和管家基因的數(shù)據(jù)趨勢，而剔除了特異性表達基因和高差異表達基因的影響。
相比較RPKM,FPKM,TPM標準化方法是除以總Read數(shù)，DESeq2標準化方法是除以一個有代表性基因的Read數(shù)，只不過這個Read數(shù)進行了變換（它除以了幾何平均Read數(shù)， reads for gene X/average for gene X）。因為更能處理存在特異性表達基因和高差異表達基因的數(shù)據(jù)。
我們按照這個DESeq2的標準化方法的思想，對圖2中的數(shù)據(jù)進行一個簡單的標準化，沒有完全按照上述的7步法，只是體會這種標準化的意思（結(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ù)基因（不計算A1CF、A2M），而樣本 2#也是同樣。
為了簡便，就用A2ML1基因了，其實結(jié)果是類似的。
原本應(yīng)該根據(jù)ln(reads for gene X/average for gene X)的值計算中位數(shù)，這里直接根據(jù)原始Read值計算，會有一定的差異。
求解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)
也就是說兩個樣本的標準化因子分別是4.16和32.66。
進行標準化變換
將原始矩陣按照樣本的不同除以各自的標準化因子，得下表：
可以發(fā)現(xiàn)，除了A2M基因有表達差異外，其他基因表達無明顯差異。
#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://statquest.org/video-index/
聯(lián)系客服
