Filter bam using SNP list in bed format with minimum mapping quality and base quality

Question

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?

Answer 1

Thanks Devon ;-) using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

a bed file 'input.bed'

and the following snippet:

private IntervalTreeMap<Boolean> 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<Boolean>();
    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 2

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.