[Bioperl-l] different results with remote-blast skript
Jonas Schaer
Jonas_Schaer at gmx.de
Sun Jun 28 10:15:18 UTC 2009
Hi again :)
please, I only have this little question:
why do I get different results with my remote::blast perl skript then on the ncbi blast homepage?
I am using blastp, the query is an amino-sequence (different results with any sequence, differences not only in number of hits but even in e-values, scores etc...), the database is 'nr'.
PLEASE help me,
thank you in advance,
Jonas
ps: my skript:
################################################################################
use Bio::Seq::SeqFactory;
use Bio::Tools::Run::RemoteBlast;
use strict;
my @blast_report;
my $prog = 'blastp';
my $db = 'nr';
my $e_val= '1e-10';
#my $e_val= '10';
my @params = ( '-prog' => $prog,
'-data' => $db,
'-expect' => $e_val,
'-readmethod' => 'SearchIO' );
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
$Bio::Tools::Run::RemoteBlast::HEADER{'GAPCOSTS'} = '11 1';
$Bio::Tools::Run::RemoteBlast::HEADER{'MAX_NUM_SEQ'} = '100';
$Bio::Tools::Run::RemoteBlast::HEADER{'EXPECT'} = '10';
$Bio::Tools::Run::RemoteBlast::HEADER{'COMPOSITION_BASED_STATISTICS'} = '1';
my $blast_seq='MGSSSVGTYHLLLVLMGAGGEQQAVQAGAEVASTEQVDGSGMAANSRGSTSGSEQPPRDSDLGLLRSLLDVAGVDRTALEVKLLALAEAGAEMPPAQDSQATAAGVVATLTSVYRQQVARAWHERDDNAFRQAHQNTAMATGPDPDDEYE';
#$v is just to turn on and off the messages
my $v = 1;
my $seqbuilder = Bio::Seq::SeqFactory->new('-type' => 'Bio::PrimarySeq');
my $seq = $seqbuilder->create(-seq =>$blast_seq, -display_id => "$blast_seq");
my $filename='temp2.out';
my $r = $factory->submit_blast($seq);
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 );
}
else
{
my $result = $rc->next_result();
$factory->save_output($filename);
$factory->remove_rid($rid);
print "\nQuery Name: ", $result->query_name(), "\n";
while ( my $hit = $result->next_hit )
{
next unless ( $v > 0);
print "\thit name is ", $hit->name, "\n";
while( my $hsp = $hit->next_hsp )
{
print "\t\tscore is ", $hsp->score, "\n";
}
}
}
}
}
@blast_report = get_file_data ($filename);
return @blast_report;
##################################################################################
More information about the Bioperl-l
mailing list