GATK的HaplotypeCaller 應(yīng)該是目前最常用的變異檢測(cè)軟件捂齐，尤其是在人類(lèi)基因組上醋粟。不過(guò)HaplotypeCaller的速度相對(duì)于其他軟件辆床，例如bcftools, freeBayes 也是最慢的乎莉，當(dāng)然這還是可以搶救一下的愚屁，只不過(guò)需要我們額外寫(xiě)一些代碼瑟由，利用--intervals參數(shù)進(jìn)行手動(dòng)并行絮重。
如下代碼僅考慮單個(gè)樣本，多個(gè)樣本的gvcf文件類(lèi)似
策略1: 按照染色體進(jìn)行運(yùn)算
最簡(jiǎn)單的方法也是最容易想到的方法歹苦，直接利把原來(lái)的任務(wù)按照染色體拆分青伤。擬南芥是5+2, 5條核基因組染色體和2條細(xì)胞器基因組染色體，那么就可以一次性投遞7個(gè)任務(wù)暂氯。
# $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: 針對(duì)染色體的NNNN區(qū)進(jìn)一步拆分
既然我們可以分成不同的染色體進(jìn)行并行運(yùn)算潮模，那能不能進(jìn)一步把染色體繼續(xù)拆分呢？這個(gè)思路也是可行痴施，畢竟真核生物的基因組大部分都不完美擎厢，也就是基因組存在著一些gap沒(méi)有被填上。gap區(qū)域肯定是沒(méi)有read能夠回貼上的辣吃，那么就可以根據(jù)這些gap進(jìn)一步的細(xì)分动遭。
這里用到了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
之后就可以根據(jù)bed文件構(gòu)建出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: 根據(jù)BAM文件特點(diǎn)進(jìn)行拆分
如果你還想追求速度厘惦，那么還可以對(duì)BAM文件進(jìn)行分析，找出在比對(duì)后結(jié)果中哪些區(qū)域是沒(méi)有read覆蓋哩簿，那么將這些區(qū)域進(jìn)行拆分宵蕉。
我們可以用bedtools genomecov 統(tǒng)計(jì)基因組上低覆蓋度區(qū)域, 選擇其中寬度大于100 bp 或1000 bp(按照需求進(jìn)行確定閾值)的區(qū)間作為分割點(diǎn)
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
分分鐘給你搞出上千個(gè)區(qū)間，你就可以瘋狂的投遞任務(wù)了节榜。
也可以簡(jiǎn)單粗暴的按照固定大小對(duì)染色體進(jìn)行劃分羡玛。
最后運(yùn)行結(jié)束之后，還需要對(duì)文件進(jìn)行合并宗苍，用到的是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"
