<p>Thanks Devon ;-) 
using samjdk : <a href="http://lindenb.github.io/jvarkit/SamJdk.html" rel="nofollow noreferrer">http://lindenb.github.io/jvarkit/SamJdk.html</a></p>

<p>a bed file 'input.bed'</p>

<p>and the following snippet:</p>

<pre><code>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;
}
</code></pre>

<p>usage:</p>

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

CC BY-SA 3.0
2018-02-15T21:55:06.853Z
2018-02-15T21:55:06.853Z
Pierre
<p>You're not going to get samtools to do that, though I wouldn't be surprised if something in <a href="https://github.com/lindenb/jvarkit" rel="nofollow noreferrer">jvarkit</a> can do it. Baring that, here's a bit of python.</p>

<pre class="lang-py prettyprint-override"><code>import pysam

bed = open("foo.bed")
bam = pysam.AlignmentFile("foo.bam")
of = pysam.AlignmentFile("filtered.bam", mode="wb", template=bam)
seen = set()  # memoize reads seen to not double print them
for line in bed:
    cols = line.strip().split("\t")
    for col in bam.pileup(cols[0], int(cols[1]), int(cols[2])):
        for b in col.pileups:
            if b.alignment.mapping_quality &lt; MIN_MAPQ:
                continue
            pos = b.query_position
            if not pos:
                continue
            if b.alignment.query_qualities[pos] &lt; MIN_PHRED:
                continue
            if b.alignment.query_name not in seen:
                seen.add(b.alignment.query_name)
                of.write(b.alignment)
of.close()
bam.close()
bed.close()
</code></pre>

<p>Or something along those lines, I haven't debugged that.</p>

CC BY-SA 3.0
2018-02-15T20:11:20.680Z
2018-02-15T20:11:20.680Z
Devon Ryan
<p>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. </p>

<p>Samtools seems to almost solve this problem, but not quite:</p>

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

<p><code>samtools view -L snp_bed -q 30 example.bam</code> 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 <code>-Q</code> option for samtools view, and the <code>-q</code> option is for mapping quality</p>

<p>Can this be achieved in samtools? Or is there another preexisting tool that can? </p>

<sam><samtools><bed>
Filter bam using SNP list in bed format with minimum mapping quality and base quality
# Filter bam using SNP list in bed format with minimum mapping quality and base quality

(Asked by: username_0 on 2018-02-15)

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?

**Tags:** sam, samtools, bed

---

## 2 answers

### Answer by username_1 on 2018-02-15. Score: 2

Thanks Devon ;-) 
using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html
a bed file 'input.bed'
and the following snippet:
private IntervalTreeMap 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();
    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->!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L->L.split("\t")).
        map(A->new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I->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< 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>=quals.length)  continue;
                if(quals[readpos]>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


---

### Answer by username_2 on 2018-02-15. Score: 1

You're not going to get samtools to do that, though I wouldn't be surprised if something in jvarkit can do it. Baring that, here's a bit of python.
import pysam

bed = open("foo.bed")
bam = pysam.AlignmentFile("foo.bam")
of = pysam.AlignmentFile("filtered.bam", mode="wb", template=bam)
seen = set()  # memoize reads seen to not double print them
for line in bed:
    cols = line.strip().split("\t")
    for col in bam.pileup(cols[0], int(cols[1]), int(cols[2])):
        for b in col.pileups:
            if b.alignment.mapping_quality < MIN_MAPQ:
                continue
            pos = b.query_position
            if not pos:
                continue
            if b.alignment.query_qualities[pos] < MIN_PHRED:
                continue
            if b.alignment.query_name not in seen:
                seen.add(b.alignment.query_name)
                of.write(b.alignment)
of.close()
bam.close()
bed.close()

Or something along those lines, I haven't debugged that.


---

