<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>

CC BY-SA 4.0
2020-05-08T15:06:34.943Z
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.
Kamil S Jaron
CC BY-SA 4.0
2020-05-09T00:22:32.027Z
I wish every bioinformatician knew this: *never* use the default biopython fastq parser on huge fastq files.
user172818
CC BY-SA 4.0
2020-05-08T14:29:02.880Z
2021-07-05T04:29:54.200Z
2021-07-05T04:29:54.200Z
James Hawley
<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>

CC BY-SA 4.0
2024-04-23T07:26:58.997Z
2024-04-23T07:27:57.667Z
2024-04-23T07:27:57.667Z
Zhuqiang Zhang
<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>

CC BY-SA 4.0
2020-05-07T22:23:54.300Z
oh... the sequences are joined.. this is nasty....i supposed you have a lot of files?
StupidWolf
<fastq><sratoolkit><fastq-dump>
How to split FASTQ reads without re-running `fastq-dump`?
# 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.
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.
# 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.


---

