How can I count the number of reads that support a variant in a bam file?

Question

I am calling variants from a human sample using bwa mem to align the reads and gatk to call the variants. I'm trying to understand why a specific variant was not called in my sample. I have checked the bam alignments in a GUI viewer and I can see that there are reads supporting the missing variant. It looks like the issue is a low allelic balance, with far more reads supporting the reference than the alternate allele but I want to get the actual numbers.

So, given a specific variant like this:

chr22   425236  C T

How can I count the number of reads in my sample.bam file that support that variant and the number that don't on Linux?

Answer 1

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:

samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

That'll give you a count of how many of each base were seen.

Answer 2

ASCIIGenome (I'm the author) has a command, filterVariantReads, designed to inspect reads having a variant at a position or range. It would go along these lines:

ASCIIGenome -fa genome.fa aln.bam

Then go to the region of interest and use:

goto chr9:4917981-4918161
filterVariantReads -r 4918011

From this:

You get:

You can also print to screen or save to file the underlying sam records with:

print > reads.sam

For several regions, the whole process can be scripted for automation.

Hope this helps!

Answer 3

The variant record in the VCF produced by GATK should include information about depth -- both total read depth (DP) and allelic depth (AD), ie depth per allele.

Keep in mind however that the caller makes its decision based on more than just allelic depth. There are some additional statistics involved, including the quality of the individual bases. Also, if any filtering was applied, your variant may have been excluded due to other statistics, like for example if there is evidence of strand bias -- eg if all the reads supporting your variant are in the same orientation, that is a classic red flag suggesting it's an artifact.