WES流程

參考鏈接：

https://mp.weixin.qq.com/s/awdjoXRYobrQAbXmAp3C0g

gatk4的使用：https://zhuanlan.zhihu.com/p/69726572?from_voters_page=true (一定要看！)

https://cloud.tencent.com/developer/article/1445600
http://www.bio-info-trainee.com/3144.html
http://www.reibang.com/p/825e7d618838

gatk官網(wǎng)：需要翻墻

堿基礦工公眾號(hào)文章--WES系列（gatk4相關(guān)）：

https://mp.weixin.qq.com/s/NXV_08mUvsJ5_RZXfwuleQ
https://mp.weixin.qq.com/s/sFnPOfXsJFHI1xwWLzyNLg
https://mp.weixin.qq.com/s/L2w6sH58F44R7n-FZkv7Ig
https://mp.weixin.qq.com/s/Sa019WuSg8fRQgkWAIG4pQ

全外顯子call snp全過程

https://zhuanlan.zhihu.com/p/45094576

0狂塘、gatk工具的使用說明

參考鏈接：https://www.plob.org/article/7005.html

gatk工作原理

（2）GATK是一個(gè)應(yīng)用于前沿科學(xué)研究的軟件薯定，不斷在更新和修正毡泻，因此刘离，在使用GATK進(jìn)行變異檢測(cè)時(shí)晕窑，最好是下載最新的版本抑片，目前的版本是2.8.1（2014-02-25）。下載網(wǎng)站：http://www.broadinstitute.org/gatk/download杨赤。

（3）在GATK使用過程中（見下面圖）敞斋，有些步驟需要用到已知變異信息截汪，對(duì)于這些已知變異，GATK只提供了人類的已知變異信息植捎，可以在GATK的FTP站點(diǎn)下載（GATK resource bundle）衙解。如果要研究的不是人類基因組，需要自行構(gòu)建已知變異焰枢，GATK提供了詳細(xì)的構(gòu)建方法丢郊。

（4）GATK在進(jìn)行BQSR和VQSR的過程中會(huì)使用到R軟件繪制一些圖，因此医咨，在運(yùn)行GATK之前最好先檢查一下是否正確安裝了R和所需要的包枫匾，所需要的包大概包括ggplot2、gplots拟淮、bitops干茉、caTools、colorspace很泊、gdata角虫、gsalib、reshape委造、RColorBrewer等戳鹅。如果畫圖時(shí)出現(xiàn)錯(cuò)誤，會(huì)提示需要安裝的包的名稱昏兆。
gatk-原始數(shù)據(jù)處理：比對(duì)+去除PCR重復(fù)+局部重比對(duì)+BQSR:https://www.plob.org/2014/04/14/7009.html
gatk-變異檢測(cè)+VQSR:https://www.plob.org/article/7023.html
gatk-變異的后續(xù)分析：https://www.plob.org/article/7035.html

1枫虏、Read Group的信息解釋

