text: user298: 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?
  
  user339: 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!
  
  user77: 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.
  
original_text: 298: <p>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. </p>
  
  <p>So, given a specific variant like this:</p>
  
  <pre><code>chr22   425236  C T
  </code></pre>
  
  <p>How can I count the number of reads in my <code>sample.bam</code> file that support that variant and the number that don't on Linux?</p>
   339: <p><a href="https://github.com/dariober/ASCIIGenome" rel="nofollow noreferrer">ASCIIGenome</a> (I'm the author) has a command, <a href="http://asciigenome.readthedocs.io/en/latest/commandHelp.html#filtervariantreads" rel="nofollow noreferrer">filterVariantReads</a>, designed to inspect reads having a variant at a position or range. It would go along these lines:</p>
  
  <pre><code>ASCIIGenome -fa genome.fa aln.bam
  </code></pre>
  
  <p>Then go to the region of interest and use:</p>
  
  <pre><code>goto chr9:4917981-4918161
  filterVariantReads -r 4918011
  </code></pre>
  
  <p>From this:</p>
  
  <p><a href="https://i.stack.imgur.com/VigzQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/VigzQ.png" alt="enter image description here"></a></p>
  
  <p>You get:</p>
  
  <p><a href="https://i.stack.imgur.com/lxoW2.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/lxoW2.png" alt="enter image description here"></a></p>
  
  <p>You can also print to screen or save to file the underlying sam records with:</p>
  
  <pre><code>print &gt; reads.sam
  </code></pre>
  
  <p>For several regions, the whole process can be scripted for automation.</p>
  
  <p>Hope this helps!</p>
   77: <p>If you don't mind a bit of manual counting, then <code>samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam</code> 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:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>That'll give you a count of how many of each base were seen.</p>
  