#!/usr/bin/perl

##############################################################
#  Written by Justin Elser 2/23/09                           #
#   Dept of BPP, Oregon State University                     #
#   elserj@science.oregonstate.edu                           #
#  Version 0.9                                               #
#   this program is meant to be part of a pipeline including #
#      interact.pl, therefore, it is written assuming the    #
#      input files are from that program (interact.pl)       #
#                                                            #
#  still to do-----                                          #
#       not all stats are computed yet                       #
##############################################################

use strict;
use warnings;

if ($#ARGV != 1) {
	print "usage :   stats_interact_Ath.pl prefix_for_in_files output_file\n";
	exit;
}
open(main_output_file, ">$ARGV[1]") || die "Error: file '$ARGV[1]' can not be opened\n";
my $out_file = "$ARGV[1]";

printf main_output_file "Num of lines\tUnique gene count\tMax Gene IDs\t# max Gene\t# min genes\n";


	open (in_file_prefix, "$ARGV[0]") || die "Error: file '$ARGV[0]' can not be opened\n";
	
	open (output_file, ">$ARGV[1]" . "_num") || die "Error: file '$ARGV[1]' can not be opened\n";

	
	# start doing the stats
	#  earlier versions of the script did this using system calls to cat, grep, uniq, etc...  
	#  implementing directly in perl has increased the performance of the program greatly
	my $count = 0;  # total number of unique interactions
	my %seen = ();  # hash that stores how many times each gene has been seen
	while (my $line = <in_file_prefix>) {
		$count++;
		my ($gene_1,$inter,$gene_2) = split ("\t",$line);
		chomp $gene_2;
		if ($gene_1 eq $gene_2) {
			$seen{"$gene_1"}++;
		}
		else {
			$seen{"$gene_1"}++;
			$seen{"$gene_2"}++;
		}
	}
	my $uniq_gene_count = keys %seen; # gives the count of unique genes
	
	# this block populates the max and min lists.  max_list holds the gene id(s) 
	#   that comes up the most.  min_list holds the gene ids that only show up once
	my $max = undef;
	my @max_list = undef;
	my @min_list = ();
	while( my ($key, $value) = each %seen ) {
		if(!defined $max) {
			$max = $value;
			@max_list = ($key);
		} elsif ($seen{$key} == $max) {
			push @max_list, $key;
		} elsif ($seen{$key} > $max) {
			$max = $value;
			@max_list = ($key);
		}
		if($value == 1) {
			push @min_list, $key;
		}
	}
	
	my $min_count = @min_list;
	
	printf main_output_file "$count\t$uniq_gene_count\t@max_list\t$max\t$min_count\n";


	
	my %hash = ();
	
	# this block will calculate the number of genes that have occurred a given number of times
	#  for example, if there are 10 genes that occur 5 times, the key of this hash is 5, value is 10
	while( my ($key, $value) = each %seen ) {
		$hash{"$value"}++;
	}
	
	foreach my $key (sort by_number keys %hash) { # see the sub description below
		my $value = $hash{$key};
		print output_file "$key\t$value\n";
	}
	

	close(in_file_prefix);
	close(output_file);
	


close(main_output_file);

# the perl sort() sorts by ASCII value, not numerically, this subroutine fixes the sort
sub by_number {
	if ($a < $b) {-1} elsif ($a > $b) {1} else {0}
}
