Re: [Seqan-dev] Random access of large FASTA file

[Johannes Dröge <johdro@mpi-inf.mpg.de>]

Hello David,
thank for your comments. I am still not confident with the design concept of memory mapped single strings in Seqan. The idea of the loop in any case is to create a StringSet which type depends on the choosen StringType. So it works this way:

1) create temporary sequences object
2) assign content from memory mapped multi-fasta file (MultiSeqFile)
3) store in StringSet which will have ownership (I guess this is done via a copy constructor)

This works fine for standard and packed string types. I would also like to have a StringSet that contains strings that are actually memory mapped from the original multi-fasta file. I thought that the assignSeq function would appropriately handle this when I use it with default-constructed memory mapped sequence object. I seems I misunderstood the design of this sequence type. Is there any way to construct such a StringSet I have in mind?

Gruß Johannes


Am Donnerstag, 7. Juli 2011 16:14:34 schrieb Weese, David:
> Hi Johannes,
> 
> I assume the value of num_records less or equal to length(db_sequences). Looking at your code it seems that you try to use a memory mapped string as a temporary variable in a large loop. Maybe not the best idea, as it would create a temporary file and deletes it in every iteration. It could be that the temporary could not be opened, you could test that with a #define SEQAN_DEBUG before including any SeqAn header.
> You should at least move all the instantiations out of the loop. Still I dont think you need a memory mapped string (seq) to store a single sequence of a multi fasta file. Also I cannot see, where you store the read sequences. It would make sense to use a single StringSet<String<..,MMap<> >, Owner<ConcatDirect<> > >  data_ that stores multiple sequences using a single memory mapped string.
> 
> HTH. If the problem still remains, please create a bug ticket with source code and example files.
> 
> Cheers,
> David
> --
> David Weese				weese@inf.fu-berlin.de
> Freie Universität Berlin		http://www.inf.fu-berlin.de/
> Institut für Informatik			Phone: +49 30 838 75246
> Takustraße 9					Algorithmic Bioinformatics
> 14195 Berlin					Room 021 
> 
> Am 06.07.2011 um 16:28 schrieb Johannes Dröge:
> 
> > Hello,
> > I am using Seqan to access a large FASTA file. In this case, I am importing the whole RefSeq DB for random access (into memory or memory-mapped). This can be quite a huge file, so I decided to go for a dynamic strategy writing a generic SequenceStorage object. It works well for
> > 
> > typedef seqan::String< seqan::Dna5 > StringType; //(default type)
> > typedef seqan::String< seqan::Dna5, seqan::Packed<> > StringType;
> > 
> > but not for
> > typedef seqan::String< seqan::Dna5, seqan::MMap<> > StringType;
> > 
> > Here is the Code that imports the data using the MMap-Trick from the HowTo and put it into a
> > 
> > StringSet< StringType > data_;
> > 
> > with an index data structure 
> > 
> > std::map< std::string, long unsigned int > id2pos_;
> > 
> > --------------------------------------------------------------------------------
> > seqan::MultiSeqFile db_sequences;
> > seqan::open( db_sequences.concat, filename.c_str(), seqan::OPEN_RDONLY );
> > seqan::split( db_sequences, seqan::Fasta() );
> > 
> > for( unsigned int i = 0; i < num_records; ++i ) {
> > 	StringType seq;
> > 	seqan::assignSeq( seq, db_sequences[i], fasta_format_ );
> > 	
> > 	std::string id;
> > 	seqan::assignSeqId( id, db_sequences[i], fasta_format_ );
> > 	id2pos_[ extractFastaCommentField( id, "gi" ) ] = seqan::assignValueById( data_, seq );
> > }
> > --------------------------------------------------------------------------------
> > 
> > 1) seqan::assignValueById() will cause a segfault at sequence number 33,924 out of 276,313 when using a StringSet with mmap strings.
> > 
> > 2) Also, I don't know how to define a StringSet using array strings.
> > 
> > 3) Using a regular Dna5 string, the how operation will take about 5 minutes. A packed string requires much longer to load. Is there any way to speed this up? I could think of a (binary) sink for a StingSet to avoid parsing and recoding every time I load the DB sequences. Is there anything like this (planned)?
> > 
> > I appreciate your help!
> > 
> > Gruß Johannes