Uniquely mapped PAIRS of reads

Question

[abascalfederico]

Hi,

For a PE-reads BAM file, I would like to keep only those pairs that map to a unique location. I know I can filter reads with multiple mappings with samtools (-bq 1), but this doesn't take into account paired information. I mean, if mate #1 maps uniquely to locus A, and mate #2 maps both to A and B, then we should be confident that #2 maps to A. Using samtools -bq 1 removes mate #2 in this case but I would like to keep it (because of the uniqueness of mate #1). Any idea how to solve this?

Thanks!

Federico

Answer 1

[Devon Ryan]

-bq 1 will not always remove multimappers, that they have a MAPQ of 0 is aligner dependent. As an example, a MAPQ of 1 can happen with multimappers produced by bowtie2.

Anyway, if you have a name-sorted file then you can write a little script with pysam to go through and filter things to whatever definition of "unique" you'd like to use.