#!/usr/bin/perl
use strict;

my $ditag_size;
my $coesive;
my $input = "";
my $output = "";
my $routput = "";
my $tcount = 0;
my $rejected = 0;
my $minimum_size = 17;
my $maximum_size = 18;
my $remove = 0;
my @mono = ();
my %hash;

my %count;

my $version = 0.35;
my $vname = "repeated_ditag_evaluator.pl";

if (!$ARGV[0] || $ARGV[0] eq "-h") {
	print "\n\t-c <integer> \t:: length of the cohesive ends[default:0]\n",
	"\t-min <integer> \t:: set minimum size of the tags [default:$minimum_size]\n",
	"\t-max <integer> \t:: set maximum size of the tags [default:$maximum_size]\n",
	"\t-R       \t:: Ignore the statistical analysis and remove all repeated ditags\n",
	"\t-i <FILE> \t:: input file name [default:STDIN]\n",
	"\t-o <FILE> \t:: output file name [default:STDOUT]\n",
	"\t-r <FILE> \t:: report file name [default:STDERR]\n",
	"\t-h      \t:: display help information\n",
	"\t-v      \t:: display current version\n",
	"\n\t$vname Version $version\n\n";
	exit;
}

if ($ARGV[0] eq "-v") {
	print ("\n\t$vname Version $version\n\n");
	exit;
}


