#!/usr/bin/perl 
# With this program, you can extract the top hit from the BLAST_Outfile. 
# How to use this program : $ perl blast_filter_v1.3pl BLAST_Outfile.txt (recomended BLAST output format option: -outfmt "7 sseqid evalue bitscore length nident pident")  
# You might add selection option below "evalue, identity, identical, score, alignment_length"

use Time::HiRes qw(gettimeofday);
my $t0 = gettimeofday;

use Getopt::Long;  
use Tie::IxHash;
my %seq;
tie %seq, 'Tie::IxHash';

# Set up defaults here if you wish
$evalue=0.0001;
$identity=;
$identical=0;
$score=0;
$alignment_length=0;


usage(), exit  if not GetOptions(
    'e|evalue=s'   => \$evalue,
    'iy|identity:i' => \$identity,
    'il|identical:i'   => \$identical,
    's|score:i'  => \$score,
    "l|alignment_length=i"	=> \$alignment_length,
	'h|help'	=> \$help,
);

#$file_to_get = "Infile.txt";
my $tmp_filename = $ARGV[0];
my $filename = $tmp_filename;

$filename=~ s/\.(.*)/_filter.txt/g;

open (FILE, $tmp_filename) or die $!;
open (OUT,">", $filename)  or die"output file error.\n";

if (!@ARGV or $help){
	
	print STDERR "\nVersion: 1.3 (2021/05/10)\n";
	print STDERR "Usage: blast_filter.pl <Infile.txt> [options]\n\n";
	print STDERR "options:   -e       INT  Minimum evalue to be calculated (default: 10)\n";
	print STDERR "           -s       INT  Minimum bit score to be calculated (default: 0)\n";
	print STDERR "           -l       INT  Minimum alignment length to be calculated (default: 0)\n";
	print STDERR "           -il      INT  Minimum identical to be calculated (default: 0)\n";
	print STDERR "           -iy      INT  Minimum % identity to be calculated (default: 0)\n";
	print STDERR "\n";
	exit;
	
}

while (<FILE>) {

    next if (/BLASTN/ or /Database/ or /Fields/ or /hits found/) ;
        
        if (/Query: ((.*))/) {

            $n=$1;
        }
        else{
           @b=split (" ",$_);  
                      
           if($b[1]<$evalue and $b[5]>$identity and $b[4]>$identical and $b[3]>$alignment_length and $b[2]>$score){  
               
               $seq{$n}.=$_ ;                              
               
               
               ($ID,$tax)=split("\#\#\#",$b[0] );
               @taxon=split(";",$tax);
               
               $K=$taxon[-7];
               $P=$taxon[-6];
               $C=$taxon[-5];
               $O=$taxon[-4];
               $F=$taxon[-3];
               $G=$taxon[-2];
               $S=$taxon[-1];

               $hash{$n}{$K}++;
               $hash2{$n}{$P}++;               
               $hash3{$n}{$C}++;
               $hash4{$n}{$O}++;                              
               $hash5{$n}{$F}++;
               $hash6{$n}{$G}++;
               $hash7{$n}{$S}++;                                                     
             
             }                        
            } 
}
close(FILE);


print OUT "#filter criteria:";

if($evalue!=10){
 print OUT "evlaue smaller than $evalue;";
}
if($score!=0){
 print OUT "bit score higher than $score;";
}
if($alignment_length!=0){
 print OUT "alignment length higher than $alignment_length;";
}
if($identical!=0){
 print OUT "identical higher than $identical;";
}             
if($identity!=0){
 print OUT "% identity higher than $identity;";
}

print OUT "\nSeqence_ID\tNumber_of_alignment\tEstimated_taxon\n"; 



foreach $f (  keys %seq){
     
    $N = scalar(split("\n", $seq{$f}));
    @seq1=split(" ",$seq{$f});
    
    
    ($ID,$tax)=split("\#\#\#",$seq1[0] );
    @seq1_taxon=split(";",$tax);
    $seq1_K=$seq1_taxon[-7];
    $seq1_P=$seq1_taxon[-6];
    $seq1_C=$seq1_taxon[-5];
    $seq1_O=$seq1_taxon[-4];
    $seq1_F=$seq1_taxon[-3];
    $seq1_G=$seq1_taxon[-2];
    $seq1_S=$seq1_taxon[-1];

    print OUT "$f\t$N";
    
            
      $k_percent=int(($hash{$f}{$seq1_K})*100/$N);
      print OUT "\t$seq1_K($k_percent)";     
           
      $p_percent=int(($hash2{$f}{$seq1_P})*100/$N);
      print OUT ";$seq1_P($p_percent)";     
    
      $c_percent=int(($hash3{$f}{$seq1_C})*100/$N);
      print OUT ";$seq1_C($c_percent)";     
    
      $o_percent=int(($hash4{$f}{$seq1_O})*100/$N);
      print OUT ";$seq1_O($o_percent)";     
      
      $fa_percent=int(($hash5{$f}{$seq1_F})*100/$N);
      print OUT ";$seq1_F($fa_percent)";     
   
      $g_percent=int(($hash6{$f}{$seq1_G})*100/$N);
      print OUT ";$seq1_G($g_percent)";     
   
      $s_percent=int(($hash7{$f}{$seq1_S})*100/$N);
      print OUT ";$seq1_S($s_percent)";     
            
    print OUT "\n";
  
}

print "done!";


my $t1 = gettimeofday;
my $elapsed = $t1-$t0;
print " time: $elapsed seconds \n";