參考鏈接：https://mp.weixin.qq.com/s/awdjoXRYobrQAbXmAp3C0g
在使用bwa進(jìn)行比對(duì)時(shí)，會(huì)有-R參數(shù)用來補(bǔ)充read group信息爬虱，這對(duì)于后續(xù)進(jìn)行call variation時(shí)必要的
read group：在sam中以@RG開頭隶债，它是用來將比對(duì)的read進(jìn)行分組的。不同的組之間測(cè)序過程被認(rèn)為是相互獨(dú)立的跑筝，這個(gè)信息對(duì)于我們后續(xù)對(duì)比對(duì)數(shù)據(jù)進(jìn)行錯(cuò)誤率分析和Mark duplicate時(shí)非常重要死讹。
1）ID，這是Read Group的分組ID曲梗，一般設(shè)置為測(cè)序的lane ID（不同lane之間的測(cè)序過程認(rèn)為是獨(dú)立的）赞警，下機(jī)數(shù)據(jù)中我們都能看到這個(gè)信息的，一般都是包含在fastq的文件名中
2）PL虏两，指的是所用的測(cè)序平臺(tái)愧旦，這個(gè)信息不要隨便寫！特別是當(dāng)我們需要使用GATK進(jìn)行后續(xù)分析的時(shí)候碘举，更是如此忘瓦！這是一個(gè)很多新手都容易忽視的一個(gè)地方，在GATK中，PL只允許被設(shè)置為：ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN這幾個(gè)信息耕皮【惩桑基本上就是目前市場(chǎng)上存在著的測(cè)序平臺(tái)，當(dāng)然凌停，如果實(shí)在不知道粱年，那么必須設(shè)置為UNKNOWN，名字方面不區(qū)分大小寫
3）SM罚拟，樣本ID台诗，同樣非常重要，有時(shí)候我們測(cè)序的數(shù)據(jù)比較多的時(shí)候赐俗，那么可能會(huì)分成多個(gè)不同的lane分布測(cè)出來拉队，這個(gè)時(shí)候SM名字就是可以用于區(qū)分這些樣本；
4）LB阻逮，測(cè)序文庫的名字粱快，這個(gè)重要性稍微低一些，主要也是為了協(xié)助區(qū)分不同的group而存在叔扼。文庫名字一般可以在下機(jī)的fq文件名中找到事哭，如果上面的lane ID足夠用于區(qū)分的話，也可以不用設(shè)置LB瓜富。
除了以上這四個(gè)之外鳍咱，還可以自定義添加其他的信息，不過如無特殊的需要与柑，對(duì)于序列比對(duì)而言谤辜，這4個(gè)就足夠了。這些信息設(shè)置好之后仅胞，在RG字符串中要用制表符（\t）將它們分開
總結(jié)：ID一般用來寫lane ID每辟，如果在測(cè)的時(shí)候一個(gè)樣本一個(gè)lane，那也可以是sample id干旧。PL必須是它指定的那幾個(gè)。SM是樣本的ID妹蔽，如果是一個(gè)樣本一個(gè)lane的話椎眯，ID=SM，如果是一個(gè)樣本多個(gè)lane的話（測(cè)序很深時(shí)）胳岂，ID是laneID,SM是樣本id编整，要做區(qū)分。LB乳丰，可以隨便設(shè)置掌测。
例子：

$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam

2、為什么比對(duì)完之后要排序（sort）产园？

FASTQ文件里面這些被測(cè)序下來的read是隨機(jī)分布于基因組上面的汞斧，第一步的比對(duì)是按照FASTQ文件的順序把read逐一定位到參考基因組上之后夜郁，隨即就輸出了，它不會(huì)也不可能在這一步里面能夠自動(dòng)識(shí)別比對(duì)位置的先后位置重排比對(duì)結(jié)果粘勒。因此竞端，比對(duì)后得到的結(jié)果文件中，每一條記錄之間位置的先后順序是亂的庙睡，我們后續(xù)去重復(fù)等步驟都需要在比對(duì)記錄按照順序從小到大排序下來才能進(jìn)行事富，所以這才是需要進(jìn)行排序的原因
[注意] 排序后如果發(fā)現(xiàn)新的BAM文件比原來的BAM文件稍微小一些，不用覺得驚訝乘陪，這是壓縮算法導(dǎo)致的結(jié)果统台，文件內(nèi)容是沒有損失的。

3啡邑、去除重復(fù)序列的原因

https://mp.weixin.qq.com/s/awdjoXRYobrQAbXmAp3C0g

什么是重復(fù)序列

