?要點
1. GATK-HaplotypeCaller簡介
2. GATK-HaplotypeCaller的基本組裝原理
3. GATK-HaplotypeCaller的安裝和使用
4. GATK-HaplotypeCaller實戰(zhàn)
hello废麻，大家好川抡，今天為大家?guī)黻P(guān)于變異檢測工具GATK-HaplotypeCaller超詳細安裝及應(yīng)用教程俘陷。
我們將持續(xù)為大家?guī)砩镝t(yī)療大數(shù)據(jù)分析一文詳解系列文章，歡迎大家關(guān)注并星標我們掉伏，可以更及時的看到我們的文章哦。
1.GATK-HaplotypeCaller簡介
基因組變異檢測是基因組學(xué)領(lǐng)域一個非常重要的問題，是遺傳性疾病溯源仆邓，物種進化等分析的前提。而目前最主流台腥、使用最廣泛的變異檢測軟件當(dāng)屬 Broad Institute 開發(fā)的 GATK(Genome Analysis ToolKit) 組件宏赘。GATK 設(shè)計之初是用于分析人類的全外顯子和全基因組數(shù)據(jù)，隨著不斷發(fā)展黎侈，現(xiàn)在也可以應(yīng)用于其他的物種察署。
GATK官網(wǎng)提供了一整套完整的變異檢測分析流程：GATK Best Practices。如下
圖示：
其中峻汉，HaplotypeCaller 是 GATK 檢測變異（SNP/INDEL）的核心模塊贴汪，主要通過單倍型的局部重組來實現(xiàn)準確的 SNP 和 INDEL 檢測脐往。GATK-HaplotypeCaller參考網(wǎng)址[1]
HaplotypeCaller 能夠從頭組裝變異活躍的區(qū)域（active region）并進行 SNP/INDEL檢測。即每當(dāng)程序遇到顯示出變異跡象的區(qū)域時扳埂，它就會丟棄現(xiàn)有的比對信息并完全重新組裝該區(qū)域中的reads业簿。因此 HaplotypeCaller 在一些傳統(tǒng)方法難以檢測的區(qū)域上會更加準確，也使得在檢測 indel 方面比 UnifiedGenotyper(gatk舊版變異檢測模塊) 等基于位置的變異檢測工具效果更好阳懂。
HaplotypeCaller 能夠處理非二倍體生物以及混合的實驗數(shù)據(jù)梅尤，還能夠正確處理可變剪切，可以適用 RNAseq數(shù)據(jù)岩调。
2. GATK-HaplotypeCaller的變異檢測的基本原理
GATK-HaplotypeCaller 模塊進行 SNP/indel 檢測的基本工作流程包含四個主要步驟：?
1) 識別活躍區(qū)域?
2) 通過重組裝活躍區(qū)域確定單體型?
3) 確定每個read的單倍型的似然值?
4) 確定基因型巷燥。
2.1 識別活躍區(qū)域
沿著參考基因組以一定的窗口滑動，統(tǒng)計比對的 mismatches, indels 和 softclips等信息計算基因組每個位置的活躍得分号枕，使用平滑算法進行處理缰揪，此處相當(dāng)于測定該區(qū)域熵值。當(dāng)熵值達到某個設(shè)定的閾值時即確定該區(qū)域為active region葱淳，用于后續(xù)組裝钝腺。
2.2 通過重組裝活躍區(qū)域確定單體型
對于每個活動區(qū)域，忽略之前的read比對結(jié)果赞厕，重新利用該區(qū)域的reads構(gòu)建一個類似 De Bruijn 的圖來組裝活躍區(qū)域并識別數(shù)據(jù)中可能存在的單倍型艳狐。然后，使用 Smith-Waterman 算法將每個單倍型與參考單倍型重新對齊坑傅，以識別潛在的變異位點僵驰。如下圖示：最優(yōu)的路徑通過構(gòu)圖的方式進行打分，得到候選的單體型路徑唁毒。
2.3 確定每個read的單倍型的似然值
對于每個活動區(qū)域蒜茴，程序使用 PairHMM 算法將每個read與每個單倍型進行成對比對, 產(chǎn)生一個單倍型似然值矩陣。然后將這些似然值邊緣化以獲得給定reads的每個潛在變異位點的等位基因可能性浆西。
2.4 確定基因型
將前面 PairHMM 步驟得到的候選單倍型的似然值應(yīng)用貝葉斯算法轉(zhuǎn)化為每個位點基因型的似然值粉私，如下圖所示：
3. GATK-HaplotypeCaller的安裝和使用
目前GATK已更新到GATK4。和GATK3相比近零，GATK4在算法上進行了優(yōu)化诺核，運行速率有所提高，而且整合了picard 軟件的功能久信。github地址[2].
3.1. 安裝
GATK軟件運行有兩種方式一種是運行通過java 調(diào)用jar包運行窖杀，一種是使用gatk腳本文件暂题，此時需要安裝Python 2.6以上版本蝌戒，不管是哪種方式纤泵，都需要機器上安裝java 8或以上版本预茄。
wget https://github.com/broadinstitute/gatk/releases/download/4.1.5.0/gatk-4.1.5.0.zip unzip gatk-4.1.5.0.zip
注意此處下載最好是直接下載gatk-**.zip包怀吻，而不要下載Source code (zip) 和 Source code (tar.gz)兩個源碼包思灌，不然還需要重新build gatk蝉揍。當(dāng)然如果你的環(huán)境無法支持那就另說了趋箩。下載完解壓即可看到可執(zhí)行程序gatk, 以及對應(yīng)的jar包，可以運行./gatk --list測試下是否可運行铆隘。
3.2. 使用說明
運行./gatk HaplotypeCaller -h查看 HaplotypeCaller 的參數(shù)細節(jié)卓舵，詳細說明到官網(wǎng)查看會更清晰一些https://gatk.broadinstitute.org/hc/en-us/articles/360040096812-HaplotypeCaller
盡管HaplotypeCaller官網(wǎng)參數(shù)非常多，但實際用上的卻不多膀钠，大部分按默認參數(shù)即最佳，這里列舉常用參數(shù)進行說明：
--input-I []? ? # BAM/SAM/CRAM file containing reads 指定輸入的bam托修、sam、cram文件--output -O? ? null? ? # File to which variants should be written 指定輸出vcf名字--reference -R? ? null? ? # Reference sequence file 指定參考序列睦刃，需要和比對時候使用的參考基因組一致--intervals -L? ? []? ? #One or more genomic intervals over which to operate 指定變異檢測的區(qū)間十酣，可以是bed文件也可以是染色體名字--emit-ref-confidence -ERC? ? NONE? ? # Mode for emitting reference confidence scores 包含以下三種變異檢測模式：? ? #1. NONE：只記錄變異位點? ? #2. BP_RESOLUTION：記錄變異和非變異位點，每個位點信息都展示????#3. GVCF：記錄變異和非變異位點耸采，其中非變異位點以block塊展示--pcr-indel-model CONSERVATIVE? ? #The PCR indel model to use設(shè)置針對PCR indel的矯正模型兴泥，包含以下幾種模式? ? #1. NONE：不會應(yīng)用專門的 PCR 錯誤模型；如果存在堿基插入/刪除質(zhì)量虾宇，它們將被使用? ? #2. HOSTILE：將應(yīng)用一個最激進的模型搓彻，該模型會犧牲真陽性以消除更多的假陽性? ? #3. AGGRESSIVE：將應(yīng)用相對激進的模型，犧牲真陽性以消除更多的假陽性? ? #4. CONSERVATIVE：將應(yīng)用一個不那么激進的模型嘱朽，該模型試圖以允許更多誤報為代價來保持較高的真陽性率
注意：
1）對于隊列分析旭贬，通常會使用-ERC GVCF參數(shù)來生成gvcf文件以支持下游多樣本聯(lián)合分析
2）對于PCR-free的樣本需要使用-pcr_indel_model NONE取消PCR indel模型。
運行命令如下：
gatk --java-options "-Xmx4g" HaplotypeCaller? \? -R Homo_sapiens_assembly38.fasta \? -I input.bam \? -O output.g.vcf.gz \? -ERC GVCF \? -L chr1
4. GATK-HaplotypeCaller實戰(zhàn)
下面我們找個NA12878樣本的bam 文件數(shù)據(jù)搪泳，具體來實踐一下吧稀轨。
4.1 數(shù)據(jù)下載
下載運行所需的數(shù)據(jù)：參考基因組及其索引，bam文件
wget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.fastawget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.dictwget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.fasta.faiwget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/hg19.chrM.bwa-mem.sort.rmdup.BQSR.bamwget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/hg19.chrM.bwa-mem.sort.rmdup.BQSR.bai
4.2 運行命令
此處我們使用Gvcf模式岸军，且只檢測chrM的變異
gatk --java-options "-Xmx4g" HaplotypeCaller? \? -R hg19.chrM.fasta \? -I hg19.chrM.bwa-mem.sort.rmdup.BQSR.bam \? -O chrM.g.vcf.gz \? -ERC GVCF \? -L chrM
4.3 輸出結(jié)果
運行完即可看到chrM.g.vcf.gz文件的生成奋刽。
此外，我們的sixoclock官網(wǎng)基于CWL (common workflow language) 對 GATK-HaplotypeCaller 軟件進行了封裝艰赞，通過我們開發(fā)的sixbox軟件可以快速進行軟件的運行佣谐。對sixbox不了解可以通過六點了官網(wǎng)了https://docs.sixoclock.net/clients/sixbox-linux.html解下。下面是具體的運行步驟如下：
1)下載cwl 源碼sixbox pull 2e932c25-8c36-4733-bb91-79a137282013或 點擊下載按鈕下載 gatk_HaplotypeCaller.cwl https://www.sixoclock.net/application/pipe/2e932c25-8c36-4733-bb91-79a137282013
2) 使用sixbox生成參數(shù)模板文件(YAML) , 并配置yaml文件
sixbox run --make-template ./HaplotypeCaller.cwl > HaplotypeCaller.job.yamlvim HaplotypeCaller.job.yaml # 編輯參數(shù)配置文件方妖，替換或設(shè)置參數(shù)以實現(xiàn)個性化分析
可以直接粘貼下方示例內(nèi)容到HaplotypeCaller.job.yaml
reference:? # type "File"? ? class: File? ? path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.fasta? ? secondaryFiles:? ? ? - class: File? ? ? ? path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.dict? ? ? - class: File? ? ? ? path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta_BwaMem_index/hg19.chrM.fasta.faioutput_vcf: chrM.g.vcf.gz? # type "string"#java_options: '-Xmx6g -XX:ParallelGCThreads=4'? # type "string"intervals:? # array of type "string" (optional)? - "chrM"input_bam:? # array of type "File"? ? class: File? ? path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/hg19.chrM.bwa-mem.sort.rmdup.BQSR.bam? ? secondaryFiles:? ? ? - class: File? ? ? ? path:? http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/hg19.chrM.bwa-mem.sort.rmdup.BQSR.baiERC: 'GVCF'? # type "string"
3) 使用sixbox運行
sixbox run ./HaplotypeCaller.cwl ./HaplotypeCaller.job.yaml#或者指定輸出目錄sixbox run --outdir /home/path ./HaplotypeCaller.cwl ./HaplotypeCaller.job.yaml
運行結(jié)果即可看到當(dāng)前目錄或者指定的輸出目錄輸出chrM.g.vcf.gz結(jié)果狭魂。
以上為我們給大家?guī)淼暮昊蚪M組裝工具Megahit基本原理知識，以及運行詳細操作過程。
如果對生物醫(yī)療健康大數(shù)據(jù)相關(guān)內(nèi)容感興趣也可以持續(xù)關(guān)注我們趁蕊。想要探索更多生物醫(yī)療大數(shù)據(jù)分析工具和軟件的介紹和使用請看六點了官網(wǎng)[3]坞生。（??點擊閱讀原文）
References
[1]GATK-HaplotypeCaller參考網(wǎng)址:https://gatk.broadinstitute.org/hc/en-us/articles/360037225632-HaplotypeCaller
[2]github地址:https://github.com/broadinstitute/gatk
[3]?https://www.sixoclock.net
推薦閱讀
?首發(fā)！生信分析--入門實踐一條龍的CWL中文教程來了
歡迎關(guān)注我們的公眾號“六點了”掷伙，原創(chuàng)生物醫(yī)療大數(shù)據(jù)分析技術(shù)文章第一時間推送是己。
六點了?
生物醫(yī)療數(shù)據(jù)處理云平臺，創(chuàng)建任柜，使用卒废，分享，團隊協(xié)作宙地。提供海量數(shù)據(jù)處理算法下載摔认、可視化編輯與配置，在線流程組合宅粥、自動生成功能参袱，以及基于解耦合的一站式生信解決方案私有云的本地化部署服務(wù)。
