上次我們整理到bwa比對后得到bam文件搂誉，下一步我們要通過GATK流程從bam文件中call variant癣丧。

一子寓、使用GATK前須知事項：

對GATK的測試主要使用的是人類全基因組和外顯子組的測序數(shù)據(jù)劲厌，而且全部是基于illumina數(shù)據(jù)格式路呜，目前還沒有提供其他格式文件（如Ion Torrent）或者實驗設計（RNA-Seq）的分析方法冠句。

GATK是一個應用于前沿科學研究的軟件轻掩，不斷在更新和修正，因此懦底，在使用GATK進行變異檢測時唇牧，最好是下載最新的版本。

在GATK使用過程中聚唐，有些步驟需要用到已知變異信息丐重，對于這些已知變異，GATK只提供了人類的已知變異信息 杆查，可以在GATK的FTP站點下載（GATK resource bundle）扮惦。如果要研究的不是人類基因組，需要自行構建已知變異亲桦，GATK提供了詳細的構建方法崖蜜。

GATK在進行BQSR和VQSR的過程中會使用到R軟件繪制一些圖，因此客峭，在運行GATK之前最好先檢查一下是否正確安裝了R和所需要的包豫领，所需要的包大概包括ggplot2、gplots舔琅、bitops等恐、caTools、colorspace搏明、gdata鼠锈、gsalib、reshape星著、RColorBrewer等购笆。如果畫圖時出現(xiàn)錯誤，會提示需要安裝的包的名稱虚循。

第一步：原始數(shù)據(jù)的處理：對原始下機fastq文件進行過濾和比對（mapping）對于Illumina下機數(shù)據(jù)推薦使用bwa進行mapping同欠。前邊已講样傍。

bwa中 -r str：定義頭文件∑趟欤‘@RG\tID:foo\tSM:bar’衫哥，如果在此步驟不進行頭文件定義，在后續(xù)GATK分析中還是需要重新增加頭文件襟锐。GATK2.0以上版本將不再支持無頭文件的變異檢測撤逢。

ID str：輸入reads集ID號；LB：read集文庫名粮坞；PL：測序平臺（illunima或solid）蚊荣；PU：測序平臺下級單位名稱（run的名稱）；SM：樣本名稱

對于最后得到的sam文件莫杈，將比對上的結(jié)果提取出來（awk即可處理）互例，即可直接用于GATK的分析。

注意：由于GATK在下游的snpcalling時筝闹，是按染色體進行callsnp的媳叨。因此，在準備原始sam文件時关顷，可以先按染色體將文件分開糊秆，這樣會提高運行速度。但是當數(shù)據(jù)量不足時议双，可能會影響后續(xù)的VQSR分析扩然，這是需要注意的。

對sam文件進行進行重新排序（reorder）

由BWA生成的sam文件時按字典式排序法進行的排序（lexicographically）進行排序的（chr10聋伦，chr11…chr19夫偶，chr1，chr20…chr22觉增，chr2兵拢，chr3…chrM，chrX逾礁，chrY）说铃，但是GATK在進行callsnp的時候是按照染色體組型（karyotypic）進行的（chrM，chr1嘹履，chr2…chr22腻扇，chrX，chrY）砾嫉，因此要對原始sam文件進行reorder幼苛。可以使用picard-tools中的ReorderSam完成焕刮。

對bam文件進行sort排序處理 一步是將sam文件中同一染色體對應的條目按照坐標順序從小到大進行排序舶沿∏奖可以使用picard-tools中SortSam完成。

第一步： Duplicates Marking

在制備文庫的過程中括荡，由于PCR擴增過程中會存在一些偏差高镐，也就是說有的序列會被過量擴增。這樣畸冲，在比對的時候嫉髓，這些過量擴增出來的完全相同的序列就會比對到基因組的相同位置。而這些過量擴增的reads并不是基因組自身固有序列邑闲，不能作為變異檢測的證據(jù)岩喷，因此，要盡量去除這些由PCR擴增所形成的duplicates监憎，這一步可以使用picard-tools來完成。

