#!/usr/bin/perl -w
#use strict;
use lib "/u/min/jack/lib";
use Bio::SearchIO;
# this program is used to process blast results against EcoliVec database and to identify contaminant ESTs.
# Jack Min @copy right, on Feb 18, 2002. Re-written in Bioperl during Bioperl Bootcamp at U. Montreal
#
# Program: bvpFinder_bio.pl  Version 1.0
#
# ######## update record ########
# Date: June 23, 04  This program is re-written using Bioperl modules as a Bioperl Bootcamp project 
#
# Criteria used:
#	identity: >= 95%;
#	Query sequence coverage: >= 80
#
# when setting up BLASTN, E value is set at 1e-10.
#############################
# Input:  blastn against EcoliVec DB
# Output: file 1: SeqID, e-value, coverage, blast_info
# 	  file 2: contaminant seqIDs
#	  file 3: good seqIDs
#	  file 4 bs ids
#	  file 5 clean EST ids (short format)
#	  file 6 all hit details
# 	  file 7 excluded hits (hits, but filtered out)
#########################################################

########################## Body of Code ######################################
# the number of argument must be 1.

if ($#ARGV != 0)
{
	print "usage: bvpFinder_bio.pl  BLASTN\n";
	print " Input files:\n";
	print "		BLASTN - output from BLASTN\n";
	print " Notes: 5 files generated automatically:\n";
	print "		(1) Details (.details)\n";
	print "		(2) contaminant EST ids (.bad)\n";
	print "		(3) clean EST ids (.good) \n";
#	print "		(4) contaminant EST ids (short format, .bs)\n";
#	print "		(5) clean EST ids (short format, .gs)\n";
	print "		(4) all hit details (.all)\n";
	print "		(5) excluded hit (.exc)\n";
	exit;
}



my $blastN = $ARGV[0];

my $f1 = $blastN.".details";
my $f2 = $blastN.".bad";
my $f3 = $blastN.".good";
# my $f4 = $blastN.".bs";
# my $f5 = $blastN.".gs";
my $f6 = $blastN.".all";
my $f7 = $blastN.".exc";

print "Please input your sequence information and other comments:\n";
my $comment = <STDIN>;
my $time = localtime();
$comment = $comment."\n".$time." Jack Min, Fungal Genomics Centre, Concordia University\n"."Contaminants are identified by bvpFinder program.\n";
$comment = $comment."bvp: bacteria, vector, and plasmid.\nIf you have questions, please do not hesitate to contact Jack at jack\@gene.concordia.ca\n";
print "In progress ...\n";
########### Parameters #########
my $query = "" ;  #query sequece: contig number or clone number
my $hitStatus = 0; # default as no hit
my $hitInfo ="";  #the first hit infor
my $hit = 0;#false
my $qLength = 0; #query length
my $aLength = 0; #aligned length
my $coverage = 0; #coverage percentage
#########################################
my $noHitCount = 0;
my $queryCount = 0;

my @detail = ();
my @badID = ();
my @goodID = ();
my @all = (); #all hit details
my @exclude = ();
my $count =0; #to indicate the first ">"
my $input = new Bio::SearchIO(-format => 'blast', -file => $blastN, -report_type => 'blastn');

while (my $result = $input->next_result)
{
	$query = $result->query_name;
	$queryCount++;
#	$print $queryCount."\n";
	$qLength = $result->query_length;
	my $hitNum = $result->num_hits;
	if ($hitNum == 0)
	{
		push (@goodID, $query."\n");
		$noHitCount++;
	}	
	else
	{
		while (my $hit = $result->next_hit)
		{
			while (my $hsp = $hit->next_hsp)
			{
				#$aLength = $hsp->query->length;?
				$aLength = $hsp->length($query);
				$expect = lc ($hsp->evalue);
				$identity = ($hsp->frac_identical)*100;
				$coverage = $aLength/$qLength * 100;
				my @tmp = split(/\./, $coverage);
				$coverage = $tmp[0];
				$hitInfo = $hit->description;
if ($query =~ /ApulSEQ3865/)
{				
print "query: ".$query."\t".$expect."\t".$coverage."\t".$aLength."\t".$hitInfo."\n";
}
				push (@all, $query."\t".$expect."\t".$coverage."\t".$hitInfo."\n");
				if ($coverage > 80 && $identity >90)
				{
					push (@detail, $query."\t".$expect."\t".$coverage."\t".$hitInfo."\n");
					push (@badID, $query."\n");
				}

				elsif ($coverage > 60 && $identity >95 && $qLength > 150)
				{
					push (@detail, $query."\t".$expect."\t".$coverage."\t".$hitInfo."\n");
					push (@badID, $query."\n");
				}
				elsif ($coverage > 45 && $identity >97 && $qLength > 200)
				{
					push (@detail, $query."\t".$expect."\t".$coverage."\t".$hitInfo."\n");
					push (@badID, $query."\n");
				}
				else
				{
					push (@goodID, $query."\n");
					push (@exclude, $query."\t".$expect."\t".$coverage."\t".$hitInfo."\n");
				}
				last;
			}# end of while hsp
			last;
		  } # end of while hit
	} # end of else
} #end of while result

open (OUT, ">$f1") || die ("Can not open file $f1 for writing!\n\n");

print OUT ($comment."\n");
print OUT ("SUMMARY OF THE REPORT:\n");
print OUT ("total number of sequences: ".$queryCount."\n");
print OUT ("number of no hit:".$noHitCount."\n");
$hitNum = $queryCount - $noHitCount;
print OUT ("number of sequences with hit: ".$hitNum."\n");
print OUT ("good sequence: ".@goodID."\n");
print OUT ("contaminants: ".@badID."\n\n");
print OUT ("query"."\t"."E-value\t"."Coverage%\t"."Hit infor\n");
print OUT (@detail);
close (OUT);


open (OUT, ">$f2") || die ("Can not open file $f2 for writing!\n\n");
print OUT (@badID);
close OUT;

open (OUT, ">$f3") || die ("Cannot open file $f3 for writing!\n");
print OUT (@goodID);
close OUT;

open (OUT, ">$f6") || die ("Can not open file $f6 for writing!\n\n");

print OUT ($comment."\n");
print OUT ("SUMMARY OF THE REPORT:\n");
print OUT ("total number of sequences: ".$queryCount."\n");
print OUT ("number of no hit:".$noHitCount."\n");
$hitNum = $queryCount - $noHitCount;
print OUT ("number of sequences with hit: ".$hitNum."\n\n");
print OUT ("query"."\t"."E-value\t"."Coverage%\t"."Hit infor\n");
print OUT (@all);
close (OUT);

open (OUT, ">$f7") || die ("Cannot open file $f7 for writing!\n");
print OUT (@exclude);
close OUT;

print "The work is done!!\n";
exit 0;
