#!/usr/bin/perl -w

use strict;
use lib "/usr2/people/khuang/bin/Genomics/Perl/modules";
use GenomicsUtils;
use Getopt::Long;
use Gene;

my $helpStr = qq{


  Description:
    
    This script parses any BLAST tab delimited output (-m 8 option).
    
    
  Details:
    
    parseBlastTab.pl <blast db> (options)

	
    <output option>:
    
           -detail: to get more info on the subject hits
	   -locus: to convert query to locus ID


           -best (hit)
	or -all

    <input option>:
    
           -i <blast results file name>
        or -d <dir> and -suffix <suffix>, to process *.suffix in dir
     
    <hit criteria>:
           -hitInfo, to just print info of hits
           -id <minimum identity percentage>, default=0
           -s <minimum score>, default=0
           -e <maximum evalue>, default=10
	   -a <minimum alignment coverage>, default=0
     
     

   Exameple:
	
        parseBlastTab.pl <blastdb> -best -i temp.blastResults -id 90 -e 1e-5 
	
	which parses temp.blast
	and returns the single best hit of each query 
	with identity >=90% AND evalue <= 1e-5

	The output is default to STDOUT

};



########
# MAIN #
########

unless (@ARGV>=3){die "$helpStr";}

my $db = $ARGV[0];
my %queries = ();


my ($best, $all, $detail, $i, $d, $suffix, $id, $s, $e , $locus, $hitInfo, $aFilter);

GetOptions( 
	    '-detail' => \$detail,
	    '-best' => \$best,
	    '-all' => \$all,
	    '-i=s' => \$i,
	    '-d=s' => \$d,
	    '-suffix=s' => \$suffix,
	    '-id=i' => \$id,
	    '-s=i' => \$s,
	    '-e=f' => \$e,
	    '-locus' => \$locus,
	    '-hitInfo' => \$hitInfo,
	    '-a=i' => \$aFilter

);

$id = 0 unless defined($id);
$s = 0 unless defined($s);
#print $e;
$e = 10 unless defined($e);
$aFilter = 0 unless defined($aFilter);
$all = 1 unless defined($best);

#print $aFilter."\t".$detail;

die "-detail is required if -a is selected" if (($aFilter)&&!($detail));

my $fastacmd = "/usr2/people/khuang/apps/blast/fastacmd";

################### Params #######################

if (($i&&$d)||($best&&$all)){die "$helpStr";}
if($d){
    die "$helpStr" unless defined($suffix);
}
if ($hitInfo){
    die "$helpStr" unless defined($detail);
}

################## Get Input Files ##################
my @blastFiles = ();

if ($i){ #if only one file
    @blastFiles = ($i);
}
elsif($d){
    opendir(DIR, $d) or die "Can't opendir $d: $!\n";
    @blastFiles = readdir(DIR);
    closedir(DIR);
}

GenomicsUtils::connect("lut");

