考慮這樣一個問題，如果要保證基因組上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了。
閱讀更多 解螺旋的礦工 的文章
