5
\$\begingroup\$

FASTQ is a notoriously bad format. This is because it uses the same @ character for the id line as it does for quality scores. Deciding what is a quality score and what is an id is a tricky endeavor with many pitfalls.

I'd like your opinion of my while 1 approach in the function read_fastq.

This works, but I'd like your tips and ideas on ways of improving it.

def read_fastq(fileH):
    """
    takes a fastq file as input
    yields idSeq, sequence and score
    for each fastq entry
    """

    #initialize the idSeq, sequence, score and index
    idSeq, sequence, score = None, None, None

    """
    main loop structure:

    An outer while loop will run until the file runs out of lines.

    If the line starts with @ and score exists, yield the id, sequence and score.
        this is where the yielding happens in our loop, it could be considered where each
        round of the loop ends

    The first id is recorded on the first iteration.

    If there is no sequence, begin an inner while loop, 
        read lines into sequence until we hit a + character, and break out of inner while loop

    if there is no score, begin an inner while loop where we will increment score until we have
         the same number of characters as sequence, this prevents interpretting a header as a score

    """

    while 1:

        line = fileH.readline()

        #break if we hit the end of the file
        if not line:
            break

        if line.startswith('@') and score:

            #before yielding the sequence, remove non-ascii characters
            sequence = legalChar(sequence)

            yield idSeq, sequence, score

            #reset to default values
            sequence = None
            score = None   
            idSeq = line.rstrip()


        elif not idSeq:
            #get our first idSeq
            idSeq = line.rstrip()
            continue

        elif not sequence:
            sequence = ""
            while not line.startswith('+'):
                sequence += line.rstrip().replace(' ', '')
                line = fileH.readline()
        elif not score:
            score = []
            #begin collecting our score, only collect as many chars as in our sequence
            while 1:
                score += line.rstrip().replace(' ', '')
                if len(score) >= len(sequence):
                    break
                else:
                    line = fileH.readline()

    #yield our final idSeq, sequence and score
    yield idSeq, sequence, score    
\$\endgroup\$
1
  • \$\begingroup\$ Thanks Jamal for adding the bioinformatics tag, I did not have enough karma for that! \$\endgroup\$ Oct 28, 2013 at 5:07

2 Answers 2

4
\$\begingroup\$

1. Bugs

All the bugs in this section are to do with handling of bad input; see §2.1 below for more about this issue.

  1. If there are no sequences in the input, your function wrongly yields an output sequence:

    >>> from StringIO import StringIO
    >>> list(read_fastq(StringIO('')))
    [(None, None, None)]
    
  2. If the input is wrongly formatted, your program can go into an infinite loop waiting for a line starting with a plus sign:

    >>> list(read_fastq(StringIO('@ID\nSEQUENCE\n')))
    ^C
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "cr32897.py", line 155, in read_fastq
        sequence += line.rstrip().replace(' ', '')
    KeyboardInterrupt
    
  3. If the input is wrongly formatted, your program can go into an infinite loop waiting for the end of the quality data:

    >>> list(read_fastq(StringIO('@ID\nSEQUENCE\n+ID\nSCORE')))
    ^C
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "cr32897.py", line 162, in read_fastq
        if len(score) >= len(sequence):
    KeyboardInterrupt
    

