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

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?

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.

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


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!