# Getting data from fastq by generator

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:
count_lines += 1



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.

• You want lines where the header endswith expression correct? Mar 6 '19 at 16:23

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.

• This is more of a comment than an answer right Devon? Mar 6 '19 at 16:14
• @d_kennetz Not really, it answers what the proper scalable method is. One could include the actual code, but I have faith in the poster's ability given a prose description. Mar 6 '19 at 16:17

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):
"""
"""
fastq_parsed = []
try:

with gzip.open(filename) as fq:
while True:
next(fq) # skip + line
next(fq) # skip qscore
except:
StopIteration # when we run out of reads, stop
return fastq_parsed

"""
Find headers that end with expression assuming the list format generated above.
"""
if (header.endswith(b'ATCACGAT')): #gzip files in bytes format, so we pass bytes object.



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):
...

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