2. Other comments on your code

  1. You don't detect or report problems in the input. You just ignore or suppress them. This leads to outright bugs, as noted in §1 above, but it also means that you don't detect other problems with the input, like these:

    • missing sequence id;
    • junk lines before the first sequence or between sequences;
    • sequence id on + line fails to match sequence id on @ line;
    • length of quality data fails to match length of sequence data;
    • quality data contains characters outside the correct range.

    You shouldn't ignore or suppress errors in the input, because you need to know about these errors in order to discover problems with your processing pipeline. What if you are being fed the wrong input? What if the input is in the wrong format (for example, FASTA instead of FASTQ)? What if there is a bug in the code that generated the input? If you don't detect errors then you will end up generating bogus output and causing trouble for whatever process is next in your processing pipeline.

    When writing code that parses input data, the code for detecting and reporting errors typically constitutes the majority of the code. See §3 below for how I would go about writing this kind of parser.

  2. The first item of each output tuple (the sequence id) starts with the @ sign, but the @ sign is not part of the id.

  3. Your claim, "Deciding what is a quality score and what is an id is a tricky endeavor with many pitfalls" is a bit exaggerated. As far as I can see there is just one pitfall, namely that you might incorrectly stop reading the quality data when you encounter a line starting with @. If there are other pitfalls, you should probably mention them in a comment in the code so that we can check that they are all avoided.

  4. It's not clear to me what the H means in the variable name fileH. Does it stand for "handle"? If that's what you mean, you should write it out for clarity: file_handle. But really there's no reason not just to use the name file here.

  5. Your code contains two docstrings, but only the first is actually available to a user via the help function. The second docstring should be a comment instead — or better still, break it up into small comments that appear next to the relevant bits of code.

  6. The actual comments could be better. There is really no need to write comments like this:

    line = fileH.readline()
    
    #break if we hit the end of the file
    if not line:
        break
    

    You can expect Python programmers to know that file.readline returns an empty string to indicate end-of-file; or if not, they can find out by running help(file.readline).

    Some of the comments are wrong, for example here:

    #get our first idSeq
    idSeq = line.rstrip()
    

    This only gets the "first" idSeq if there was no previous line starting with an @.

  7. For an infinite loop, you write:

    while 1:
    

    but it's conventional to write:

    while True:
    
  8. In Python, a file object is also an iterator that yields the lines in the file, so instead of

    line = fileH.readline()
    
    #break if we hit the end of the file
    if not line:
        break
    

    you can write:

    line = next(fileH)
    

    This raises StopIteration when there are no more lines.

  9. You build up the sequence string by repeated addition:

    sequence = ""
    while not line.startswith('+'):
        sequence += line.rstrip().replace(' ', '')
        line = fileH.readline()
    

    In Python this is an anti-pattern that results in quadratic space and time consumption. The efficient way to build up a string from components is to put the components in a list and then call str.join:

    sequence_lines = []
    while not line.startswith('+'):
        sequence_lines.append(line.rstrip().replace(' ', ''))
        line = fileH.readline()
    sequence = ''.join(sequence_lines)
    
  10. You build up the score in the form of a list rather than a string:

    score = []
    while 1:
        score += line.rstrip().replace(' ', '')
        ...
    

    which means that your results come out like this:

    ('@ID', 'GATTTGGG', ['!', "'", "'", '*', '(', '(', '(', '('])
    

    Is that really what you want? I would have expected score to come out as a string, like this:

    ('ID', 'GATTTGGG', "!''*((((")
    

3. Revised code

Here's some revised code that fixes the problems above, and is much more robust in its detection and reporting of errors in the input.

class Error(Exception):
    pass

class Line(str):
    """A line of text with associated filename and line number."""
    def error(self, message):
        """Return an error relating to this line."""
        return Error("{0}({1}): {2}\n{3}"
                     .format(self.filename, self.lineno, message, self))

class Lines(object):
    """Lines(filename, iterator) wraps 'iterator' so that it yields Line
    objects, with line numbers starting from 1. 'filename' is used in
    error messages.

    """
    def __init__(self, filename, iterator):
        self.filename = filename
        self.lines = enumerate(iterator, start=1)

    def __iter__(self):
        return self

    def __next__(self):
        lineno, s = next(self.lines)
        line = Line(s)
        line.filename = self.filename
        line.lineno = lineno
        return line

    # For compatibility with Python 2.
    next = __next__

