text: user240: 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.
  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)
  
  ```
  
  
  user240: 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.
  # 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.
  
original_text: 240: <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>
   240: <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>
  