在线一区久视频,久伊人免费在线视频,精品天堂一

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è)專門為了定制一套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)

$\text{cpm}_i=\frac{\text{read count of Gene}_i}{\text{total reads}/10^6}$

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è)縮放因子：

$d_i=\frac{S_{baseline}}{S_i}$

然后基于該縮放因子對(duì)特定的sample中的每個(gè)基因的read count進(jìn)行標(biāo)準(zhǔn)化（縮放）：

$x_i'=d_ix_i$

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=\log_2\frac{\text{sample}_i}{ref}$

LFC過(guò)大或過(guò)小都表示具有偏倚性，LFC越大表示reads數(shù)在sample_i中越高，即偏向sample_i；LFC越小表示reads數(shù)在ref中越高，即偏向ref

衡量reads數(shù)的量：read的幾何平均數(shù) (read geometric mean, RGM)

RGM越大表示基因reads越多，RGM越小表示基因reads越少

結(jié)合偏倚度、read數(shù)畫出散點(diǎn)圖：

偏倚度小、表達(dá)量居中的那些基因落在圖中的紅線附近

由參考代表基因集計(jì)算樣本的標(biāo)準(zhǔn)化因子：

對(duì)這些代表基因集計(jì)算加權(quán)平均數(shù)：

$\frac{\sum_i^n LFC_i\times \text{ReadCount}_i}{\sum_i^n \text{ReadCount}_i}$

該加權(quán)平均數(shù)就能代表該樣本的標(biāo)準(zhǔn)化因子，只是經(jīng)過(guò)了log變換，所以需要恢復(fù)為它的正值：

$\text{ScalingFactor}=2^{\text{weight-average}}$

2. 為什么說(shuō)FPKM和RPKM都錯(cuò)了？

2.1. FPKM和RPKM分別是什么

FPKM和RPKM分別是什么

RPKM是Reads Per Kilobase per Million的縮寫，它的計(jì)算方程非常簡(jiǎn)單：

$RPKM=\frac{10^6 \times n_r}{L \times N}$

FPKM是Fregments Per Kilobase per Million的縮寫，它的計(jì)算與RPKM極為類似，如下：

$FPKM=\frac{10^6 \times n_f}{L \times N}$

與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)該是：

$r_g=\frac{\text{mRNA}_g}{\text{mRNA}_{total}}$

則一個(gè)樣本中基因表達(dá)豐度的均值為

$r_{mean}=\frac{1}{g_{total}} \sum_g^G r_g = \frac{1}{g_{total}} \frac{\sum_g^G \text{mRNA}_g}{\text{mRNA}_{total}}$

而

$\sum_g^G \text{mRNA}_g=\text{mRNA}_{total}$

所以

$r_{mean}=\frac{1}{g_{total}}$

這個(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)度分別如下：

$r_{mean}$ 值是:

$r_{mean}=\frac 15=0.2$

如果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à)）：

$\frac{n_f}{L}$

第二部分的統(tǒng)計(jì)量是測(cè)序獲得的總有效Fregment數(shù)量的百萬(wàn)分之一：

$\frac{N}{10^6}$

這么一拆，就可以看出這個(gè)公式的問(wèn)題了：邏輯上根本說(shuō)不通嘛！

尤其是第二部分（N/10⁶），本來(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/10⁶并不能描述樣本的轉(zhuǎn)錄本總數(shù)量。N/10⁶的大小其實(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)了，稱之為TPM —— Transcripts Per Million，它的計(jì)算是：

$TPM=\frac{\frac{n_r \times read_l}{g_l}}{T}=\frac{n_r \times read_l \times 10^6}{g_l \times T}$