<p>Hi everyone,</p>

<p>I have a number of Human Heart RNASeq samples. I have generated a mpileup file on each of them using a bed file containing a list of positions for which I want the mpileup to be made. I have included the -B option to turn off the BAQ computation in mpileup.</p>

<p>Single sample mpileup:</p>

<pre><code>samtools mpileup -uf hg19.fa -B -l snps_chr17.bed sample1.bam &gt; sample1.mpileup
</code></pre>

<p>Multiple sample mpileup (all samples):</p>

<pre><code>samtools mpileup -uf hg19.fa -B -l snps_chr17.bed -b bam_files_list.bam &gt; all_samples.mpileup
</code></pre>

<p>Now, I want to output the <strong><em>frequency of nucleotides A, T, G and C at each of those positions</em></strong>. There is a perl script <a href='https://github.com/riverlee/pileup2base/blob/master/pileup2baseindel.pl'>pileup2baseindel</a> which generates <em>exactly</em> the output that I want, only that it is incorrect (inconsistent with what is seen in IGV). This program takes a single/multi-sample <a href='https://github.com/riverlee/pileup2base/blob/master/test.mpileup'>mpileup file</a> and produces the following <a href='https://github.com/riverlee/pileup2base/blob/master/sample1.txt'>output</a>. But like I said, the output is incorrect.</p>

<p>I have searched a lot but can't seem to find any program that does the same. Does anyone know any tool/script that would give me such an output with/without taking the mpileup file as an input. Suggestions on alternatives/options are welcomed.</p>

<p>[UPDATE] There was a bug in the perl script and I have updated the code, it works fine now. I have contacted the author and as soon as the author responds, I will send it to him so that he can update it. I would still like to get suggestions and try out other programs. </p>

<p>[UPDATE] Thanks everyone who helped me with this problem. Each of your suggestions worked but I accepted <a href='/u/117/'>Chris Miller</a>'s answer recommending the use of <a href='https://github.com/genome/bam-readcount'>bam-readcount</a> because of its easy of use and the elaborate output that it generates. It directly takes, as input, a <a href='http://samtools.sourceforge.net/SAM1.pdf'>bam</a> file as well as a list of positions in bed format for which you want the nucleotide frequency. So you don't need to create an intermediate mpileup file - as in <a href='https://github.com/riverlee/pileup2base'>pileup2baseindel</a> or run the program against positions that you are not interested in - as in <a href='https://github.com/alimanfoo/pysamstats'>pysamstats</a>.</p>

<p>Thanks!</p>

Another option is just to run [bam-readcount][1], giving it the bams and specified positions. It's more straightforward than creating an mpileup and parsing that ugly string.

 [1]: https://github.com/genome/bam-readcount
