#! /usr/local/bin/perl -w

if ($#ARGV < 2)	# If you're given < 2 arguments
{
	print STDERR "\nUSAGE:  $0 BLAT_file sorted_BLAT_results\n";
	print STDERR "Sorts a BLAT file by choosing the hit with the best score for each query.\n";
	print STDERR "Score = 3 * matches - mismatches\n";
	print STDERR "In the case of >1 hits with the same score, all are selected.\n";
	exit(0);
}

$blat_output = $ARGV[0];
$sorted_BLAT_results = $ARGV[1];

$best_score = $last_query = 0;
$lineNum = 0;

# Open temp file (that will be made non-redundant at the end)
open (SORTED_BLAT, ">$sorted_BLAT_results.tmp") || die "cannot open $sorted_BLAT_results for writing: $!";

open (RAW_BLAT, $blat_output) || die "cannot open $blat_output for reading: $!";
while(<RAW_BLAT>)
{
	chomp($_);
	
	if (/^\d/)	# Skip header lines
	{
		@fields = split(/\t/, $_);
		$line = $_;
		$query = $fields[9];
		$score = 3 * $fields[0] - $fields[1];

		# If we're looking at data for the same query sequence again
		if ($query eq $last_query)
		{
			if ($score > $best_score)
			{
				@best_fields = @fields;
				$best_score = $score;

				# Make format like usual BLAT "psl" output
				$best_line = $line;
			}
			elsif ($score == $best_score)	# Two hits with the same score
			{
				$best_line .= "\n$line";	# Add this hit to the other
			}
		}
		else	# If we're looking at a query sequence for the first time
		{
			if ($best_line)	# Have a best lines for the previous query
			{
				# Print the best line(s) for the previous query
			
				print SORTED_BLAT "$best_line\n";
			}
			@best_fields = @fields;
			$best_score = $score;
			$best_line = $line;
		}
		$last_query = $query;
	}
	else
	{
		if ($lineNum < 5)	# header lines
		{
			if (! $lineNum)	# first header line: Should say "psLayout version 3"
			{
				print SORTED_BLAT "$_\tSorted BLAT output - best hit(s) for each query sequence\n";
			}
			else
			{
				print SORTED_BLAT "$_\n";
			}
		}
	}
	$lineNum++;
}

# Get the last one
print SORTED_BLAT "$best_line\n";

close (SORTED_BLAT);

# Optional: Start with $sorted_BLAT_results.tmp and print non-redundant lines to $sorted_BLAT_results

makeNonRedundant();

#########################  Subroutines  #########################

sub makeNonRedundant
{
	open (SORTED_BLAT, "$sorted_BLAT_results.tmp") || die "cannot open $sorted_BLAT_results.tmp for reading: $!";	
	open (FINAL_SORTED_BLAT, ">$sorted_BLAT_results") || die "cannot open $sorted_BLAT_results for writing: $!";	
	while (<SORTED_BLAT>)
	{		
		if (! $gotThisLine{$_})
		{
			print FINAL_SORTED_BLAT $_;
		}
		else
		{
			# print STDERR "Already got this hit:\n$_";
			$numRemoved++;
		}
		$gotThisLine{$_} = 1;
	}
	
	print STDERR "Sorted and removed $numRemoved redundancies\nFinal non-redundant sorted output is $sorted_BLAT_results\n";
	
	unlink("$sorted_BLAT_results.tmp");
}
