extract ALL mapped paired-end reads from BAM file

Question

[romain.derelle]

Hi,

I performed a mapping of paired-end reads using BWA mem.

I would like now to extract from the BAM file all paired-end reads for which at least one of the read is mapped (the mate read can be unmapped) and export both paired reads in the fastq format.

My commands:

samtools view -b -F 4 ali.bam > ali_filtered.bam
bedtools bamtofastq -i ali_filtered.bam -fq reads_1.fastq -fq2 reads_2.fastq

The problem I have is that I have a lot of error messages while running the bedtools command line:

*****WARNING: Query A00904:12:HGCJFDRXX:1:1101:15022:1235 is marked as paired, but its mate does not occur next to it in your BAM file. 
Skipping.

I'm now wondering whether my samtools command line is correct or if I'm missing all cases in which 1 of the 2 reads is unmapped.

Any help would be highly appreciated.

Thanks

[iraun]

I think -F 12 might be what you are looking for. You can use Explain samflags for the explanation :).

Answer 1

[lieven.sterck]

I would personally add the sorting on name in the samtools cmdline. This will put similarly named reads next to each other in the (sorted) BAM file.

samtools sort -n 

(the -n is important to sort on names otherwise the default is on left-most coordinate/position)

[Carlo Yague]

Yes, sorting on name is mandatory before using bamtofastq as explained in bedtools manual.