[原]測序量估計

 

考慮到cnblog不適合基因組領域這種類型的文章，進過多番折騰，終於用jekyll+github搭了個獨立博http://www.stbioinf.com/，如今博客已經搬遷至這裏，很是歡迎你們訂閱！　　

 

考慮這樣一個問題，「若是要保證基因組上95%的區域其覆蓋深度在30x以上的話，那麼最低的平均測序深度應該是多少？」。

　　關於測序量的估計，對於作生物信息的人來說應算是屢見不鮮了，多數時候咱們都能直接根據以往項目的經驗來得到，或是說的更具體些，在變異檢測中通常要有25x以上的覆蓋度才能獲得一個比較靠譜的結果，因而以此爲目的給出測序量的估計值；固然少數狀況下也會有直接拍腦殼拍出一個值來的瘋狂行爲，不過嘛，雖然說是拍腦殼，但也不是隨便拍的，拍腦殼的背景靠的是身後豐富的經驗。相對更好一些的估計方式就是直接模擬數據，不過老是用模擬數據仍是讓人以爲麻煩，最好是能不用花多少時間，也不用作不少的計算就能脫口給出。我想在這裏說一下這種狀況下個人解法。固然了並不必定徹底準確，僅做交流，歡迎各位看官拍磚。

　　閒話說完，回到上面的問題，在不經過數據模擬也不「拍腦殼」的狀況下，要如何才能估算出一個合理的值呢？其實在做出任何推斷以前咱們都應當要先有一個合理的前提假設，或者說是理論依據來做爲後續分析的基礎。咱們都知道短序列測序的一個特色是，在理論狀況下位點被覆蓋到的深度符合泊松分佈（測序沒什麼問題的話，實際的情形也相差很少），但實際上在這種狀況下用正態分佈來考慮也是合理的，做爲一個估計值，偏差也是可以接受的，這是咱們的基礎。之因此想用正態分佈來考慮，是由於正態分佈有許多方便於計算的性質。其中一個頗有用的法則，就是「68-95-99」法則，意思就是距離均值一個標準差的區域圍起來的面積大約是整體的68%，2個標準差的區域範圍的面積是整體的95%，3個標準差區域範圍佔到了整體的99%，若是你本身想要驗證這一法則也並不困難，只需作些積分就能算出來，但這裏就不作計算了。以下圖，均值用μ表示，標準差用σ表示。

　　

　　如今事情就很簡單了，從圖中咱們能夠看出，只要30x深度的位置在-2σ如下，那麼就能達到理論的要求。要獲得這一結果，問題就只剩下一個了，此時咱們只須要知道測序深度分佈的標準差就能粗略估計出此時咱們所須要的最低平均測序深度。雖然這個標準差跟許多因素有關，這裏以illumina公司的Hiseq系列測序儀爲例子，依照以往基因組重測序的經驗，σ約等於10x。那麼，簡單算一下，此刻，理論上咱們只須要測50x就可使得基因組上有97.7%的區域其覆蓋深度在30x以上了，注意這裏不是95%了，由於咱們的區域其實是[-2σ, +∞)，而不是[-2σ,+2σ]! 再除掉一些邊邊角角的偏差，50x這個值在這裏應當是合理的了。

　　以上計算都是以正態分佈爲基礎而作出的估計。固然了，若是必定要用泊松分佈去推算也能夠，只是運算起來會麻煩不少。此外，若是是不一樣系列或是不一樣公司的測序儀，σ就不必定是10了。