GATK HP is missing high coverage variant sites

[goutham.atla]

I am trying to use GATK HP ( v3.5 ) to call SNPs in amplicon seq data of a small genome of around 600bp. As shown in the attachment, the variants are not called between location 50 and 60, despite high coverage across many samples. ( There are total 96 samples ).

https://www.dropbox.com/s/a4311lbti2sn9ze/seq.png?dl=0

The base qualities are above 30 and mapping quality is also 60 ( bwa mem ). I also did not remove duplicates ( not marked as well ) as its amplicon seq data.

The command I used was ( multisample SNP calling ) :

java -Xmx50g -jar GenomeAnalysisTK.jar -nct 2 -R<in.fasta> -T HaplotypeCaller -I merged_samples.bam -o gatk_out_raw_snps_indels.vcf --min_base_quality_score 30

Its the same case even with base quality 20.

Thanks in advance.

[vivek_]

Have you checked the default filtering parameters for SNP clusters? That looks like a lot of SNPs within a 20 base region, which also makes me question if the mapping quality of the reads spanning that region are really good.

[goutham.atla]

Yes. This is a viral amplicon data, which I would expect many variants. The mapping quality is 60, as assigned by the bwa-mem and base qualities are above 30. As I did not mark duplicates, I would expect it to consider all the reads. But is there any problem with very high depth ? Like more than 5k reads per base ? Should I randomly downsample few reads ?

[vivek_]

May be and there is a lot of discussion/confusion on haplotypecaller's downsampling technique as well. In general haplotypecaller has been very dicey in case of multisample calling for me even in human data.

[goutham.atla]

Yeah. My initial plan was to use freebayes, but it was failing due to memory issues. Erik garrison is looking in to the issue with freebayes.