[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");
}
}