目錄
1. 標準化
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. 為什么說FPKM和RPKM都錯了柜候?
2.1. FPKM和RPKM分別是什么
2.2. 什么樣才算好的統(tǒng)計量
2.3. FPKM和RPKM犯的錯
2.4. TPM是一個合適的選擇
1. 標準化
目前肌括,轉(zhuǎn)錄組測序(RNA-seq)分析是非常成熟的研究手段枫浙,有眾多分析工具和方法供大家使用朦拖,其中我磁,對基因或轉(zhuǎn)錄本的讀段數(shù)目(read count)進行歸一化是一個非常重要的分析過程孽文,如何對基因區(qū)域進行準確的定量和歸一化,是大家十分關心的核心問題之一夺艰。
無疑芋哭,轉(zhuǎn)錄組測序雙端數(shù)據(jù)分析中,目前FPKM是最常用的歸一化方法郁副,那FPKM歸一化方法是最準確的嗎减牺?隨著生物信息分析技術(shù)的快速發(fā)展,F(xiàn)PKM或許已經(jīng)是“明日黃花”存谎。
歸一化的基本背景
總的來說拔疚,傳統(tǒng)的轉(zhuǎn)錄組定量方法是相對定量,一個基因的定量結(jié)果很大程度上會受到基因的長度和測序深度的影響既荚≈墒В基因長度越長、測序深度越高恰聘,得到該基因的read counts就越多句各,相對表達水平也越高。所以晴叨,
在進行下游分析的時候凿宾,例如聚類、主成分分析篙螟,如果不進行數(shù)據(jù)歸一化直接使用原始read count菌湃,簡直就是耍流氓。
因此遍略,表達量歸一化的精確計算需要同時考慮基因長度惧所、測序深度等信息。
用總reads進行均一化可能最簡單绪杏,其基于以下兩個基本假設:
1下愈、絕大多數(shù)的gene表達量不變;
2蕾久、高表達量的gene表達量不發(fā)生改變势似;
但在轉(zhuǎn)錄組中,通常一小部分極高豐度基因往往會貢獻很多reads,如果這些“位高權(quán)重”的基因還是差異表達的履因,則會影響所有其它基因分配到的reads數(shù)障簿,而且,兩個樣本總mRNA量完全相同的前提假設也過于理想了栅迄。那如何比較呢站故,各個方家使出渾身解數(shù),有用中位數(shù)的毅舆,有用75分位數(shù)的西篓,有用幾何平均數(shù)的,有用TMM(trimmed mean of Mvalues)的等等憋活,總之要
找一個更穩(wěn)定的參考值岂津。
除FPKM/RPKM外,標準化方法有:
1.1. House-keeping gene(s)
矯正的思路很簡單悦即,就是在變化的樣本中尋找不變的量
那么在不同RNA-seq樣本中吮成,那些是不變的量呢?一個很容易想到的就是管家基因(House-keeping gene(s))
那么 Human 常用的 House-keeping gene 怎么確定辜梳?
目前大家用的比較多的一個human housekeeping gene list 來源于下面這篇文章赁豆,是2013年發(fā)表在 Cell系列的 Trends in Genetics 部分的一篇文章
1.2. spike-in
使用Housekeeping gene的辦法來進行相對定量,這種辦法在一定程度上能夠解決我們遇到的問題冗美。但其實這種辦法有一個非常強的先驗假設:housekeeping gene的表達量不怎么發(fā)生變化。其實housekeeping gene list有幾千個析二,這幾千個基因有一定程度上的變化是有可能的粉洼。
spike-in方法:在RNA-Seq建庫的過程中摻入一些預先知道序列信息以及序列絕對數(shù)量的內(nèi)參。這樣在進行RNA-Seq測序的時候就可以通過不同樣本之間內(nèi)參(spike-in)的量來做一條標準曲線叶摄,就可以非常準確地對不同樣本之間的表達量進行矯正属韧。
spike-in是已知的其他物種基因組,可以對基因表達數(shù)據(jù)進行絕對值的矯正
個人認為蛤吓,F(xiàn)PKM/RPKM/TPM? 都是基因的相對表達值
DESeq2會對測序深度進行矯正宵喂,將普遍高表達的樣本A看作是測序深度過高所導致,從而影響差異基因的篩選
RNA的spike-in是一組序列已知的RNA轉(zhuǎn)錄本会傲,目前普遍使用的是ERCC( The External RNA Control Consortium)搞出來的一組RNA序列锅棕。當然,也可以使用序列相似度較低的物種的序列作為spike-in淌山。如兩種常見的酵母裸燎,S.pombe和S. cerevisiae,他們的序列可以作為彼此的spike-in進行后續(xù)矯正泼疑。
通過在樣品制備過程中德绿,混入指定數(shù)量的spike-in,我們就可以知道不同樣本中的基因絕對比表達值。
如等細胞數(shù)的樣本A和樣本B移稳,在每個樣本中蕴纳,我加入了等量的spike-in。最后分析發(fā)現(xiàn)个粱,spike-in占樣本A的1%古毛,但是占樣本B的5%。這表明樣本A的RNA表達量也許普遍比樣本B的表達量高五倍左右几蜻。
下面是一個簡單的草圖喇潘,希望可以幫助理解,左邊是細胞梭稚,右邊是RNA
比較常用的spike-in類型:ERCC Control RNA
ERCC = External RNA Controls Consortium
ERCC就是一個專門為了定制一套spike-in RNA而成立的組織颖低,這個組織早在2003年的時候就已經(jīng)宣告成立。主要的工作就是設計了一套非常好用的spike-in RNA弧烤,方便microarray忱屑,以及RNA-Seq進行內(nèi)參定量
1.3. CPM
CPM(count-per-million)
1.4. TCS (Total Count Scaling)
簡單來說,就是找出多個樣本中l(wèi)ibrary size為中位數(shù)的樣本暇昂,作為參考樣本莺戒,將所有的樣本的library size按比例縮放到參考樣本的水平
選擇一個library size為中位數(shù)的sample,以它為baseline急波,計算出其它sample對于baseline的normalization factor从铲,即一個縮放因子:
然后基于該縮放因子對特定的sample中的每個基因的read count進行標準化(縮放):
1.5. Quantile
簡單來說,就是排序后求平均澄暮,然后再回序
在R里面名段,推薦用preprocessCore 包來做quantile normalization,不需要自己造輪子啦泣懊!
但是需要明白什么時候該用 quantile normalization伸辟,什么時候不應該用,就復雜很多了
1.6. Median of Ratio (DESeq2)
該方法基于的假設是馍刮,即使處在不同條件下的不同個樣本信夫,大多數(shù)基因的表達是不存在差異的,實際存在差異的基因只占很小的部分卡啰,?那么我們?只需要將這些穩(wěn)定的部分找出來静稻,作為標準化的內(nèi)參,依據(jù)內(nèi)參算出各個樣本的標準化因子
參考?
(1)對每個基因計算幾何平均數(shù)碎乃,得到一個假設的參考樣本(pseudo-reference sample)
genesampleAsampleBpseudo-reference sample
EF2A1489906sqrt(1489 * 906) = 1161.5
ABCD12213sqrt(22 * 13) = 17.7
…………
(2)對每個樣本的每個基因對于參考樣本計算Fold Change
genesampleAsampleBpseudo-reference sampleratio of sampleA/refratio of sampleB/ref
EF2A14899061161.51489/1161.5 = 1.28906/1161.5 = 0.78
ABCD1221316.922/16.9 = 1.3013/16.9 = 0.77
MEFV793410570.2793/570.2 = 1.39410/570.2 = 0.72
BAG1764256.576/56.5 = 1.3542/56.5 = 0.74
MOV105211196883.7521/883.7 = 0.5901196/883.7 = 1.35
………………
(3)獲取每個樣本中Fold Change的中位數(shù)姊扔,我們就得到了非DE基因代表的Fold Change,該基因就是我們選擇的該樣本的內(nèi)參基因梅誓,它的Fold Change就是該樣本的標準化因子
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ù)基因的表達是不存在差異的佛南。
它與DESeq2的不同之處在于對內(nèi)參的選擇上:
DESeq2選擇一個內(nèi)參基因,它的 Ratio/Fold-Change 就是標準化因子嵌言。
edgeR?選擇一組內(nèi)參基因集合嗅回,它們對標準化因子的計算均有貢獻:加權(quán)平均。
(1)移除所有未表達基因
(2)從眾多樣本中找出一個數(shù)據(jù)趨勢較為平均的樣本作為參考樣本
對所有樣本求總Read數(shù)摧茴;
各樣本除以各自的總Read數(shù)绵载,得到修正Read數(shù);
求出各自樣本修正Read數(shù)的Q3值(第3個四分位數(shù))苛白;
所有的Q3值求平均娃豹,與平均Q3相差最小的樣本即是參考樣本。
(3)找出每個樣本中的代表基因集购裙,參考這些代表基因集的fold change懂版,計算出該樣本的標準化因子
尋找樣本的代表基因集:依據(jù)基因的偏倚程度和Reads數(shù)大小選出——偏倚程度小、reads數(shù)居中的基因
衡量偏倚度的量:LFC (log fold change)
LFC過大或過小都表示具有偏倚性躏率,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ù)畫出散點圖:
偏倚度小夯到、表達量居中的那些基因落在圖中的紅線附近
由參考代表基因集計算樣本的標準化因子:
對這些代表基因集計算加權(quán)平均數(shù):
該加權(quán)平均數(shù)就能代表該樣本的標準化因子嚷缭,只是經(jīng)過了log變換,所以需要恢復為它的正值:
2. 為什么說FPKM和RPKM都錯了耍贾?
2.1. FPKM和RPKM分別是什么
在RNA-Seq中峭状,最簡單的定量基因表達量(gene expression)的方法就是將RNA-Seq數(shù)據(jù)比對到相應的參考序列上時,會有比對到各個基因的read數(shù)量逼争,稱為raw read counts。
但是如果要比較不同樣本中基因的表達量劝赔,光有raw counts是遠遠不夠的誓焦,因為raw cread counts受到很多因素的影響,如目標基因的轉(zhuǎn)錄本長度(transcript length)着帽、總的有效比對的read數(shù)量(即測序深度 sequencing depth)以及測序的偏差(sequencing bias)等等杂伟,這些因素是如何影響raw read counts的后面會有解釋。
那么為了將不同樣本的基因表達量歸一化到一個能夠量化比較的標準上仍翰,科學家們采取的措施是將raw counts同時除以目標基因的外顯子長度之和(也就是目標基因轉(zhuǎn)錄本長度)和總的有效比對的read總數(shù)赫粥。這就是RPKM的定義
? RPKM = (10^6 * nr)?/ (L * N)
其中 nr?代表比對至目標基因的read數(shù)量;
L代表目標基因的外顯子長度之和除以1000予借,單位是Kb越平,不是bp频蛔;
N是總的有效比對至基因組的reads數(shù)量。
注意這里的?nr:在single-end測序中秦叛,一個read就是一個read晦溪。而在pair-end測序中,若一對paired-read 都比對上了挣跋,當做兩個read三圆;若只有一個read比對上,另一個未比對上避咆,當做一個read計算舟肉。
類似的,F(xiàn)PKM的定義如下?
?FPKM = (10^6 * nf)?/ (L * N)
?其中?nf?代表比對至目標基因的fragment數(shù)量查库;
L代表目標基因的外顯子長度之和除以1000路媚,單位是Kb,不是bp膨报;
N是總的有效比對至基因組的fragment數(shù)量磷籍。
?注意這里的?nf:在single-end測序中,F(xiàn)PKM將read當做fragment計算现柠,此時FPKM和RPKM是相同的院领。而在pair-end測序? 中,? 若一堆paired-read 都比對上了够吩,當做一個fragment比然。
RPKM是Reads Per Kilobase per Million的縮寫,它的計算方程非常簡單:
FPKM是Fregments Per Kilobase per Million的縮寫周循,它的計算與RPKM極為類似强法,如下:
與RPKM唯一的區(qū)別為:F是fragments,R是reads湾笛,
如果是PE(Pair-end)測序饮怯,每個fragments會有兩個reads,F(xiàn)PKM只計算兩個reads能比對到同一個轉(zhuǎn)錄本的fragments數(shù)量嚎研,而RPKM計算的是可以比對到轉(zhuǎn)錄本的reads數(shù)量而不管PE的兩個reads是否能比對到同一個轉(zhuǎn)錄本上蓖墅。
如果是SE(single-end)測序,那么FPKM和RPKM計算的結(jié)果將是一致的临扮。
這兩個量的計算方式的目的是為了解決計算RNA-seq轉(zhuǎn)錄本豐度時的兩個bias:
1论矾、相同表達豐度的轉(zhuǎn)錄本,往往會由于其基因長度上的差異杆勇,導致測序獲得的Read(Fregment)數(shù)不同贪壳。總的來說蚜退,越長的轉(zhuǎn)錄本闰靴,測得的Read(Fregment)數(shù)越多彪笼;
2、由測序文庫的不同大小而引來的差異传黄。即同一個轉(zhuǎn)錄本杰扫,其測序深度越深,通過測序獲得的Read(Fregment)數(shù)就越多膘掰。
2.2. 什么樣才算好的統(tǒng)計量
首先章姓,RNA轉(zhuǎn)錄本的表達豐度 是什么?识埋?
事實上凡伊,對于任何一個制備好的樣本,它上面任何一個基因的表達量(或者說豐度)窒舟,都將已是一個客觀存在的值系忙,這個值是不管你改變了多少測序環(huán)境都不會變的。而且總共有多少個不同的基因在表達惠豺,實際上也已經(jīng)是客觀定好了的银还。一旦我們開始以這樣一種“先知”的形式來理解的時候,有趣的事情就開始出現(xiàn)了洁墙。
對于樣本X蛹疯,其有一個基因G被轉(zhuǎn)錄了mRNA_g次,同時樣本X中所有基因的轉(zhuǎn)錄總次數(shù)假定是mRNA_total, 那么正確描述基因g轉(zhuǎn)錄豐度的值應該是:
則一個樣本中基因表達豐度的均值為
而
而
所以
這個期望值竟然和測序狀態(tài)無關热监!僅僅由樣本中基因的總數(shù)所決定的
也就是說捺弦,對于同一個物種,不管它的樣本是哪種組織(正常的或病變的)孝扛,也不管有多少個不同的樣本列吼,只要它們都擁有相同數(shù)量的基因,那么它們的r_mean都將是一致的苦始!?
由于上面的結(jié)果是在理論情況下推導出來的寞钥,實際上我們無法直接計算這個r,那么我們可以嘗試通過其他方法來近似估計r陌选,只要這些近似統(tǒng)計量可以隱式地包含這一恒等關系即可凑耻。
(那么)現(xiàn)在,我們回過頭來看看FPKM和RPKM的計算式柠贤,就會發(fā)現(xiàn)它們根本做不到。
2.3. FPKM和RPKM犯的錯
實際數(shù)據(jù)來證明
假定有兩個來自同一個個體但不同組織的樣本(比如 同一個病人的 肝組織 和 脾組織 )X和Y类缤,這個個體只有5個基因臼勉,分別為A、B餐弱、C宴霸、D和E它們的長度分別如下:
值是:
如果FPKM或RPKM是一個合適的統(tǒng)計量的話囱晴,那么,樣本X和Y的平均FPKM(或RPKM)應該相等瓢谢。
以下這個表格列出的分別是樣本X和Y在這5個基因分別比對上的fregment數(shù)和各自總的fregment數(shù)量:
可以算出:
樣本X在這5個基因上的FPKM均值FPKM_mean = 5,680;
樣本Y在這5個基因上的FPKM均值FPKM_mean = 161,840
很明顯畸写,它們根本不同,而且差距相當大
究竟為什么會有如此之大的差異氓扛?
可以從其公式上找到答案
首先枯芬,我們可以把FPKM的計算式拆分成如下兩部分。
第一部分的統(tǒng)計量是對一個基因轉(zhuǎn)錄本數(shù)量的一個等價描述(雖然嚴格來講也沒那么等價):
第二部分的統(tǒng)計量是測序獲得的總有效Fregment數(shù)量的百萬分之一:
這么一拆采郎,就可以看出這個公式的問題了:邏輯上根本說不通嘛千所!
尤其是第二部分(N/106),本來式子的第一部分是為了描述一個基因的轉(zhuǎn)錄本數(shù)量蒜埋,那么正常來講淫痰,第二部分就應該是樣本的轉(zhuǎn)錄本總數(shù)量(或至少是其總數(shù)量的等價描述)才能形成合理的比例關系,而且可以看出來FPKM/RPMK是有此意的整份,這本來就是這個統(tǒng)計量的目的待错。
但是N/106并不能描述樣本的轉(zhuǎn)錄本總數(shù)量。
N/106的大小其實是由RNA-seq測序深度所決定的烈评,并且是一個和總轉(zhuǎn)錄本數(shù)量無直接線性關系的統(tǒng)計量——N與總轉(zhuǎn)錄本數(shù)量之間的關系還受轉(zhuǎn)錄本的長度分布所決定火俄,而這個分布往往在不同樣本中是有差異的!
比如础倍,雖然有些基因烛占,有效比對到它們的Fregment數(shù)是相等的,但總的來講長度越長的基因沟启,其被轉(zhuǎn)錄的次數(shù)就越少忆家。也就是說,N必須將各個被轉(zhuǎn)錄的基因的長度考慮進去才能正確描述總體的轉(zhuǎn)錄本數(shù)德迹!而FPKM(RPKM)顯然沒有做到這一點芽卿,這便是FPKM(RPKM)出錯的內(nèi)在原因。
FPKM和RPKM通過同時除以L(轉(zhuǎn)錄本長度)和N(有效比對的Read(Fragment)總數(shù))的辦法胳搞,最終將不同樣本(或者同個樣本在不同條件下)的轉(zhuǎn)錄本豐度歸一化到一個能夠進行量化比較的標準上卸例。
以上一切看起來都很合理
但是!<∫恪筷转!
既然說了測序獲得的read(fragment)受到基因長度的影響,RPKM和FPKM計算中也去除了目標基因長度的影響悬而,
但是除以N時沒有考慮到這個影響呜舒,N是總的有效比對的read(fragment),它同樣會受到各個轉(zhuǎn)錄基因長度(distribution of transcript lengths)的影響笨奠。
所以FPKM/RPKM是不準確的袭蝗。那么有沒有一個統(tǒng)計量能解決這個問題呢唤殴?有!那就是TPM
2.4. TPM是一個合適的選擇
這個統(tǒng)計量在2012年所發(fā)表的一篇討論RPKM的文章(RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.)中就被提出來了到腥,稱之為TPM —— Transcripts Per Million朵逝,它的計算是:
TPMi= {( nr/Lr?)*10^6 } / sum( nr/Lr+……..+ nm/Lm?)
nr:mapping到目標基因上的read數(shù);
Lr:目標基因的外顯子長度的總和乡范。
在一個樣本中一個基因的TPM:先對每個基因的read數(shù)用基因的長度進行校正配名,
之后再用校正后的這個基因read數(shù) (nr/Lr)? ?與 校正后的這個樣本的所有校正后的read數(shù)(sum( nr/Lr+……..+ nm/Lm?))求商。
沒錯篓足!TPM不是除以有效比對的read總數(shù)段誊,而是除以經(jīng)過基因長度歸一化后的有效比對的read總數(shù),即歸一化后的測序深度栈拖。
因此连舍,TPM在計算不同樣本的基因表達量比較時,是更加準確的統(tǒng)計量涩哟。
其中索赏,readl是比對至基因g的平均read長度,gl是基因g的外顯子長度之和(這里無需將其除以1000了)贴彼。
在不考慮比對的剪切的情況下潜腻,readl這個值往往都是一個固定值(如100bp或者150bp等),因此我們也可以將readl統(tǒng)一約掉器仗,那么分子就會蛻變成RPKM計算式的第一部分融涣,但留著會更合理。
這樣子精钮,整個統(tǒng)計量就很好理解了威鹿,分子是某個基因g的轉(zhuǎn)錄本數(shù)(的等價描述),分母則為樣本中總轉(zhuǎn)錄本的數(shù)量轨香,兩者的比值TPM——便是基因g的轉(zhuǎn)錄本豐度忽你!
簡單計算之后我們就可以發(fā)現(xiàn)?TPM的均值是一個獨立于樣本之外的恒定值,它等于:
這個值剛好是r_mean的一百萬倍臂容,滿足等價描述的關系科雳。
TPM值就是RPKM的百分比!E肌糟秘!
TPM的優(yōu)點是什么呢?很明顯球散,所有基因的TPM值加起來肯定是1M尿赚,因為百分比的總和就是1嘛,與樣本無關,各個樣本都可以保證TPM庫是一樣的吼畏,這樣比較更有意義!`业啤泻蚊!
最后還是貼上公式吧浅妆!
既然FPKM/RPKM是錯的拌阴,那為什么大家還在用,而且還真找到了(能被實驗所驗證)有價值的結(jié)果呢涕侈?
羹奉。而且我們都知道2008那篇關于RPKM的文章用實驗結(jié)果證明了秒旋,RPKM是一個合適的統(tǒng)計量,符合qPCR的驗證結(jié)果诀拭。
但歸根到底迁筛,眼見未必為實,實驗是表象的耕挨,我們更應該從其本質(zhì)意義和原理上去考慮细卧。FPKM/RPKM之所以看起來會是一個合適的值,我想主要原因有二:
其一筒占,它們和TPM之間存在一定的正比關系贪庙。這在通過它們各自的數(shù)學計算式就可以看出來(以RPKM的計算為例):
而且在同一個樣本內(nèi)部由于T,N 和 read_l實際上都是定值翰苫,因此同個樣本內(nèi)止邮,的RPKM和TPM是可以恒等轉(zhuǎn)換的,只是在 樣本 與 樣本之間就不行奏窑,因為它們之間的轉(zhuǎn)換因子(百分比)大小不一了导披!
如以上例子,對于樣本X良哲,TPM轉(zhuǎn)換到RPKM的轉(zhuǎn)換因子(百分比)為:0.0284盛卡,但樣本Y的轉(zhuǎn)換因子(百分比)為:0.8092。
而由于這個標準的改變筑凫,會導致其原本所要描述的“轉(zhuǎn)錄本豐度”變得不可比較滑沧。
然,這其實不是最根本的原因巍实,更本質(zhì)的原因是滓技,這個轉(zhuǎn)換會對本來已經(jīng)正確標準化了的結(jié)果,再次做了一次無意義的不等變換棚潦,最終導致了結(jié)果不可解釋令漂。
如何理解呢,后文會有補充,這里先簡單說一下:這個數(shù)學轉(zhuǎn)換式子僅是告訴了我們這樣子來計算是可行的叠必,但是在RNA-seq的實際應用場景中荚孵,它其實是無生物(或物理)意義的;
其二纬朝,實驗驗證的精度是有限的收叶,常用的qPCR也只能給出定性的比較結(jié)果,而且實驗驗證也未必總能成功共苛。
5. 總結(jié)
現(xiàn)在回過頭來總結(jié)一下判没。事實上,F(xiàn)PKM/RPKM最大的問題就在于其無意義性隅茎。
我們所要表達出來的任何統(tǒng)計量澄峰,它的變化都應該要能對應到其物理或生物過程中的變化,如果做不到這一點辟犀,那么這個統(tǒng)計量往往都是無意義的俏竞,用它得到的結(jié)果就算看起來符合預期也只不過是數(shù)值上的巧合,本質(zhì)上是不可解釋的踪蹬。
FPKM/RPKM的分母(N/106)并不具有任何形式的生物意義胞此,它所能表達出來的這個量,只能代表測序深度的變化跃捣,而無法作為表達生物過程的量漱牵,比如無法代表(等價代表)樣本中轉(zhuǎn)錄本的總量。
一個統(tǒng)計量該如何計算疚漆,說到底都只是一個“術(shù)”的問題酣胀,而我們應該盡可能在接近其本質(zhì)意義的地方去確定。
FPKM/RPKM 和 TPM 存在一定的正比關系娶聘,因此我們在使用FPKM/RPK時闻镶,有些時候確實也能獲得可以被實驗所驗證的“好”結(jié)果,但其實它是一個橡皮筋丸升,它的單位刻度是會隨著樣本的不同而改變的铆农。
到頭來,樣本之間的差異比較實際上也只是在不同的標準下進行的狡耻,這樣的比較就算得到了所謂的“好”結(jié)果墩剖,那又有什么意義,根本就是個錯誤的東東夷狰。想想就是由于這種統(tǒng)計量岭皂,我們一定已經(jīng)獲得了許多的假陽性結(jié)果,同時也肯定錯過了許多本來真正有意義的差異沼头,真是彎路走盡也不知爷绘,而且還浪費了大堆的心情和時間书劝。
文章:《A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis.?Briefings in Bioinformatics.10.1093/bib/bbs046.》對7種主要的RNA-seq標準化方法(但不包含本文提到的TPM)做了一個詳細的比較,它用實際結(jié)果進行比較(不同于本文所用的數(shù)學方式)土至。也得出了RPKM這個統(tǒng)計量是應該被摒棄的結(jié)論购对,因為它所描述出來的結(jié)果是最不合理的,其實所有類似于RPKM(包括FPKM)的統(tǒng)計量在描述轉(zhuǎn)錄本豐度的時候都應該被摒棄陶因。
TPM的缺陷
轉(zhuǎn)錄組數(shù)據(jù)定量歸一化方法有很多洞斯,經(jīng)典的RPKM/FPKM因其本身固有的缺陷,越來越多的學者采用TPM這一冉冉升起的新星坑赡,大有取而代之的勢頭。
其實么抗,不管TPM毅否、RPKM還是FPKM都是相對定量的歸一化方法。
定量的前提需要樣本的表達量變化比較穩(wěn)定蝇刀,不能出現(xiàn)整體的上調(diào)或下調(diào)螟加,或者個別基因表達量發(fā)生劇烈變化,否則定量歸一化方法可能會失效吞琐。
另外捆探,傳統(tǒng)轉(zhuǎn)錄組測序在信息分析過程中通常不會去除 duplicate reads,
因為根本不知道這些reads是多個表達拷貝的結(jié)果站粟,還是文庫構(gòu)建中PCR duplication產(chǎn)生的結(jié)果黍图。
為了在源頭實現(xiàn)精確定量,可以在reads中追加序列唯一的UMI(Unique Melocular Identifier)分子標簽奴烙,這樣攜帶相同UMI標簽的reads認為是duplicate reads助被,保留一條質(zhì)量值最高的read即可,從而實現(xiàn)較為準確的絕對定量切诀。
6.如何實現(xiàn)絕對定量揩环。
轉(zhuǎn)錄組測序的終極目的是基于表達量來發(fā)掘背后的生物學問題,問題是表達量真的準確嗎幅虑?
序列偏好丰滑、cDNA反轉(zhuǎn)錄、文庫PCR擴增倒庵、測序擴增等都會增加解讀數(shù)據(jù)的難度褒墨。如何解釋常規(guī)轉(zhuǎn)錄組數(shù)據(jù)面需要解決的問題比較多,不僅僅是定量這一個方面哄芜。
忽如一夜春風來貌亭,最近各個科服大廠都在討論轉(zhuǎn)錄組UMI定量的事情。UMI正如火如荼的使用在單細胞轉(zhuǎn)錄組的研究中认臊,同時整合barcode圃庭、UMI信息對單細胞數(shù)據(jù)進行解讀。
早在2012年,關于digital轉(zhuǎn)錄組UMI定量的文章就已發(fā)表剧腻,作者系統(tǒng)的討論了UMI或barcode序列的設計思路拘央、性能驗證等工作∈樵冢總之灰伟,UMI定量更加準確、測序序列可以相互校正從而提高序列準確性儒旬,更重要的是對于低拷貝轉(zhuǎn)錄本的定量也更為準確栏账。
建庫定量示意圖如下:
基于二代測序免不了進行轉(zhuǎn)錄本組裝,組裝過程可能引入組裝錯誤或剪切體的丟失栈源。
而三代測序所測即所得的特點則不存在上述問題的困擾挡爵,與UMI/barcode相結(jié)合不失為一種更高效的思路,以市面上比較流行的PacBio三代測序平臺為例甚垦,在克服轉(zhuǎn)錄本產(chǎn)出低茶鹃、片段選擇等問題后其轉(zhuǎn)錄本準確定量則水到渠成。
(1)孟浩巍《生物信息學100個基礎問題 —— 第38題 當轉(zhuǎn)錄組普遍變化時RNA-Seq怎么進行分析(1)艰亮?》
(2)孟浩巍《生物信息學100個基礎問題 —— 第38題 當轉(zhuǎn)錄組普遍變化時RNA-Seq怎么進行分析(2)闭翩?》
(3)【生信菜鳥團】quantile normalization到底對數(shù)據(jù)做了什么?
(5)生信菜鳥團:StatQuest生物統(tǒng)計學專題 - library normalization進階之edgeR的標準化方法
(6)http://bib.oxfordjournals.org/content/early/2012/09/15/bib.bbs046.full
http://www.ncbi.nlm.nih.gov/pubmed/22872506