foreach (@blastFiles){
    my $file = $_;
    if (defined($d)){ 
	next unless /\.$suffix/ ;
	$file = $d."/".$file;
    }
    
    ############ Read each input files#################

    my %best = ();
    my $name = "cog.rpsblast" if $file =~ /cog.rpsblast/;
    $name = "pfam.rpsblast" if $file =~ /pfam.rpsblast/;
    $name = "smart.rpsblast" if $file =~ /smart.rpsblast/;
    $name = "kegg.blastp" if $file =~ /kegg|ec/;
    $name = "phage" if $file =~ /phage/;
    $name = "transposase" if $file =~ /transposase/;
    $name = "parsed" unless defined($name);
    open AUTOANNO, ">>AutoAnno.$name.txt" or die "AutoAnno.$name.txt: $!\n";
    open HITINFO, ">>HitInfo.$name.txt" or die "HitInfo.$name.txt: $!\n";
    open BLAST, "$file" or die "$file: $!\n";
    while (<BLAST>){
	s/\n//;
	my ($query, $subject, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score) = split /\t/;
	my $geneObj = Gene::new(locusId=>$query) or die "Can't create a gene object for $query";
	my $geneL = length($geneObj->protein());

	my $subjectID = $subject;
	my $subDesc = $subject;
	my $cdd = "";

	my $database = "";
	$database = "PFAM" if $db =~ /Pfam/i;
	$database = "SMART" if $db =~ /Smart/i;
	$database = "COG" if $db =~ /Cog/i;
	$database = "KEGG" if $db =~ /kegg/i;
	$database = "PHAGE" if $db =~ /phage/i;
	$database = "TRANSPOSASE" if $db =~ /transposase/i;

	next if $score < $s || $identity < $id || $evalue > $e;
	my $sSeq = "";
	my $sLength = "";
	if ($detail){

	    my @temp = `$fastacmd -s "$subject" -d $db`;
	    ( $subject = shift @temp ) =~ s/>|\n//g;
	    $subject =~ s/lcl\|//;
	    my @temp2 = ();
	    if ($database eq "KEGG"){
		($subjectID, @temp2) = split / /, $subject;
		$subDesc = join ",", @temp2;
	    }
	    elsif ( ($database eq "PHAGE")||($database eq "TRANSPOSASE") ){
		$subDesc = $subject;
		$subject = $subjectID;

	    }
	    else{
		($subjectID, @temp2) = split /,/, $subject;
		($cdd, $subjectID) = split / /, $subjectID; 
		if ($subjectID =~ /COG(\d+)/){
		    ($subjectID = $1) =~ s/^0*//;
		}
		$cdd =~ s/gnl\|CDD\|//;
		$subDesc = join ",", @temp2;
	    }
	    chomp @temp;
	    $sLength = length(join "", @temp);
	}
	
	if ($locus){
	    
	    $query =~ s/gi\|(\d+)\|*.*/$1/;
	    my $tempR = GenomicsUtils::query(qq{SELECT locusId FROM Synonym WHERE name="$query"});
	    $query = $tempR->[0][0];
	   
	}

	if ($all){
	    if ($hitInfo){
		#print $query."\t".$geneL."==\n";
		#print $sLength."====\n";
		print HITINFO join "\t", $subjectID, $cdd, $sLength, $subDesc."\n";     
		if ($aFilter){
		    if ( (100*($alignment/$geneL) >= $aFilter) && (100*($alignment/$sLength) >= $aFilter) ){
			print AUTOANNO join "\t", $database, $query, $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n";
		    }
		}
		else{    
		    
		    print AUTOANNO join "\t", $database, $query, "1", $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n";
		}
	    }
	}

	elsif ($best){
	    $best{$query}++;
	    
	    if ($hitInfo){
		print HITINFO join "\t", $subjectID, $cdd, $sLength, $subDesc."\n" if $best{$query} == 1;     
		if ($aFilter){
		    if ( (100*($alignment/length($geneL)) >= $aFilter)&&( 100*($alignment/$sLength) >= $aFilter )){
			print AUTOANNO join "\t", $database, $query, $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n" if $best{$query}==1;
		    }
		}
		else{    
		    
		    print AUTOANNO join "\t", $database, $query, "1", $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n" if $best{$query}==1;
		}
	    }
	   
	    else{ 
		if ($aFilter){
		    if ( (100*($alignment/length($geneL)) >= $aFilter)&&( 100*($alignment/$sLength) >= $aFilter )){
			print AUTOANNO join "\t", $database, $query, $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n" if $best{$query}==1;;
		    }
		}
		else {
		    print AUTOANNO join "\t", $database, $query, "1", $subjectID, $identity, $alignment, $mismatch, $gap, $qBegin, $qEnd, $sBegin, $sEnd, $evalue, $score."\n" if $best{$query}==1;
		    
		}
	    }
	    
	}

    }
    close BLAST;
    


}


GenomicsUtils::disconnect();
exit;