去重復的過程是給這些序列設置一個flag以標志它們婶溯，方便GATK的識別鲸阔。還可以設置 REMOVE_DUPLICATES=true 來丟棄duplicated序列。對于是否選擇標記或者刪除迄委，對結(jié)果應該沒有什么影響褐筛，GATK官方流程里面給出的例子是僅做標記不刪除。這里定義的重復序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置叙身，那么就認為這樣的reads是由PCR擴增而來渔扎，就會被GATK標記。

最主要目的就是盡量減小文庫構建時引入文庫的PCR bias

標記后對結(jié)果文件生成索引文件信轿。

GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatkref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fastasnp=/public/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gzindel=/public/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gzforsamplein{7E5239.L1,7E5240,7E5241.L1}doecho $sample$GATK --java-options"-Xmx20G -Djava.io.tmpdir=./"MarkDuplicates \? ? -I $sample.bam \? ? -O ${sample}_marked.bam \? ? -M $sample.metrics \1>${sample}_log.mark2>&1samtools index ${sample}_marked_fixed.bam

第二步：Local realignment around indels

這一步的目的就是將比對到indel附近的reads進行局部重新比對晃痴，將比對的錯誤率降到最低。一般來說财忽，絕大部分需要進行重新比對的基因組區(qū)域倘核，都是因為插入/缺失的存在，因為在indel附近的比對會出現(xiàn)大量的堿基錯配即彪，這些堿基的錯配很容易被誤認為SNP紧唱。還有，在比對過程中隶校，比對算法對于每一條read的處理都是獨立的漏益，不可能同時把多條reads與參考基因組比對來排錯。因此深胳，即使有一些reads能夠正確的比對到indel绰疤，但那些恰恰比對到indel開始或者結(jié)束位置的read也會有很高的比對錯誤率，這都是需要重新比對的舞终。

主要分為兩步：

通過運行RealignerTargetCreator來確定要進行重新比對的區(qū)域峦睡。

java-jarGenomeAnalysisTK.jar-Rhg19.fa-TRealignerTargetCreator-Ihg19.reorder.sort.addhead.dedup_04.bam-ohg19.dedup.realn_06.intervals-knownMills_and_1000G_gold_standard.indels.hg38.vcf-known1000G_phase1.indels.hg38.vcf參數(shù)說明：-R： 參考基因組翎苫；-T： 選擇的GATK工具；-I：? 輸入上一步所得bam文件榨了；-o： 輸出的需要重新比對的基因組區(qū)域結(jié)果煎谍；-maxInterval:允許進行重新比對的基因組區(qū)域的最大值，不能太大龙屉，太大耗費會很長時間呐粘，默認值500；-known:已知的可靠的indel位點转捕，指定已知的可靠的indel位點作岖，重比對將主要圍繞這些位點進行，對于人類基因組數(shù)據(jù)而言五芝，可以直接指定GATKresourcebundle里面的indel文件（必須是vcf文件）痘儡。

通過運行IndelRealigner在這些區(qū)域內(nèi)進行重新比對

-R hg19.fa-T IndelRealigner-targetIntervals hg19.dedup.realn_06.intervals-I hg19.reorder.sort.addhead.dedup_04.bam-o hg19.dedup.realn_07.bam-known Mills_and_1000G_gold_standard.indels.hg38.vcf-known1000G_phase1.indels.hg38.vcf

注意：1. 第一步和第二步中使用的輸入文件（bam文件）、參考基因組和已知indel文件必須是相同的文件枢步。

2. 當在相同的基因組區(qū)域發(fā)現(xiàn)多個indel存在時沉删，這個工具會從其中選擇一個最有可能存在比對錯誤的indel進行重新比對，剩余的其他indel不予考慮醉途。

第三步：.Base quality score recalibration

這一步是對bam文件里reads的堿基質(zhì)量值進行重新校正矾瑰，使最后輸出的bam文件中reads中堿基的質(zhì)量值能夠更加接近真實的與參考基因組之間錯配的概率。這一步適用于多種數(shù)據(jù)類型隘擎，包括illunima殴穴、solid、454货葬、CG等數(shù)據(jù)格式采幌。在GATK2.0以上版本中還可以對indel的質(zhì)量值進行校正，這一步對indel calling非常有幫助震桶。

