Getting data from fastq by generator

Question

I have a task in a training that I have to read and filter the 'good' reads of big fastq files. I downsampled, got the code working, saving in a python dictionary. But turns out the original files are huge and I rewrite the code to give a generator. It did work for the down-sampled sample. But I was wondering if it's a good idea to get out all the data and filtering in a dictionary. Does anybody here a better idea? I am asking because I am learning python on my own and I am not sure if it makes sense.

I got some ideas from a code in Biostar:

import sys
import gzip

filename = sys.argv[1]

def parsing_fastq_files(filename):

    with gzip.open(filename, "rb") as infile:
        count_lines = 0
            for line in infile:
            line = line.decode()
            if count_lines % 4 == 0:
                ids = line[1:].strip()
                yield ids
            if count_lines == 1:
                reads = line.rstrip()
                yield reads
        count_lines += 1

total_reads = parsing_fastq_files(filename)
print(next(total_reads))
print(next(total_reads))

I now need to figure out to get the data filtered by using 'if value.endswith('expression'):' but if I use a dict for example, but that's my doubt because of the number of keys and values.

Answer 1

Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.

Answer 2

I do not have an example of fastqs with 'expression' as the last string in the header, so instead I will filter based off of a specific barcode used for demultiplexing. I have a gzipped fastq of 250 reads in size. Based on Devon's answer, I have first generated a fastq parser which stores only (header, seq) as a list of tuples, and then looks inside the list of tuples for a specific index. Some of the indexes vary by 1 nucleotide so those will be left out. I can prove this by printing the length of the 2 lists.

Here is the code:

import sys
import gzip

fastq = sys.argv[1]

def parser(filename):
    """
    Generate a list of tuples (header, read)
    """
    fastq_parsed = []
    try:

        with gzip.open(filename) as fq:
            header = next(fq)
            read = next(fq)
            fastq_parsed.append((header[1:-1], read[:-1])) # we don't want new-line chars or @ sign in reads
            while True:
                next(fq) # skip + line
                next(fq) # skip qscore
                header = next(fq) # store next read header
                read = next(fq) # store next read seq
                fastq_parsed.append((header[1:-1], read[:-1]))
    except:
        StopIteration # when we run out of reads, stop
    return fastq_parsed

def expression_finder(read_list):
    """
    Find headers that end with expression assuming the list format generated above.
    """
    filtered_reads = []
    for header, read in read_list:
        if (header.endswith(b'ATCACGAT')): #gzip files in bytes format, so we pass bytes object.
            filtered_reads.append((header.decode("utf-8"), read.decode("utf-8"))) # list items decoded to utf-8
    return filtered_reads


read_list = parser(fastq)
print(len(read_list))
subset_reads = expression_finder(read_list)
print(len(subset_reads))

Which returns:

(base) [dkennetz@nodecn203  test]$ python3.6 fastq_parser.py subset.fastq.gz
250
173

We can see that of the 250 reads, 173 have the correct index. I have hardcoded the filter in the second function, but you could also pass it as an argument to the function, like:

def expression_finder(read_list, filter):
    ...
        if header.endswith(filter):
            ...

subset_reads = expression_finder(read_list, b'expression')

which would also work. But include the b'' in your string syntax because the first function still passes a bytes object.