#!/usr/bin/perl -w              

$main_dirname = "";

$title = $ARGV[0];

use Bio::AlignIO;
use Bio::SeqIO;
use Bio::SearchIO;
use Bio::Seq;

$in_dirname = $main_dirname . "blast_out/";
$out_dirname = $main_dirname ."blast_result/";

$bl_seq_name = $in_dirname . $title . "_blast.txt";
$unknown_filename = $out_dirname . $title. "_unknown.xls";
$pai1_filename = $out_dirname . $title . "_pai1_match.xls";
$reject_blast_filename = $out_dirname . $title. "_blast_reject.xls";
$display_filename = $out_dirname . $title. "_display.txt";
$count_filename = $out_dirname . "master_count_from_complete_blast.csv";
$seq_pos_filename = $out_dirname . $title . "_start_end_blast.xls";

if (-e $display_filename) {`rm $display_filename`}
if (-e $unknown_filename) {`rm $unknown_filename`}
if (-e $pai1_filename) {`rm $pai1_filename`}
if (-e $reject_blast_filename) {`rm $reject_blast_filename`}

open(PAI,">$pai1_filename") || die "cannot open " . $pai1_filename;
open(REJECT_B,">$reject_blast_filename") || die "cannot open " . $reject_blast_filename;
print REJECT_B join(",","Seq Id","Query String","Consenses","Ref Start","Ref End","Seq Len","Percent","Seq Start","Seq End","Strand","# Gaps") . "\n";
open(DISPLAY,">>$display_filename") || die "cannot open " . $display_filename;
open(UNKNOWN,">$unknown_filename") || die "cannot open " . $unknown_filename;

$display = 0;
$seq_cnt = 0; 
$unknown_cnt = 0;
$pai_cnt = 0;
$blast_reject_cnt = 0;
$seq_cnt = 0;

$per_homo_value = 80;
$query_hit_length = 90;
undef(%seq_start);
undef(%seq_end);
undef(%seq_len);

$start_time = time();
$st_time = $start_time;
## use seq_Id for bl2seq format, use desc for blastxml format
print DISPLAY "start searchIO " . $bl_seq_name . " start " . $st_time . "\n";
$seq_aln = Bio::SearchIO->new(-format => 'blast', -file   => $bl_seq_name);    
while ($result = $seq_aln->next_result) {
    $seq_hits = $result->num_hits;
    $seq_id = $result->query_name;
    $seq_id =~ s/\@//;
    @tmp_id = split(/\s+/,$seq_id);
    $seq_id = $tmp_id[0];
#print "hits " . $result->num_hits . "\n";
    $seq_cnt ++;
    if ($seq_cnt % 100000 == 0) {
	$t_time = time();
#	print "time " . ($t_time - $st_time) . "\n";
#	print "Seq Cnt " . $seq_cnt . " file " . $bl_seq_name . "\n";
#	print "time\t" . ($t_time - $st_time) . "\n";
#	print "Seq Cnt\t" . $seq_cnt . "\tfile\t" . $bl_seq_name . "\n";
#	print "pai_cnt\t" . $pai_cnt . "\tunknown_cnt\t" . $unknown_cnt . "\tbl_reject_cnt\t" . $blast_reject_cnt . "\n";
	print DISPLAY "time\t" . ($t_time - $st_time) . "\n";
	print DISPLAY "Seq Cnt\t" . $seq_cnt . "\tfile\t" . $bl_seq_name . "\n";
	print DISPLAY"pai_cnt\t" . $pai_cnt . "\tunknown_cnt\t" . $unknown_cnt . "\tbl_reject_cnt\t" . $blast_reject_cnt . "\n";
#	$b=<STDIN>;
	$st_time = $t_time;
    }
    if ($seq_hits) {
	if ($seq_hits > 1) {
	    print DISPLAY "Multiple hits, " . $seq_hits . " in sequence Id "  . $desc . "\n";
#	    print "Multiple hits in " . $desc . "\n";
	}
	$hit = $result->next_hit; 
	$aln_p = $hit->next_hsp;
	$seq_ref_id = $hit->name;
	$query_string = $aln_p->query_string;
	$seq_length = $aln_p->length('query');
	$hit_length = $aln_p->length('hit');
	$seq_percent = $aln_p->percent_identity;
	$seq_start = $aln_p->start('query');
	$seq_end = $aln_p->end('query');
	$seq_gaps = $aln_p->gaps;
	$seq_ref_start = $aln_p->start('hit');
	$seq_ref_end = $aln_p->end('hit');
	$seq_consensus_string = $aln_p->hit_string;
	$strand = $aln_p->strand('hit');
	if ($display) {
	    print "seq_ref Id " . $seq_ref_id . "\n";
	    print "seq Id " . $seq_id . "\n";
	    print "length " . $seq_length . "\n";
	    print "hit_length " . $hit_length . "\n";
	    print "percent_identity " . $seq_percent . "\n";
	    print "consensus string " .$seq_consensus_string. "\n";
	    print "query_string " . $query_string . "\n";
	    print "gaps1 " . $seq_gaps . "\n";
	    print "start1 " . $seq_start . "\n";
	    print "end1 " . $seq_end . "\n";
	    print "start2 " . $seq_ref_start . "\n";
	    print "end2 " . $seq_ref_end . "\n";
	    print "seq_cnt " . $seq_cnt . "\n";
	    print "desc " . $desc . "\n";
	    print "per_homo_value " . $per_homo_value . "\n";
	    print "hit length " . $query_hit_length . "\n";
#		$b = <STDIN>;
	}
	if (exists($seq_len{$seq_length})) {
	    $seq_len{$seq_length}[0] ++;
	} else {
	    $seq_len{$seq_length}[0] = 1;
	}
	if (exists($seq_start{$seq_ref_start})) {
	    $seq_start{$seq_ref_start}[0] ++;
	} else {
	    $seq_start{$seq_ref_start}[0] = 1;
	}
	if (exists($seq_end{$seq_ref_end})) {
	    $seq_end{$seq_ref_end}[0] ++;
	} else {
	    $seq_end{$seq_ref_end}[0] = 1;
	}
	if ($seq_percent > $per_homo_value && $hit_length > $query_hit_length) {
	    if (exists($seq_len{$seq_length})) {
		$seq_len{$seq_length}[1] ++;
	    } else {
		$seq_len{$seq_length}[1] = 1;
	    }
	    if (exists($seq_start{$seq_ref_start})) {
		$seq_start{$seq_ref_start}[1] ++;
	    } else {
		$seq_start{$seq_ref_start}[1] = 1;
	    }
	    if (exists($seq_end{$seq_ref_end})) {
		$seq_end{$seq_ref_end}[1] ++;
	    } else {
		$seq_end{$seq_ref_end}[1] = 1;
	    }
		print PAI join("\t",$seq_id,$query_string,$seq_consensus_string,$seq_ref_start,$seq_ref_end,$seq_length,$seq_percent,$seq_start,$seq_end,$strand,$seq_gaps) . "\n";
	    $pai_cnt ++;
	} else {
	    if (exists($seq_len{$seq_length})) {
		$seq_len{$seq_length}[2] ++;
	    } else {
		$seq_len{$seq_length}[2] = 1;
	    }
	    if (exists($seq_start{$seq_ref_start})) {
		$seq_start{$seq_ref_start}[2] ++;
	    } else {
		$seq_start{$seq_ref_start}[2] = 1;
	    }
	    if (exists($seq_end{$seq_ref_end})) {
		$seq_end{$seq_ref_end}[2] ++;
	    } else {
		$seq_end{$seq_ref_end}[2] = 1;
	    } 
	    print REJECT_B  join("\t",$seq_id,$query_string,$seq_consensus_string,$seq_ref_start,$seq_ref_end,$seq_length,$seq_percent,$seq_start,$seq_end,$strand,$seq_gaps) . "\n";
	    $blast_reject_cnt ++;
	}
    } else { ##seq_hits 
#	print UNKNOWN $desc . "\n";
	print UNKNOWN $seq_id . "\n";
	$unknown_cnt ++;
    } ##seq_hits
}## while hits = result->next_hits

