I am trying to retrieve uniquely mapped reads from a BAM file (paired-end) mapped with BWA MEM, with the scanBam function from Rsamtools:
# setting isSecondaryAlignment=FALSE param <- ScanBamParam(flag=scanBamFlag(isUnmappedQuery=FALSE, isSecondaryAlignment=FALSE), what=c("qname", "pos", "qwidth", "rname")) # scanning BAM readsfile <- "test.bam" aln <- scanBam(readsfile, param=param)[]
The problem is that there remain some secondary alignments:
# testing whether some read IDs are present more than twice (because it is paired-end data, we expect them to be present exactly 2 times if only uniquely mapped reads remain): any(table(aln$qname) > 2) # TRUE
I read that I could also use the "mapqFilter" in ScanBamParam, but I also read that it might not be the most reliable way to proceed. Also, I am sure which MAPQ filter to use with BWA MEM...
The following is working with samtools (outside of R), which I found in this post:
samtools view -F 2048 -bo filtered.bam original.bam
Is there a way to do the same think with Rsamtools, i.e. using the appropriate SAM flag (2048 corresponds to any secondary alignment).
Either that, or someone can advice a robust alternative solution to this problem (or point me to something I might be doing wrong)?