一咨堤、問題
在按照GATK流程做堿基質(zhì)量重校正時(shí)呆馁，需要用到gatk-bundle中下載的一些已知位點(diǎn)vcf文件秸抚，如：
- resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf
- resources_broad_hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf
- resources_broad_hg38_v0_Homo_sapiens_assembly38.known_indels.vcf
但是在運(yùn)行下述腳本時(shí)凄诞，提示--known-sites這里存在錯(cuò)誤斗锭。
ref_genome_fasta_file=~/database/ref_genome/human/hg38/gatk-bundle/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
dbsnp=~/database/ref_genome/human/hg38/gatk-bundle/resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf
Mills_1000G_indel=~/database/ref_genome/human/hg38/gatk-bundle/resources_broad_hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
assembly38_indel=~/database/ref_genome/human/hg38/gatk-bundle/resources_broad_hg38_v0_Homo_sapiens_assembly38.known_indels.vcf.gz
#### BQSR
echo "Start ${i} BQSR for ${id}." `date`
$GATK --java-options "-Xmx10G -Djava.io.tmpdir=./${dir}/tmp" BaseRecalibrator \
-R ${GENOME} \
-I ./${dir}/${id}_markdup.bam \
--known-sites ${dbsnp} \
--known-sites ${Mills_1000G_indel} \
--known-sites ${assembly38_indel} \
-O ./${dir}/${id}_BaseRecal.table \
1>./${dir}/${id}_BaseRecal.log 2>&1
echo -e "BQSR for ${id} has finised at:" `date` "\n"
二、解決方法
參考該帖后偿渡，發(fā)現(xiàn)原因是因?yàn)間atk-bundle下載的vcf文件是二進(jìn)制的臼寄，需要先用bcftools處理一下就可以正確運(yùn)行BQSR了。
bcftools view -O z -o tmp.vcf.gz resources_broad_hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf
mv tmp.vcf.gz resources_broad_hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
bcftools index -t resources_broad_hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
bcftools view -O z -o tmp.vcf.gz resources_broad_hg38_v0_Homo_sapiens_assembly38.known_indels.vcf
mv tmp.vcf.gz resources_broad_hg38_v0_Homo_sapiens_assembly38.known_indels.vcf.gz
bcftools index -t resources_broad_hg38_v0_Homo_sapiens_assembly38.known_indels.vcf.gz
