目錄
- 1. 標(biāo)準(zhǔn)化
- 1.1. House-keeping gene(s)
- 1.2. spike-in
- 1.3. CPM
- 1.4. TCS
- 1.5. Quantile
- 1.6. Median of Ration
- 1.7. TMM
- 2. 為什么說(shuō)FPKM和RPKM都錯(cuò)了絮短？
- 2.1. FPKM和RPKM分別是什么
- 2.2. 什么樣才算好的統(tǒng)計(jì)量
- 2.3. FPKM和RPKM犯的錯(cuò)
- 2.4. TPM是一個(gè)合適的選擇
1. 標(biāo)準(zhǔn)化
由于不同文庫(kù)測(cè)序深度不同穴墅，比較前當(dāng)然要進(jìn)行均一化！用總reads進(jìn)行均一化可能最簡(jiǎn)單忆家，其基于以下兩個(gè)基本假設(shè)：
- 絕大多數(shù)的gene表達(dá)量不變；
- 高表達(dá)量的gene表達(dá)量不發(fā)生改變顾彰；
但在轉(zhuǎn)錄組中缤沦，通常一小部分極高豐度基因往往會(huì)貢獻(xiàn)很多reads，如果這些“位高權(quán)重”的基因還是差異表達(dá)的疗涉，則會(huì)影響所有其它基因分配到的reads數(shù)拿霉，而且，兩個(gè)樣本總mRNA量完全相同的前提假設(shè)也過(guò)于理想了咱扣。那如何比較呢绽淘，各個(gè)方家使出渾身解數(shù)，有用中位數(shù)的闹伪，有用75分位數(shù)的沪铭，有用幾何平均數(shù)的，有用TMM(trimmed mean of Mvalues)的等等偏瓤，總之要找一個(gè)更穩(wěn)定的參考值杀怠。
1.1. House-keeping gene(s)
矯正的思路很簡(jiǎn)單，就是在變化的樣本中尋找不變的量
那么在不同RNA-seq樣本中厅克，那些是不變的量呢赔退？一個(gè)很容易想到的就是管家基因 (House-keeping gene(s))
那么 Human 常用的 House-keeping gene 怎么確定？
目前大家用的比較多的一個(gè)human housekeeping gene list 來(lái)源于下面這篇文章证舟，是2013年發(fā)表在 Cell系列的 Trends in Genetics 部分的一篇文章
1.2. spike-in
使用Housekeeping gene的辦法來(lái)進(jìn)行相對(duì)定量硕旗，這種辦法在一定程度上能夠解決我們遇到的問(wèn)題。但其實(shí)這種辦法有一個(gè)非常強(qiáng)的先驗(yàn)假設(shè)：housekeeping gene的表達(dá)量不怎么發(fā)生變化褪储。其實(shí)housekeeping gene list有幾千個(gè)卵渴，這幾千個(gè)基因有一定程度上的變化是有可能的
spike-in方法：在RNA-Seq建庫(kù)的過(guò)程中摻入一些預(yù)先知道序列信息以及序列絕對(duì)數(shù)量的內(nèi)參。這樣在進(jìn)行RNA-Seq測(cè)序的時(shí)候就可以通過(guò)不同樣本之間內(nèi)參（spike-in）的量來(lái)做一條標(biāo)準(zhǔn)曲線鲤竹，就可以非常準(zhǔn)確地對(duì)不同樣本之間的表達(dá)量進(jìn)行矯正
比較常用的spike-in類型：ERCC Control RNA
ERCC = External RNA Controls Consortium
ERCC就是一個(gè)專門(mén)為了定制一套spike-in RNA而成立的組織浪读，這個(gè)組織早在2003年的時(shí)候就已經(jīng)宣告成立昔榴。主要的工作就是設(shè)計(jì)了一套非常好用的spike-in RNA，方便microarray碘橘，以及RNA-Seq進(jìn)行內(nèi)參定量
1.3. CPM
CPM(count-per-million)
1.4. TCS (Total Count Scaling)
簡(jiǎn)單來(lái)說(shuō)互订，就是找出多個(gè)樣本中l(wèi)ibrary size為中位數(shù)的樣本，作為參考樣本痘拆，將所有的樣本的library size按比例縮放到參考樣本的水平
選擇一個(gè)library size為中位數(shù)的sample仰禽，以它為baseline，計(jì)算出其它sample對(duì)于baseline的normalization factor纺蛆，即一個(gè)縮放因子：
然后基于該縮放因子對(duì)特定的sample中的每個(gè)基因的read count進(jìn)行標(biāo)準(zhǔn)化（縮放）：
1.5. Quantile
簡(jiǎn)單來(lái)說(shuō)吐葵，就是排序后求平均，然后再回序
在R里面桥氏，推薦用preprocessCore 包來(lái)做quantile normalization温峭，不需要自己造輪子啦！
但是需要明白什么時(shí)候該用quantile normalization字支，什么時(shí)候不應(yīng)該用凤藏，就復(fù)雜很多了
1.6. Median of Ratio (DESeq2)
該方法基于的假設(shè)是，即使處在不同條件下的不同個(gè)樣本堕伪，大多數(shù)基因的表達(dá)是不存在差異的揖庄，實(shí)際存在差異的基因只占很小的部分那么我們只需要將這些穩(wěn)定的部分找出來(lái)，作為標(biāo)準(zhǔn)化的內(nèi)參欠雌，依據(jù)內(nèi)參算出各個(gè)樣本的標(biāo)準(zhǔn)化因子
（1）對(duì)每個(gè)基因計(jì)算幾何平均數(shù)蹄梢，得到一個(gè)假設(shè)的參考樣本(pseudo-reference sample)
| gene | sampleA | sampleB | pseudo-reference sample |
|---|---|---|---|
| EF2A | 1489 | 906 | sqrt(1489 * 906) = 1161.5 |
| ABCD1 | 22 | 13 | sqrt(22 * 13) = 17.7 |
| … | … | … | … |
（2）對(duì)每個(gè)樣本的每個(gè)基因?qū)τ趨⒖紭颖居?jì)算Fold Change
| gene | sampleA | sampleB | pseudo-reference sample | ratio of sampleA/ref | ratio of sampleB/ref |
|---|---|---|---|---|---|
| EF2A | 1489 | 906 | 1161.5 | 1489/1161.5 = 1.28 | 906/1161.5 = 0.78 |
| ABCD1 | 22 | 13 | 16.9 | 22/16.9 = 1.30 | 13/16.9 = 0.77 |
| MEFV | 793 | 410 | 570.2 | 793/570.2 = 1.39 | 410/570.2 = 0.72 |
| BAG1 | 76 | 42 | 56.5 | 76/56.5 = 1.35 | 42/56.5 = 0.74 |
| MOV10 | 521 | 1196 | 883.7 | 521/883.7 = 0.590 | 1196/883.7 = 1.35 |
| … | … | … | … | … | … |
（3）獲取每個(gè)樣本中Fold Change的中位數(shù)，我們就得到了非DE基因代表的Fold Change桨昙，該基因就是我們選擇的該樣本的內(nèi)參基因检号，它的Fold Change就是該樣本的標(biāo)準(zhǔn)化因子
normalization_factor_sampleA <- median(c(1.28, 1.3, 1.39, 1.35, 0.59))
normalization_factor_sampleB <- median(c(0.78, 0.77, 0.72, 0.74, 1.35))
1.7. TMM (Trimmed Mean of M value, edgeR)
該方法的思想與DESeq2的Median of Ratio相同，假設(shè)前提都是：大多數(shù)基因的表達(dá)是不存在差異的
它與DESeq2的不同之處在于對(duì)內(nèi)參的選擇上：
DESeq2選擇一個(gè)內(nèi)參基因蛙酪，它的Ratio/Fold-Change就是標(biāo)準(zhǔn)化因子
edgeR選擇一組內(nèi)參基因集合，它們對(duì)標(biāo)準(zhǔn)化因子的計(jì)算均有貢獻(xiàn)：加權(quán)平均
（1）移除所有未表達(dá)基因
（2）從眾多樣本中找出一個(gè)數(shù)據(jù)趨勢(shì)較為平均的樣本作為參考樣本
對(duì)所有樣本求總Read數(shù)翘盖；
各樣本除以各自的總Read數(shù)桂塞，得到修正Read數(shù)；
求出各自樣本修正Read數(shù)的Q3值（第3個(gè)四分位數(shù)）馍驯；
所有的Q3值求平均阁危，與平均Q3相差最小的樣本即是參考樣本。
（3）找出每個(gè)樣本中的代表基因集汰瘫，參考這些代表基因集的fold change狂打，計(jì)算出該樣本的標(biāo)準(zhǔn)化因子
尋找樣本的代表基因集：依據(jù)基因的偏倚程度和Reads數(shù)大小選出——偏倚程度小、reads數(shù)居中的基因
- 衡量偏倚度的量：LFC (log fold change)
LFC過(guò)大或過(guò)小都表示具有偏倚性混弥，LFC越大表示reads數(shù)在samplei中越高趴乡，即偏向samplei；LFC越小表示reads數(shù)在ref中越高，即偏向ref
- 衡量reads數(shù)的量：read的幾何平均數(shù) (read geometric mean, RGM)
RGM越大表示基因reads越多晾捏，RGM越小表示基因reads越少
結(jié)合偏倚度蒿涎、read數(shù)畫(huà)出散點(diǎn)圖：
偏倚度小、表達(dá)量居中的那些基因落在圖中的紅線附近
由參考代表基因集計(jì)算樣本的標(biāo)準(zhǔn)化因子：
對(duì)這些代表基因集計(jì)算加權(quán)平均數(shù)：
該加權(quán)平均數(shù)就能代表該樣本的標(biāo)準(zhǔn)化因子惦辛，只是經(jīng)過(guò)了log變換劳秋，所以需要恢復(fù)為它的正值：
2. 為什么說(shuō)FPKM和RPKM都錯(cuò)了？
2.1. FPKM和RPKM分別是什么
FPKM和RPKM分別是什么
RPKM是Reads Per Kilobase per Million的縮寫(xiě)胖齐，它的計(jì)算方程非常簡(jiǎn)單：
FPKM是Fregments Per Kilobase per Million的縮寫(xiě)玻淑，它的計(jì)算與RPKM極為類似，如下：
與RPKM唯一的區(qū)別為：F是fragments呀伙，R是reads补履，如果是PE(Pair-end)測(cè)序，每個(gè)fragments會(huì)有兩個(gè)reads区匠，F(xiàn)PKM只計(jì)算兩個(gè)reads能比對(duì)到同一個(gè)轉(zhuǎn)錄本的fragments數(shù)量干像，而RPKM計(jì)算的是可以比對(duì)到轉(zhuǎn)錄本的reads數(shù)量而不管PE的兩個(gè)reads是否能比對(duì)到同一個(gè)轉(zhuǎn)錄本上。如果是SE(single-end)測(cè)序驰弄，那么FPKM和RPKM計(jì)算的結(jié)果將是一致的麻汰。
這兩個(gè)量的計(jì)算方式的目的是為了解決計(jì)算RNA-seq轉(zhuǎn)錄本豐度時(shí)的兩個(gè)bias：
- 相同表達(dá)豐度的轉(zhuǎn)錄本，往往會(huì)由于其基因長(zhǎng)度上的差異戚篙，導(dǎo)致測(cè)序獲得的Read（Fregment）數(shù)不同五鲫。總的來(lái)說(shuō)岔擂，越長(zhǎng)的轉(zhuǎn)錄本位喂，測(cè)得的Read（Fregment）數(shù)越多；
- 由測(cè)序文庫(kù)的不同大小而引來(lái)的差異乱灵。即同一個(gè)轉(zhuǎn)錄本塑崖，其測(cè)序深度越深，通過(guò)測(cè)序獲得的Read（Fregment）數(shù)就越多痛倚。
2.2. 什么樣才算好的統(tǒng)計(jì)量
首先规婆，到底什么是RNA轉(zhuǎn)錄本的表達(dá)豐度這個(gè)問(wèn)題
對(duì)于樣本X，其有一個(gè)基因g被轉(zhuǎn)錄了mRNA_g次蝉稳，同時(shí)樣本X中所有基因的轉(zhuǎn)錄總次數(shù)假定是mRNA_total, 那么正確描述基因g轉(zhuǎn)錄豐度的值應(yīng)該是：
則一個(gè)樣本中基因表達(dá)豐度的均值為
而
所以
這個(gè)期望值竟然和測(cè)序狀態(tài)無(wú)關(guān)抒蚜！僅僅由樣本中基因的總數(shù)所決定的
也就是說(shuō)，對(duì)于同一個(gè)物種耘戚，不管它的樣本是哪種組織（正常的或病變的）嗡髓，也不管有多少個(gè)不同的樣本，只要它們都擁有相同數(shù)量的基因收津，那么它們的r_mean都將是一致的
由于上面的結(jié)果是在理論情況下推導(dǎo)出來(lái)的饿这，實(shí)際上我們無(wú)法直接計(jì)算這個(gè)r浊伙，那么我們可以嘗試通過(guò)其他方法來(lái)近似估計(jì)r，只要這些近似統(tǒng)計(jì)量可以隱式地包含這一恒等關(guān)系即可
2.3. FPKM和RPKM犯的錯(cuò)
實(shí)際數(shù)據(jù)來(lái)證明
假定有兩個(gè)來(lái)自同一個(gè)個(gè)體不同組織的樣本X和Y蛹稍，這個(gè)個(gè)體只有5個(gè)基因吧黄，分別為A、B唆姐、C拗慨、D和E它們的長(zhǎng)度分別如下：
值是:
如果FPKM或RPKM是一個(gè)合適的統(tǒng)計(jì)量的話，那么奉芦，樣本X和Y的平均FPKM（或RPKM）應(yīng)該相等缩麸。
以下這個(gè)表格列出的分別是樣本X和Y在這5個(gè)基因分別比對(duì)上的fregment數(shù)和各自總的fregment數(shù)量：
可以算出：樣本X在這5個(gè)基因上的FPKM均值FPKM_mean = 5,680;樣本Y在這5個(gè)基因上的FPKM均值FPKM_mean = 161,840
很明顯善延，它們根本不同氢烘，而且差距相當(dāng)大
究竟為什么會(huì)有如此之大的差異执桌？
可以從其公式上找到答案
首先，我們可以把FPKM的計(jì)算式拆分成如下兩部分先巴。
第一部分的統(tǒng)計(jì)量是對(duì)一個(gè)基因轉(zhuǎn)錄本數(shù)量的一個(gè)等價(jià)描述（雖然嚴(yán)格來(lái)講也沒(méi)那么等價(jià)）：
第二部分的統(tǒng)計(jì)量是測(cè)序獲得的總有效Fregment數(shù)量的百萬(wàn)分之一：
這么一拆其爵，就可以看出這個(gè)公式的問(wèn)題了：邏輯上根本說(shuō)不通嘛！
尤其是第二部分（N/106）伸蚯，本來(lái)式子的第一部分是為了描述一個(gè)基因的轉(zhuǎn)錄本數(shù)量摩渺，那么正常來(lái)講，第二部分就應(yīng)該是樣本的轉(zhuǎn)錄本總數(shù)量（或至少是其總數(shù)量的等價(jià)描述）才能形成合理的比例關(guān)系剂邮，而且可以看出來(lái)FPKM/RPMK是有此意的摇幻，這本來(lái)就是這個(gè)統(tǒng)計(jì)量的目的。
但是N/106并不能描述樣本的轉(zhuǎn)錄本總數(shù)量挥萌。N/106的大小其實(shí)是由RNA-seq測(cè)序深度所決定的绰姻，并且是一個(gè)和總轉(zhuǎn)錄本數(shù)量無(wú)直接線性關(guān)系的統(tǒng)計(jì)量——N與總轉(zhuǎn)錄本數(shù)量之間的關(guān)系還受轉(zhuǎn)錄本的長(zhǎng)度分布所決定，而這個(gè)分布往往在不同樣本中是有差異的引瀑！
2.4. TPM是一個(gè)合適的選擇
這個(gè)統(tǒng)計(jì)量在2012年所發(fā)表的一篇討論RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就被提出來(lái)了狂芋，稱之為T(mén)PM —— Transcripts Per Million，它的計(jì)算是：
簡(jiǎn)單計(jì)算之后我們就可以發(fā)現(xiàn)TPM的均值是一個(gè)獨(dú)立于樣本之外的恒定值憨栽，它等于：
這個(gè)值剛好是r_mean的一百萬(wàn)倍银酗，滿足等價(jià)描述的關(guān)系。
參考資料：
(1) 孟浩巍《生物信息學(xué)100個(gè)基礎(chǔ)問(wèn)題 —— 第38題 當(dāng)轉(zhuǎn)錄組普遍變化時(shí)RNA-Seq怎么進(jìn)行分析(1)徒像？》
(2) 孟浩巍《生物信息學(xué)100個(gè)基礎(chǔ)問(wèn)題 —— 第38題 當(dāng)轉(zhuǎn)錄組普遍變化時(shí)RNA-Seq怎么進(jìn)行分析(2)？》
(3) 【生信菜鳥(niǎo)團(tuán)】quantile normalization到底對(duì)數(shù)據(jù)做了什么蛙讥？
