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
filterVariantReads -r 4918011

From this:

enter image description here

You get:

enter image description here

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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.