[Bioperl-l] Parsing Blast Output
Brian Osborne
brian_osborne@cognia.com
Tue, 11 Jun 2002 09:52:22 -0400
Mat and Anthony,
I'd recently modified the example code in bptutorial.pl so that it runs,
I've inserted it below. Note there's no call to next_result(). I'd also
modified the example code in RemoteBlast.pm, it should run as well but I'm
guessing these fixes are only in version bioperl-live and version 1.0.1.
Perhaps you're using 1.0?
Brian O.
$run_remoteblast = sub {
print "\nBeginning run_remoteblast example... \n";
eval { require Bio::Tools::Run::RemoteBlast; };
if ( $@ ) {
print STDERR "Cannot load Bio::Tools::Run::RemoteBlast\n";
print STDERR "Cannot run run_remoteblast example:\n$@\n";
} else {
my (@params, $remote_blast_object, $blast_file, $r, $rc,
$database);
$database = 'ecoli';
@params = ('-prog' => 'blastp',
'-data' => $database,
'-expect' => '1e-10',
'-readmethod' => 'BPlite' );
$remote_blast_object = Bio::Tools::Run::RemoteBlast->new(@params);
$blast_file = Bio::Root::IO->catfile("t","data","ecolitst.fa");
$r = $remote_blast_object->submit_blast( $blast_file);
print "submitted Blast job\n";
while ( my @rids = $remote_blast_object->each_rid ) {
foreach my $rid ( @rids ) {
$rc = $remote_blast_object->retrieve_blast($rid);
"retrieving results...\n";
if( !ref($rc) ) { # $rc not a reference => either error
# or job not yet finished
if( $rc < 0 ) {
$remote_blast_object->remove_rid($rid);
print "Error return code for BlastID code $rid ... \n";
}
sleep 5;
} else {
$remote_blast_object->remove_rid($rid);
while ( my $sbjct = $rc->nextSbjct ) {
print "sbjct name is ", $sbjct->name, "\n";
while ( my $hsp = $sbjct->nextHSP ) {
print "score is ", $hsp->score, "\n";
}
}
}
}
}
}
return 1;
} ;
-----Original Message-----
From: bioperl-l-admin@bioperl.org [mailto:bioperl-l-admin@bioperl.org]On
Behalf Of Wiepert, Mathieu
Sent: Tuesday, June 11, 2002 9:12 AM
To: 'bioperl-l@bioperl.org'
Subject: RE: [Bioperl-l] Parsing Blast Output
I got a bounce message, if this shows up twice, I apologize.
Hi,
I see that you copied the example from the RemoteBlast.pm, and it isn't
working. I modified your code, and got an error from BPLite myself:
Can't locate object method "next_result" via package "Bio::Tools::BPlite" at
testremoteblast.pl line 32
Which is true (according to the documentation, I didn't look at the code for
the module). But, the docs say basically say that your line of code:
my $rc = $factory->retrieve_blast($rid);
is already returning the BPLite object, so you don't need the line:
my $result = $rc->next_result (you already have it in $rc)
I modified the code to work (there were other errors in the script, the POD
should be fixed perhaps, to have working code). My command line was
(because I don't have 1.0 bioperl installed):
perl -I ~/bioperl_latest/lib/site_perl/5.6.0 testremoteblast.pl
#!/usr/bin/perl -w
use Bio::Tools::Run::RemoteBlast;
use strict;
my $v = 1;
my $prog = 'blastn';
my $db = 'nr';
my $e_val= '1e-10';
my @params = ( '-prog' => $prog,
'-data' => $db,
'-expect' => $e_val );
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
$v = 1;
my $str = Bio::SeqIO->new(-file=>'comt.cdna.fa' , '-format' => 'fasta' );
my $input = $str->next_seq();
# Blast a sequence against a database:
my $r = $factory->submit_blast($input);
print STDERR "waiting..." if( $v > 0 );
while ( my @rids = $factory->each_rid ) {
foreach my $rid ( @rids ) {
my $rc = $factory->retrieve_blast($rid);
if( !ref($rc) ) {
if( $rc < 0 )
$factory->remove_rid($rid);
}
print STDERR "." if ( $v > 0 );
sleep 5;
} else {
$factory->remove_rid($rid);
# my $result = $rc->nextSbjct;
print "db is ", $rc->database(), "\n";
print "query is ", $rc->query(), "\n";
my $count = 0;
while( my $hit = $rc->nextSbjct )
$count++;
next unless ( $v > 0);
print "hit name is ", $hit->name, "\n";
while( my $hsp = $hit->nextHSP ) {
print "score is ", $hsp->score, "\n";
}
}
}
}
}
If there are errors in my code, I apologize, I am not the best at this
either :-)
-Mat
_______________________________________________
Bioperl-l mailing list
Bioperl-l@bioperl.org
http://bioperl.org/mailman/listinfo/bioperl-l