Problem with extracting reverse-reverse read pairs

Question

[Wang Cong]

Hi, I am trying to extract reverse-reverse read pairs from a bam file (all reads are paired and mapped).

I used the following command:

samtools view --expr 'flag.paired && !flag.unmap && !flag.munmap && flag.reverse && flag.mreverse' Pooled_valid.bam -o Pooled_valid.RR.bam

However, I found that all the reads are extracted (the input and output have the exactly same number of reads). May I ask if there is a problem with my command? Thanks

[Pierre Lindenbaum]

Please create a Minimal, Reproducible Example to enable others to exactly reproduce your issue on their machines.

Answer 1

[Carlo Yague]

The command doesnt' seem wrong to me but depending on the version of samtools, it might not work. It is arguably more robust to use the -f and -F options along with the numerical flags as it works with all samtools versions:

samtools view -f 49 -F 12 Pooled_valid.bam -o Pooled_valid.RR.bam

-f 49 means read paired (0x1), read reverse strand (0x10), mate reverse strand (0x20)

-F 12 mean NOT read unmapped (0x4) NOT mate unmapped (0x8)*

[jkbonfield]

Agreed the -f and -F options will be much faster, as the expression syntax is an interpretted string.

Note you can use -f and -F with strings, eg -f PAIRED,REVERSE,MREVERSE which may be more descriptive if you don't like the bit fields. Also the samtools flags command is a handy guide to the numeric and string names available.