[Bioperl-l] Parsing Blast output for nohit sequences
Marc Logghe
Marc.Logghe@devgen.com
Mon, 5 Nov 2001 08:18:36 +0100
Hi,
I had the same problem a while ago, at least when using Bio::Tools::Blast.
You should try Bio::Tools::BPLite, this 'no hits' bug appears not to show up
there.
Hope this helps,
Marc
> -----Original Message-----
> From: Sammons, Scott [mailto:zno6@cdc.gov]
> Sent: Wednesday, October 31, 2001 6:04 PM
> To: 'bioperl-l@bioperl.org'
> Subject: [Bioperl-l] Parsing Blast output for nohit sequences
>
>
> 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,
> );
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l@bioperl.org
> http://bioperl.org/mailman/listinfo/bioperl-l
>