$outline = "Start Pos\tTotal\tFnd Start\tReject Start\n";
foreach $key ( sort {$a<=>$b} keys %seq_start) {
    if (!$seq_start{$key}[0]) {
	$seq_start{$key}[0] = 0;
    }
    if (!$seq_start{$key}[1]) {
	$seq_start{$key}[1] = 0;
    }
    if (!$seq_start{$key}[2]) {
	$seq_start{$key}[2] = 0;
    }
    $outline .= join("\t",$key,$seq_start{$key}[0],$seq_start{$key}[1],$seq_start{$key}[2]) . "\n";
}
$outline .= "\nEnd Pos\tTotal\tFnd End\tReject End\n";
foreach $key( sort {$a<=>$b} keys %seq_end) {
    if (!$seq_end{$key}[0]) {
	$seq_end{$key}[0] = 0;
    }
    if (!$seq_end{$key}[1]) {
	$seq_end{$key}[1] = 0;
    }
    if (!$seq_end{$key}[2]) {
	$seq_end{$key}[2] = 0;
    }
    $outline .= join("\t",$key,$seq_end{$key}[0],$seq_end{$key}[1],$seq_end{$key}[2]) . "\n";
}
$outline .= "\nSeq Length\tTotal\tFnd Len\tReject Len\n";
foreach $key( sort {$a<=>$b} keys %seq_len) {
    if (!$seq_len{$key}[0]) {
	$seq_len{$key}[0] = 0;
    }
    if (!$seq_len{$key}[1]) {
	$seq_len{$key}[1] = 0;
    }
    if (!$seq_len{$key}[2]) {
	$seq_len{$key}[2] = 0;
    }
    $outline .= join("\t",$key,$seq_len{$key}[0],$seq_len{$key}[1],$seq_len{$key}[2]) . "\n";
}

open(OUT,">$seq_pos_filename") || die "cannot open " . $seq_pos_filename;
print OUT $outline;
close(OUT);

close(PAI);
close(UNKNOWN);
close(REJECT_B);
$end_time = time();
#print "Total Time " . ($end_time - $start_time) . "\n";
print DISPLAY "\ntotal time\t" . ($end_time - $start_time) . "\n";
print DISPLAY "Seq Cnt\t" . $seq_cnt . "\n";
print DISPLAY "PAI Cnt\t" . $pai_cnt . "\tUnknown:\t". $unknown_cnt . "\tblast reject\t" . $blast_reject_cnt . "\n";
close(DISPLAY);
open(OUT,">>$count_filename") || die "cannot open " . $count_filename;
print OUT join("\t",$pai_cnt,$unknown_cnt,$blast_reject_cnt,$seq_cnt,($end_time - $start_time),$bl_seq_name). "\n";
close(OUT);


exit 0;
