今天，繼續把發在「解螺旋技術交流圈」的部分主題整理出來，分享給你。
1
請問對於同一份BAM文件使用samtools depth和用samtools mpileup跑出來的位點的depth有何差異？
你會注意到這個差異，應該是由於你所用的是Pair-End（PE）測序的數據吧，如果是SE數據，差異其實很小。對於PE測序數據主要有兩個地方的差異：
（1）第一個差異，對於PE數據，mpileup默認會把不正常比對的PE Read（比如read1和read2的比對位置彼此間的距離超過插入片段長度的波動範圍或者read1與read2有一條沒有比對上）先排除掉再做計算，但samtools depth則不會，depth默認不做任何過濾，只要比上就算。這也是我們會看到samtools depth計算的覆蓋深度往往都高於mpileup的最主要原因。如果要讓兩者一致，可以在mpileup中加上 -A 參數，強制留下不正常的PE比對結果即可；
（2）它們之間的第二個差異是，在默認情況下，mpileup還會過濾掉測序質量值低於13的鹼基，depth默認不過濾。
雖然調整一下參數就可以保證兩者一樣。但我並不建議這麼做，雖說mpileup這裡得到的是高質量的覆蓋深度，但是說到底它和samtools depth的目的還是不同的。
此外，如果要更好地計算比對數據的覆蓋深度和覆蓋度的話，samtools depth雖然能夠勝任，但是功能還是比較單一，而且由於每個位點都會輸出，導致結果文件總是很巨大，我還是比較推薦使用bedtools2來完成，如下圖，它的功能和輸出形式要更加豐富。
bedtools2計算基因組覆蓋度的不同模式
2
為什麼WES的數據無法使用VQSR進行變異質控？
其實不只是WES，還包括很多小panel的數據，如果樣本量比較少的話基本都無法使用VQSR進行變異的質控。其原因就在VQSR的原理上。
VQSR的核心原理是利用機器學習算法構造一個區分「好」變異和「壞」變異的分類器。這個分類器在GATK中是通過GMM模型來構造的，它在構造的時候並不是盲目地使用所有數據來進行構造，而是挑出和已知的變異集合Overlap的位點（通常是HapMap數據集）——並分配相應的可信度權重來進行訓練。
基於群體遺傳的原理，這些已知且被嚴格驗證的變異（如HapMap數據）會被認為是更加靠譜的變異，因此在初始化的時候先把它們當作是「好」的——也就是正確的變異。這個初始變異集很重要，然後利用這些好變異訓練一個區分好變異的GMM，接著對全部數據進行打分，再把評分最低的那些拿出來，構成一個最不像正確變異的集合，用來構造一個區分壞變異的GMM，用來專門識別壞變異。最後同時用好和壞的GMM再一次同時對變異進行打分，看每個變異更像誰，就能夠評判出這個變異可信的質量值了。越靠近好的GMM，質量就越高，這就是VQSR過濾的大致原理（如下圖）。
VQSR區分好變異和壞變異的分類器
為了得到理想好的結果，VQSR在進行模型訓練的時候就有一個最低可用位點數目的要求——通常是好和壞變異可供訓練的數目必須超過5000個，如果Overlap位點太少，是無法用於訓練一個合適的模型的，這對於全基因組來說是沒任何問題的，但外顯子區域加起來也就差不多50Mb左右，長度不大，單個樣本裡面包含的變異數目大約30K-40K。這些位點本來就不多，它們和已知高質量變異集Overlap的就更少了,最終就導致達不到模型訓練的最低要求。所以單個樣本的WES（或者樣本數量較少的WES）都無法使用VQSR進行質控，小Panel的測序數據也是同理。
但隨著樣本數目的增加，群體中會有更多的變異也在這些外顯子區域中被發現，從而增大了這個可用的訓練集合，直到滿足了最低訓練要求，按照經驗，通常是30個樣本（隨著捕獲區域的差別，會略有差異），這也是為什麼對於WES數據而言，GATK會提到至少需要30個樣本才能進行VQSR的原因。
3
樣本量重要，還是測序深度重要？
我認為是樣本量遠比測序深度重要。只要有足夠多的樣本，我們甚至可以用很低的測序深度（比如1x）獲得這些樣本中每個人準確的genotype和群體的遺傳頻譜。這是為什麼？
其中一個核心原因是人類這個物種具有單一祖先起源，這也是一個重要的前提假設。但同時我想強調一點，這裡的「單一」並不是特指只有一個個體，而是指形成這個群體（比如說現代人，甚至就只是中國的漢族人）的祖先歸結起來只有為數不多的若干個部落。在這種情況下，人群多樣性的源頭實際上就主要來自這些部落之間的基因交流和融合。至於什麼是基因交流，大家可以自行腦補。
另一個核心原因是時間不夠。人類其實是一個很年輕的群體，特別是現代智人（我們這一波），遺傳的分化歷史很短，按照目前估算大約是10萬年前才開始。而群體出現遺傳差異的動力主要有兩個：（1）基因組自身的突變和重組；（2）生殖細胞在形成配子過程中發生的重組。但基因組突變和重組的速率都是很低的，大概只有10^-8次方左右。也就是說一個人因為突變所帶來的遺傳差異，積累起來大約是30-100個。這個只是序列上的突變（主要是點突變），重組雖然有所不同——它是大範圍序列的交換，影響的範圍很大，但是一般不認為它直接帶來序列突變。我們可以理解為它帶來的是突變在整個群體中的擴散和分配。
然而，10萬年的時間，差不多只有5000代人，這個數字放在物種遺傳的歷史上是很短暫的一瞬，這個時間跨度不足以引起整個群體的多樣性爆發。對於東亞人來說則更少，目前發表過的研究表明，東亞人的歷史更短，大概起源於6萬年前，所以你會在千人基因組項目中看到東亞人（特別是漢族人）內部的分化差異極小。最終歸結起來，人類這個群體中單倍體的組合數目是非常有限的。
所以如果要揭示一個特定群體的遺傳圖譜，我們大可不必對全體樣本都進行高深度測序，只需要把其中一部分人進行深測獲得較高質量的變異集合，然後其他樣本則直接使用低深度測序（甚至是定製的晶片測序，不過我更偏向於選擇低深度全基因組測序），再結合連鎖不平衡遺傳定律，我們就完全有能力推斷那些沒被充分覆蓋的區域中的具體基因型，千人基因組和冰島人就是這樣的一個例子。
GATK的HaplotypeCaller算法實際上也是利用這樣的原理實現了更加準確的變異檢測的。在變異檢測時，GATK會利用所有樣本的數據，預先構造出這個群體的Haplotype組合（這應該也是HaplotypeCaller這個名字的由來），以及這個組合中各個單體型在群體中的後驗概率，然後再依據每個樣本自己的比對數據，通過貝葉斯原理計算出各個樣本在每個位點上的基因型和各自基因型的後驗概率。如果參與分析的樣本足夠多，那麼理論上它就能夠構建出更加準確的Haplotype組合，然後反過來就會提升各個樣本的變異檢測結果。
4
怎麼通過LD衰減距離去看群體的一個遺傳多樣性呢？
LD本身反應的是一個物種基因組上發生過的重組情況。基因組的重組在每一代都會發生，如果一個群體越古老，那麼可以預期它基因組中發生過重組的次數就越多，那麼相應的它的LD長度就會越短，從而這個族群的遺傳多樣性就越高。比如在現代人類中，遺傳多樣性最高的是非洲人，他們歷史最久遠，而我們東亞黃種人，多樣性則是最低的。如果我們要通過基因晶片對非洲人的某些特徵進行全基因組關聯分析，那麼理論上適合這個群體的晶片密度要比我們黃種人的高。
5
生物信息工程師可以分為多少種類型？
總的來說包含三個大的分類導向：
第一類，技術導向，目標是開發更好的算法，思考如何利用數理和計算機等方面的知識提供更好的工具和平臺。幫助解決組學問題，比如編寫比對算法、組裝算法、變異檢測算法、質控程序等，當然也包括編寫生產級別的數據分析流程（如標準化WGS流程），這一類型的生信工程師解決的是生產工具的問題。
第二類，數據導向/問題導向，或者叫「業務」導向——這裡的業務包括科學研究和商業應用。主要是解決生物和組學問題、遺傳諮詢等，如癌症研究、群體遺傳學等。這類人更多的是工具的使用者，他們會根據具體的「業務」需要組合最合適的算法和工具來解決問題，這一類人需要較深的生物和基因遺傳學知識背景。同時，必須對自己所在的領域有一個完整的認識，知道在什麼場景下需要什麼數據，應用什麼算法，使用什麼數理知識和什麼工具，才能更好地解決問題——其實這一類人也是真正知道該做什麼分析流程的人。
關於這一類生信工程師，或者應該稱為「基因組學專家」更加合適，他們包含很多方面，比如群體遺傳學、動植物基因組學、進化、腫瘤研究、醫學基因檢測、消費級基因檢測、遺傳諮詢等。他/她們通常是依據「業務」目標，運用相應的技術手段和工具（包括WGS、WES、RNAseq、甲基化測序、相關組學分析方法等）解決達成目標道路上的問題。這裡每一個都可以再進一步展開，總的來說，這個類型是工具的使用方，具體組學問題的解決者。
上面這兩類看起來各有特點，掌握的知識點各有側重，但其實並不能割裂，真正做得好的人，都是兩類通吃的（可能只是兩強相較，某一類更突出）。因為能深刻理解生物問題和組學問題的人，才能創造出真正合適的工具和流程。
第三類，資源和人導向，或者叫「Boss」/PI導向。這些人由於各自成長經歷的不同，可能已經和上面的情況有所出入了（很難說會全都懂），他們中有些可能更擅長於去找資源，搭橋，做連接。他們更多的不是解決具體問題，而是儘可能地提出好問題，發現好方向，並為提供解決這些問題創造環境和條件。這一類人其實往往也是第一類和第二類人發展在後面的一個方向。
※ ※ ※
在這裡你可以結識到全國優秀的基因組學和生物信息學專家，同時可以分享你的經驗、見解和思考，有問題也可以向我提問和圈裡的星友們提問。