Filter bam using SNP list in bed format with minimum mapping quality and base quality I have a bam file and a bed file that defines a list of SNPs. I would like to filter the bam file to contain only those reads with a minimum mapping quality that overlap at least one SNP with a minimum base quality. 

Samtools seems to almost solve this problem, but not quite:

samtools depth -b snp_list -q 20 -Q 30 example.bam allows me to count the number of alignments meeting my desired criteria, but not to produce a bam with these alignments.

samtools view -L snp_bed -q 30 example.bam allows me to filter the alignments overlapping any SNP with minimum mapping quality, but not with a minimum base quality at the SNP sites. There is no -Q option for samtools view, and the -q option is for mapping quality

Can this be achieved in samtools? Or is there another preexisting tool that can?
Thanks Devon ;-) 
using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

a bed file 'input.bed'

and the following snippet:

private IntervalTreeMap&lt;Boolean&gt; treeMap=null;
private final int BASE_QUALITY= 20;
@Override
public Object apply(final SAMRecord record) {
if(record.getReadUnmappedFlag() || record.getCigar()==null) return null;
if(treeMap==null) {
    treeMap = new  IntervalTreeMap&lt;Boolean&gt;();
    final String bedStr;
    try {
        bedStr = IOUtil.slurp(new java.io.File("input.bed"));
        }
    catch(java.io.IOException err)
        {
        throw new RuntimeIOException(err);
        }
    Arrays.stream(bedStr.split("[\n]")).
        filter(L-&gt;!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L-&gt;L.split("\t")).
        map(A-&gt;new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I-&gt;treeMap.put(I,true));
    }


byte quals[]= record.getBaseQualities();
int refpos= record.getStart();
for(CigarElement ce:record.getCigar()) {
    CigarOperator op= ce.getOperator();
    if(op.consumesReferenceBases())
        {
        if(op.consumesReadBases()) {
            for(int i=0;i&lt; ce.getLength();i++)
                {
                Interval t = new Interval(record.getContig(),refpos+i,refpos+i);
                if(!treeMap.containsOverlapping(t)) continue;
                int readpos = record.getReadPositionAtReferencePosition(refpos+i);
                if(readpos==0 ) continue;
                readpos--;
                if(readpos&gt;=quals.length)  continue;
                if(quals[readpos]&gt;BASE_QUALITY) return true;
                }
            }

        refpos+=ce.getLength();
        }

    }
return null;
}


usage:

samtools view -bu -L input.bed input.bam |\
   java -jar dist/samjdk.jar --body -f snippet.code