quickly getting paired-end fastq of all unmapped reads from BAM file

Question

[dp]

I've aligned my paired end reads with BWA and have a BAM file. Now I want to create 2 fastq files for all the paired-end reads that didn't correctly map.

Is the following command correct, if not can you suggest how to fix it? samtools fastq -1 filtered_reads1.fq -2 filtered_reads2.fq -n -G 0x2 -@ <n_threads>

(I've gone through some answers on the topic, and they've managed to confuse...)

[GenoMax]

Using reformat.sh from BBMap suite (make sure samtools is available in your $PATH):

reformat.sh in=your.bam out=stdout.sam unmappedonly=t | reformat.sh in=stdin.sam out1=R1.fq.gz out2=R2.fq.gz

[dp]

Thanks! This will keep only those with either pair unmapped? And it will keep both reads in the pair in that case? Should I be handling singletons somehow?

[GenoMax]

Do you need singletons? If you do then look at pairedonly= option. Turn it off or on based on your needs in first part of reformat.

[dp]

Thanks! I'm not sure? If the original read files are paired then I guess I won't have any singletons?

[GenoMax]

If unmapped reads are included in your BAM's then you should have no singletons (even if only one member of a pair is mapped).

[dp]

@genomax thanks for the answer! Another related question: If I aligned the reads to two references to create bam1 and bam2 and I want to keep all the reads that are unmapped in bam1 or mapped in bam2 (put another way: either they can be unmapped in bam1 or, if they are mapped in bam1, then they also need to be mapped in bam2), is there an easy way to do that?

[GenoMax]

Depending on how many you have you could use a different program from BBMap suite. filterbyname.sh (run the program without any options to see in-line help) to identify ones that are in both files (mapped or unmapped).

[dp]

Thanks for pointing it out!