舉例說明植榕，在reads堿基質(zhì)量值被校正之前，我們要保留質(zhì)量值在Q25以上的堿基尼夺，但是實際上質(zhì)量值在Q25的這些堿基的錯誤率在1%尊残，也就是說質(zhì)量值只有Q20，這樣就會對后續(xù)的變異檢測的可信度造成影響淤堵。還有寝衫，在邊合成邊測序的測序過程中，在reads末端堿基的錯誤率往往要比起始部位更高拐邪。另外慰毅，AC的質(zhì)量值往往要低于TG。BQSR的就是要對這些質(zhì)量值進行校正扎阶。

利用工具BaseRecalibrator汹胃，根據(jù)一些known sites婶芭，生成一個校正質(zhì)量值所需要的數(shù)據(jù)文件，GATK網(wǎng)站以“.grp”為后綴命名着饥。

-R$ref\? ? -I${sample}_marked_fixed.bam? \? ? --known-sites$snp\? ? --known-sites$indel\? ? -O${sample}_recal.table \? ? 1>${sample}_log.recal 2>&1

第四步：ApplyBQSR

利用已知snp數(shù)據(jù)庫犀农，通過機器學習方法調(diào)整堿基質(zhì)量分數(shù)

-R$ref\? ? -I${sample}_marked_fixed.bam? \? ? -bqsr${sample}_recal.table \? ? -O${sample}_bqsr.bam \? ? 1>${sample}_log.ApplyBQSR? 2>&1

這時候，GATK流程文件準備工作結(jié)束宰掉，完成了variants calling所需要的所有準備呵哨，生成了用于下一步變異檢測的bam文件。

第五步：Variant Calling

GATK在這一步里面提供了兩個工具進行變異檢測——UnifiedGenotyper和HaplotypeCaller

HaplotypeCaller不支持Reduce之后的bam文件轨奄，因此孟害，當選擇使用HaplotypeCaller進行變異檢測時，不需要進行Reduce reads

-R$ref\? ? -I${sample}_bqsr.bam \? ? ? --dbsnp$snp\? ? ? -O${sample}_raw.vcf \? ? ? 1>${sample}_log.HC 2>&1

對輸入的bam文件中的所有樣本進行變異檢測挪拟，最后生成一個vcf文件挨务，vcf文件中會包含所有樣本的變異位點和基因型信息。但是現(xiàn)在所得到的結(jié)果是最原始的玉组、沒有經(jīng)過任何過濾和校正的Variants集合谎柄。這一步產(chǎn)生的變異位點會有很高的假陽性，尤其是indel球切，因此，必須要進行進一步的篩選過濾绒障。這一步還可以指定對基因組的某一區(qū)域進行變異檢測吨凑，只需要增加一個參數(shù) -L：target_interval.list，target_interval.list 格式是bed格式文件户辱。

bed文件格式鸵钝，三列必須

TCGAGA#對應bed文件的坐標應為#chrome start endchr1? ? ? ? ? ? 0? ? 5

注意：GATK進行變異檢測的時候，是按照染色體排序順序進行的（先call chr1庐镐，然后chr2恩商，然后chr3…最后chrY），并非多條染色體并行檢測的必逆，因此怠堪，如果數(shù)據(jù)量比較大的話，建議分染色體分別進行名眉，對性染色體的變異檢測可以同常染色體方法粟矿。

大多數(shù)參數(shù)的默認值可以滿足大多數(shù)研究的需求，因此损拢，在做變異檢測過程中陌粹，如果對參數(shù)意義不是很明確，不建議修改福压。

第六步： 對原始變異檢測結(jié)果進行過濾（hard filter and VQSR）

這一步的目的就是對上一步call出來的變異位點進行過濾掏秩，去掉不可信的位點或舞。這一步可以有兩種方法，一種是通過GATK的VariantFiltration蒙幻，另一種是通過GATK的VQSR（變異位點質(zhì)量值重新校正）進行過濾映凳。

GATK是推薦使用VQSR的，但使用VQSR數(shù)據(jù)量一定要達到要求杆煞，數(shù)據(jù)量太小無法使用高斯模型魏宽。在使用VAQR時，indel和snp要分別進行决乎。

VQSR原理介紹：