#
#first read command line to collect arguments
#
while (my $arg=shift(@ARGV)) {
	#
	#output file
	#
	if ($arg eq "-o") {
		if ($arg = shift(@ARGV)) {
			$output = $arg;
		}
		else {
			print STDERR "OUTPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}

	#
	#input file
	#
	elsif ($arg eq "-i") {
		if ($arg = shift(@ARGV)) {
			$input = $arg;
		}
		else {
			print STDERR "INPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}

	#
	#output file
	#
	if ($arg eq "-r") {
		if ($arg = shift(@ARGV)) {
			$routput = $arg;
		}
		else {
			print STDERR "REPORT OUTPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}

	#
	#sets coesive end analisys 
	#
	elsif ($arg eq "-c") {
		if ($arg = shift(@ARGV)) {
			$coesive = $arg;
			if ($coesive < 0) {
				print STDERR "COESIVE END SIZE INVALID!\n";
				exit;
			}
		}
		else {
			$coesive = 0;
		}
	}
	#
	#set the minimum size of the tag
	#
	elsif ($arg eq "-min") {
		if ($arg = shift(@ARGV)) {
			$minimum_size = $arg;
		}
		else {
			print STDERR "TAG MINIMUM SIZE NOT SPECIFIED!\n";
			exit;
		}
	}
	#
	#set the maximum size of the tag
	#
	elsif ($arg eq "-max") {
		if ($arg = shift(@ARGV)) {
			$maximum_size = $arg;
		}
		else {
			print STDERR "TAG MAXIMUM SIZE NOT SPECIFIED!\n";
			exit;
		}
	}
	#
	#Remove all repeated ditags
	#
	elsif ($arg eq "-R") {
		$remove = 1;
	}
}


my %ditag_count = ();
my %tag_count = ();

if ($input) {
	open INPUT, "<$input" or die "COULD NOT OPEN $input : $!\n";
	$input = "<INPUT>";
}
else {
	$input = "<>";
}

if ($output) {
	open OUTPUT, ">$output" or die "COULD NOT OPEN $output : $!\n";
}

if ($routput) {
	open ROUTPUT, ">$routput" or die "COULD NOT OPEN $routput : $!\n";
}

if ($remove) {
	my %uniq;
	while (my $ditag = eval($input)) {
		if (!defined $uniq{$ditag}) {
			&printout ($ditag, "\n");
			$uniq{$ditag} = 1;
		}
	}
	exit;
}

while (my $ditag = eval($input)) {

	chomp ($ditag);
	&clean ($ditag);
	#
	#get the ditag size
	#
	my $dt_size = length ($ditag);
	my $tag1;
	my $tag2;

	my $cut_th = &max ($dt_size + $coesive - $maximum_size, $minimum_size);
	$tag1 = substr ($ditag, 0, $cut_th);
	$tag1 = &cleanX ($tag1);
	$tag2 = &reverseComplement (substr ($ditag, -$cut_th));
	$tag2 = &cleanX ($tag2);

	if (!$tag1 && $tag2) {
		push (@mono, $tag2);
		next;
	}
	elsif ($tag1 && !$tag2) {
		push (@mono, $tag1);
		next;
	}
	elsif (!$tag1 && !$tag2) {
		next;
	}

	$hash{$ditag} = "$tag1#$tag2";

	my $tail1 = &getTail ($tag1);
	my $tail2 = &getTail ($tag2);

	if (defined $count{$tail1}) {
		$count{$tail1}++;
	}
	else {
		$count{$tail1} = 1;
	}

	if (defined $count{$tail2}) {
		$count{$tail2}++;
	}
	else {
		$count{$tail2} = 1;
	}
	$tcount += 2;

	if (exists($ditag_count{$ditag})) {
		$ditag_count{$ditag}++;
	}
	else {
		$ditag_count{$ditag} = 1;
	}

	if (exists($tag_count{$tag1})) {
		$tag_count{$tag1} += 1;
	}
	else {
		$tag_count{$tag1} = 1;
	}

	if (exists($tag_count{$tag2})) {
		$tag_count{$tag2} += 1;
	}
	else {
		$tag_count{$tag2} = 1;
	}
}
close INPUT;


my $uniq = 0;
foreach my $ditag (keys(%ditag_count)) {

	#&printerr ("ditag = $ditag\n");

	$uniq++;

	(my $tag1, my $tag2) = split (/#/, $hash{$ditag});

	my $tail1 = &getTail ($tag1);
	my $tail2 = &getTail ($tag2);

	my $freq_tag1 = $tag_count{$tag1}/$count{$tail1};
	my $freq_tag2 = $tag_count{$tag2}/$count{$tail2}; 

#	my $expect = $freq_tag1 * $freq_tag2 * $count{$tail1} * $count{$tail2};

	my $p = ($tag_count{$tag1} / $count{$tail1}) * ($tag_count{$tag2} / $count{$tail2});

	my $var = $tcount * $p * (1 - $p);

	#
	#expect = prob of selecting 2 tags  times the numer of ditags in the bins
	#
	my $expect = $p * ($count{$tail1} + $count{$tail2}) / 2;

	if ($tag1 eq $tag2) {
		$expect /= 2;
	}

	#if ($tail1  eq $tail2){
	#	$expect /= 2;
	#}


	my $xc = (1.96 * $var) / sqrt ($count{$tail1} + $count{$tail2}) + $expect; 

	#
	#Print the data, rebuild the Ditag to be print correctly
	#
	&printerr ("Ditag: ", $tag1, &reverseComplement (substr ($tag2, 0, -2)) ,"= $ditag_count{$ditag}\n",
		"\tTag1: $tag1 => $tag_count{$tag1};\tp = ",$tag_count{$tag1}/$count{$tail1}, "\n",
		"\tTag2: $tag2 => $tag_count{$tag2};\tp = ",$tag_count{$tag2}/$count{$tail2}, "\n\n",
		"\t\tp = $p\n",
		"\t\tvar = $var\n",
		"\t\tX\\ = ", $expect, "\n",
		"\tXc = $xc\n");
	my $margin = $xc - $ditag_count{$ditag};

	my $result;
	if ($ditag_count{$ditag} == 1) {
		$result = "ACCEPTED by Small Sample";
	}
	elsif (($xc >= $ditag_count{$ditag})) {
		$result = "ACCEPTED by $margin";
	}
	else {
		$result = "**REJECTED** by $margin";
		$rejected++;
	}

	&printerr ("\n\t\t$result\n\n");
	my $extra = int (&max (0, -$margin));
	for (my $i = 0; $i < $ditag_count{$ditag} - $extra; $i++) {
		&printout ($ditag, "\n"); 
	}
}

foreach my $mono (@mono) {
	&printout ($mono);
}
&printerr ("$rejected REJECTED of ", $uniq, "\n\t", 100 * $rejected / $uniq, "\n");

sub getTail {
	my $tag = shift();
	my $tail = 'X';
	if ($coesive) {
		$tail = substr ($tag, -$coesive);
	}
	return $tail;
}
#
#return the reverse complement of the DNA sequence
#
sub reverseComplement {
	my $string = shift();
	my $base;
	my @split;
	@split = split (//, $string);
	$string = "";
	while (@split) {
		$base = pop (@split);
		if ($base eq "a") {
			$string .= "t";
		}
		elsif ($base eq "t") {
			$string .= "a";
		}
		elsif ($base eq "c") {
			$string .= "g";
		}
		else {
			$string .= "c";
		}
	}

	return $string;
}



sub printout {
	if ($output) {
		while (my $p = shift()) {
			print OUTPUT $p;
		}
	}
	else {
		while (my $p = shift()) {
			print $p;
		}
	}
}


sub printerr {
	if ($routput) {
		while (my $p = shift()) {
			print ROUTPUT $p;
		}
	}
	else {
		while (my $p = shift()) {
			print STDERR $p;
		}
	}
}


sub max {
	my $a = shift;
	my $b = shift;
	if ($a > $b) {
		return $a;
	}
	else {
		return $b;
	}
}


sub clean {
	my $string = shift;

	$string =~ s/^\s+//;
	$string =~ s/\s+$//;

	return $string;
}

sub cleanX {
	my $string = shift;

	$string =~ s/^X+//;
	$string =~ s/X+$//;

	return $string;
}

