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區(qū)進一步拆分

既然我們可以分成不同的染色體進行并行運算，那能不能進一步把染色體繼續(xù)拆分呢？這個思路也是可行，畢竟真核生物的基因組大部分都不完美，也就是基因組存在著一些gap沒有被填上。gap區(qū)域肯定是沒有read能夠回貼上的，那么就可以根據這些gap進一步的細分。

GAP區(qū)域

這里用到了seqkit locate 定位染色體上NNN， 然后用bedtools complement 得到后續(xù)輸入的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文件進行分析，找出在比對后結果中哪些區(qū)域是沒有read覆蓋，那么將這些區(qū)域進行拆分。

低覆蓋區(qū)間

我們可以用bedtools genomecov 統(tǒng)計基因組上低覆蓋度區(qū)域, 選擇其中寬度大于100 bp 或1000 bp(按照需求進行確定閾值)的區(qū)間作為分割點

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

分分鐘給你搞出上千個區(qū)間，你就可以瘋狂的投遞任務了。

也可以簡單粗暴的按照固定大小對染色體進行劃分。

最后運行結束之后，還需要對文件進行合并，用到的是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"