Variant calling using GATK

Question

[evelyn]

I am using GATK for variant calling. It works and give results but the number of SNPs found is significantly very low as compared to samtools. Is this fine or there is an error?

Answer 1

[Ram]

Short version: It is fine and expected, as they work differently. samtools is more sensitive compared to GATK's callers. More reading:

• Difference Between Samtools And Gatk Algorithms
• https://gatkforums.broadinstitute.org/gatk/discussion/2159/very-different-snp-calling-by-gatk-and-samtools

[evelyn]

I understand there will be difference in the number of SNPs but the vcf file size by gatk is 10,000 times less than the vcf file size by samtools. Is this large difference normal?

[Ram]

File size is not an indicator in most cases, but 10,000 times smaller does sound a little suspect. Can you check the GATK logs and ensure it ran to completion?

[evelyn]

Yes, it ran to completion and didn't give any error.

[Ram]

Just to be sure, they're both the same format, correct? Variants can be in VCF/BCF formats and either of those formats can be compressed/uncompressed. Does running htsfile <filename> give the same output on both? Plus, the one from GATK is not a gVCF file, correct?

If you could give us the exact commands you used, that would help a lot.

[evelyn]

I am using this code to call SNPs for DNA seq paired end files aligned using bowtie2:

Using GATK:

module load gatk/4.1.0.0 picard/2.9.2
picard CreateSequenceDictionary R=ref.fa O=ref.dict
java -jar picard.jar ValidateSamFile \
I=bowtie.sorted.bam \
MODE=SUMMARY
java -jar picard.jar AddOrReplaceReadGroups \
I=bowtie.sorted.bam \
O=bowtie_output.bam \
RGID=1 \
RGLB=lib1 \
RGPL=illumina \
RGPU=unit1 \
RGSM=20
java -jar picard.jar ValidateSamFile \
I=bowtie_output.bam \
MODE=SUMMARY
samtools index bowtie_output.bam
gatk --java-options "-Xmx4G" HaplotypeCaller -R ref.fa -I bowtie_output.bam -O bowtie_gatk.vcf

Using samtools:

bcftools mpileup -f ref.fa bowtie.sorted.bam > bowtie_samtools.vcf