4
$\begingroup$

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?

$\endgroup$

2 Answers 2

2
$\begingroup$

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 
$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.