#!/usr/bin/perl -w

=head1 NAME 

fastm9_to_table  - turn FASTA -m 9 output into NCBI -m 9 tabular output

=head1 SYNOPSIS

 fastm9_to_table [-e evaluefilter] [-b bitscorefilter] [--header] [-o outfile] inputfile1 inputfile2 ... 

=head1 DESCRIPTION

Comand line options:
  -e/--evalue   evalue    -- filter by evalue
  -b/--bitscore bitscore  -- filter by bitscore
  --header                -- boolean flag to print column header 
  -o/--out                -- optional outputfile to write data, 
                             otherwise will write to STDOUT
  -h/--help               -- show this documentation

Not technically a SearchIO script as this doesn't use any Bioperl
components but is a useful and fast.  The output is tabular output
with the standard NCBI -m9 columns.

 queryname
 hit name
 percent identity
 alignment length
 number mismatches 
 number gaps
 query start  (if on rev-strand start > end)
 query end 
 hit start (if on rev-strand start > end)
 hit end 
 evalue
 bit score

Additionally 3 more columns are provided
 fasta score
 sw-score
 percent similar
 query length
 hit length
 query gaps
 hit gaps

=head1 AUTHOR - Jason Stajich

Jason Stajich jason_at_bioperl-dot-org

=cut

use strict;
use Getopt::Long;
my $hitsection = 0;
my %data;

my ($evalue,$bitscore,$header,$outfile) = ( 10, 0);
GetOptions(
	   'b|bitscore|bits:f'   => \$bitscore,
	   'e|evalue:f'          => \$evalue,
	   'header'              => \$header,
	   'o|out|outfile:s'     => \$outfile,
	   'h|help'              => sub { exec('perldoc',$0); exit; }
	   );

my $outfh;
if( $outfile ) { 
    open($outfh, ">$outfile") || die("$outfile: $!");
} else { 
    $outfh = \*STDOUT; 
}
	   
# query start -- an0
# query en    -- ax0
# hit start   -- an1
# hit end     -- ax1

my @fields = qw(qname hname percid alen mmcount gapcount 
		qstart qend hstart hend evalue score bits fs sw-score 
		percsim qlen hlen qgap hgap);
print $outfh "#",uc(join("", map{ sprintf("%-10s",$_) } @fields)), "\n" if $header;

while(<>) {
    my $linestr = $_;
    if( /^\s*\d+>>>(\S+).+/ ) {	
	$data{'qname'} = $1;
	if( /\-?\s+(\d+)\s+(aa|nt)\s+$/ ){
	    $data{'qlen'} = $1;
	}
    } elsif( $hitsection && (/^>>>\Q$data{'qname'}/ || /^>>>/) ) {	
	$hitsection = 0;
    } elsif( /^The best scores are:/ ) {	       
	$hitsection = 1;
    } elsif( /^\s+$/ ) {
    } elsif( $hitsection ) {
	if( s/^(\S+)\s+(.+)\(\s*(\d+)\)\s+// ) {
	    my ($hit, $desc,$hitlen) = ($1,$2,$3);
	    my ($dir) = ( s/^\[(r|f)\]\s+// );
	    my @line = split(/\s+/,$_);
	    $data{'hname'}    = $hit;
	    $data{'hlen'}     = $hitlen;
	    $data{'score'}    = shift @line;
	    $data{'bits'}     = shift @line;
	    $data{'evalue'}   = shift @line;
	    $data{'percid'}   = shift @line;
	    
	    $data{'percsim'}  = shift @line;
	    $data{'sw-score'} = shift @line;
	    $data{'alen'}     = shift @line;
	    $data{'qstart'}   = shift @line;
	    $data{'qend'}     = shift @line;
	    $data{'pn0'}      = shift @line; # pn0
	    $data{'px0'}      = shift @line; # px0
	    $data{'hstart'}   = shift @line; # an1 
	    $data{'hend'}     = shift @line; # ax1
	    $data{'pn1'}      = shift @line; # pn1
	    $data{'px1'}      = shift @line; # px1
	    # query + hit gaps
	    $data{'qgap'}     = shift @line;
	    $data{'hgap'}     = shift @line;
	    $data{'gapcount'} = $data{'qgap'} + $data{'hgap'};
	    $data{'fs'}       = shift @line;
	    
	    $data{'mmcount'} = $data{'alen'} - ( int($data{'percid'} * $data{'alen'}) + $data{'gapcount'});
	    #next if( $data{'evalue'} > $evalue || 
	#	     $data{'bits'} < $bitscore );
	    
	    for ( $data{'percid'}, $data{'percsim'} ) {
		$_ = sprintf("%.2f",$_*100);
	    }
	    print $outfh join( "\t",map { $data{$_} } @fields),"\n";
	} else { 
	    # print STDERR "unrecognized line \n$linestr";
	}
    } else { 
	# warn("skipping a line like this: $_");
    }
}