bin/bp_search2table

   1 #!/usr/bin/perl
   2
   3 =head1 NAME
   4
   5 bp_search2table - turn SearchIO parseable reports into tab delimited format like NCBI's -m 9
   6
   7 =head1 SYNOPSIS
   8
   9   bp_search2table -f fasta -i file.FASTA -o output.table
  10
  11 =head1 DESCRIPTION
  12
  13 Turn SearchIO reports into a tabular format like NCBI's -m 9 output.
  14
  15 =head1 FEEDBACK
  16
  17 =head2 Mailing Lists
  18
  19 User feedback is an integral part of the evolution of this and other
  20 Bioperl modules. Send your comments and suggestions preferably to
  21 the Bioperl mailing list.  Your participation is much appreciated.
  22
  23   bioperl-l@bioperl.org                  - General discussion
  24   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  25
  26 =head2 Reporting Bugs
  27
  28 Report bugs to the Bioperl bug tracking system to help us keep track
  29 of the bugs and their resolution. Bug reports can be submitted via
  30 email or the web:
  31
  32   https://github.com/bioperl/bioperl-live/issues
  33
  34 =head1 AUTHOR
  35
  36   Jason Stajich jason_at_bioperl-dot-org
  37
  38 =cut
  39
  40 use strict;
  41 use warnings;
  42 use Bio::SearchIO;
  43 use Getopt::Long;
  44
  45 my ($format, $file,$output) = ('blast');
  46
  47 GetOptions(
  48            'f|format:s'   => \$format,
  49            'i|input:s'    => \$file,
  50            'o|output:s'   => \$output);
  51
  52 if( @ARGV ) {
  53     $file = shift;
  54 }
  55
  56 my $in = Bio::SearchIO->new(-format => $format,
  57                             -file   => $file);
  58 my $out;
  59 if( $output ) {
  60     open $out, '>', $output or die "Could not write file '$output': $!\n";
  61 } else {
  62     $out = \*STDOUT;
  63 }
  64
  65 while( my $r = $in->next_result ) {
  66     while( my $hit = $r->next_hit ) {
  67         while( my $hsp = $hit->next_hsp ) {
  68             my $mismatchcount = $hsp->length('total') -
  69                 ($hsp->num_conserved + $hsp->gaps('total'));
  70             print $out join("\t", ( $r->query_name,
  71                                     $hit->name,
  72                                     sprintf("%.2f",$hsp->percent_identity),
  73                                     $hsp->length('total'),
  74                                     $mismatchcount,
  75                                     $hsp->gaps('total'),
  76                                     # flip start/end on rev strand
  77                                     $hsp->query->strand < 0 ?
  78                                     ( $hsp->query->end,
  79                                       $hsp->query->start ) :
  80                                     ( $hsp->query->start,
  81                                       $hsp->query->end ),
  82                                     $hsp->hit->strand < 0 ?
  83                                     ( $hsp->hit->end,
  84                                       $hsp->hit->start ) :
  85                                     ( $hsp->hit->start,
  86                                       $hsp->hit->end ),
  87
  88                                     $hsp->evalue,
  89                                     # chance this to $hsp->sw_score
  90                                     # if you would rather have that
  91                                     # it will only work for FASTA parsing though!
  92                                     $hsp->bits)),"\n";
  93         }
  94     }
  95 }