1
$\begingroup$

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)

```
$\endgroup$
5
  • 2
    $\begingroup$ My sense is that running fastq-dump takes about the same or even less time than trying out other methods with the possibility of running into errors and other issues. You could also use fasterq-dump which is faster but writes unzipped files and then use pigz or bgzip to compress your output. $\endgroup$ – Phoenix Mu May 7 '20 at 21:52
  • 1
    $\begingroup$ oh... the sequences are joined.. this is nasty....i supposed you have a lot of files? $\endgroup$ – StupidWolf May 7 '20 at 22:23
  • $\begingroup$ the time consuming part is the IO, not the splitting or processing.. so chances are, by the time you finish that, you are better off redoing fastq-dump $\endgroup$ – StupidWolf May 7 '20 at 22:29
  • $\begingroup$ I would also guess that this is IO problem and practically all working solutions will have about the same speed. $\endgroup$ – Kamil S Jaron May 8 '20 at 12:08
  • $\begingroup$ The IO for the python code was definitely the issue. It was iterating at ~ 3600 reads/s, which, for a FASTQ with 500M reads would take ~1.6 d to complete. $\endgroup$ – James Hawley May 8 '20 at 14:18
1
$\begingroup$

My comment above mentions the slow IO for the Python code, 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

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.

$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$ – Kamil S Jaron May 8 '20 at 15:06
  • $\begingroup$ 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 $\endgroup$ – James Hawley May 8 '20 at 15:17
  • 1
    $\begingroup$ I wish every bioinformatician knew this: never use the default biopython fastq parser on huge fastq files. $\endgroup$ – user172818 May 9 '20 at 0:22
  • $\begingroup$ I love how how extensive that 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 $\endgroup$ – James Hawley May 9 '20 at 14:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.