今天，繼續(xù)把發(fā)在“解螺旋技術(shù)交流圈”的部分主題整理出來(lái)，分享給你。
1. 請(qǐng)問(wèn)對(duì)于同一份BAM文件使用samtools depth和用samtools mpileup跑出來(lái)的位點(diǎn)的depth有何差異？
你會(huì)注意到這個(gè)差異，應(yīng)該是由于你所用的是Pair-End（PE）測(cè)序的數(shù)據(jù)吧，如果是SE數(shù)據(jù)，差異其實(shí)很小。對(duì)于PE測(cè)序數(shù)據(jù)主要有兩個(gè)地方的差異：
（1）第一個(gè)差異，對(duì)于PE數(shù)據(jù)，mpileup默認(rèn)會(huì)把不正常比對(duì)的PE Read（比如read1和read2的比對(duì)位置彼此間的距離超過(guò)插入片段長(zhǎng)度的波動(dòng)范圍或者read1與read2有一條沒(méi)有比對(duì)上）先排除掉再做計(jì)算，但samtools depth則不會(huì)，depth默認(rèn)不做任何過(guò)濾，只要比上就算。這也是我們會(huì)看到samtools depth計(jì)算的覆蓋深度往往都高于mpileup的最主要原因。如果要讓兩者一致，可以在mpileup中加上 -A 參數(shù)，強(qiáng)制留下不正常的PE比對(duì)結(jié)果即可；
（2）它們之間的第二個(gè)差異是，在默認(rèn)情況下，mpileup還會(huì)過(guò)濾掉測(cè)序質(zhì)量值低于13的堿基，depth默認(rèn)不過(guò)濾。
雖然調(diào)整一下參數(shù)就可以保證兩者一樣。但我并不建議這么做，雖說(shuō)mpileup這里得到的是高質(zhì)量的覆蓋深度，但是說(shuō)到底它和samtools depth的目的還是不同的。
此外，如果要更好地計(jì)算比對(duì)數(shù)據(jù)的覆蓋深度和覆蓋度的話，samtools depth雖然能夠勝任，但是功能還是比較單一，而且由于每個(gè)位點(diǎn)都會(huì)輸出，導(dǎo)致結(jié)果文件總是很巨大，我還是比較推薦使用bedtools2來(lái)完成，如下圖，它的功能和輸出形式要更加豐富。
2. 為什么WES的數(shù)據(jù)無(wú)法使用VQSR進(jìn)行變異質(zhì)控？
其實(shí)不只是WES，還包括很多小panel的數(shù)據(jù)，如果樣本量比較少的話基本都無(wú)法使用VQSR進(jìn)行變異的質(zhì)控。其原因就在VQSR的原理上。
VQSR的核心原理是利用機(jī)器學(xué)習(xí)算法構(gòu)造一個(gè)區(qū)分“好”變異和“壞”變異的分類器。這個(gè)分類器在GATK中是通過(guò)GMM模型來(lái)構(gòu)造的，它在構(gòu)造的時(shí)候并不是盲目地使用所有數(shù)據(jù)來(lái)進(jìn)行構(gòu)造，而是挑出和已知的變異集合Overlap的位點(diǎn)（通常是HapMap數(shù)據(jù)集）——并分配相應(yīng)的可信度權(quán)重來(lái)進(jìn)行訓(xùn)練。
基于群體遺傳的原理，這些已知且被嚴(yán)格驗(yàn)證的變異（如HapMap數(shù)據(jù)）會(huì)被認(rèn)為是更加靠譜的變異，因此在初始化的時(shí)候先把它們當(dāng)作是“好”的——也就是正確的變異。這個(gè)初始變異集很重要，然后利用這些好變異訓(xùn)練一個(gè)區(qū)分好變異的GMM，接著對(duì)全部數(shù)據(jù)進(jìn)行打分，再把評(píng)分最低的那些拿出來(lái)，構(gòu)成一個(gè)最不像正確變異的集合，用來(lái)構(gòu)造一個(gè)區(qū)分壞變異的GMM，用來(lái)專門識(shí)別壞變異。最后同時(shí)用好和壞的GMM再一次同時(shí)對(duì)變異進(jìn)行打分，看每個(gè)變異更像誰(shuí)，就能夠評(píng)判出這個(gè)變異可信的質(zhì)量值了。越靠近好的GMM，質(zhì)量就越高，這就是VQSR過(guò)濾的大致原理（如下圖）。
為了得到理想好的結(jié)果，VQSR在進(jìn)行模型訓(xùn)練的時(shí)候就有一個(gè)最低可用位點(diǎn)數(shù)目的要求——通常是好和壞變異可供訓(xùn)練的數(shù)目必須超過(guò)5000個(gè)，如果Overlap位點(diǎn)太少，是無(wú)法用于訓(xùn)練一個(gè)合適的模型的，這對(duì)于全基因組來(lái)說(shuō)是沒(méi)任何問(wèn)題的，但外顯子區(qū)域加起來(lái)也就差不多50Mb左右，長(zhǎng)度不大，單個(gè)樣本里面包含的變異數(shù)目大約30K-40K。這些位點(diǎn)本來(lái)就不多，它們和已知高質(zhì)量變異集Overlap的就更少了,最終就導(dǎo)致達(dá)不到模型訓(xùn)練的最低要求。所以單個(gè)樣本的WES（或者樣本數(shù)量較少的WES）都無(wú)法使用VQSR進(jìn)行質(zhì)控，小Panel的測(cè)序數(shù)據(jù)也是同理。
但隨著樣本數(shù)目的增加，群體中會(huì)有更多的變異也在這些外顯子區(qū)域中被發(fā)現(xiàn)，從而增大了這個(gè)可用的訓(xùn)練集合，直到滿足了最低訓(xùn)練要求，按照經(jīng)驗(yàn)，通常是30個(gè)樣本（隨著捕獲區(qū)域的差別，會(huì)略有差異），這也是為什么對(duì)于WES數(shù)據(jù)而言，GATK會(huì)提到至少需要30個(gè)樣本才能進(jìn)行VQSR的原因。
3. 樣本量重要，還是測(cè)序深度重要？
我認(rèn)為是樣本量遠(yuǎn)比測(cè)序深度重要。只要有足夠多的樣本，我們甚至可以用很低的測(cè)序深度（比如1x）獲得這些樣本中每個(gè)人準(zhǔn)確的genotype和群體的遺傳頻譜。這是為什么？
其中一個(gè)核心原因是人類這個(gè)物種具有單一祖先起源，這也是一個(gè)重要的前提假設(shè)。但同時(shí)我想強(qiáng)調(diào)一點(diǎn)，這里的“單一”并不是特指只有一個(gè)個(gè)體，而是指形成這個(gè)群體（比如說(shuō)現(xiàn)代人，甚至就只是中國(guó)的漢族人）的祖先歸結(jié)起來(lái)只有為數(shù)不多的若干個(gè)部落。在這種情況下，人群多樣性的源頭實(shí)際上就主要來(lái)自這些部落之間的基因交流和融合。
另一個(gè)核心原因是時(shí)間不夠。人類其實(shí)是一個(gè)很年輕的群體，特別是現(xiàn)代智人（我們這一波），遺傳的分化歷史很短，按照目前估算大約是10萬(wàn)年前才開(kāi)始。而群體出現(xiàn)遺傳差異的動(dòng)力主要有兩個(gè)：（1）基因組自身的突變和重組；（2）生殖細(xì)胞在形成配子過(guò)程中發(fā)生的重組。但基因組突變和重組的速率都是很低的，大概只有10^-8次方左右。也就是說(shuō)一個(gè)人因?yàn)橥蛔兯鶐?lái)的遺傳差異，積累起來(lái)大約是30-100個(gè)。這個(gè)只是序列上的突變（主要是點(diǎn)突變），重組雖然有所不同——它是大范圍序列的交換，影響的范圍很大，但是一般不認(rèn)為它直接帶來(lái)序列突變。我們可以理解為它帶來(lái)的是突變?cè)谡麄€(gè)群體中的擴(kuò)散和分配。
然而，10萬(wàn)年的時(shí)間，差不多只有5000代人，這個(gè)數(shù)字放在物種遺傳的歷史上是很短暫的一瞬，這個(gè)時(shí)間跨度不足以引起整個(gè)群體的多樣性爆發(fā)。對(duì)于東亞人來(lái)說(shuō)則更少，目前發(fā)表過(guò)的研究表明，東亞人的歷史更短，大概起源于6萬(wàn)年前，所以你會(huì)在千人基因組項(xiàng)目中看到東亞人（特別是漢族人）內(nèi)部的分化差異極小。最終歸結(jié)起來(lái)，人類這個(gè)群體中單倍體的組合數(shù)目是非常有限的。
所以如果要揭示一個(gè)特定群體的遺傳圖譜，我們大可不必對(duì)全體樣本都進(jìn)行高深度測(cè)序，只需要把其中一部分人進(jìn)行深測(cè)獲得較高質(zhì)量的變異集合，然后其他樣本則直接使用低深度測(cè)序（甚至是定制的芯片測(cè)序，不過(guò)我更偏向于選擇低深度全基因組測(cè)序），再結(jié)合連鎖不平衡遺傳定律，我們就完全有能力推斷那些沒(méi)被充分覆蓋的區(qū)域中的具體基因型，千人基因組和冰島人就是這樣的一個(gè)例子。
GATK的HaplotypeCaller算法實(shí)際上也是利用這樣的原理實(shí)現(xiàn)了更加準(zhǔn)確的變異檢測(cè)的。在變異檢測(cè)時(shí)，GATK會(huì)利用所有樣本的數(shù)據(jù)，預(yù)先構(gòu)造出這個(gè)群體的Haplotype組合（這應(yīng)該也是HaplotypeCaller這個(gè)名字的由來(lái)），以及這個(gè)組合中各個(gè)單體型在群體中的后驗(yàn)概率，然后再依據(jù)每個(gè)樣本自己的比對(duì)數(shù)據(jù)，通過(guò)貝葉斯原理計(jì)算出各個(gè)樣本在每個(gè)位點(diǎn)上的基因型和各自基因型的后驗(yàn)概率。如果參與分析的樣本足夠多，那么理論上它就能夠構(gòu)建出更加準(zhǔn)確的Haplotype組合，然后反過(guò)來(lái)就會(huì)提升各個(gè)樣本的變異檢測(cè)結(jié)果。
4. 怎么通過(guò)LD衰減距離去看群體的一個(gè)遺傳多樣性呢？
LD本身反應(yīng)的是一個(gè)物種基因組上發(fā)生過(guò)的重組情況。基因組的重組在每一代都會(huì)發(fā)生，如果一個(gè)群體越古老，那么可以預(yù)期它基因組中發(fā)生過(guò)重組的次數(shù)就越多，那么相應(yīng)的它的LD長(zhǎng)度就會(huì)越短，從而這個(gè)族群的遺傳多樣性就越高。比如在現(xiàn)代人類中，遺傳多樣性最高的是非洲人，他們歷史最久遠(yuǎn)，而我們東亞黃種人，多樣性則是最低的。如果我們要通過(guò)基因芯片對(duì)非洲人的某些特征進(jìn)行全基因組關(guān)聯(lián)分析，那么理論上適合這個(gè)群體的芯片密度要比我們黃種人的高。
5. 生物信息工程師可以分為多少種類型？
總的來(lái)說(shuō)包含三個(gè)大的分類導(dǎo)向：
第一類，技術(shù)導(dǎo)向，目標(biāo)是開(kāi)發(fā)更好的算法，思考如何利用數(shù)理和計(jì)算機(jī)等方面的知識(shí)提供更好的工具和平臺(tái)。幫助解決組學(xué)問(wèn)題，比如編寫(xiě)比對(duì)算法、組裝算法、變異檢測(cè)算法、質(zhì)控程序等，當(dāng)然也包括編寫(xiě)生產(chǎn)級(jí)別的數(shù)據(jù)分析流程（如標(biāo)準(zhǔn)化WGS流程），這一類型的生信工程師解決的是生產(chǎn)工具的問(wèn)題。
第二類，數(shù)據(jù)導(dǎo)向/問(wèn)題導(dǎo)向，或者叫“業(yè)務(wù)”導(dǎo)向——這里的業(yè)務(wù)包括科學(xué)研究和商業(yè)應(yīng)用。主要是解決生物和組學(xué)問(wèn)題、遺傳咨詢等，如癌癥研究、群體遺傳學(xué)等。這類人更多的是工具的使用者，他們會(huì)根據(jù)具體的“業(yè)務(wù)”需要組合最合適的算法和工具來(lái)解決問(wèn)題，這一類人需要較深的生物和基因遺傳學(xué)知識(shí)背景。同時(shí)，必須對(duì)自己所在的領(lǐng)域有一個(gè)完整的認(rèn)識(shí)，知道在什么場(chǎng)景下需要什么數(shù)據(jù)，應(yīng)用什么算法，使用什么數(shù)理知識(shí)和什么工具，才能更好地解決問(wèn)題——其實(shí)這一類人也是真正知道該做什么分析流程的人。
關(guān)于這一類生信工程師，或者應(yīng)該稱為“基因組學(xué)專家”更加合適，他們包含很多方面，比如群體遺傳學(xué)、動(dòng)植物基因組學(xué)、進(jìn)化、腫瘤研究、醫(yī)學(xué)基因檢測(cè)、消費(fèi)級(jí)基因檢測(cè)、遺傳咨詢等。他/她們通常是依據(jù)“業(yè)務(wù)”目標(biāo)，運(yùn)用相應(yīng)的技術(shù)手段和工具（包括WGS、WES、RNAseq、甲基化測(cè)序、相關(guān)組學(xué)分析方法等）解決達(dá)成目標(biāo)道路上的問(wèn)題。這里每一個(gè)都可以再進(jìn)一步展開(kāi)，總的來(lái)說(shuō)，這個(gè)類型是工具的使用方，具體組學(xué)問(wèn)題的解決者。
上面這兩類看起來(lái)各有特點(diǎn)，掌握的知識(shí)點(diǎn)各有側(cè)重，但其實(shí)并不能割裂，真正做得好的人，都是兩類通吃的（可能只是兩強(qiáng)相較，某一類更突出）。因?yàn)槟苌羁汤斫馍飭?wèn)題和組學(xué)問(wèn)題的人，才能創(chuàng)造出真正合適的工具和流程。
第三類，資源和人導(dǎo)向，或者叫“Boss”/PI導(dǎo)向。這些人由于各自成長(zhǎng)經(jīng)歷的不同，可能已經(jīng)和上面的情況有所出入了（很難說(shuō)會(huì)全都懂），他們中有些可能更擅長(zhǎng)于去找資源，搭橋，做連接。他們更多的不是解決具體問(wèn)題，而是盡可能地提出好問(wèn)題，發(fā)現(xiàn)好方向，并為提供解決這些問(wèn)題創(chuàng)造環(huán)境和條件。這一類人其實(shí)往往也是第一類和第二類人發(fā)展在后面的一個(gè)方向。
技術(shù)交流圈往期精華
我的微信公眾號(hào)：解螺旋的礦工 歡迎關(guān)注更及時(shí)了解更多信息。
這是知識(shí)星球：『解螺旋技術(shù)交流圈』，是一個(gè)我與讀者朋友們的私人朋友圈。我有9年前沿而完整的生物信息學(xué)、NGS領(lǐng)域的工作經(jīng)歷，在該領(lǐng)域發(fā)有多篇Nature級(jí)別的科學(xué)文章，我也希望借助這個(gè)知識(shí)星球把自己的一些微薄經(jīng)驗(yàn)分享給更多對(duì)組學(xué)感興趣的伙伴們。
自從星球正式運(yùn)行以來(lái)，已經(jīng)過(guò)去了6個(gè)月，星球的成員也已經(jīng)超過(guò)220人了。所分享的主題超過(guò)了500個(gè)，回答的問(wèn)題超過(guò)了140個(gè)，精華70個(gè)。我在知識(shí)星球上留下的文字估計(jì)也已經(jīng)超過(guò)10萬(wàn)字，加上大家的就更多了，相信接下來(lái)星球的內(nèi)容一定還會(huì)不斷豐富。另外，上周獲得了知識(shí)星球官方評(píng)選的“最優(yōu)質(zhì)星球”優(yōu)秀獎(jiǎng)。
這是知識(shí)星球上 第一個(gè)真正與基因組學(xué)和生物信息學(xué)強(qiáng)相關(guān)的圈子。我希望能夠借此營(yíng)造一個(gè)高質(zhì)量的組學(xué)知識(shí)圈和人脈圈，通過(guò)提問(wèn)、彼此分享、交流經(jīng)驗(yàn)、心得等，彼此更好地學(xué)習(xí)生信知識(shí)，提升基因組數(shù)據(jù)分析和解讀的能力。
在這里你可以結(jié)識(shí)到全國(guó)優(yōu)秀的基因組學(xué)和生物信息學(xué)專家，同時(shí)可以分享你的經(jīng)驗(yàn)、見(jiàn)解和思考，有問(wèn)題也可以向我提問(wèn)和圈里的星友們提問(wèn)。
知識(shí)星球邀請(qǐng)鏈接：「解螺旋技術(shù)交流圈」
