biopython/biopython — pull request #219 by mmokrejs (2013-08-24 12:38:04+00:00)
Title: Do not break on SFF file created by broken sff-dump-2.1.[7-9] which do not zero-out some fields.

The sra-tools email archives are not publicly accessible so I rather attach here the whole email thread, explaining what is the issue with SFF files. The files are usable and therefore we should drop the ValueError from biopython. Wish this makes it into final 1.62.

Date: Sat, 17 Mar 2012 20:01:20 -0500
To: Martin Mokrejs
From: SRA Toolkit Development Group
Subject: Re: [Sra-tools] SffIO
CC: Peter Cock ,
NLM/NCBI List sra-tools
Content-Type: text/plain; charset="us-ascii"; format=flowed

Hi Martin,

The 2.1.6 sra-dump should work fine for you, and the 2.1.10 release
is scheduled to appear in a few days (week of the 19th), so it's
likely either approach will work fine.

The padding issue occurs on nearly every spot. Technically, it occurs
on all but the longest spot emitted so far due to reuse of an output
buffer. So if spot 1 had 200 bases, while spots 2 and 3 had
non-decreasing lengths, the padding will be zeroed. As soon as any
spot is emitted that is shorter than the longest emitted so far, the
padding will contain left-over data in the buffer.

It's worth repeating that the significant data are 100% correct. It's
fine to insist on zeroing out the padding, and we are issuing an
immediate fix because we also feel it's important, but it cannot be
confused with containing corrupt data - only a format violation.

So in summary, you should expect nearly but not all rows to have non
zero padding. sff-dump.2.1.6 will work fine for you until you get the
next release.

SRA files are in a different category than simple flat files like SFF
and do not have a "format" per se. They are described by a schema,
and are container units for a massive distributed database storage
technology used to implement the SRA. They are indexed and
compressed, and the format is directly executable. The only purpose
of exporting them to formats such as SFF is to use tools that do not
yet read SRA directly. If you were to look at an SRA file using, say
vdb-dump, you would be able to see a number of columns in the table.
We record any data present in the SFF related to clips -
{adapter/quality}:{left/right}, but typically only find
clip-quality-right. But we store the bases, phred score qualities,
signal and position mapping in every run originating in SFF format.
We also record a read description and in some cases the linker used,
which is omitted in SFF.

All of this is to say that SRA files do not have any singular format,
but rather a description of storage given in schema. They can be
converted into many other formats if desired, or read directly. Since
they are always compressed and remain indexed and directly
addressable, SRA storage is the most compact way of storing
biological data, as well as being fully spec'd with an easy to use
API.

Regarding the tool you're interested in, it would be trivial for you
to write such a tool that loads SFF into SRA and then uses the SRA
API to read trimmed data and rewrite it in any format you wish. Let
us know if you'd like help.

Thanks,

SRA Toolkit Team

