1
0
Entering edit mode
5.5 years ago
abascalfederico ★ 1.2k

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

paired end reads samtools bam file • 1.8k views
2
Entering edit mode
5.5 years ago

-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.