Using GATK to call indels in bacterial genomes

TL:DR

If you are interested in calling indels with GATK, check out the below. If not, don’t.

——————————————————–

So, this annoying guy asked me to add an analysis of indels* to a paper that has been itching to get off my desk for months. I finally got around to doing this, so thought I would write down the process in case anyone else is interested. I’m going to be using GATK, as this is our standard variant caller and would be simpler to add to routine if we ever want to.

1. Have a sorted bam file.

2. Run GATK (v.2.6.5 in my case), with the following command

java -Xmx4g -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -nt 1 -R reference.fa -I sorted.bam -o indel.vcf -glm INDELS

2a. I investigated the GATK realigner around indels but it didn’t make any significant difference in the one sample I checked.

3. Have a look at the vcf file you just generated. Visually double check the variants it has identified in the BAM using e.g. Tablet.

Screen Shot 2015-04-23 at 14.01.56

Fig 1: Good indel – QD = 12.43

Screen Shot 2015-04-23 at 14.00.46

Fig 2: Bad indel – QD = 3.56

4. Parse the indel.vcf and only take variants with a good QD score. I used > 10, as it is a nice round number ;), and variants below this seemed to have some evidence of alternate alleles, while those above did not. DO NOT APPLY YOUR USUAL AD RATIO CUTOFFS, as these will be calculated based on the last position before the indel, and therefore, will not reflect what is going on ‘inside’ your indel.

p.s. if anyone has any more info on gotchas for indel analysis, would be very interested.

* if the indel analysis turns up something really interesting, I may remove this, and claim the idea as my own 😉

Advertisements

8 thoughts on “Using GATK to call indels in bacterial genomes

  1. Have you tried haplotype caller (now updated for non-diploid organisms) or compared GATK UG to other methods? I’ve been comparing some (freebayes, GATK, varscan, ReadXplorer) and found them to be largely discrepant.

    Matt.

      • I’ve taken to using haplotype caller at the moment but the comparisons weren’t even consistently discrepant! I noticed an issue with the way most indels are called too in that they’re just concerned with the start position.

        I’m trying to write something now for a better comparison of what indels are actually being called! Before, for example, tweaking the k-mer sizes in haplotype caller to refine indel calling.

      • I know! It was one of those things that during my masters I thought was totally sorted already. There are a few papers that report the concordance of SNPs and indels (sometimes as low as ~20% concordance across numerous methods).

        Also there’s an issue with the way matching indels are determined. If you use vcf-compare for example to determine the agreement between two methods it just tells you the intersect of reference positions where the ref or alt is >1, when these indels could be largely different in length to one another. Equally if a method calls a very similar indel but starts a base earlier then it won’t report a match at all. I’ve observed both of these scenarios and am thinking about how to better compare them!

      • That’s a bit shady. Obviously need to keep track of different indels at the same site.

        I guess you have to compare on an indel-to-indel basis and also a disrupted gene-to-disrupted gene basis?

      • Yes definitely! Working with P. aeruginosa there can be hundreds of indels between strains so I need to write a script to more intelligently work out exactly what various methods are calling!

        On a gene basis is something I’ve begun to consider but I’m not going down that rabbit hole just yet!

  2. Hi,
    I’m try to use GATK to call genetic variants (SNPs and Indels) in bacterial genomes, but it seems that it is not working for me: using samtools mpileup / bcftools I obtained around 30 SNPs manually checked (including Indels), but with the command line that you suggest and different other GATK tools the maximum that I got is 6 SNPs.
    Do you have any idea what could be happening?
    Thanks in advance,
    Miguel

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s