[Bioperl-l] Problems running blast
Barry Moore
barry.moore at genetics.utah.edu
Fri Feb 6 01:02:52 EST 2004
Hong-
You don't have to make your sequence into a fasta file. Have a look at
the documentation for the submit_blast method of the
Bio::Tools::Run::RemoteBlast module where it tells you that the input
can be a sequence object, a reference to an array of sequence objects,
or the filename of a fasta file. If your script already has your
sequence as any of the Bioperl sequence objects, then you are ready to
go. If your script has your sequence as a simple string, it is quite
easy to convert that to a PrimarySeq object which you can then submit to
BLAST. The following script (adapted from the module documentation)
suggests one way of converting a string to a PrimaySeq object and
submitting it to BLAST. See the example code in the Synopsis section of
the RemoteBlast module documentation mentioned above for examples of how
to submit a sequence object, or a fasta file to BLAST.
Barry
----------------------------------------------------------------------------------------------
#!/usr/bin/perl
use strict;
use warnings;
use Bio::PrimarySeq;
use Bio::Tools::Run::RemoteBlast;
#Your sequence as a string
my $sequence_string =
"atggagagcagaggcccactggctacctcgcgcctgctgctgttgctgctgttgctacta";
#Initialize string as new sequence
my $seq = new Bio::PrimarySeq(-seq => $sequence_string,
-display_id => "Your_favorite_gene");
#Build the BLAST factory
my $BLAST_factory = Bio::Tools::Run::RemoteBlast->new('-prog' =>
'blastn',
'-data' => 'nr',
'-expect' => .001,
'-readmethod' =>
'SearchIO' );
#Submit the sequence object to NCBI's BLAST server
my $job = $BLAST_factory->submit_blast($seq);
print STDERR "Blasting sequence ";
#Load the RIDs returned for the BLAST job submitted (in this case only one)
while ( my @rids = $BLAST_factory->each_rid ) {
#Iterate over RIDs
foreach my $rid ( @rids ) {
#Hit the server for a result on RID
my $blast_results = $BLAST_factory->retrieve_blast($rid);
#Was a result returned?
if( !ref($blast_results) ) {
#If so and it returned an error remove that RID from the stack
if ($blast_results < 0) {
$BLAST_factory->remove_rid($rid);
}
print STDERR "."; #Keep staring at the dots
sleep 5; #Plays nice with the servers
}
#If a result was returned and it isn't an error, then pass it to a
variable...
else {
my $result = $blast_results->next_result();
$BLAST_factory->remove_rid($rid); #...and remove it's RID from the
stack.
#Check the result for a hit...
my $hit = $result->next_hit;
if (ref($hit)) {
my $hsp = $hit->next_hsp;
#...collect some values from the result, hit and hsp objects and do
#something with them.
my $q_name = $result->query_name();
my $h_name = $hit->name;
my $evalue = $hsp->evalue();
print "\nQuery name: $q_name\nHit name: $h_name\nLowest
e-value: $evalue\n";
}
}
}
}
Hong Ching Lee wrote:
>Hey everyone,
>
>I have a question about whether i can run remote blast using just a string
>or whether i have to make it into a fasta format file. Can anyone help me
>with this?
>
>Thank You,
>Hong
>_______________________________________________
>Bioperl-l mailing list
>Bioperl-l at portal.open-bio.org
>http://portal.open-bio.org/mailman/listinfo/bioperl-l
>
More information about the Bioperl-l
mailing list