Answers.Body: <p>My comment above mentions the slow IO for the Python code (it was iterating at ~ 3600 reads/s, which, for a FASTQ with 500M reads would take ~1.6 d to complete), so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using <code>zcat | awk</code> with a <code>switch-case</code> on the line number to stream the data.</p>
  <pre class="lang-sh prettyprint-override"><code># first pass to generate _R1 FASTQ
  zcat {input} | awk '{switch (NR % 4) {case 1: print <span class="math-container">$1; break; case 2: print substr($</span>1, 1, 76); break; case 3: print &quot;+&quot;; break; case 0: print substr($1, 1, 76)} }' | gzip -nc &gt; {output}_R1.fastq.gz
  
  # second pass to generate _R2 FASTQ
  zcat {input} | awk '{switch (NR % 4) {case 1: print <span class="math-container">$1; break; case 2: print substr($</span>1, 77, 76); break; case 3: print &quot;+&quot;; break; case 0: print substr($1, 77, 76)} }' | gzip -nc &gt; {output}_R2.fastq.gz
  </code></pre>
  <p>I played around with trying to pipe each half of the combined read into its own file, but it was difficult to set that up and generalize it for the multiple FASTQs that I wanted to fix. I didn't feel like over-engineering something, and this was straightforward enough. It worked much quicker than I expected.</p>
  <p>All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution.
  If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.</p>
  <p>But I will probably make use of <code>fasterq-dump</code> next time, and make sure I specify to split the read mates.</p>
  <p>I love how how extensive the default biopython fastq parser library is, but the slow IO for large sequencing files makes it unusable in a lot of cases, unfortunately. I wonder if there's a good way to use C bindings to make this faster, or whether you just need to use a different package/language.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-08T15:06:34.943Z
Answers.Text: I was thinking of two runs of awk, but I concluded that it can not be faster as you have to run it twice. Thank you for proving my intuition wrong.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-09T00:22:32.027Z
Answers.Text: I wish every bioinformatician knew this: *never* use the default biopython fastq parser on huge fastq files.
Answers.UserDisplayName: user172818
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-08T14:29:02.880Z
Answers.LastActivityDate: 2021-07-05T04:29:54.200Z
Answers.LastEditDate: 2021-07-05T04:29:54.200Z
Answers.UserDisplayName: James Hawley
Answers.Body: <p>As reads have equal length, so all 152-bp reads can be split into 76-bp read pairs accordingly. We can get Read1 by trimming the tail 76 bp, and get Read2 by trimming the heading 76 bp. No need to re-download or re-run fastq-dump.</p>
  <p>Specifically, use seqtk to do the job:</p>
  <pre><code>seqtk trimfq -e 76 merged_sra_file.fastq.gz | gzip -c &gt; Read1.fq.gz
  seqtk trimfq -b 76 merged_sra_file.fastq.gz | gzip -c &gt; Read2.fq.gz
  </code></pre>
  <p>I've tried and it runs very fast.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-04-23T07:26:58.997Z
