How to split FASTQ reads without re-running `fastq-dump`?

Question

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)

```

Answer 1

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.

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