首先什么是重復(fù)序列饺谬，重復(fù)序列是在進(jìn)行PCR擴(kuò)增時(shí)，由同一個(gè)DNA分子產(chǎn)生了很多的相同的拷貝谣拣。重復(fù)序列的存在會(huì)導(dǎo)致對(duì)于變異的判斷產(chǎn)生錯(cuò)誤募寨，主要有以下幾點(diǎn)：
1）DNA在打斷的時(shí)候會(huì)發(fā)生一些變異，而PCR會(huì)擴(kuò)大這個(gè)信號(hào)森缠，導(dǎo)致假陽性的出現(xiàn)拔鹰。
2）PCR過程會(huì)引入新的變異，這些變異越早發(fā)生贵涵，那其在后續(xù)的擴(kuò)增中錯(cuò)誤的拷貝會(huì)越多列肢，導(dǎo)致假陽性
3）PCR本身存在序列偏好性，如果存在真實(shí)的變異后宾茂，PCR產(chǎn)生了偏好性瓷马，如對(duì)reference序列擴(kuò)增偏向強(qiáng)烈，那變異的堿基信息會(huì)減少跨晴，導(dǎo)致假陰性欧聘，反之，導(dǎo)致假陽性端盆。
4）目前使用的主流工具怀骤，GATK、Samtools焕妙、Platpus等這種利用貝葉斯原理的變異檢測(cè)算法都是認(rèn)為所用的序列數(shù)據(jù)都不是重復(fù)序列（即將它們和其他序列一視同仁地進(jìn)行變異的判斷蒋伦，所以帶來誤導(dǎo)），因此必須要進(jìn)行標(biāo)記（去除）或者使用PCR-Free的測(cè)序方案
其次是如何識(shí)別或去除重復(fù)序列焚鹊，既然PCR擴(kuò)增是把同一段DNA序列復(fù)制出很多份痕届，那么這些序列在經(jīng)過比對(duì)之后它們一定會(huì)定位到基因組上相同的位置，比對(duì)的信息看起來也將是一樣的！于是研叫，我們就可以根據(jù)這個(gè)特點(diǎn)找到這些重復(fù)序列了锤窑！事實(shí)上，現(xiàn)有的工具包括Samtools和Picard中去除重復(fù)序列的算法也的確是這么做的蓝撇。不同的地方在于果复，samtools的rmdup是直接將這些重復(fù)序列從比對(duì)BAM文件中刪除掉，而Picard的MarkDuplicates默認(rèn)情況則只是在BAM的FLAG信息中標(biāo)記出來渤昌，而不是刪除虽抄，因此這些重復(fù)序列依然會(huì)被留在文件中，只是我們可以在變異檢測(cè)的時(shí)候識(shí)別到它們独柑，并進(jìn)行忽略迈窟。
3.1、為什么有的數(shù)據(jù)處理中需要merge bam忌栅，有的不需要车酣？
在進(jìn)行了PCR去重之后，隨后要進(jìn)行局部區(qū)域重比對(duì)索绪，但是對(duì)于測(cè)序量較大的樣本湖员，會(huì)產(chǎn)生多個(gè)bam文件，因此在這種情況下需要使用samtools將多個(gè)bam merge起來構(gòu)成該樣本的完整的bam

samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

4瑞驱、局部區(qū)域重比對(duì)
局部區(qū)域重比對(duì)也叫Indel 局部區(qū)域重比對(duì)娘摔，其目的是將BWA比對(duì)過程中所發(fā)現(xiàn)有潛在序列插入或者序列刪除（insertion和deletion，簡稱Indel）的區(qū)域進(jìn)行重新校正唤反。這個(gè)過程往往還會(huì)把一些已知的Indel區(qū)域一并作為重比對(duì)的區(qū)域凳寺，那為什么要做這種矯正呢？
其根本原因來自于參考基因組的序列特點(diǎn)和BWA這類比對(duì)算法本身彤侍，注意這里不是針對(duì)BWA肠缨，而是針對(duì)所有的這類比對(duì)算法，包括bowtie等盏阶。這類在全局搜索最優(yōu)匹配的算法在存在Indel的區(qū)域及其附近的比對(duì)情況往往不是很準(zhǔn)確晒奕，特別是當(dāng)一些存在長Indel、重復(fù)性序列的區(qū)域或者存在長串單一堿基（比如般哼，一長串的TTTT或者AAAAA等）的區(qū)域中更是如此吴汪。
另一個(gè)重要的原因是在這些比對(duì)算法中，對(duì)堿基錯(cuò)配和開gap的容忍度是不同的蒸眠。具體體現(xiàn)在罰分矩陣的偏向上，例如杆融，在read比對(duì)時(shí)楞卡，如果發(fā)現(xiàn)堿基錯(cuò)配和開gap都可以的話，它們會(huì)更偏向于錯(cuò)配。但是這種偏向錯(cuò)配的方式蒋腮，有時(shí)候卻還會(huì)反過來引起錯(cuò)誤的開gap淘捡！這就會(huì)導(dǎo)致基因組上原本應(yīng)該是一個(gè)長度比較大的Indel的地方，被錯(cuò)誤地切割成多個(gè)錯(cuò)配和短indel的混合集池摧，這必然會(huì)讓我們檢測(cè)到很多錯(cuò)誤的變異焦除。而且，這種情況還會(huì)隨著所比對(duì)的read長度的增長（比如三代測(cè)序的Read作彤，通常都有幾十kbp）而變得越加嚴(yán)重膘魄。
因此，我們需要有一種算法來對(duì)這些區(qū)域進(jìn)行局部的序列重比對(duì)竭讳。這個(gè)算法通常就是大名鼎鼎的Smith-Waterman算法创葡，它非常適合于這類場(chǎng)景，可以極其有效地實(shí)現(xiàn)對(duì)全局比對(duì)結(jié)果的校正和調(diào)整绢慢，最大程度低地降低由全局比對(duì)算法的不足而帶來的錯(cuò)誤灿渴。而且GATK的局部重比對(duì)模塊，除了應(yīng)用這個(gè)算法之外胰舆，還會(huì)對(duì)這個(gè)區(qū)域中的read進(jìn)行一次局部組裝骚露，把它們連接成為長度更大的序列，這樣能夠更進(jìn)一步提高局部重比對(duì)的準(zhǔn)確性缚窿。
下圖給大家展示一個(gè)序列重比對(duì)之前和之后的結(jié)果棘幸，其中灰色的橫條指的是read，空白黑線指的是deletion滨攻，有顏色的堿基指的是錯(cuò)配堿基够话。

