pick up some unigene sequences based on a id set

[Lisa]

Hello,

I am wondering if anyone would like to share some scripts to pick up sequences from fasta format dataset (Unigene) based on a group of id. Thanks a lot.


Lisa

[simonandrews]

Rather than write scripts to do this I'd normally create a local database using the tools in EMBOSS (dbxfasta for a fasta database). Once you've indexed it then you can do batch retrieval of sequences really easily using the seqret program.

If you'd rather handle this in a script then there are existing indexing schemes available in most languages. For Perl you could use a simple Fasta index via something like Bio::Index::Fasta or go the more formal route to create a Bio:: DB::Flat database

[Lisa]

Thanks a lot. That would definitely help!
Lisa

[maubp]

If you want to using scripting, for a Python option you might try Biopython's Bio.SeqIO.index() functionality which can be used on a FASTA file or many other formats. See:

http://news.open-bio.org/news/2009/09/biopython-seqio-index/

[Lisa]

I downloaded a Unigene set fron NCBI, which is in fasta format as the following,

>gnl|UG|Les#S50866699 Solanum lycopersicum cDNA, clone: LEFL1002BB12, HTC in leaf /gb=AK319185 /gi=225321
324 /ug=Les.1 /len=3282
CAACTCCCAATTC......

I created a index with the following sample code ( eg, test1.pl) based the data set with the command as

perl test1.pl myindex Les.seq.uniq

# Complete code for making an index for several
# fasta files ----------test1.pl
use Bio::Index::Fasta;
use strict;
my $Index_File_Name = shift;
my $inx = Bio::Index::Fasta->new(
'-filename' => $Index_File_Name,
'-write_flag' => 1);
$inx->make_index(@ARGV);


Then, I tried the following sample code ( test2.pl) to retrieve a sequence with the command as

perl test2.pl myindex UG:Les#S50866699

# Print out several sequences present in the index
# in Fasta format --------test2.pl
use Bio::Index::Fasta;
use strict;
my $Index_File_Name = shift;
my $inx = Bio::Index::Fasta->new('-filename' => $Index_File_Name);
my $out = Bio::SeqIO->new('-format' => 'Fasta','-fh' => *STDOUT);
foreach my $id (@ARGV) {
my $seq = $inx->fetch($id); # Returns Bio::Seq object
$out->write_seq($seq);
}

Finally, I got this error message


------------- EXCEPTION -------------
MSG: Did not provide a valid Bio::PrimarySeqI object
STACK Bio::SeqIO::fasta::write_seq /usr/lib/perl5/site_perl/5.8.8/Bio/SeqIO/fasta.pm:178
STACK toplevel retrieve_seq.pl:12
-------------------------------------


What is the wrong I made? Please help! Thanks a lot!

[maubp]

I'd try asking on the BioPerl mailing list - they're a helpful bunch