2
$\begingroup$

I wrote some throwaway code for ingesting paired-end reads from a Fastq file into memory.

def ingest(stream1, stream2):
    while True:
        record1, record2 = list(), list()
        for _ in range(4):
            record1.append(next(stream1))
            record2.append(next(stream2))
        yield record1, record2

I had initially put the for loop in a try/Except block, anticipating I'd need to explicitly break out of the endless while loop. However, the code seems to work correctly even after I removed it. I even "corrupted" the final read in one of the paired files, and the ingest function performs reasonably (the last read pair is never yielded).

My question is...why does this code work, since I never explicitly break out of the loop?

Here's an example of how the code is invoked.

>>> with open('seqR1.fastq', 'r') as fh1, open('seqR2.fastq', 'r') as fh2:
...     counter = 0
...     for r1, r2 in interleave(fh1, fh2):
...         counter += 1
>>> counter
5000
>>> r1
['@read5000/1\n', 'AGGTCTGGGGACAGGTGTTTTACAGGCTCCTTTAGGGGTGGCAAGTCTCCCGGAAACATGATTAGTGAAGACACGTTAACAGCGAAACTGACAGAGAACCTCAACATTCTCTGGTAGGAAAAGCCAGGATCCTGGCAGAGGAAGCAACGATTCGAGCAGGGAACTGTGCAGGACAGAGAACAATGCGGCCCCAATTTTGATCCAGGGGTAGTTAGCTTAGCTAAGGAGTTCCCCAGTCCCTAAGTGGAGAGGGGATCCAGAACCCTACAATTTATTATTATATCTTAGCGGATATATCTAT\n', '+\n', 'CCCCCGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@GGGGGGGFAGGGCGGGGGFGGGGGGGGGCGGGGGDGGGGFGGGGGGGGGCGGGGGGGGGGGGGGGCGGGGGGGGGGGGGGGGG@G<GGGGGGGGGGGGFGGFGFF:DGGFFCGGGGCGGFEGGGGGGEGGGFB:GGGDG=FGGGG,GAGGGDGGFGGGCD3GFG+G?G*FGGG<GCGGGGC@CD<FFGF*:+G=FE<;GG,8C7GF=;9D=G<C5C66G6C>7)FD*6G7D)9;*C:)4D5+9>@G*=0):)*;1*504\n']
>>> r2
['@read5000/2\n', 'TGTGGGGTTCTGGATCCCCTCTCCACTTAGGGACTGGGGAACTCCTTAGCTAAGCTAATTACCCCTGGATCAAAATTGGGGCCGCATTGTTCTCTGTCCTGCACAGTTCCCTGCTCGAATCGTTGCTTCCTCTGCCAGGATCCTGGCTTTTCCTACCAGAGAATGTTGAGGTTCTCTGTCATTTTCGCTGCTAATGTGTCTTCACTAATCGTGTTTCCGGGAGACTTGCCACCCCTAAAGGAGCTTGTTGAACACTTGTCACCAGACCTAAGTCCTTAGGGATGAACCTGTGTCTCGGCAC\n', '+\n', 'CCCCBGC9FGCD-GGGFGGGEFEFF@CCGFGGGGFF+DG6,F69+FGG+FGF<FCC;GGEBFGGEDECFGGG6GDG7BFG9GCG9G9GDGGGG<GDFFFGEGC,<GC=GFG:,:A:EE+GD<:E,FGF=FCAGEFF++G,G,<EGGFGFFA>*FG9EGFCFCGC,B>43GFF@E@2>9ED@,568D=A;C+=5@E/F53D>,6F*9@+,E,C=*0+*@56<*9@:/F,1,:*)*14*0CC=*1)*0*,*)0;+)*))/**)**;0;//9)1+)10.()5/)070*)7*C**3/)(2*****\n']
$\endgroup$
4
$\begingroup$

Your code only runs prior to Python 3.7, now it raises RuntimeError: generator raised StopIteration. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via from __future__ import generator_stop.

Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:

it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

You should not write code relying on the bug

$\endgroup$
2
$\begingroup$

EDIT: Chris_Rands's answer is more authoritative and clarifies that this behavior is essentially a bug. (I had only tested my example on Python 3.6.8 and hadn't seen PEP 479.)

When next runs out of lines, it raises StopIteration, and apparently that stops the while loop similarly to calling break without you needing to catch the exception. I couldn't see a definitive source on that but you can see it in action if you catch and re-raise the exception:

def ingest2(stream1, stream2):
    while True:
        try:
            record1, record2 = list(), list()
            for _ in range(4):
                record1.append(next(stream1))
                record2.append(next(stream2))
            yield record1, record2
        except StopIteration:
            print("StopIteration raised")
            raise

You can also try to catch it outside the loop, but it doesn't get re-raised at that point. (If you try raising one yourself in a while loop, though, it doesn't catch it in the same way, so there must be some extra Python magic going on. I can't see any difference between the two exception objects so I'd be curious what's going on there.) Neat example. You could also make a pair of FASTQ parsers using Biopython's SeqIO.parse for more complex/permanent situations if needed.

$\endgroup$

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.