image.png

相信大家都可以明顯地看到在序列重比對(duì)之前，在這個(gè)區(qū)域的比對(duì)數(shù)據(jù)是多么的糟糕光绕，如果就這樣進(jìn)行變異檢測(cè)女嘲，那么一定會(huì)得到很多假的結(jié)果。而在經(jīng)過局部重比對(duì)之后诞帐，這個(gè)區(qū)域就變得非常清晰而分明欣尼，它原本發(fā)生的就只是一個(gè)比較長的序列刪除（deletion）事件，但在原始的比對(duì)結(jié)果中卻被錯(cuò)誤地用堿基錯(cuò)配和短的Indel所代替停蕉。
在gatk中愕鼓，該過程包含了兩個(gè)步驟
第一步，RealignerTargetCreator 慧起，目的是定位出所有需要進(jìn)行序列重比對(duì)的目標(biāo)區(qū)域（如下圖）菇晃；
第二步，IndelRealigner蚓挤，對(duì)所有在第一步中找到的目標(biāo)區(qū)域運(yùn)用算法進(jìn)行序列重比對(duì)磺送，最后得到捋順了的新結(jié)果驻子。

以上這兩個(gè)步驟是缺一不可的，順序也是固定的估灿。而且崇呵，需要指出的是，這里的-R參數(shù)輸入的human.fasta不是BWA比對(duì)中的索引文件前綴馅袁，而是參考基因組序列（FASTA格式）文件域慷，下同。

另外汗销，在重比對(duì)步驟中犹褒，我們還看到了兩個(gè)陌生的VCF文件，分別是：1000G_phase1.indels.b37.vcf和Mills_and_1000G_gold_standard.indels.b37.vcf大溜。這兩個(gè)文件來自于千人基因組和Mills項(xiàng)目化漆，里面記錄了那些項(xiàng)目中檢測(cè)到的人群Indel區(qū)域。我上面其實(shí)也提到了钦奋，候選的重比對(duì)區(qū)除了要在樣本自身的比對(duì)結(jié)果中尋找之外座云，還應(yīng)該把人群中已知的Indel區(qū)域也包含進(jìn)來，而這兩個(gè)是我們?cè)谥乇葘?duì)過程中最常用到的付材。這些文件你可以很方便地在GATK bundle ftp中下載（ftp://ftp.broadinstitute.org/bundle/）朦拖，注意一定要選擇和你的參考基因組對(duì)應(yīng)的版本，我們這里用的是b37版本厌衔。
不過在使用GATK的HaplotypeCaller模塊的時(shí)候璧帝，不需要進(jìn)行重比對(duì)，因?yàn)樵撍惴ㄔ谠O(shè)計(jì)的時(shí)候就已經(jīng)考慮了重比對(duì)富寿，而使用其它工具的時(shí)候睬隶，重比對(duì)是必須的。