這個模型是根據(jù)已有的真實變異位點（人類基因組一般使用HapMap3中的位點队询，以及這些位點在Omni 2.5M SNP芯片中出現(xiàn)的多態(tài)位點）來訓練，最后得到一個訓練好的能夠很好的評估真?zhèn)蔚腻e誤評估模型构诚，可以叫他適應性錯誤評估模型蚌斩。這個適應性的錯誤評估模型可以應用到call出來的原始變異集合中已知的變異位點和新發(fā)現(xiàn)的變異位點，進而去評估每一個變異位點發(fā)生錯誤的概率范嘱，最終會給出一個得分送膳。這個得分最后會被寫入vcf文件的INFO信息里，叫做VQSLOD丑蛤，就是在訓練好的混合高斯模型下叠聋，一個位點是真實的概率比上這個位點可能是假陽性的概率的log odds ratio（對數(shù)差異比），因此受裹，可以定性的認為碌补，這個值越大就越好。

VQSR主要分兩個步驟棉饶，這兩個步驟會使用兩個不同的工具：VariantRecalibrator和ApplyRecalibration厦章。

VariantRecalibrator：通過大量的高質(zhì)量的已知變異集合的各個注釋（包括很多種，后面介紹）的值來創(chuàng)建一個高斯混合模型照藻，然后用于評估所有的變異位點袜啃。這個文件最后將生成一個recalibration文件。

原理簡單介紹： 這個模型首先要拿到真實變異數(shù)據(jù)集和上一步驟中得到的原始變異數(shù)據(jù)集的交集幸缕，然后對這些SNP值相對于具體注釋信息的分布情況進行模擬群发，將這些變異位點進行聚類，最后根據(jù)聚類結(jié)果賦予所有變異位點相應的VQSLOD值发乔。越接近聚類核心的變異位點得到的VQSLOD值越高也物。

ApplyRecalibration：這一步將模型的各個參數(shù)應用于原始vcf文件中的每一個變異位點，這時列疗，每一個變異位點的注釋信息列中都會出現(xiàn)一個VQSLOD值滑蚯，然后模型會根據(jù)這個值對變異位點進行過濾，過濾后的信息會寫在vcf文件的filter一列中。

原理簡單介紹：在VariantRecalibrator這一步中告材，每個變異位點已經(jīng)得到了一個VQSLOD值了坤次，同時，這些LOD值在訓練集里也進行了排序斥赋。當你在這一步中設置一個tranche sensitivity 的閾值（這個閾值一般是一個百分數(shù)缰猴，如設置成99%），那么疤剑，如果LOD值從大到小排序的話滑绒，這個程序就會認為在這個訓練集中，LOD值在前99%的是可信的隘膘，當這個值低于這個閾值疑故，就認為是錯誤的。最后弯菊，程序就會用這個標準來過濾上一步call出來的原始變異集合纵势。如果LOD值超過這個閾值，在filter那一列就會顯示PASS管钳，如果低于這個值就會被過濾掉钦铁，但是這些位點仍然會顯示在結(jié)果里面，只不過會在filter那一列標示出他所屬于的tranche sensitivity 的名稱才漆。在設置tranche sensitivity 的閾值時牛曹，要兼顧敏感度和質(zhì)量值。

tranche值的設定

前面提到了醇滥，這個值得設定是用來在后續(xù)的ApplyRecalibration中如何根據(jù)這個閾值來過濾變異位點的黎比，也就是說，如果這個值設定的比較高的話腺办，那么最后留下來的變異位點就會多焰手，但同時假陽性的位點也會相應增加糟描；如果設定的低的話怀喉，雖然假陽性會減少，但是會丟失很多真實的位點船响。因此躬拢，跟選擇注釋時一樣，可以run兩遍VariantRecalibrator见间，第一遍的時候多寫幾個閾值聊闯，第一遍跑完之后看結(jié)果，看那個閾值好米诉，選擇一個最好的閾值菱蔬，再run一遍VariantRecalibrator。至于說怎么區(qū)分好壞，有幾個標準：

看結(jié)果中已知變異位點與新發(fā)現(xiàn)變異位點之間的比例拴泌，這個比例不要太大魏身，因為大多數(shù)新發(fā)現(xiàn)的變異都是假陽性，如果太多的話蚪腐，可能假陽性的比例就比較大箭昵；

