examples/searchio/hspwriter.pl

   1 #!/usr/bin/env perl
   2
   3 # Demonstrates the use of a SearchIO Blast parser and a SearchWriterI object
   4 # for producing tab-delimited output of HSP data from a Blast report
   5 # input stream.
   6 #
   7 # Each row in the output represents data for a single HSP.
   8 #
   9 # This parser represents a new and improved version of Bio::Tools::Blast.
  10 #
  11 # Usage:
  12 #   STDIN:  stream containing one or more BLAST or PSI-BLAST reports.
  13 #   STDOUT: none, but generates an output file "hspwriter.out"
  14 #           containing tab-delimited data on a per-HSP basis.
  15 #   STDERR: Progress info and any errors.
  16 #
  17 # In this example, we create a SearchIO parser that screens out hits
  18 # based on expect (or P) scores and a default HSPTableWriter. This writer
  19 # provides the same functionality as the original Bio::Tools::Blast::table2()
  20 # function (i.e., a tab-delimited summary of each hit per row).
  21 # HSPTableWriter, however, is customizable so you can specify just the columns
  22 # you want to have in the output table.
  23 #
  24 # For more documentation about the writer, including
  25 # a complete list of columns, execute:
  26 #   perldoc Bio::SearchIO::Writer::HSPTableWriter.
  27 #
  28 # For more documentation about working with Blast result objects,
  29 # see docs for these modules:
  30 #   Bio::Search::Result::BlastResult
  31 #   Bio::Search::Iteration::IterationI
  32 #   Bio::Search::Hit::BlastHit
  33 #   Bio::Search::HSP::BlastHSP
  34 #
  35 # For more documentation about the Blast parser, see docs for
  36 #   Bio::SearchIO
  37 #
  38 # Author: Steve Chervitz <sac@bioperl.org>
  39
  40 use strict;
  41 use lib '../../';
  42
  43 use Bio::SearchIO;
  44 use Bio::SearchIO::Writer::HSPTableWriter;
  45
  46 # These are the columns that will be in the output table of BLAST results.
  47 my @columns = qw(
  48                  query_name
  49                  query_length
  50                  hit_name
  51                  hit_length
  52                  rank
  53                  expect
  54                  frac_identical_query
  55                  length_aln_query
  56                  gaps_total
  57                  strand_query
  58                  strand_hit
  59                 );
  60
  61
  62 print STDERR "\nUsing SearchIO->new()\n";
  63
  64 # Note that all parameters for the $in, $out, and $writer objects are optional.
  65 # Default in = STDIN; Default out = STDOUT; Default writer = all columns
  66 # In this example, we're reading from STDIN and  writing to a STDOUT
  67 my $in     = Bio::SearchIO->new( -format => 'blast',
  68                                  -fh => \*ARGV
  69                                );
  70 my $writer = Bio::SearchIO::Writer::HSPTableWriter->new( -columns => \@columns );
  71 my $out    = Bio::SearchIO->new( -format => 'blast',
  72                                  -writer => $writer,
  73                                  -file   => ">hspwriter.out" );
  74
  75 while ( my $result = $in->next_result() ) {
  76   printf STDERR "\nReport %d: $result\n", $in->result_count;
  77
  78   if( $result->hits ) {
  79     $out->write_result($result, ($in->result_count - 1 ? 0 : 1) );
  80   }
  81   else {
  82     print STDERR "Hitless Blast Report: $result ";
  83     print STDERR ($result->no_hits_found ? "\n" : "(filtered)\n");
  84   }
  85
  86   ## For a simple progress monitor, uncomment this line:
  87   #print STDERR "."; print STDERR "\n" if $in->result_count % 50 == 0;
  88 }
  89
  90 printf STDERR "\n%d Blast report(s) processed.\n", $in->result_count;
  91 printf STDERR "Output sent to file: %s\n",  $out->file if $out->file;
  92