Answers.LastActivityDate: 2024-04-23T07:27:57.667Z
Answers.LastEditDate: 2024-04-23T07:27:57.667Z
Answers.UserDisplayName: Zhuqiang Zhang
Body: <p>I downloaded a few FASTQs from the SRA that were paired-end 76 bp.
  When I look at the FASTQs, I get something like this:</p>
  
  <pre><code>@SRR7702334.1 1 length=152
  AAGATNTTTTTCAGTATCAATCAGANGNNANANNCNNCTGCCCCAAATCCTGCAAAGAGCATGCAACGTCCTTAGAAGACAACGTCTTGCTCTGTCTCACCCAGGCTGGAGCATAGGGGCATGATCGATCATGACAAGGAGTTTACATTTTG
  </code></pre>
  
  <p>I forgot to add the <code>--split-3</code> option to my <code>fastq-dump</code> command.
  Is there an easy and fast way, without re-downloading the files, to split the reads in half and write them to two output FASTQs?</p>
  
  <p>I tried the following using Biopython, and it'll work, but it's going to take a long time to finish.</p>
  
  <pre class="lang-py prettyprint-override"><code>from os.path import basename
  from Bio import SeqIO
  from Bio.Seq import Seq
  from Bio.SeqRecord import SeqRecord
  from Bio.Alphabet import IUPAC
  from gzip import open as gzopen
  import argparse
  from tqdm import tqdm
  
  
  def split_read(read, pos=152 // 2):
      """
      Split read in half to produce 2 reads
  
      Parameters
      ----------
      read_id : Bio.SeqRecord.SeqRecord
          Read to be split
      """
      return {
          # R1/R2 as keys
          "R" + str(i + 1): SeqRecord(
              # first and second halves of sequence string
              seq=Seq(str(read.seq[(pos * i):(pos * (i + 1))]), IUPAC.ambiguous_dna),
              id=read.id,
              description="",
              letter_annotations={"phred_quality": read.letter_annotations["phred_quality"][(pos * i):(pos * (i + 1))]}
          ) for i in range(2)
      }
  
  def split_fastq(fastq, prefix="split"):
      """
      Split reads from a FASTQ
  
      Parameters
      ----------
      fastq : str
          Path to input FASTQ file
      prefix : str
          Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)
      """
      # check for file format
      if fastq.endswith(".gz"):
          fq_handle = gzopen(fastq, "rt")
      else:
          fq_handle = open(fastq, "r")
      out_handles = {r: gzopen(prefix + "_" + r + ".fastq.gz", "wt") for r in ["R1", "R2"]}
      # load reads for random access
      records = SeqIO.parse(fq_handle, "fastq")
      for read in tqdm(records):
          # split read into each mate and fix header information
          mates = split_read(read)
          # write to each output file
          for m in mates:
              out_handles[m].write(mates[m].format("fastq"))
      fq_handle.close()
      for m in mates:
          out_handle[m].close()
  
  if __name__ == "__main__":
      # parse command line arguments
      PARSER = argparse.ArgumentParser(description="Split a FASTQ into 2 FASTQs by splitting each read in half")
      PARSER.add_argument(
          "fq",
          type=str,
          help="Input FASTQ to split"
      )
      PARSER.add_argument(
          "-p",
          "--prefix",
          type=str,
          help="Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)",
          default="split"
      )
      ARGS = PARSER.parse_args()
      # run splitting function
      split_fastq(ARGS.fq, ARGS.prefix)
  
  <span class="math-container">```</span>
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-05-07T22:23:54.300Z
Comments.Text: oh... the sequences are joined.. this is nasty....i supposed you have a lot of files?
Comments.UserDisplayName: StupidWolf
Tags: <fastq><sratoolkit><fastq-dump>
Title: How to split FASTQ reads without re-running `fastq-dump`?
ThreadText: # How to split FASTQ reads without re-running `fastq-dump`?
  
  (Asked by: username_0 on 2020-05-07)
  
  I downloaded a few FASTQs from the SRA that were paired-end 76 bp.
  When I look at the FASTQs, I get something like this:
  
  ```
  @SRR7702334.1 1 length=152
  AAGATNTTTTTCAGTATCAATCAGANGNNANANNCNNCTGCCCCAAATCCTGCAAAGAGCATGCAACGTCCTTAGAAGACAACGTCTTGCTCTGTCTCACCCAGGCTGGAGCATAGGGGCATGATCGATCATGACAAGGAGTTTACATTTTG
  
  ```
  
  I forgot to add the `--split-3` option to my `fastq-dump` command.
  Is there an easy and fast way, without re-downloading the files, to split the reads in half and write them to two output FASTQs?
  
  I tried the following using Biopython, and it'll work, but it's going to take a long time to finish.
  
  ```py
  from os.path import basename
  from Bio import SeqIO
  from Bio.Seq import Seq
  from Bio.SeqRecord import SeqRecord
  from Bio.Alphabet import IUPAC
  from gzip import open as gzopen
  import argparse
  from tqdm import tqdm
  
  def split_read(read, pos=152 // 2):
      """
      Split read in half to produce 2 reads
  
      Parameters
      ----------
      read_id : Bio.SeqRecord.SeqRecord
          Read to be split
      """
      return {
          # R1/R2 as keys
          "R" + str(i + 1): SeqRecord(
              # first and second halves of sequence string
              seq=Seq(str(read.seq[(pos * i):(pos * (i + 1))]), IUPAC.ambiguous_dna),
              id=read.id,
              description="",
              letter_annotations={"phred_quality": read.letter_annotations["phred_quality"][(pos * i):(pos * (i + 1))]}
          ) for i in range(2)
      }
  
  def split_fastq(fastq, prefix="split"):
      """
      Split reads from a FASTQ
  
      Parameters
      ----------
      fastq : str
          Path to input FASTQ file
      prefix : str
          Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)
      """
      # check for file format
      if fastq.endswith(".gz"):
          fq_handle = gzopen(fastq, "rt")
      else:
          fq_handle = open(fastq, "r")
      out_handles = {r: gzopen(prefix + "_" + r + ".fastq.gz", "wt") for r in ["R1", "R2"]}
      # load reads for random access
      records = SeqIO.parse(fq_handle, "fastq")
      for read in tqdm(records):
          # split read into each mate and fix header information
          mates = split_read(read)
          # write to each output file
          for m in mates:
              out_handles[m].write(mates[m].format("fastq"))
      fq_handle.close()
      for m in mates:
          out_handle[m].close()
  
  if __name__ == "__main__":
      # parse command line arguments
      PARSER = argparse.ArgumentParser(description="Split a FASTQ into 2 FASTQs by splitting each read in half")
      PARSER.add_argument(
          "fq",
          type=str,
          help="Input FASTQ to split"
      )
      PARSER.add_argument(
          "-p",
          "--prefix",
          type=str,
          help="Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)",
          default="split"
      )
      ARGS = PARSER.parse_args()
      # run splitting function
      split_fastq(ARGS.fq, ARGS.prefix)
  
  ```
  
  ```
  
  **Tags:** fastq, sratoolkit, fastq-dump
  
  > **Comments:**
  > **username_1**: oh... the sequences are joined.. this is nasty....i supposed you have a lot of files?  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2020-05-08. Score: 2 | Accepted answer
  
  My comment above mentions the slow IO for the Python code (it was iterating at ~ 3600 reads/s, which, for a FASTQ with 500M reads would take ~1.6 d to complete), so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using `zcat | awk` with a `switch-case` on the line number to stream the data.
  
  ```sh
  # first pass to generate _R1 FASTQ
  zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 1, 76); break; case 3: print "+"; break; case 0: print substr($1, 1, 76)} }' | gzip -nc > {output}_R1.fastq.gz
  
  # second pass to generate _R2 FASTQ
  zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 77, 76); break; case 3: print "+"; break; case 0: print substr($1, 77, 76)} }' | gzip -nc > {output}_R2.fastq.gz
  
  ```
  
  I played around with trying to pipe each half of the combined read into its own file, but it was difficult to set that up and generalize it for the multiple FASTQs that I wanted to fix. I didn't feel like over-engineering something, and this was straightforward enough. It worked much quicker than I expected.
  
  All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution.
  If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.
  
  But I will probably make use of `fasterq-dump` next time, and make sure I specify to split the read mates.
  
  I love how how extensive the default biopython fastq parser library is, but the slow IO for large sequencing files makes it unusable in a lot of cases, unfortunately. I wonder if there's a good way to use C bindings to make this faster, or whether you just need to use a different package/language.
  
  > **Comments:**
  > **username_2**: I was thinking of two runs of awk, but I concluded that it can not be faster as you have to run it twice. Thank you for proving my intuition wrong.  
  > **username_3**: I wish every bioinformatician knew this: *never* use the default biopython fastq parser on huge fastq files.  
  
  ---
  
  ### Answer by username_4 on 2024-04-23. Score: 2
  
  As reads have equal length, so all 152-bp reads can be split into 76-bp read pairs accordingly. We can get Read1 by trimming the tail 76 bp, and get Read2 by trimming the heading 76 bp. No need to re-download or re-run fastq-dump.
  
  Specifically, use seqtk to do the job:
  
  ```
  seqtk trimfq -e 76 merged_sra_file.fastq.gz | gzip -c > Read1.fq.gz
  seqtk trimfq -b 76 merged_sra_file.fastq.gz | gzip -c > Read2.fq.gz
  
  ```
  
  I've tried and it runs very fast.
  
  
  ---
  
  