There is a much simpler solution than all the ones showed here, which uses GATK's UnifiedGenotyper to do the base calling in all the wanted sites with the aid of the EMIT_ALL_SITES capability. the wanted sites should be previously defined in a bed file with the following format:

chrXXX  PositionOfInterest-1  PositionOfInterest

the final command would be something like this:

java -jar GenomeAnalysisTK.jar \
  -T UnifiedGenotyper \
  -R human_hg19.fa \
  -D dbsnp_135.hg19.vcf \
  -I input.bam \
  -o out.vcf \
  --intervals positions.bed \
  --output_mode EMIT_ALL_SITES

Entering the dbsnp file is completely optional, but it will output the rs code if a particular position has it, which is something very useful for us. I've just tested this command very recently because we are about to start using a 27 SNPs Sequenom multplex as quality control and sample identification for all the samples to be run in our sequencer, and it takes only a couple of seconds to run on an entire 10-15GB BAM exome file.

Ok here is an ugly (global variables...) C program that will output the information you need. You need to walk over the CIGAR string for each read. Please, not that I'm not sure about to handle the cigar operations involving a 'CLIP', feel free to correct my code.

Compilation:

gcc prog.c -I ./samtools-0.1.18 -L samtools-0.1.18 -lbam -lz

Execute:

$ ./a.out  toy.bam ref:7-10 
7   r001[1] T
8   r001[2] T
9   r001[3] A
9   r002[4] A
9   r003[6] A
10  r001[4] G
10  r002[5] G
10  r003[7] N
11  r001[5] A
11  r002[6] A
11  r003[8] N
12  r001[6] T
12  r002[7] T
12  r003[9] N
13  r001[7] A
13  r002[8] A
13  r003[10]    N
14  r001[8] A
14  r002[9] A
14  r003[11]    N
15  r001[13]    G
15  r002[14]    T
16  r001[14]    A
16  r002[15]    A
17  r001[15]    T
17  r002[16]    A
18  r001[16]    A
18  r002[17]    A
20  r001[17]    C
21  r001[18]    T
22  r001[19]    G

Source:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "bam.h"
#include "sam.h"

static bamFile in=0;
static bam_header_t* header=0;
static bam_index_t *idx=0;

// callback for bam_fetch()
static int fetch_func(const bam1_t *b, void *data)
    {
    bam_plbuf_t *buf = (bam_plbuf_t*)data;
    bam_plbuf_push(b, buf);
    return 0;
    }

static int  scan_all_genome_func(uint32_t tid, uint32_t pos, int depth, const bam_pileup1_t *pl, void *data)
        {    
        int i;   
        for( i=0;i<depth ;++i)
        {
        const bam_pileup1_t *p=&pl[i];
        bam1_t* aln=p->b;
        uint32_t *cigar = bam1_cigar(aln);
        //const uint8_t* seq=bam1_seq(aln);
        int read_seq_index=0;
        uint32_t genomic_pos=aln->core.pos;
        uint32_t  cigaridx;
        for (cigaridx = 0; cigaridx < aln->core.n_cigar; ++cigaridx)
            {
            int k;
            int op = cigar[cigaridx] & BAM_CIGAR_MASK;
            int n_op= cigar[cigaridx] >> BAM_CIGAR_SHIFT;
            for(k=0;k< n_op;++k)
            {
            switch(op)
                {
                case BAM_CMATCH:
                {
                if(genomic_pos==pos)
                    {
                    char base=bam_nt16_rev_table[bam1_seqi(bam1_seq(aln),read_seq_index)];
                    fprintf(stdout,"%d\t%s[%d]\t%c\n",
                        genomic_pos+1,
                        bam1_qname(aln),
                        read_seq_index+1,
                        base

                        );
                    }
                ++read_seq_index;
                ++genomic_pos;
                break;
                }
                //deletion from the reference
                case BAM_CDEL:
                {
                ++genomic_pos;
                break;
                }
                //insertion to the reference
                case BAM_CHARD_CLIP://?
                case BAM_CPAD://?
                case BAM_CSOFT_CLIP:
                case BAM_CINS:
                {
                ++read_seq_index;
                break;
                }
                default:
                {
                fprintf(stderr,"#cigar::op not handled: %d\n",op) ;
                break;
                }
                }
            }
            }
        }
        return 0;
        }

int main(int argc,char** argv)
    {
    int ref,beg,end;
    bam_plbuf_t* buf;
    in= bam_open(argv[1],"rb");
    if(in==NULL) exit(-1); 

    header= bam_header_read(in);
    if(header==NULL) exit(-1);
    bam_init_header_hash(header);
    idx = bam_index_load(argv[1]);
    if(idx==NULL) exit(-1);

    bam_parse_region(header, argv[2], &ref, &beg, &end); // parse the region
    if (ref < 0) exit(-1);
    buf= bam_plbuf_init( scan_all_genome_func,NULL);
    if(buf==NULL) exit(-1);
    bam_fetch(in, idx, ref, beg, end, buf, fetch_func);
    bam_plbuf_push(0, buf); // finalize pileup
    bam_plbuf_destroy(buf);
    bam_header_destroy(header);
    bam_index_destroy(idx);
    bam_close(in);
    return EXIT_SUCCESS;
    }

I know this post is really old, but in case someone is still looking, there is a tool called bam-readcount which should do exactly what you want.

Simply run bam-readcount your_bam -l site_list where site_list contains the regions of interest and has the format chr:start-stop (1 based and tab seperated)

As with almost every life problem, this is one that can [surprisingly] be handled by mpileup:

cat rsBED/rs2857845.bed 
chr4    3026385 3026386

bcftools mpileup --regions-file rsBED/rs2857845.bed \
  --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
  input.bam \
  --fasta-ref refdata-GRCh38-2.1.0/fasta/genome.fa

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  CH01848
chr4    3026386 .   A   T,<*>   0   .   DP=20;ADF=6,6,0;ADR=2,4,0;AD=8,10,0;I16=6,2,6,4,295,11299,365,13851,480,28800,600,36000,183,4439,203,4917;QS=0.44697,0.55303,0;VDB=0.558487;SGB=-0.670168;RPB=0.906029;MQB=1;MQSB=1;BQB=0.991158;MQ0F=0 PL:DP:SP:ADF:ADR:AD 224,0,174,248,204,255:18:2:6,6,0:2,4,0:8,10,0

Let's pipe it into bcftools query to tidy it up nicely:

bcftools mpileup --regions-file rsBED/rs2857845.bed \
  --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
  input.bam \
  --fasta-ref refdata-GRCh38-2.1.0/fasta/genome.fa |\
  bcftools query --format '%CHROM\t%POS\t%REF\t%ALT[\t%AD]\n'

chr4    3026386 A   T,<*>   8,10,0

We don't even need the [really] biased Buffalo New York reference genome - just add --no-reference to bcftools mpileup

bcftools mpileup --regions-file rsBED/rs2857845.bed \
  --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
  input.bam \
  --no-reference |\
  bcftools query --format '%CHROM\t%POS\t%ALT[\t%AD]\n'

chr4    3026386 T,A,<*> 0,11,8,0

Be also weary of filter flags, such as:

-q, --min-MQ INT        skip alignments with mapQ smaller than INT [0]
-Q, --min-BQ INT        skip bases with baseQ/BAQ smaller than INT [13]

Kevin

I know this post is really old, but in case someone is still looking, there is a tool called bam-readcount which should do exactly what you want.

Simply run bam-readcount your_bam -l site_list where site_list contains the regions of interest and has the format chr:start-stop (1 based and tab seperated)