「GATK 4」如何提高HaplotyperCaller的效率

GATK的HaplotypeCaller 应该是目前最常用的变异检测软件，尤其是在人类基因组上。不过HaplotypeCaller的速度相对于其他软件，例如bcftools, freeBayes 也是最慢的，当然这还是可以抢救一下的，只不过需要我们额外写一些代码，利用--intervals参数进行手动并行。

如下代码仅考虑单个样本，多个样本的gvcf文件类似

策略1: 按照染色体进行运算

最简单的方法也是最容易想到的方法，直接利把原来的任务按照染色体拆分。拟南芥是5+2, 5条核基因组染色体和2条细胞器基因组染色体，那么就可以一次性投递7个任务。

# $SM指的是你样本名(不包括名字后缀)
# $REF指的是你的参考基因组名字
chroms=($(grep '>' $REF |sed 's/>//' | tr '\n' ' '))
for chr in ${chroms[@]}
do
        if [ ! -f ${SM}.${chr}.vcf.gz ]; then
        gatk HaplotypeCaller -R $REF -I ${SM}_mkdup.bam --genotyping-mode DISCOVERY \
                --intervals ${chr} -stand-call-conf 30 --sample-ploidy 2 \
                -O ${SM}.${chr}.vcf.gz &
        fi
done && wait

策略2: 针对染色体的NNNN区进一步拆分

既然我们可以分成不同的染色体进行并行运算，那能不能进一步把染色体继续拆分呢？这个思路也是可行，毕竟真核生物的基因组大部分都不完美，也就是基因组存在着一些gap没有被填上。gap区域肯定是没有read能够回贴上的，那么就可以根据这些gap进一步的细分。

GAP区域

这里用到了seqkit locate 定位染色体上NNN， 然后用bedtools complement 得到后续输入的bed文件。

seqkit locate -j 20 -Pp '"N{10,}"' /data/reference/genome/TAIR10/Athaliana.fa | tail -n+2 | awk '{print $1"\t"$5"\t"$6}' > ath_n_pos.bed
bioawk -c fastx '{print $name"\t"length($seq)}' /data/reference/genome/TAIR10/Athaliana.fa > tair10.txt
bedtools complement -i ath_n_pos.bed -g tair10.txt  > tair10_interval.bed

之后就可以根据bed文件构建出interval

BED="tair10_interval.bed"
chroms=($(awk '{print $1":"$2+1"-"$3}' $BED | tr '\n' ' '))
for chr in ${chroms[@]}
do
        if [ ! -f ${SM}.${chr}.vcf.gz ]; then
        gatk HaplotypeCaller -R $REF -I ${SM}_mkdup.bam --genotyping-mode DISCOVERY \
                --intervals ${chr} -stand-call-conf 30 --sample-ploidy 2 \
                -O ${SM}.${chr}.vcf.gz &
        fi
done && wait

策略3: 根据BAM文件特点进行拆分

如果你还想追求速度，那么还可以对BAM文件进行分析，找出在比对后结果中哪些区域是没有read覆盖，那么将这些区域进行拆分。

低覆盖区间

我们可以用bedtools genomecov 统计基因组上低覆盖度区域, 选择其中宽度大于100 bp 或1000 bp(按照需求进行确定阈值)的区间作为分割点

bedtools genomecov -bga -ibam ${SM}.bam -g tair10.txt | awk '$4 < 3' | bedtools merge -i - | awk '{if ($3-$2 > 1000) print $0}' > result.bed 
bedtools complement -i  result.bed -g tair10.txt  > tair10_interval.bed

分分钟给你搞出上千个区间，你就可以疯狂的投递任务了。

也可以简单粗暴的按照固定大小对染色体进行划分。

最后运行结束之后，还需要对文件进行合并，用到的是MergeVcfs。

## merge chr.vcf
merge_vcfs=""
for chr in ${chroms[@]}; do
        merge_vcfs=${merge_vcfs}" -I ${SM}.${chr}.vcf.gz"
done && gatk MergeVcfs ${merge_vcfs} -O ${SM}.HC.vcf.gz && echo "Vcfs haved been successfully merged"