[Bioperl-l] Parsing FASTA files into PrimarySeq objects
Robert Buels
rmb32 at cornell.edu
Wed Feb 17 17:52:00 UTC 2010
OK, well are you guys going to roll them back and do a branch? Broken
trunk is not good.
R
Jason Stajich wrote:
> hmm yes. Let's roll back those Florent.
>
> I guess we need to also discuss the concept of branches again here for
> major changes like this. The speed fix is helpful for short read but
> probably not for all systems so just having the ability to provide the
> right type when you use your scripts focused on short-reads might be the
> best way to go forward and preserve backwards compatibility.
>
> -jason
>
> Chris Fields wrote:
>> Yes, except that tests are now failing on live. Bio::PrimarySeq isn't
>> a Bio::SeqI, so it's a bit of an API change in a lot of places
>> expecting a Bio::SeqI. We should make this optional for the time
>> being, until we discuss it more on list.
>>
>> chris
>>
>> On Feb 16, 2010, at 10:57 AM, Jason Stajich wrote:
>>
>>
>>> Great! thanks for committing it.
>>>
>>> Florent Angly wrote:
>>>
>>>> Thanks for your reply Jason. It's very valuable since you coded the
>>>> modules in question!
>>>>
>>>> Now that I understand the purpose of these modules, I did a little
>>>> bit of documentation clarification. I also figured out that I should
>>>> implement the 'type' method in Bio::Seq::SeqFastaSpeedFactory and
>>>> made it support creating Bio::PrimarySeq. This is all committed to SVN.
>>>>
>>>> When sequences are large, you are right that it does not make much
>>>> difference in time or memory resource to use one sequence type over
>>>> the other. However, that's different for short sequences. I took a
>>>> FASTA file weighting 39.3 MB and containing 360,000 sequences of 106
>>>> bp on average.
>>>>
>>>> CASE 1:
>>>> Requested factory Bio::Seq::SeqFactory
>>>> Requested sequence type is Bio::Seq
>>>> Memory used: 329 MB
>>>> Time spent: 1m8.522s
>>>>
>>>> CASE 2:
>>>> Requested factory Bio::Seq::SeqFactory
>>>> Requested sequence type is Bio::PrimarySeq
>>>> Memory used: 246 MB
>>>> Time spent: 0m52.663s
>>>>
>>>> CASE 3:
>>>> Requested factory Bio::Seq::SeqFastaSpeedFactory
>>>> Requested sequence type is Bio::Seq
>>>> Memory used: 298 MB
>>>> Time spent: 0m31.292s
>>>>
>>>> CASE 4:
>>>> Requested factory Bio::Seq::SeqFastaSpeedFactory
>>>> Requested sequence type is Bio::PrimarySeq
>>>> Memory used: 231 MB
>>>> Time spent: 0m28.500s
>>>>
>>>> Clearly, when only simple sequence attributes are needed or FASTA
>>>> files are used, it's fast and memory efficient to use the
>>>> SeqFastaSpeedFactory with PrimarySeq.
>>>>
>>>> Florent
>>>>
>>>>
>>>> On 16/02/10 10:03, Jason Stajich wrote:
>>>>
>>>>> I don't think that aspect of the documentation vs interface was
>>>>> ever implemented - the interface object doesn't specify a type
>>>>> method or init argument even though the documentation says so. Not
>>>>> really sure why not, this was ages ago unfortunately.
>>>>>
>>>>> This particular factory definitely assumes you are building
>>>>> Bio::Seq objects - you can try and subclass and build your own to
>>>>> see if it makes much of a difference in speed/memory - I would
>>>>> posit it won't make a significant difference but be interested to
>>>>> see what you find.
>>>>>
>>>>> Just make your own Bio::Seq::PrimarySeqFactory object based on
>>>>> Bio::Seq::FastaSpeedFactory and simplify that code so it doesn't
>>>>> build Bio::Seq object wrapper around the Bio::PrimarySeq and do
>>>>> some perf tests so we'll know if it makes any difference here.
>>>>>
>>>>> -jason
>>>>> Florent Angly wrote:
>>>>>
>>>>>> Hi all,
>>>>>>
>>>>>> I am trying to reduce memory usage and speedup reading FASTA files
>>>>>> using the facilities provided by BioPerl.
>>>>>>
>>>>>> The first thing I noticed is that when using Bio::SeqIO::fasta,
>>>>>> the objects returned are Bio::Seq, not Bio::PrimarySeq objects.
>>>>>> Bio::PrimarySeq sequences are lighter than Bio::Seq sequences, so
>>>>>> it's what I need. See code below:
>>>>>>
>>>>>>> use warnings;
>>>>>>> use strict;
>>>>>>> use Data::Dumper;
>>>>>>> use Bio::SeqIO;
>>>>>>> my $in = Bio::SeqIO->new(-fh=>\*DATA);
>>>>>>> my $seqfactory = $in->sequence_factory; #
>>>>>>> Bio::Factory::ObjectBuilderI
>>>>>>> print "The factory is a ".ref($seqfactory)."\n";
>>>>>>> $seqfactory->type('Bio::PrimarySeq'); # gives an error
>>>>>>> my $seq = $in->next_seq;
>>>>>>> print Dumper($seq);
>>>>>>> __END__
>>>>>>>
>>>>>>>> seq1 a small test sequence q
>>>>>>>>
>>>>>>> ACGTACGACTACGACTAGCGCCATCAGC
>>>>>>>
>>>>>> It returns:
>>>>>>
>>>>>>> $VAR1 = bless( {
>>>>>>> 'primary_id' => 'seq1',
>>>>>>> 'primary_seq' => bless( {
>>>>>>> 'display_id' => 'seq1',
>>>>>>> 'primary_id' => 'seq1',
>>>>>>> 'desc' => 'a small
>>>>>>> test sequence',
>>>>>>> 'seq' =>
>>>>>>> 'ACGTACGACTACGACTAGCGCCATCAGC',
>>>>>>> 'alphabet' => 'dna'
>>>>>>> }, 'Bio::PrimarySeq' )
>>>>>>> }, 'Bio::Seq' );
>>>>>>>
>>>>>> Actually, we have a Bio::Seq containing a Bio::PrimarySeq. I
>>>>>> really only need the Bio::PrimarySeq. Looking at the documentation
>>>>>> for Bio::SeqIO I found that I could in theory adjust the sequence
>>>>>> factory and sequence builder to my liking... I tried this:
>>>>>>
>>>>>>> use warnings;
>>>>>>> use strict;
>>>>>>> use Data::Dumper;
>>>>>>> use Bio::SeqIO;
>>>>>>> my $in = Bio::SeqIO->new(-fh=>\*DATA);
>>>>>>> my $seqfactory = $in->sequence_factory; #
>>>>>>> Bio::Factory::ObjectBuilderI
>>>>>>> print "The factory is a ".ref($seqfactory)."\n";
>>>>>>> $seqfactory->type('Bio::PrimarySeq'); # gives an error
>>>>>>> my $seq = $in->next_seq;
>>>>>>> print Dumper($seq);
>>>>>>> __END__
>>>>>>>
>>>>>>>> seq1 a small test sequence
>>>>>>>>
>>>>>>> ACGTACGACTACGACTAGCGCCATCAGC
>>>>>>>
>>>>>> This returns:
>>>>>>
>>>>>>> The factory is a Bio::Seq::SeqFastaSpeedFactory
>>>>>>> Can't locate object method "type" via package
>>>>>>> "Bio::Seq::SeqFastaSpeedFactory" at ./seqbuilder_test_3.pl line
>>>>>>> 12,<DATA> line 1.
>>>>>>>
>>>>>> According to Bio::Seq::FastaSpeedFactory's documentation:
>>>>>>
>>>>>>> If you want the factory to create Bio::Seq objects instead
>>>>>>> of the default Bio::PrimarySeq objects, use the -type parameter
>>>>>>>
>>>>>> So, PrimarySeq should be the default type, even though in my case
>>>>>> it seems not to be. Second, I can't seem to use the -type method
>>>>>> to change what the return type is... It errors.
>>>>>>
>>>>>> Any ideas??? Thanks,
>>>>>>
>>>>>> Florent
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> Bioperl-l mailing list
>>>>>> Bioperl-l at lists.open-bio.org
>>>>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>>>>>>
>>> _______________________________________________
>>> Bioperl-l mailing list
>>> Bioperl-l at lists.open-bio.org
>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>>>
>>
>>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>
--
Robert Buels
Bioinformatics Analyst, Sol Genomics Network
Boyce Thompson Institute for Plant Research
Tower Rd
Ithaca, NY 14853
Tel: 503-889-8539
rmb32 at cornell.edu
http://www.sgn.cornell.edu
More information about the Bioperl-l
mailing list