[Bioperl-l] Writing blast report with HitTableWriter

Vilanova,David,LAUSANNE,NRC/BS david.vilanova@rdls.nestle.com
Mon, 24 Jun 2002 18:00:02 +0200


Hello guys,
I've been playing with a script i found on the exemple section. I try to
customize it but i have a problem.
Basically all the fields are equal to 0 in the output after the parsing of
the blast report.
Does anyone know why ???
Thanks,
David

Output below:
---	-------	-------
	0		0	0	0.0e+00	0.00	0	0

	0		0	0	0.0e+00	0.00	0	0

	0		0	0	0.0e+00	0.00	0	0


Script:

use Bio::SeqIO; 
use Bio::Tools::Run::StandAloneBlast;
use Bio::SearchIO;
use Bio::SearchIO::Writer::HitTableWriter;

#Check STDIN
die "Usage: perl script.pl input_file_to_blast output_file\n" unless @ARGV
eq '2';

($file,$output) = @ARGV;
$raw_blast_output = "raw_blast.".$output;

printf ("Preparing files....\n");

open (BLAST_REPORT,">$output") or die "cannot open $output for writing\n";


#Define blast parameters.
my $blast = Bio::Tools::Run::StandAloneBlast->new(
						  database => 'nrdb',
						  program => 'blastx',
						  '_READMETHOD' => 'Blast',
						  # to automatically create
a
						  # Bio::Tools::SearchIO
object rather
						  # than a
Bio::Tools::BPlite one
						  
						 );

# Here I add more blast parameters
$blast->a(8); #nb of cpu
$blast->e(0.0000000001); #e-10
#$blast->m(7); #generate XML output
$blast->o($raw_blast_output);

#*********open file containing the sequences to blast**********

my $Seq_in = Bio::SeqIO->new ( -file   => $file,   
			       -format => 'fasta');		  


$seq = $Seq_in->next_seq(); #Read sequence
$blast_out = $blast->blastall($seq); #Run blast



# These are the columns that will be in the output table of BLAST results.
my @columns = qw(
		 query_name
		 query_length
                 hit_name
                 hit_length
		 num_hsps
                 expect
                 frac_identical_query
                 length_aln_query
                 gaps_total
                 strand_query
                 strand_hit
		);





my $in     =  Bio::SearchIO->new( -format => 'blast',
				  -hit_filter => $filt_func,
				  -file => $raw_blast_output,
				 );

my $writer = Bio::SearchIO::Writer::HitTableWriter->new( -columns =>
\@columns, 
						       );
my $out    = Bio::SearchIO->new(  -writer => $writer,
				  -file   => ">$output" );
my $hit_count = 0;

while ( my $blast_result = $in->next_result() ) {
  #printf STDERR "\nReport %d: $blast\n", $in->report_count;
  
  if( $blast_result->hits ) {
    $hit_count++;
    $out->write_result($blast_result, $hit_count==1 );
  }
  else {
    print STDERR "Hitless Blast Report: ";
    print STDERR ($blast->no_hits_found ? "\n" : "(filtered)\n");
  }
}