Question: Variant calling using GATK
1
gravatar for evelyn
27 days ago by
evelyn30
evelyn30 wrote:

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?

snp • 159 views
ADD COMMENTlink modified 27 days ago • written 27 days ago by evelyn30
0
gravatar for RamRS
27 days ago by
RamRS22k
Houston, TX
RamRS22k wrote:

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

ADD COMMENTlink modified 27 days ago • written 27 days ago by RamRS22k

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?

ADD REPLYlink written 25 days ago by evelyn30

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?

ADD REPLYlink written 24 days ago by RamRS22k

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

ADD REPLYlink written 24 days ago by evelyn30

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.

ADD REPLYlink modified 24 days ago • written 24 days ago by RamRS22k

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
ADD REPLYlink modified 21 days ago • written 21 days ago by evelyn30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 931 users visited in the last hour