text: user4222: 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.
  
  user77: 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.
  
  user3669: 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.
  
original_text: 4222: <p>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.</p>
  
  <p>I got some ideas from a code in Biostar:</p>
  
  <pre><code>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))
  </code></pre>
  
  <p>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.</p>
   77: <p>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.</p>
   3669: <p>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.</p>
  
  <p>Here is the code:</p>
  
  <pre><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))
  </code></pre>
  
  <p>Which returns:</p>
  
  <pre><code>(base) [dkennetz@nodecn203  test]$ python3.6 fastq_parser.py subset.fastq.gz
  250
  173
  </code></pre>
  
  <p>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:</p>
  
  <pre><code>def expression_finder(read_list, filter):
      ...
          if header.endswith(filter):
              ...
  
  subset_reads = expression_finder(read_list, b'expression')
  </code></pre>
  
  <p>which would also work. But include the <code>b''</code> in your string syntax because the first function still passes a bytes object.</p>
  