http://qs321.pair.com?node_id=1204589

TJCooper has asked for the wisdom of the Perl Monks concerning the following question:

I currently load a FASTA formatted file, which contains DNA sequence separated by headers denoted by '>' i.e.:

>header1 AGACATGATCGACTGGACACAATTTACGTAGCTG >header2 AAATACTAGGGCAACACACACACACACCACACAC >header3 AGTTAGTCCAGTAGTTTACAGTACGACTGATCGT

via bioPerl (Bio::SeqIO), and store the entire FASTA file into memory as a hash:

my $fasta = $ARGV[0]; my $seqio = Bio::SeqIO->new(-file => $fasta); while(my $seqobj = $seqio->next_seq) { my $id = $seqobj->display_id; my $seq = $seqobj->seq; $sequences{$id} = $seq; } open my $IN2, '<', $ARGV[1] or die "$!"; while (<$IN2>) { # Subsample the FASTA sequences }

Or I could code it manually (which is likely quicker given the amount of overhead in Bio::SeqIO).

Based on a second input file, I then subsample the FASTA sequences and do some processing.

Each line of $IN2 includes the header from its associated FASTA sequence.

This works fine, and on modern computers the memory usage isn't too much of a problem even with larger genomes (Human etc.). Nevertheless, i'd like to take advantage of the fact that my second input file ($IN2) is always sorted by header and thus, for example, the first 5000 lines of it only need the FASTA sequence from >header1. This could potentially save me having to load the entire FASTA file into memory and instead only load the sequence that is currently being used.

I'm however having a tough time trying to figure out how to actually do this so any suggestions or code would be greatly appreciated.