一.RNA-seq數(shù)據(jù)分析的基礎(chǔ)知識
我們平常拿到的差異分析數(shù)據(jù)都長這樣: 橫軸是基因,縱軸是基因在每個樣本中測到的read數(shù),也就是count數(shù)。 想知道RNA-seq的每個基因count值是如何分布的?我們來畫個柱狀圖: ![]() 縱軸是基因,橫軸是count數(shù),這個圖揭示了某個樣本中所有基因的count數(shù)分布。 上圖體現(xiàn)了RNA-seq 基因count數(shù)分布的3個基本特點
對于這種類型的數(shù)據(jù),很明顯正態(tài)分布模型就不適用了(注意,這個與芯片數(shù)據(jù)不同哦),那我們應(yīng)該用什么模型去擬合RNA-seq gene的count數(shù)的分布呢? 泊松分布,二項分布都在一定程度上能擬合這種偏態(tài)分布,但是目前的觀點認為,負二項分布在擬合RNA-seq的count數(shù)分布上表現(xiàn)最佳,因為RNA-seq counts數(shù)還具有這樣的分布規(guī)律: ![]() 隨著基因表達量的增加,基因count數(shù)分布的方差快速增加。低表達基因中count數(shù)分布的方差大小不一(學(xué)名:具有方差異質(zhì)性)。泊松分布和二項分布表示,這么奇怪的分布,我干不來,找別人吧!哼!但是負二項分布對于這種分布擬合度較高,它的公式長這樣:(給有需要的人吧) ![]()
別急,哈佛大佬速速來解答疑問: 首先,為啥要那么多重復(fù)呢?因為重復(fù)會讓差異分析更加精確,可信!基因表達差異也有能因為一些無關(guān)的原因,實驗污染,不明確的技術(shù)偏差等等。(如上圖),我們基因差異分析的目的找到的差異是我們實驗組和對照組之間真正有意義的生物學(xué)差異。因此,足夠的生物學(xué)重復(fù)(3組以上)是必要滴 ![]() 舉個栗子:geneA平均表達在處理后是處理前的兩倍,但是這么大的組內(nèi)離散真的會有顯著的組間差異嗎?當(dāng)然不行。而且這種離散度的趨勢只有3個以上的重復(fù)才能充分體現(xiàn),所以不能偷懶哦。 還不信?上圖! ![]() https:///10.1093/bioinformatics/btt688 這篇2014發(fā)表在bioinformatics的文章表明,生物學(xué)重復(fù)比增加測序深度更能提高差異分析的精確度和可信度。 二. DEseq2 測序流程詳解 講了半天鋪墊,想必大家已經(jīng)迫不及待了,來,上正菜!: DEseq2測序的流程 ![]()
(1)DEseq2的基因標(biāo)準(zhǔn)化原理:DEseq有自己的一套count標(biāo)準(zhǔn)化程序:其實TPM之類的標(biāo)準(zhǔn)化方法雖然解決了基因長度和測序深度的影響的問題,但還是不能解決一個問題:那就是測序文庫組成不同造成的差異 這種差異的來源是一個基因被敲減了,完全沒表達了,因而影響到了其他基因。 ![]() DEseq2使用提高中位表達基因的辦法解決了這個問題。想知道具體如何解決的?生信菜鳥團的這篇文章可以完美滿足你的需求(向大佬膜拜):https://mp.weixin.qq.com/s/fw0muJwF-cz2ki2LkflpVw (2) 初始步驟:評估樣本之間整體的相似性 請回答下列問題:哪些樣本與別的樣本相似,哪些不同?樣本間的區(qū)別來自哪里?這種差異與我們需要找的是否一致? ![]() 使用聚類和PCA來得到答案 至于PCA和聚類的原理?你同樣可以在生信菜鳥團的專欄里得到答案:https://mp.weixin.qq.com/s/h0pZ0_ZK9BK-K1XI6ZKDzg --PCA的https://mp.weixin.qq.com/s/1rOVTuAxwjgTDThTThtViw --聚類的
![]() 差異分析分為多個部分: 1.計算離散度 2.擬合并壓縮基因的分布使之更適合建模 3.建模并進行統(tǒng)計學(xué)檢驗 1計算size factors 使用size factors對reads進行標(biāo)準(zhǔn)化(就是我上面說的那個原理,剛剛只是說了原理,這個是在軟件中的運行步驟)
![]() DEseq2需要處理這種表達量和方差的關(guān)系。我們不希望差異基因都是方差大的高表達基因。所以,DEseq2使用dispersion代替variation, 來表達基因的表達和方差,公式為 Dispersion is calculated by Var = μ + α*μ^2, where α = dispersion, Var = variance, and μ = mean, giving the relationship 因此dispersion與variation還有mean expression關(guān)系為: ![]() 離散度(dispersion)與表達量成負相關(guān),與方差成正相關(guān) Dispersion可以反應(yīng)同一平均表達量下基因的離散程度。 ![]() 如何計算呢?只有3-6個重復(fù)難以計算出可信的基因dispersion,所以DEseq2假設(shè)有著相同表達的基因有大致相同的統(tǒng)計分布。所以share the information across gene ,用多個大致相同表達的基因的方差,一起算出幾個基因的dispersion
![]() 紅色就是擬合曲線,擬合基因表達與dispersion的關(guān)系 黑色是基因點
![]() 壓縮的程度由: (1)gene離線的距離 (2)樣本量 決定 Shrink方法對降低差異分析中的假陽性率至關(guān)重要。 因為只有更靠近曲線才能更好的擬合模型。但是過高的gene不會被壓縮,因為他們可能是因為技術(shù)或生物學(xué)原因產(chǎn)生的離群點。(圖上被圈起來的小點)
這是一個表示你的數(shù)據(jù)有沒有很好的擬合模型的方法。如果點在曲線周圍,說明擬合度較高,如果分散,說明擬合度較低。 plotDispEsts() ![]() 這是差的 ![]() 如果像這樣下降又增加。高表達對應(yīng)的應(yīng)該是高dispersion這說明高表達基因?qū)?yīng)的variation降低了,提示可能有離群樣本或者污染。
三.解讀結(jié)果 ![]() 1.p-value Adj-p是經(jīng)過統(tǒng)計學(xué)檢驗后經(jīng)過FDR調(diào)整過的結(jié)果 啥?什么是FDR,那你應(yīng)該看看這篇文章:https://mp.weixin.qq.com/s/dDi7I8LcWSl40JmbkSqCfQ 為啥P-value會出現(xiàn)NA呢? 可能因為在差異分析之前被篩掉了,這樣搞更能提高差異分析的效能,DEseq2不會物理移走gene,但是會出現(xiàn)NA,可能出現(xiàn)NA的情況有: (1)gene在所有樣本中都是0 (2)gene中有一個樣本出現(xiàn)離群—cook檢驗 (3)gene有著極低的表達—對應(yīng)著high dispersion,只要你提高p值閾值,應(yīng)該不會出現(xiàn)這種情況。也有可能被independent filting給干掉。但只有adj p會受到影響 2.log2Foldchange
超過P值的樣本很多,有沒有更嚴(yán)格的指標(biāo)呢? Log2FC:(需要注意,這里差異的倍數(shù)取Log2了哦) log2 (normalized_counts_group1 / normalized_counts_group2) 但是log2FC也有問題,它只考慮了表達值,沒考慮組內(nèi)差異,這該怎么辦呢? DEseq軟件已經(jīng)幫你想好了: 為了更精確地計算log2FC,DEseq2對 (1)low counts (2)high dispersion 的基因進行了零壓縮,在全基因表達的范圍內(nèi)。 ![]() 例子:綠色和紫色是兩個在C57小鼠中表達的gene,紫色gene方差更大。在壓縮之前是實線數(shù)據(jù),在壓縮后是虛線數(shù)據(jù)。明顯,經(jīng)過壓縮,方差的數(shù)據(jù)LFC減少。因此,雖然兩個數(shù)據(jù)標(biāo)準(zhǔn)化后的counts相同,但LFC不同。 注意!壓縮不會改變顯著差異的基因數(shù),只會調(diào)整logFC的大小。 |
|