Answers.Body: <p>Not sure what is the reason you are realigning (older version of the genome? different mapper?) but probably you want to have read pairs even if only one in the pair was mapped to the region of interest.</p>
  <p>Assuming you do have a FASTQ file:</p>
  <pre><code>samtools view aligment.bam 21:31,831,502-31,896,094 \
  | cut -f 1 | sort | uniq &gt; selected_read_names.txt
  
  seqkit grep -f selected_read_names.txt seqs.r1.fq.gz -o result.r1.fq.gz
  seqkit grep -f selected_read_names.txt seqs.r2.fq.gz -o result.r2.fq.gz
  
  </code></pre>
  <p>Certainly not the fastest route but you will not drop imho perfectly usable reads.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-05-14T09:05:20.510Z
Answers.LastActivityDate: 2024-05-14T09:05:20.510Z
Answers.UserDisplayName: darked89
Answers.Body: <p>As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of <code>samtools collate</code>.</p>
  
  <pre><code>… | samtools fastq -F 0x900 -@ 48 \
      -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -
  </code></pre>
  
  <p>Your <code>samtools fastq</code> command is not doing anything to siphon off these singleton reads. You should add <code>-s /dev/null</code> to get rid of them.</p>
  
  <p>See the discussion in <a href="https://github.com/samtools/samtools/issues/874" rel="noreferrer">samtools/samtools#874</a> (especially the part starting <a href="https://github.com/samtools/samtools/issues/874#issuecomment-399403984" rel="noreferrer">here</a>). I have not looked at how the <a href="https://github.com/samtools/samtools/commit/f69e678bc16ebb93b6f5b16cd23a4ef7af0f7319" rel="noreferrer">recent changes from jkbonfield</a>, which will land in an upcoming samtools version (after 1.9), might affect this.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-09T14:14:46.997Z
Answers.Text: Oooh, that explains why I previously thought this worked! My script used to include `-s /dev/null` but then I removed it because I thought it had no effect. The documentation of these different outputs is somewhat confusing.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-09T14:07:58.503Z
Answers.LastActivityDate: 2019-07-09T14:07:58.503Z
Answers.UserDisplayName: John Marshall
Body: <p>I thought I had figured this one out. But apparently not.</p>
  <p>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.</p>
  <p>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):</p>
  <pre class="lang-sh prettyprint-override"><code>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 -
  </code></pre>
  <p>(Apart from the region filtering, this command is essentially the same that <a href="https://twitter.com/lh3lh3/status/1132756684789768202" rel="nofollow noreferrer">word of god</a> tells us to use.)</p>
  <p>The resulting files are not properly paired:</p>
  <pre><code>⟩⟩⟩ diff &lt;(gunzip -c reads_R1.fastq.gz | sed -n -e 1~4p -e 10q) \
  1 ⟩   &lt;(gunzip -c reads_R2.fastq.gz | sed -n -e 1~4p -e 10q)
  1,2d0
  &lt; @A00346:17:H52H7DSXX:3:2142:32316:24518
  &lt; @A00346:17:H52H7DSXX:3:1153:6705:26412
  3a2,3
  &gt; @A00346:17:H52H7DSXX:3:1173:9118:3129
  &gt; @A00346:17:H52H7DSXX:3:1109:29414:11757
  </code></pre>
  <p>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.</p>
  <p>And running <code>bwa mem</code> on the files gives me the error:</p>
  <pre><code>[mem_sam_pe] paired reads have different names: […]
  </code></pre>
  <p>How can I ensure that only properly paired reads are included in the FASTQ files?</p>
  <p>Some reads fall in the target region but their pairs fall outside it. I’d naively expect <code>samtools fastq</code> to do the right thing here, and actually output <em>pairs of reads</em>, or no read at all. Because resolving this efficiently separately is pretty hard.</p>
  
Tags: <sam><samtools><fastq>
Title: Convert BAM to properly paired FASTQ files
ThreadText: # Convert BAM to properly paired FASTQ files
  
  (Asked by: username_0 on 2019-07-09)
  
  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):
  
  ```sh
  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](https://twitter.com/lh3lh3/status/1132756684789768202) 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?
  
  Some reads fall in the target region but their pairs fall outside it. I’d naively expect `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.
  
  **Tags:** sam, samtools, fastq
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2019-07-09. Score: 8 | Accepted answer
  
  As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of `samtools collate`.
  
  ```
  … | samtools fastq -F 0x900 -@ 48 \
      -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -
  
  ```
  
  Your `samtools fastq` command is not doing anything to siphon off these singleton reads. You should add `-s /dev/null` to get rid of them.
  
  See the discussion in [samtools/samtools#874](https://github.com/samtools/samtools/issues/874) (especially the part starting [here](https://github.com/samtools/samtools/issues/874#issuecomment-399403984)). I have not looked at how the [recent changes from jkbonfield](https://github.com/samtools/samtools/commit/f69e678bc16ebb93b6f5b16cd23a4ef7af0f7319), which will land in an upcoming samtools version (after 1.9), might affect this.
  
  > **Comments:**
  > **username_0**: Oooh, that explains why I previously thought this worked! My script used to include `-s /dev/null` but then I removed it because I thought it had no effect. The documentation of these different outputs is somewhat confusing.  
  
  ---
  
  ### Answer by username_1 on 2024-05-14. Score: 0
  
  Not sure what is the reason you are realigning (older version of the genome? different mapper?) but probably you want to have read pairs even if only one in the pair was mapped to the region of interest.
  
  Assuming you do have a FASTQ file:
  
  ```
  samtools view aligment.bam 21:31,831,502-31,896,094 \
  | cut -f 1 | sort | uniq > selected_read_names.txt
  
  seqkit grep -f selected_read_names.txt seqs.r1.fq.gz -o result.r1.fq.gz
  seqkit grep -f selected_read_names.txt seqs.r2.fq.gz -o result.r2.fq.gz
  
  ```
  
  Certainly not the fastest route but you will not drop imho perfectly usable reads.
  
  
  ---
  
  