[Bioperl-l] Parsing Blast output for nohit sequences
Sammons, Scott
zno6@cdc.gov
Wed, 31 Oct 2001 12:03:50 -0500
Greetings:
I am a relatively new bioperler. I am trying to parse a large batch of
blast
searches and build a fasta file of the subjects that have no hits. The code
below works fine when I want to do something with the hits. But by the time
process_blast gets to the data, the nohits are already parsed out. I can't
seem to get to the blast parse while the nohits are still there. Any help
would
be appreciated.
-Scott Sammons
Centers for Disease Control and Prevention
1600 Clifton Rd, NE
MS G36
Atlanta, GA 30333
404-639-3560 (Voice)
404-639-1331 (FAX)
ssammons@cdc.gov
==============
use Bio::Index::Fasta;
use Bio::SeqIO;
use Bio::Tools::Blast qw(:obj);
my ($inx, $seqobj, $outputfile, $blastobj, $params);
# parse blast output and create a fasta file from subject sequences with
nohits
# First, we create a index from the fasta file that was used to create the
blast output
# call this with cat blastoutfile | nohit_blast2fasta.pl fastafile
sub process_blast {
my $hit;
my $blast_object = shift;
if ($blast_object->num_hits('total') == 0) {
$seqobj = $inx->fetch($blast_object->name);
printf ">%s %s\n",$seqobj->display_id,$seqobj->desc;
printf "%s\n", $seqobj->seq;
}
$blast_object->destroy;
}
$inx = Bio::Index::Fasta->new(
-filename => "/tmp/t$$.idx",
-write_flag => 1);
$inx->make_index(@ARGV);
$Blast->parse( -parse => 1,
-exec_func => \&process_blast,
);