def read_fastq(filename, iterator):
    """Read FASTQ data from 'iterator' (which may be a file object or any
    other iterator that yields strings) and generate tuples (sequence
    name, sequence data, quality data). 'filename' is used in error
    messages.

    """
    # This implementation follows the FASTQ specification given here:
    # <http://nar.oxfordjournals.org/content/38/6/1767.full>
    import re
    at_seqname_re = re.compile(r'@(.+)$')
    sequence_re = re.compile(r'[!-*,-~]*$')
    plus_seqname_re = re.compile(r'\+(.*)$')
    quality_re = re.compile(r'[!-~]*$')

    lines = Lines(filename, iterator)
    for line in lines:
        # First line of block is @<seqname>.
        m = at_seqname_re.match(line)
        if not m:
            raise line.error("Expected @<seqname> but found:")
        seqname = m.group(1)
        try:
            # One or more lines of sequence data.
            sequence = []
            for line in lines:
                m = sequence_re.match(line)
                if not m:
                    break
                sequence.append(m.group(0))
            if not sequence:
                raise line.error("Expected <sequence> but found:")

            # The line following the sequence data consists of a plus
            # sign and an optional sequence name (if supplied, it must
            # match the sequence name from the start of the block).
            m = plus_seqname_re.match(line)
            if not m:
                raise line.error("Expected +[<seqname>] but found:")
            if m.group(1) not in ['', seqname]:
                raise line.error("Expected +{} but found:".format(seqname))

            # One or more lines of quality data, containing the same
            # number of characters as the sequence data.
            quality = []
            n = sum(map(len, sequence))
            while n > 0:
                line = next(lines)
                m = quality_re.match(line)
                if not m:
                    raise line.error("Expected <quality> but found:")
                n -= len(m.group(0))
                if n < 0:
                    raise line.error("<quality> is longer than <sequence>:")
                quality.append(m.group(0))

            yield seqname, ''.join(sequence), ''.join(quality)

        except StopIteration:
            raise line.error("End of input before sequence was complete:")
\$\endgroup\$
5
  • \$\begingroup\$ Wow, Thank you Gareth for the incredibly detailed response. I appreciate the thouroughness you went to to provide this. \$\endgroup\$ Oct 22, 2013 at 21:07
  • \$\begingroup\$ Gareth, before I was calling the function like this: for lineNum,(idSeq, sequence, score) in enumerate(read_fastq(options.inFile)): Now I am having trouble with the iterator. Can you show how I would, instantiate the iterator and call this function? \$\endgroup\$ Oct 23, 2013 at 19:10
  • \$\begingroup\$ Have you tried help(read_fastq)? \$\endgroup\$ Oct 24, 2013 at 8:21
  • \$\begingroup\$ Ok, so your code and response is awesome, I am really amazed at how well done this is, thank you. One thing I am currently working on is printing warnings (printing to stderr) for specific problems. For instance, if theres an empty sequence line, just print to stderr, yield an empty string and continue. How would you recommend to best implement that? Thanks! \$\endgroup\$ Oct 28, 2013 at 5:04
  • \$\begingroup\$ So with some small modifications to your code, I have a warning system in place for empty sequences or quality scores, but left in the errors for missing sequence name or + line. Here is my current implementation of the warning: def warn(self, message): """yield a warning relating to this line""" sys.stderr.write("Warning: {0}({1}): {2}\n{3}".format(self.filename, self.lineNum, message, self)) \$\endgroup\$ Oct 28, 2013 at 6:01
1
\$\begingroup\$

I have not tried to debug your code, but it seems pretty typical for a simple ad hoc parser - a state machine, implemented with if-statements and state variables (those you give None at the beginning).

I do not know if you are writing one-time script or a part of larger system, but if you want to make it in a more elegant way, then some python library could be used, which will allow you to describe the grammar of the FASTQ in an declarative way.

For example, pyparsing or lepl could be used.

The grammar of FASTQ is in fact very simple (judging by what is written on the Wikipedia page), so the grammar would be quite simple too. It may even be you will find official grammar for the FASTQ, which will make maintanace of the code very simple and bug-free too with more confidence.

See also https://stackoverflow.com/questions/223866/elegant-structured-text-file-parsing

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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