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 Bio::SearchIO;
  27
  28 my $in = Bio::SearchIO->new( -format => 'psiblast',
  29                              -fh => \*ARGV,
  30                              -signif => 0.1,
  31                              -verbose => 0 );
  32 my @hitless_reports = ();
  33
  34 while ( my $blast = $in->next_result() ) {
  35
  36     if( $blast->hits ) {
  37       while( my $feature = $blast->next_feature() ) {
  38         print "Feature from ", $feature->start, " to ", $feature->end, "\n";
  39       }
  40     }
  41     else {
  42       push @hitless_reports, $blast;
  43     }
  44 }
  45
  46 printf STDERR "\n%d Blast report(s) processed.\n", $in->result_count;
  47 printf STDERR "\n%d reports had no hits:\n", scalar(@hitless_reports);
  48
  49 foreach my $blast (@hitless_reports) {
  50     print STDERR "No hits for query ", $blast->query_name;
  51     print STDERR ($blast->no_hits_found ? "\n" : "(filtered)\n")
  52 ;
  53 }
  54
  55