I thought I had figured this one out. But apparently not.
I need to convert a BAM file of paired-end alignments to two FASTQ files of paired reads to realign them, with a twist: I only want reads that fall within a defined region.
Currently I am using the following command (the region is an example; it also fails on different regions; at any rate, this is on hs37d5 but presumably the issue would exist for other references too):
samtools view -b -u -@ 48 alignments.bam 21:31,831,502-31,896,094 \
| samtools collate -n 128 -u -O -@ 48 - /tmp/bam-to-fastq- \
| samtools fastq -F 0x900 -@ 48 \
-0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -
(Apart from the region filtering, this command is essentially the same that word of god tells us to use.)
The resulting files are not properly paired:
⟩⟩⟩ diff <(gunzip -c reads_R1.fastq.gz | sed -n -e 1~4p -e 10q) \
1 ⟩ <(gunzip -c reads_R2.fastq.gz | sed -n -e 1~4p -e 10q)
1,2d0
< @A00346:17:H52H7DSXX:3:2142:32316:24518
< @A00346:17:H52H7DSXX:3:1153:6705:26412
3a2,3
> @A00346:17:H52H7DSXX:3:1173:9118:3129
> @A00346:17:H52H7DSXX:3:1109:29414:11757
As far as I can see from performing a spot check, these are all reads which are marked as paired in the BAM file, but for which only one mate falls within the filtered region.
And running bwa mem
on the files gives me the error:
[mem_sam_pe] paired reads have different names: […]
How can I ensure that only properly paired reads are included in the FASTQ files?
samtools fastq
to do the right thing here, and actually output pairs of reads, or no read at all. Because resolving this efficiently separately is pretty hard. $\endgroup$samtools collate
should take care of this issue: If both reads of a pair are present, they will be adjacent; so no fancy data structure is necessary. In principle this should even be possible using a simple additional pass usingsed
but I’m currently failing at that. $\endgroup$