?要點(diǎn)
1. GATK-HaplotypeCaller簡(jiǎn)介
2. GATK-HaplotypeCaller的基本組裝原理
3. GATK-HaplotypeCaller的安裝和使用
4. GATK-HaplotypeCaller實(shí)戰(zhàn)
hello，大家好，今天為大家?guī)?lái)關(guān)于變異檢測(cè)工具GATK-HaplotypeCaller超詳細(xì)安裝及應(yīng)用教程。
我們將持續(xù)為大家?guī)?lái)生物醫(yī)療大數(shù)據(jù)分析一文詳解系列文章，歡迎大家關(guān)注并星標(biāo)我們，可以更及時(shí)的看到我們的文章哦。
1.GATK-HaplotypeCaller簡(jiǎn)介
基因組變異檢測(cè)是基因組學(xué)領(lǐng)域一個(gè)非常重要的問(wèn)題，是遺傳性疾病溯源，物種進(jìn)化等分析的前提。而目前最主流、使用最廣泛的變異檢測(cè)軟件當(dāng)屬 Broad Institute 開發(fā)的 GATK(Genome Analysis ToolKit) 組件。GATK 設(shè)計(jì)之初是用于分析人類的全外顯子和全基因組數(shù)據(jù)，隨著不斷發(fā)展，現(xiàn)在也可以應(yīng)用于其他的物種。
GATK官網(wǎng)提供了一整套完整的變異檢測(cè)分析流程：GATK Best Practices。如下
圖示：
其中，HaplotypeCaller 是 GATK 檢測(cè)變異（SNP/INDEL）的核心模塊，主要通過(guò)單倍型的局部重組來(lái)實(shí)現(xiàn)準(zhǔn)確的 SNP 和 INDEL 檢測(cè)。GATK-HaplotypeCaller參考網(wǎng)址[1]
HaplotypeCaller 能夠從頭組裝變異活躍的區(qū)域（active region）并進(jìn)行 SNP/INDEL檢測(cè)。即每當(dāng)程序遇到顯示出變異跡象的區(qū)域時(shí)，它就會(huì)丟棄現(xiàn)有的比對(duì)信息并完全重新組裝該區(qū)域中的reads。因此 HaplotypeCaller 在一些傳統(tǒng)方法難以檢測(cè)的區(qū)域上會(huì)更加準(zhǔn)確，也使得在檢測(cè) indel 方面比 UnifiedGenotyper(gatk舊版變異檢測(cè)模塊) 等基于位置的變異檢測(cè)工具效果更好。
HaplotypeCaller 能夠處理非二倍體生物以及混合的實(shí)驗(yàn)數(shù)據(jù)，還能夠正確處理可變剪切，可以適用 RNAseq數(shù)據(jù)。
2. GATK-HaplotypeCaller的變異檢測(cè)的基本原理
GATK-HaplotypeCaller 模塊進(jìn)行 SNP/indel 檢測(cè)的基本工作流程包含四個(gè)主要步驟：?
1) 識(shí)別活躍區(qū)域?
2) 通過(guò)重組裝活躍區(qū)域確定單體型?
3) 確定每個(gè)read的單倍型的似然值?
4) 確定基因型。
2.1 識(shí)別活躍區(qū)域
沿著參考基因組以一定的窗口滑動(dòng)，統(tǒng)計(jì)比對(duì)的 mismatches, indels 和 softclips等信息計(jì)算基因組每個(gè)位置的活躍得分，使用平滑算法進(jìn)行處理，此處相當(dāng)于測(cè)定該區(qū)域熵值。當(dāng)熵值達(dá)到某個(gè)設(shè)定的閾值時(shí)即確定該區(qū)域?yàn)閍ctive region，用于后續(xù)組裝。
2.2 通過(guò)重組裝活躍區(qū)域確定單體型
對(duì)于每個(gè)活動(dòng)區(qū)域，忽略之前的read比對(duì)結(jié)果，重新利用該區(qū)域的reads構(gòu)建一個(gè)類似 De Bruijn 的圖來(lái)組裝活躍區(qū)域并識(shí)別數(shù)據(jù)中可能存在的單倍型。然后，使用 Smith-Waterman 算法將每個(gè)單倍型與參考單倍型重新對(duì)齊，以識(shí)別潛在的變異位點(diǎn)。如下圖示：最優(yōu)的路徑通過(guò)構(gòu)圖的方式進(jìn)行打分，得到候選的單體型路徑。
2.3 確定每個(gè)read的單倍型的似然值
對(duì)于每個(gè)活動(dòng)區(qū)域，程序使用 PairHMM 算法將每個(gè)read與每個(gè)單倍型進(jìn)行成對(duì)比對(duì), 產(chǎn)生一個(gè)單倍型似然值矩陣。然后將這些似然值邊緣化以獲得給定reads的每個(gè)潛在變異位點(diǎn)的等位基因可能性。
2.4 確定基因型
將前面 PairHMM 步驟得到的候選單倍型的似然值應(yīng)用貝葉斯算法轉(zhuǎn)化為每個(gè)位點(diǎn)基因型的似然值，如下圖所示：
3. GATK-HaplotypeCaller的安裝和使用
目前GATK已更新到GATK4。和GATK3相比，GATK4在算法上進(jìn)行了優(yōu)化，運(yùn)行速率有所提高，而且整合了picard 軟件的功能。github地址[2].
3.1. 安裝
GATK軟件運(yùn)行有兩種方式一種是運(yùn)行通過(guò)java 調(diào)用jar包運(yùn)行，一種是使用gatk腳本文件，此時(shí)需要安裝Python 2.6以上版本，不管是哪種方式，都需要機(jī)器上安裝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)兩個(gè)源碼包，不然還需要重新build gatk。當(dāng)然如果你的環(huán)境無(wú)法支持那就另說(shuō)了。下載完解壓即可看到可執(zhí)行程序gatk, 以及對(duì)應(yīng)的jar包，可以運(yùn)行./gatk --list測(cè)試下是否可運(yùn)行。
3.2. 使用說(shuō)明
運(yùn)行./gatk HaplotypeCaller -h查看 HaplotypeCaller 的參數(shù)細(xì)節(jié)，詳細(xì)說(shuō)明到官網(wǎng)查看會(huì)更清晰一些https://gatk.broadinstitute.org/hc/en-us/articles/360040096812-HaplotypeCaller
盡管HaplotypeCaller官網(wǎng)參數(shù)非常多，但實(shí)際用上的卻不多，大部分按默認(rèn)參數(shù)即最佳，這里列舉常用參數(shù)進(jìn)行說(shuō)明：
--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 指定參考序列，需要和比對(duì)時(shí)候使用的參考基因組一致--intervals -L? ? []? ? #One or more genomic intervals over which to operate 指定變異檢測(cè)的區(qū)間，可以是bed文件也可以是染色體名字--emit-ref-confidence -ERC? ? NONE? ? # Mode for emitting reference confidence scores 包含以下三種變異檢測(cè)模式：? ? #1. NONE：只記錄變異位點(diǎn)? ? #2. BP_RESOLUTION：記錄變異和非變異位點(diǎn)，每個(gè)位點(diǎn)信息都展示????#3. GVCF：記錄變異和非變異位點(diǎn)，其中非變異位點(diǎn)以block塊展示--pcr-indel-model CONSERVATIVE? ? #The PCR indel model to use設(shè)置針對(duì)PCR indel的矯正模型，包含以下幾種模式? ? #1. NONE：不會(huì)應(yīng)用專門的 PCR 錯(cuò)誤模型；如果存在堿基插入/刪除質(zhì)量，它們將被使用? ? #2. HOSTILE：將應(yīng)用一個(gè)最激進(jìn)的模型，該模型會(huì)犧牲真陽(yáng)性以消除更多的假陽(yáng)性? ? #3. AGGRESSIVE：將應(yīng)用相對(duì)激進(jìn)的模型，犧牲真陽(yáng)性以消除更多的假陽(yáng)性? ? #4. CONSERVATIVE：將應(yīng)用一個(gè)不那么激進(jìn)的模型，該模型試圖以允許更多誤報(bào)為代價(jià)來(lái)保持較高的真陽(yáng)性率
注意：
1）對(duì)于隊(duì)列分析，通常會(huì)使用-ERC GVCF參數(shù)來(lái)生成gvcf文件以支持下游多樣本聯(lián)合分析
2）對(duì)于PCR-free的樣本需要使用-pcr_indel_model NONE取消PCR indel模型。
運(yùn)行命令如下：
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實(shí)戰(zhàn)
下面我們找個(gè)NA12878樣本的bam 文件數(shù)據(jù)，具體來(lái)實(shí)踐一下吧。
4.1 數(shù)據(jù)下載
下載運(yùn)行所需的數(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 運(yùn)行命令
此處我們使用Gvcf模式，且只檢測(cè)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é)果
運(yùn)行完即可看到chrM.g.vcf.gz文件的生成。
此外，我們的sixoclock官網(wǎng)基于CWL (common workflow language) 對(duì) GATK-HaplotypeCaller 軟件進(jìn)行了封裝，通過(guò)我們開發(fā)的sixbox軟件可以快速進(jìn)行軟件的運(yùn)行。對(duì)sixbox不了解可以通過(guò)六點(diǎn)了官網(wǎng)了https://docs.sixoclock.net/clients/sixbox-linux.html解下。下面是具體的運(yùn)行步驟如下：
1)下載cwl 源碼sixbox pull 2e932c25-8c36-4733-bb91-79a137282013或 點(diǎn)擊下載按鈕下載 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ù)以實(shí)現(xiàn)個(gè)性化分析
可以直接粘貼下方示例內(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運(yùn)行
sixbox run ./HaplotypeCaller.cwl ./HaplotypeCaller.job.yaml#或者指定輸出目錄sixbox run --outdir /home/path ./HaplotypeCaller.cwl ./HaplotypeCaller.job.yaml
運(yùn)行結(jié)果即可看到當(dāng)前目錄或者指定的輸出目錄輸出chrM.g.vcf.gz結(jié)果。
以上為我們給大家?guī)?lái)的宏基因組組裝工具M(jìn)egahit基本原理知識(shí)，以及運(yùn)行詳細(xì)操作過(guò)程。
如果對(duì)生物醫(yī)療健康大數(shù)據(jù)相關(guān)內(nèi)容感興趣也可以持續(xù)關(guān)注我們。想要探索更多生物醫(yī)療大數(shù)據(jù)分析工具和軟件的介紹和使用請(qǐng)看六點(diǎn)了官網(wǎng)[3]。（??點(diǎn)擊閱讀原文）
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
推薦閱讀
?一文詳解基因組denovo組裝原理和實(shí)戰(zhàn)
?首發(fā)！生信分析--入門實(shí)踐一條龍的CWL中文教程來(lái)了
歡迎關(guān)注我們的公眾號(hào)“六點(diǎn)了”，原創(chuàng)生物醫(yī)療大數(shù)據(jù)分析技術(shù)文章第一時(shí)間推送。
六點(diǎn)了?
生物醫(yī)療數(shù)據(jù)處理云平臺(tái)，創(chuàng)建，使用，分享，團(tuán)隊(duì)協(xié)作。提供海量數(shù)據(jù)處理算法下載、可視化編輯與配置，在線流程組合、自動(dòng)生成功能，以及基于解耦合的一站式生信解決方案私有云的本地化部署服務(wù)。
