#!/usr/bin/perl

# this program is used to process blast results against EcoliVec database and to identify contaminant ESTs.
# Jack Min @copy right, on Feb 18, 2002.
#
# Program: bvpFinder.pl  Version 1.0
#
# ######## update record ########
#
# 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
#########################################################

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

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



$blastN = $ARGV[0];

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

print "Please input your sequence information and other comments:\n";
$comment = <STDIN>;
$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";
$midLine = 0; #the line below the "Query:" to show identity and similarity

########### Parameters #########
$query = "" ;  #query sequece: contig number or clone number
$hitStatus = 0; # default as no hit
$hitInfo ="";  #the first hit infor
$frame =""; #used for translation
$head = 0; #false
$hit = 0;#false
$qLength = 0; #query length
$aLength = 0; #aligned length
$coverage = 0; #coverage percentage
#########################################
$noHitCount = 0;
$queryCount = 0;

@detail = ();
@badID = ();
@goodID = ();
@bsID = (); #short format of ID
@gsID = (); #short format of ID
@all = (); #all hit details
@exclude = ();
$count =0; #to indicate the first ">"

open (IN, $blastN) || die ("Cannot open file $blastN\n\n");
$offSet = tell(IN);

while ($record = get_next_record(IN))
{
	@recordA = split(/\n/, $record);
	foreach $BlastContent (@recordA) #foreach start
	{#foreachloop #1
	chomp ($BlastContent);

	if ($BlastContent =~ /^Query= /)

	{
		@array1 = split (" ", $BlastContent);
		$hFlag = 1;
		$query = $array1[1];
		@qA = split(/\_/, $query);
		$sQuery = $qA[0];
                $queryCount++;
		print $queryCount."\n";
  		next;
	}

	elsif ($BlastContent =~ /No hits found/)
	{
	        push (@goodID, $query."\n"); #store the no hit query
             	push (@gsID, $sQuery."\n");
		$noHitCount++;
		$query = "";
                next;
	}

	elsif ($BlastContent =~ /letters\)/)
	{
		@tmp = split(/\(/, $BlastContent);
		@tmp1 = split(/ /, $tmp[1]);
		$qLength = $tmp1[0];
	#	print $query."\t",$qLength."\n";
                next;
	}
       elsif ($BlastContent =~ /^Sequences/)
	{
		$hit = 1; # there is a hit at least
		next;
	}



       elsif ($BlastContent =~ /^>/)
       {
         if ($hit ==1)
	 {
             $firstHit = 1;
             $hitInfo = $BlastContent;
             $position1 = index ($BlastContent, " ");
	     $sName = substr ($BlastContent, $position1+1);
             $sName = lc($sName);
	  }
        }

        else
        { # else - OUT1
	    if ($hit ==1 && $firstHit ==1)
	    {# if$hit and $firstHit both true)

		if ($BlastContent =~ /Expect/)
		{
                        $hitInfo = $hitInfo.$BlastContent;
			@expectLine = split(/,/, $BlastContent);
			$eString = $expectLine[1];
			@subString = split (/ /, $eString);
			$expect = $subString[3];
			$expect = lc($expect);
		}

		elsif ($BlastContent =~ /Identities =/)
                {
			@identity = split(/[()]/, $BlastContent);
                        $identity = $identity[1] * 1;
			@tL = split(/\//, $BlastContent);
			@tL2 = split(/ /, $tL[0]);
			$aLength = pop (@tL2);
			#print "aLength: ".$aLength."\n";
			$coverage = $aLength/$qLength * 100;
			@t = split(/\./, $coverage);
			$coverage = $t[0];
			#print "coverage: ".$coverage."\n";
##########for debug
#print $BlastContent."\n";
#print "Query: ".$query."\t query length:".$qLength."\t alignment length:".$aLength."\t"."coverage: ".$coverage."\n";

			push (@all, $sQuery."\t".$expect."\t".$coverage."\t".$hitInfo.$BlastContent."\n");

			if ($coverage > 80 && $identity >90)
			{
			push (@detail, $sQuery."\t".$expect."\t".$coverage."\t".$hitInfo.$BlastContent."\n");
			push (@badID, $query."\n");
			push (@bsID, $sQuery."\n");
			$hit = 0;
			$firstHit = 0;
			$hitInfo = '';
			next;
			}

			elsif ($coverage > 60 && $identity >95 && $qLength > 150)
			{
			push (@detail, $sQuery."\t".$expect."\t".$coverage."\t".$hitInfo.$BlastContent."\n");
			push (@badID, $query."\n");
			push (@bsID, $sQuery."\n");
			$hit = 0;
			$firstHit = 0;
			$hitInfo = '';
			next;
			}
			elsif ($coverage > 40 && $identity >97 && $qLength > 200)
			{
			push (@detail, $sQuery."\t".$expect."\t".$coverage."\t".$hitInfo.$BlastContent."\n");
			push (@badID, $query."\n");
			push (@bsID, $sQuery."\n");
			$hit = 0;
			$firstHit = 0;
			$hitInfo = '';
			next;
			}
			else
			{
				push (@goodID, $query."\n");
				push (@gsID, $sQuery."\n");
				push (@exclude, $sQuery."\t".$expect."\t".$coverage."\t".$hitInfo.$BlastContent."\n");
				$firstHit = 0;
                        	$hit = 0;
				$hitInfo = 0;
			}
		}

  	 	else
		{
			$hitInfo = $hitInfo.$BlastContent;
		}
	  } # end of #both $hit and $firstHit are true
	} # end of Out-1
	} #end of for loop
$offSet = tell(IN);
}#end of while

close (IN);

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, ">$f4") || die ("Cannot open file $f4 for writing!\n");
print OUT (@bsID);
close OUT;

open (OUT, ">$f5") || die ("Cannot open file $f5 for writing!\n");
print OUT (@gsID);
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;

################
sub get_next_record
{
	my($fh) = @_;
	use strict;
	my($offSet);
	my($record);
	my($seperator) = $/;
	$/ = "BLASTN";
	$record =<$fh>;
	$/ =$seperator;
	return $record;
}