看保留的變異數(shù)目，這個就要根據(jù)具體的需求進行選擇了回季。

看TI/TV值家制，對于人類全基因組，這個值應該在2.15左右泡一，對于外顯子組颤殴，這個值應該在3.2左右，不要太小或太大瘾杭，越接近這個數(shù)值越好诅病，這個值如果太小，說明可能存在比較多的假陽性粥烁。

千人中所選擇的tranche值是99

注意：Indel不支持tranche值的選擇贤笆，另外，一部分注釋類型在做indel的校正時也不支持讨阻，具體信息可以詳查GATK網(wǎng)站芥永。

當數(shù)據(jù)量太小時，可能高斯模型不會運行钝吮，因為變異位點數(shù)滿足不了模型的統(tǒng)計需求埋涧。這時候可以通過降低--maxGaussian的值，讓程序運行奇瘦。這個值表示的是程序?qū)⒆儺愇稽c分成的最大的組數(shù)棘催，降低這個值讓程序把變異位點聚類到更少的組里面，使每個組中的變異位點數(shù)增加來滿足統(tǒng)計需求耳标，但是這樣做降低程序分辨真?zhèn)蔚哪芰Υ及印Ｒ虼耍谶\行程序的時候次坡，要對各方面進行權衡呼猪。

過濾后生成的vcf文件的格式說明，即每一列所代表的的內(nèi)容砸琅，可參考下面的網(wǎng)站宋距，有詳細的說明

http://www.broadinstitute.org/gatk/guide/article?id=1268

第七步：初步注釋分析利用GATK中的VariantEval來完成，或者用annovar注釋

ANNOVAR是由王凱編寫的一個注釋軟件症脂，可以對SNP和indel進行注釋谚赎，也可以進行變異的過濾篩選淫僻。

ANNOVAR能夠利用最新的數(shù)據(jù)來分析各種基因組中 的遺傳變異。主要包含三種不同的注釋方法壶唤，Gene-based Annotation（基于基因的注釋）嘁傀、Region-based Annotation（基于區(qū)域的注釋）、Filter-based Annotation（基于篩選的注釋）视粮。ANNOVAR由Perl編寫细办。

優(yōu)點：提供多個數(shù)據(jù)可直接下載、支持多種格式蕾殴、注釋直觀笑撞；缺點：沒有數(shù)據(jù)庫的物種無法注釋。

│? annotate_variation.pl#主程序钓觉，功能包括下載數(shù)據(jù)庫茴肥，三種不同的注釋│? coding_change.pl#可用來推斷蛋白質(zhì)序列│? convert2annovar.pl#將多種格式轉(zhuǎn)為.avinput的程序│? retrieve_seq_from_fasta.pl#用于自行建立其他物種的轉(zhuǎn)錄本│? table_annovar.pl#注釋程序，可一次性完成三種類型的注釋│? variants_reduction.pl#可用來更靈活地定制過濾注釋流程│├─example#存放示例文件│└─humandb#人類注釋數(shù)據(jù)庫

ANNOVAR下載數(shù)據(jù)庫

命令示例

# -buildver 表示version# -downdb 下載數(shù)據(jù)庫的指令# -webfrom annovar 從annovar提供的鏡像下載荡灾，不加此參數(shù)將尋找數(shù)據(jù)庫本身的源# humandb/ 存放于humandb/目錄下

輸出的csv文件將包含輸入的5列主要信息以及各個數(shù)據(jù)庫里的注釋瓤狐，此外，table_annoval.pl可以直接對vcf文件進行注釋（不需要轉(zhuǎn)換格式）批幌，注釋的內(nèi)容將會放在vcf文件的“INFO”那一欄础锐。

說明：本文主要說明wes分析流程的理論，代碼僅作示例荧缘。

參考文獻

1.http://www.reibang.com/p/49d035b121b8

2.https://blog.csdn.net/zhu_si_tao/article/details/53321374

3.https://www.plob.org/article/9976.html

作者：鳳凰_0949

鏈接：http://www.reibang.com/p/f9ac13d3c938

來源：簡書

簡書著作權歸作者所有皆警，任何形式的轉(zhuǎn)載都請聯(lián)系作者獲得授權并注明出處。