FASTAid - index a FASTA sequence database for fast random-access sequence retrieval

    use FASTAid qw( create_index retrieve_entry );

    my $FASTA_file = 'lots_of_seqs.fa';

    # index the file of FASTA seqs
    create_index($FASTA_file);

    # retrieve one or more FASTA seqs from the file
    my $seq_array_ref = retrieve_entry($FASTA_file, 'NM_00204', 'AA23456');

    foreach my $seq ( @{$seq_array_ref} ) {

        ...do something with each FASTA sequence...
    }

 Jarret Glasscock C<< <glasscock_cpan@mac.com> >>
 Dave Messina C<< <dave-pause@davemessina.net> >>

 This software was developed at the Genome Sequencing Center at Washington
 University, St. Louis, MO.

 Copyright (C) 2004-6 Glasscock, Messina. All Rights Reserved.

 This software is provided "as is" without warranty of any kind.

package FASTAid;

use version; $VERSION = qv('0.0.4');

# PRAGMAS
use strict;
use warnings;

# INCLUDES
use Carp;

sub create_index {
    my ($fasta) = @_;
    my $index = $fasta . '.fec';

    my ( %data, $begin, $id );
    open( DB, $fasta )
        or croak( qq{could not open FASTA file $fasta:\n}, qq{$!\n} );

    # record offsets of records in perl database
    while (<DB>) {
        if ( /^>(\S+)/) {
            $id    = $1;
            $begin = tell(DB) - length($_);
			if ( defined $data{$id} ) {
	            croak "There is already an entry ID = $id\n";
	        }
			else { $data{$id} = "$begin" }
	    }
    }
    close DB;

	# test for empty index
	return 0 if (scalar (keys %data) == 0);

    # Write out index
    open( OUT, ">$index" )
        or croak( qq{Cannot write FASTAid index $index:\n}, qq{$!\n} );
    foreach my $key ( sort { $a cmp $b } keys %data ) {
        print OUT $key, " ", $data{$key}, "\n";
    }
    close OUT;

    return 1;
}

sub retrieve_entry {
    my ( $fasta, @ids ) = @_;

	# make sure we have IDs to retrieve
	croak "Must supply at least one ID" if scalar @ids == 0;

    my @seqs;    # where we store the FASTAs we're returning

    # Get the INDEX into an array so we can easily work with it
    my ( @IND1, @IND2 );
    my $index = $fasta . '.fec';

    open( INDEX, $index )
        or croak( qq{cannot open FASTAid index $index:\n}, qq{$!\n} );
    while (<INDEX>) {
        if ( $_ =~ /^(\S+)\s+(\S+)$/ ) {
            push @IND1, $1;
            push @IND2, $2;
        }
    }
    close INDEX;

    # Get the FASTA file into an index before each $id is looked for
    open( DB, "$fasta" )
        or croak( qq{cannot open FASTA file $fasta:\n}, qq{$!\n} );

    foreach my $id (@ids) {
        my ( $low, $high ) = ( 0, @IND1 - 1 );
        my $seq;    # the FASTA we're returning

        my $try;
        TRY:
        while ( $low <= $high ) {
            $try = int( ( $low + $high ) / 2 );    # middle
            $low = $try + 1, next TRY if $IND1[$try] lt $id;
            $high = $try - 1, next TRY if $IND1[$try] gt $id;
            last;                                  # found it!!!
        }

        # translate this back to the corresponding value
        # in the byte index portion of lookup
        my $trans_try = $IND2[$try];

        # go to the offset in the FASTA file
        if ( $trans_try >= 0 && $IND1[$try] eq $id ) {
            seek DB, $trans_try, 0;    # position the file pointer there
            my $def = <DB>;            # defline is first line
            $seq = $def;               # put the defline into our
                                       # return string

            ENTRY:
            while (<DB>) {             # save the other lines
                last ENTRY if $_ =~ /^>/;
                $seq .= $_;
            }

            # add the seq to the array we're returning
            push @seqs, $seq;
        }
        else {
            carp "Entry ID = $id not found!\n";
        }
    }

    # send back an array of FASTA sequences
    return \@seqs;
}

1;