Another option, for querying any ad hoc range (chrN:X-Y):

$ echo -e "chrN\tX\tY" | bedops --element-of 1 <(bam2bed < reads.bam) - > overlapping_reads.bed

To query any given set of genomic ranges in variants.bed:

$ bedops --element-of 1 <(bam2bed < reads.bam) variants.bed > overlapping_reads.bed

The bam2bed tool will preserve all fields, which may be useful.

If variants are in VCF, then a second bash subshell to call vcf2bed can convert them to BED, e.g.:

$ bedops --element-of 1 <(bam2bed < reads.bam) <(vcf2bed < variants.vcf) > overlapping_reads.bed

I would use my tool sam2tsv to extract the name of the reads having (or not) the mutation.

For a specific variant, you can use bamql

For instance, if you want to extract all reads from a bam file :

• which support the alternate 'A' in snp : chr17:29827429 G>A

  bamql -f yourfile.bam 'chr(17) & nt_exact(29827429,A)' -o A.bam
  

• which support the reference 'G' in snp : chr17:29827429 G>A

  bamql -f yourfile.bam 'chr(17) & nt_exact(29827429,G)' -o B.bam
  

I guess you can loop over variant file and do this for each variant. You can also filter with mate information. See documentation Here

In the off chance that you really do need to have the full alignments in SAM format in separate files, then you can code something up in either python with pysam or C with HTSlib. In both cases you can use the pileup functions to just grab all alignments that overlap a given position (or set of them) and then iterate over the results. The pileup functions will tell you which base of each alignment overlaps a given position, so you can check what genotype it supports and treat it accordingly.

Edit: But really, if you think you want the full alignments in separate files then you should probably rethink what you're doing. The two answers posted before mine are vastly simpler and should suffice for most normal needs.

Thanks for the input, guys. I've realized that bam-readcount ultimately is the most appropriate tool for my purposes.