examples/searchio/psiblast_features.pl

   1 #!/usr/bin/perl
   2
   3 # Example usage of a SearchIO::psiblast parser of traditional format Blast
   4 # and PSI-Blast reports.
   5 # Illustrates how to grab a set of SeqFeatures from a Blast report.
   6 # This parser represents a new and improved version of Bio/Tools/Blast.pm.
   7 #
   8 # Usage:
   9 #   STDIN:  stream containing one or more BLAST or PSI-BLAST reports.
  10 #   STDOUT: feature start, end data
  11 #   STDERR: Processing info, such as the number of reports processed
  12 #           and the number of hitless reports.
  13 #
  14 # For more documentation about working with Blast result objects,
  15 # see to documentation for these modules:
  16 #   Bio::Search::Result::BlastResult
  17 #   Bio::Search::Hit::BlastHit
  18 #   Bio::Search::HSP::BlastHSP
  19 #
  20 # For more documentation about the PSI-Blast parser, see docs for
  21 #   Bio::SearchIO::psiblast
  22 #
  23 # Author: Steve Chervitz <sac@bioperl.org>
  24
  25 use strict;
  26 use lib '../../';
  27 use Bio::SearchIO;
  28
  29 my $in = Bio::SearchIO->new( -format => 'psiblast',
  30                              -fh => \*ARGV,
  31                              -signif => 0.1,
  32                              -verbose => 0 );
  33 my @hitless_reports = ();
  34
  35 while ( my $blast = $in->next_result() ) {
  36
  37     if( $blast->hits ) {
  38       while( my $feature = $blast->next_feature() ) {
  39         print "Feature from ", $feature->start, " to ", $feature->end, "\n";
  40       }
  41     }
  42     else {
  43       push @hitless_reports, $blast;
  44     }
  45 }
  46
  47 printf STDERR "\n%d Blast report(s) processed.\n", $in->result_count;
  48 printf STDERR "\n%d reports had no hits:\n", scalar(@hitless_reports);
  49
  50 foreach my $blast (@hitless_reports) {
  51     print STDERR "No hits for query ", $blast->query_name;
  52     print STDERR ($blast->no_hits_found ? "\n" : "(filtered)\n")
  53 ;
  54 }
  55
  56