GATK之HaplotypeCaller

GATK的主要功能其實就是識別變異位點，其餘功能都是錦上添花。因此這一次學習GATK尋找變異位點的工具。
在GATK的文檔中，與變異位點識別相關的有9個工具，分別是：

Name	Summary
ApplyRecalibration	Apply a score cutoff to filter variants based on a recalibration table
CalculateGenotypePosteriors	Calculate genotype posterior likelihoods given panel data
GATKPaperGenotyper	Simple Bayesian genotyper used in the original GATK paper
GenotypeGVCFs	Perform joint genotyping on gVCF files produced by HaplotypeCaller
HaplotypeCaller	Call germline SNPs and indels via local re-assembly of haplotypes
MuTect2	Call somatic SNPs and indels via local re-assembly of haplotypes
RegenotypeVariants	Regenotypes the variants from a VCF containing PLs or GLs.
UnifiedGenotyper	Call SNPs and indels on a per-locus basis
VariantRecalibrator	Build a recalibration model to score variant quality for filtering purposes

稍微看了一下GATK最佳實踐文檔，發現目前用的最多的是HaplotypeCaller，準確率高，可是比較吃配置，因此運行時間會比較久。不過因爲HaplotypeCaller的工做原理，直接省去了BQSR和indel realignment步驟，因此對於一個variant calling流程而言，能夠直接比對，去重複後運行HaplotypeCaller。

簡介

功能： 經過局部單倍型重組裝尋找體細胞的SNP和indel
分類： 變異位點探索功能

HaplotypeCaller，簡稱HC，能過經過對活躍區域（也就是與參考基因組不一樣處較多的區域）局部重組裝，同時尋找SNP和INDEL。換句話說，當HC看到一個地方好活躍呀，他就無論以前的比對結果了，直接對這個地方進行從新組裝，因此就比傳統的基於位置（position-based)的工具，如前任UnifiedGenotyper好用多了。

HC能夠處理非二倍體物種和混池實驗數據，可是不推薦用於癌症或者體細胞。
HC還能夠正確處理可變剪切，因此也能用於RNA-Seq數據。

HC有一個GVCF模式，產生各個樣本的中間基因組文件(gVCF),而後用於多樣本的聯合基因定型(joint genotyping)，效率很高呢。

工做原理

HaplotypeCaller的核心操做就是四步：

1. 尋找活躍區域，就是和參考基因組不一樣部分較多的區域

2. 經過對該區域進行局部重組裝，肯定單倍型（haplotypes）。就是這一步能夠省去indel realignment

3. 在給定的read數據下，計算單倍型的可能性。

4. 分配樣本的基因型

    

 

使用說明

輸入：BAM文件
輸出：原始的VCF文件或者gVCF文件，未過濾SNP和INDEL。通常而言，在進行下游分析，須要對這些數據要麼手動要麼VQSR過濾。

案例:
DNA-seq variant-calling

 

  java -jar GenomeAnalysisTK.jar \
     -R reference.fasta \
     -T HaplotypeCaller \
     -I sample1.bam [-I sample2.bam ...] \
     [--dbsnp dbSNP.vcf] \
     [-stand_call_conf 30] \
     [-L targets.interval_list] \
     -o output.raw.snps.indels.vcf

RNA-seq variant-calling

 

  java -jar GenomeAnalysisTK.jar \
     -R reference.fasta \
     -T HaplotypeCaller \
     -I sample1.bam \
     [--dbsnp dbSNP.vcf] \
     -stand_call_conf 20 \
     -o output.raw.snps.indels.vcf

其餘感受比較使用的參數：

參數名	默認值	概要
—genotyping_mode/-gt_mode	DISCOVERY	如何處理識別的等位基因
—min_base_quality_score/ -mbq	10	最低鹼基質量
—sample_ploidy/-ploidy	2	樣本是幾倍體？
—standard_min_confidence_threshold_for_calling/-stand_call_conf	10.0	肯定variant的最低phred-scaled 置信閾值
—annotation/-A	-	variant註釋內容
—pcr_indel_model/-pcrModel	CONSERVATIVE	PCR indel模式

注：對於我作mapping-by-sequencing而言，須要結果有ref和alt鹼基的支持數，因此選項-A必定要跟上StrandAlleleCountsBySample。

這是我用於MBS的流程的GATK部分

 

##############################
# Vaiant Discovery with GATK #
##############################mkdir -p ${wd_dir}/analysis/variant_gatk
recal_files=$(ls *recal_reads.bam)
for file in $recal_files
do
    java -Xmx16g -jar $gatk -T HaplotypeCaller \
    -R $reference  -I $file --genotyping_mode DISCOVERY \
    --standard_min_confidence_threshold_for_calling 10 \
    -A StrandAlleleCountsBySample \
    -o ../variant_gatk/${file%%.*}_raw_variants.vcf
done