At 6:49 AM -0400 3/17/12, Martin Mokrejs wrote:
>Hi "SRA",
> would you mind investigating how often/when the error occurs? The sff
>file written by biopython (until it exited) contained 1-2 entries as far as I
>had a quick look. Is there something you can learn from that? Not me probably.
>;)
> I also vote for more sanity checks.
> I use biopython to update the clip_adapter_right fields so I need python
>to be able to parse the input SFF and write out a corrected file. Do you
>recommend me to use sff-dump-2.1.6 or rather wait for 2.1.10 binaries to come
>out?
> I somewhat did not bother to learn how is the SFF data packed into
>.sra files
>and hence do not understand whether the SFF data in them are those original
>SFF files produced by Roche just packed into or whether you re-wrote them into
>.sra structures. Would you enlighten me please? ;)
>
> BTW, would it be much work for you to write something like sfffile
>from Roche
>what could read SFF file and output a new one? I am specifically interested in
>an option to write new SFF files with trimmed reads based on the
>clip_adapter_right
>field (in the input file). The field would be left empty of course
>in the output
>SFF file. One could save quite some space using this. No, I am not
>talking about
>clip_qual fields nor do I need a tool which produces qual+adap
>clipped SFF files
>although that would be also handy, and follow specs fully. ;)
>
>Best,
>Martin
>
>SRA Toolkit Development Group wrote:
>> Hi Martin & Peter:
>>
>> It looks like there may be two bugs here. The first bug is in the
>>recent releases of sff-dump - it produces correct data, but does
>>not zero-out padding of the rows as it did previously. This not
>>only does not conform with SFF specs, but produces the result that
>>there can be no reliable checksum of the file. There was no leakage
>>of information, however, and no loss or corruption of actual data
>>which is why sffinfo worked without problem.
>>
>> The second bug would appear to be attaching any significance to
>>data values within the padding areas of SFF on the part of SffIO.
>>The format specifies they should be zeroed, but they contain no
>>data and should not be read or given any significance.
>>
>> You may continue to use the SFFs you have with total confidence in
>>their content, and the upcoming release of 2.1.10 will contain a
>>fix that zeros out padding.
>>
>> Thank you very much,
>>
>> SRA Toolkit Team
>>
>>
>> On Mar 16, 2012, at 7:56 AM, Martin Mokrejs wrote:
>>
>>> Dear SRA team,
>>> I think I found a bug in sff-dump-2.1.9 at least. I have some
>>>sff files which cannot be parsed by biopython. It seems this is a
>>>broken SFF file written by sff-dump-2.1.9 (contrary to 2.1.6).
>>>Please read further below.
>>>
>>> Peter Cock wrote:
>>>> On Thu, Mar 15, 2012 at 11:50 PM, Martin Mokrejs
>>>> wrote:
>>>>> Hi Peter,
>>>>> have you seens SffIO crashing here?
>>>>>
>>>>> _wrote1 = SeqIO.write(_fixed_sff_records, prefix +
>>>>>".with_annotated_adapters.sff", "sff")
>>>>> File
>>>>>"/usr/lib64/python2.7/site-packages/Bio/SeqIO/__init__.py", line
>>>>>421, in write
>>>>> count = writer_class(handle).write_file(sequences)
>>>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>>>line 887, in write_file
>>>>> for record in records:
>>>>> File "myprog.py", line 3470, in fix_sff_records
>>>>> for _record in SeqIO.parse(prefix + ".sff", "sff"):
>>>>> File
>>>>>"/usr/lib64/python2.7/site-packages/Bio/SeqIO/__init__.py", line
>>>>>541, in parse
>>>>> for r in i:
>>>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>>>line 802, in SffIterator
>>>>> trim)
>>>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>>>line 573, in _sff_read_seq_record
>>>>> % padding)
>>>>> ValueError: Post quality 6 byte padding region contained data
>>>>>
>>>>> It happens for ERR016587.sff file extracted from ERR016587.sra. ;)
>>>>>
>>>>>
>>>>> For ERR016602.sff I get:
>>>>>
>>>>> ...
>>>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>>>line 573, in _sff_read_seq_record
>>>>> % padding)
>>>>> ValueError: Post quality 7 byte padding region contained data
>>>>>
>>>>
>>>> That's new. If my code is right, the SFF file is corrupt. How big is it?
>>>
>>> I wen to check the files and re-created the SFF and FASTA files
>>>from the "broken" SFF and also, re-fetched .sra from SRA FTP
>>>server. Then, I extracted the SFF file again from the new SRA file
>>>using both, sff-dump 2.1.6 and 2.1.9.
>>>
>>>
>>> $ md5sum ERR016587*.fasta ERR016587*.sff ERR016587*.sra
>>> f7807c1ea3d7b681d9ad2d7b70d8323e ERR016587-from-bad-sff.fasta
>>> f7807c1ea3d7b681d9ad2d7b70d8323e ERR016587-from-re-extracted-sra.fasta
>>> 57223d8ae31eed29fe9e8f60cc71a02a ERR016587-bad.sff
>>> 72ebe1c4b44ef54731b40e3915401eec ERR016587-re-extracted-2.1.6.sff
>>> cd4da7b5f4e72c3870eb812298410b36 ERR016587-maybe-bad.sra
>>> cd4da7b5f4e72c3870eb812298410b36 ERR016587.sra
>>> $ ls -la ERR016587*.fasta ERR016587*.sff ERR016587*.sra
>>> -rw------- 1 mmokrejs users 6736514448 Mar 1 02:38 ERR016587-bad.sff
>>> -rw-r--r-- 1 mmokrejs users 1188503102 Mar 12 13:11
>>>ERR016587-from-bad-sff.fasta
>>> -rw------- 1 mmokrejs users 1188503102 Mar 16 11:50
>>>ERR016587-from-re-extracted-sra.fasta
>>> -r-------- 1 mmokrejs users 2502489652 Jan 19 22:09 ERR016587-maybe-bad.sra
>>> -rw------- 1 mmokrejs users 6736514448 Mar 16 10:32
>>>ERR016587-re-extracted-2.1.6.sff
>>> -rw------- 1 mmokrejs users 2502489652 Mar 16 11:32 ERR016587.sra
>>> $ ~/sratoolkit.2.1.9-centos_linux64/sff-dump ERR016587.sra
>>> Written 1237359 spots for ERR016587.sra
> >> Written 1237359 spots total
>>> $ mv ERR016587.sff ERR016587-from-new-sra-with-2.1.9.sff
>>> $ md5sum ERR016587-from-new-sra-with-2.1.9.sff
>>> 57223d8ae31eed29fe9e8f60cc71a02a ERR016587-from-new-sra-with-2.1.9.sff
>>> $ /usr/bin/sff-dump.2.1.6 ERR016587.sra
>>> Written 1237359 spots for ERR016587.sra
>>> Written 1237359 spots total
>>> $ md5sum ERR016587.sff
>>> 72ebe1c4b44ef54731b40e3915401eec ERR016587.sff
>>> $
>>>
>>>
>>> The checksum of the newly re-extracted SFF file using
>>>sff-dump-2.1.9 is same like the of the broken ERR016587-bad.sff
>>>file. Therefore, I conclude this is a bug with sff-dump-2.1.9.
>>>Interestingly, Roche's sffinfo can extract the sequences without
>>>any error. Probably it is broken in some other data structure.
>>>
>>> Contrary to that, sff-dump-2.1.6 extracts a SFF file with
>>>different checksum. Alhough I did not re-test that using biopython
>>>I bet sff-dump in between 2.1.6 to 2.1.9 got broken. Unless you
>>>say that sff-dump can produce SFF files with different checksums
>>>for some reasons, I believe the file I parsed using biopython was
>>>created by sff-dump-2.1.9.
>>>
>>> Anyways, here is the proof.
>>>
>>> $ python
>>> Python 2.7.2 (default, Feb 29 2012, 04:38:36)
>>> [GCC 4.5.3] on linux2
>>> Type "help", "copyright", "credits" or "license" for more information.
>>>>>> from Bio import SeqIO
>>>>>> def get_generator(filename):
>>> ... for _record in SeqIO.parse(filename, "sff"):
>>> ... yield _record
>>> ...
>>>>>> _sff_records = get_generator("ERR016587-bad.sff")
>>>>>> _wrote1 = SeqIO.write(_sff_records,
>>>>>>"ERR016587-from-biopython.sff", "sff")
>>> Traceback (most recent call last):
>>> File "", line 1, in
>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/__init__.py",
>>>line 412, in write
>>> count = writer_class(handle).write_file(sequences)
>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>line 834, in write_file
>>> for record in records:
>>> File "", line 2, in get_generator
>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/__init__.py",
>>>line 532, in parse
>>> for r in i:
>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>line 749, in SffIterator
>>> trim)
>>> File "/usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py",
>>>line 566, in _sff_read_seq_record
>>> % padding)
>>> ValueError: Post quality 6 byte padding region contained data
>>>>>>
>>>
>>> Now when I parse and write the SFF file written by sff-dump-2.1.6:
>>>
>>>>>> _sff_records = get_generator("ERR016587.sff")
>>>>>> _wrote1 = SeqIO.write(_sff_records,
>>>>>>"ERR016587-from-biopython.sff", "sff")
>>> /usr/lib64/python2.7/site-packages/Bio/SeqIO/SffIO.py:1002:
>>>UserWarning: Read GINDECP06F3I2S has file offset 4244900184, which
>>>is too large to store in the Roche SFF index structure. No index
>>>block will be recorded.
>>> "index block will be recorded." % (name, offset))
>>>>>> _wrote1
>>> 1237359
>>>>>>
>>> [1]+ Stopped python
>>> $ sffinfo ERR016587-from-biopython.sff | head
>>> Common Header:
>>> Magic Number: 0x2E736666
>>> Version: 0001
>>> Index Offset: 0
>>> Index Length: 0
>>> # of Reads: 1237359
>>> Header Length: 1440
>>> Key Length: 4
>>> # of Flows: 1400
>>> Flowgram Code: 1
>>> $ ls -la ERR016587-from-biopython.sff
>>> -rw------- 1 mmokrejs users 6736514448 Mar 16 12:35
>>>ERR016587-from-biopython.sff
>>> $ md5sum ERR016587*.sff
>>> 57223d8ae31eed29fe9e8f60cc71a02a ERR016587-bad.sff
>>> 72ebe1c4b44ef54731b40e3915401eec ERR016587-from-biopython.sff
>>> 57223d8ae31eed29fe9e8f60cc71a02a ERR016587-from-new-sra-with-2.1.9.sff
>>> 72ebe1c4b44ef54731b40e3915401eec ERR016587-re-extracted-2.1.6.sff
>>> 72ebe1c4b44ef54731b40e3915401eec ERR016587.sff
>>> $
>>>
>>> So, biopython writes exactly same SFF file back (so its parsing
>>>and writing works).
>>> Biopython chokes on file created by sff-dump-2.1.9 while not on
>>>the one from sff-dump-2.1.6.
>>>
>>> The "UserWarning:" is explained here by Peter:
>>>

>>> That's not an error, its a warning. Your SFF file is over 4 GB, which
>>> is the limit for the Roche index. Therefore biopython won't
>>>include an index
>>> at teh ond of the SFF file.
> >> You'll see the same if you combine to large SFF files with the Roche
>>> tools.
>>>

>>>
>>> Not sure if this is what is causing sff-dump-2.1.9 to break.
>>>
>>> Good luck,
>>> Martin
>>> P.S.: I got similar errors with few other SFF files. I suppose I
>>>also created them
>>> by sff-dump-2.1.9 but could study these more would you not be
>>>able to reproduce this.
>>>
>>> _______________________________________________
>>> Sra-tools mailing list
>>> Sra-tools@ncbi.nlm.nih.gov
>>> http://www.ncbi.nlm.nih.gov/mailman/listinfo/sra-tools
>>
>>
>
>_______________________________________________
>Sra-tools mailing list
>Sra-tools@ncbi.nlm.nih.gov
>http://www.ncbi.nlm.nih.gov/mailman/listinfo/sra-tools