#!/usr/bin/perl
#  Originally written by Justin Elser 2/2/09
#  Version 0.95 - completed 2/19/09
#  Version 0.96 - still working on
#  This file takes the two files that define the network and the orthologies and
#       outputs a predicted interactome
#
#   updated 2/19/09 to filter out duplicates A -> B = B -> A in the interactions
#
#  0.96 will add the ability to use a separate (optional) paralog list
use strict;
use warnings;

# print usage if other than three arguments to the script are found
if ($#ARGV < 2 || $#ARGV > 3) {
 print "usage: interact.pl Network_file Orthology_file Paralog_file(optional) Output_file\n";
 exit;
}

my $network_file_name = $ARGV[0];
my $ortho_file_name = $ARGV[1];

my $para_file_name = "";
my $out_file_name = "";
my $para = "";

if ($#ARGV ==3) {
	$para_file_name = $ARGV[2];
	$out_file_name = $ARGV[3];
	$para=1; # true
}else{
	$out_file_name = $ARGV[2];
	$para=0; #false
}

my %para_hash;

# read in the paralogs first so it is not in the loop, but the hash is available to the loop
if ($para) {
	
	open(paralog_file, "$para_file_name") || die "Error: file '$para_file_name' can not be opened\n";
	
	my $min_ident_para = 0; # set the threshold for paralog identities
	while(<paralog_file>) {
		my $entry = $_;
		# strip off newline characters
		$entry =~ s/\r//g;
		$entry =~ s/\n//g;
		# split the columns into 3 separate variables
		my ($gene_1, $gene_2, $ident) = split("\t", $entry);
		# skip if $ident is empty or Nan
		next if(!defined($ident));
		next if($ident =~ /D/);
		next if($ident eq "");
		next if($ident =~ /\%/);
		# change gene ids to all caps
		$gene_1 =~ tr/a-z/A-Z/;
		$gene_2 =~ tr/a-z/A-Z/;
		
		# put the genes in ASCII order to help remove dupes
		if($gene_1 gt $gene_2) {
			my $temp = $gene_1;
			$gene_1 = $gene_2;
			$gene_2 = $temp;
		}
		
		# only use paralogs with high confidence values
		if($ident>=$min_ident_para) {
			# most genes will show up many times, so push values for hash key
			if(defined($para_hash{$gene_1})) {
				my $oldhash = $para_hash{$gene_1};
				$oldhash->{'gene2'} = "$oldhash->{'gene2'}\t$gene_2";;
				$para_hash{$gene_1} = $oldhash;
				#print "$oldhash\n";
			}else{
				my %hash;
				$hash{'gene1'} = $gene_1;
				$hash{'gene2'} = $gene_2;
				$para_hash{$gene_1} = \%hash;
			}
		}
		
	}

	close(paralog_file);

}

#while (my $key = each %para_hash) {
#	my $value = $para_hash{$key}->{'gene2'};
#	print "$key\t$value\n";
#}
#exit;



# Set the minimum % identity to be used.  Higher values mean higher confidence in 
#  orthology
#my $min_ident=33;
# loop through minimum identity, decreasing by 5
for(my $min_ident=100; $min_ident>=0; $min_ident-=50) {

# Two files are to be read as input.  This one describes the interactions (network map)
open(network_file,"$network_file_name") || die "Error: file '$network_file_name' can not be opened\n";

# import the network map entries into a hash table
my %map_hash;
my @map_array;

while(<network_file>){
	my $entry = $_;
	# strip off newline characters
	$entry =~ s/\r//g;
	$entry =~ s/\n//g;
	# split the 3 columns into separate variables
	my ($map_gene_id, $interaction_type, $interaction_gene) = split("\t", $entry);
	# change gene ids to all caps
	$map_gene_id =~ tr/a-z/A-Z/;
	$interaction_gene =~ tr/a-z/A-Z/;
	
	
	# some genes have multiple interactions
	if(defined($map_hash{$map_gene_id})){
		my $oldhash = $map_hash{$map_gene_id};
		$oldhash->{'type'} = "$oldhash->{'type'}\t$interaction_type";;
		$oldhash->{'inter'} = "$oldhash->{'inter'}\t$interaction_gene";;
		$map_hash{$map_gene_id} = $oldhash;
	}else{
		my %hash;
		$hash{'map'} = $map_gene_id;
		$hash{'type'} = $interaction_type;
		$hash{'inter'} = $interaction_gene;
		$map_hash{$map_gene_id} = \%hash;
		push @map_array, $map_gene_id;
	}
}

# close file
close(network_file);


# This file gives the orthologs
open(ortho_file,"<$ortho_file_name") || die "Error: file '$ortho_file_name' can not be opened\n";

my %ortho_hash;
my @ortho_array;

# import the ortholog entries into a hash table
while(<ortho_file>){
	my $id1 = $_;
	# strip off newline characters
	$id1 =~ s/\r//g;
	$id1 =~ s/\n//g;
	# split the 3 columns into separate variables
	my ($gene_id, $ortho_id, $ident_id) = split("\t", $id1);
	# strip off gene id suffix and convert to all caps
	$gene_id =~ s/\-TAIR\-G//g;
	$gene_id =~ tr/a-z/A-Z/;
	$ortho_id =~ tr/a-z/A-Z/;
	# check to make sure the %identity is above the minimum defined above
	if($ident_id =~ /\d/ && $ident_id>=$min_ident){
		# if gene_id already found, add the new ortholog to hash
		if(defined($ortho_hash{$gene_id})){
			my $oldhash = $ortho_hash{$gene_id};
			$oldhash->{'ortho'} = "$oldhash->{'ortho'}\t$ortho_id";;
			$ortho_hash{$gene_id} = $oldhash;
		# if new gene_id, create hash element
		}else{
			my %hash;
			$hash{'ortho'} = $ortho_id;
			$hash{'gene'} = $gene_id;
			$ortho_hash{$gene_id} = \%hash;
			push @ortho_array, $gene_id;
		}
	}
}
#  close file
close(ortho_file);

#  open file for output
open(output_file,">$out_file_name" . "_$min_ident" . ".sif") || die "Error: file '$out_file_name' can not be opened for writing\n";
my $out_file = "$out_file_name" . "_$min_ident" . ".sif";


foreach my $id_entry (@map_array) {
	my $inter_type=$map_hash{$id_entry}->{'type'};
	my $gene_map=$map_hash{$id_entry}->{'inter'};
	#print "$id_entry\t$gene_map\n";
	
	# create array of interaction types
	my @inter_array = split("\t", $inter_type);
		
	# create array of interaction genes
	my @gene_map_array = split("\t", $gene_map);
	
	# check to see if ortholog exists for reference gene
	if(exists($ortho_hash{$id_entry})) {
		# create array of orthologs for reference gene
		my @gene_ortho_array = split("\t", $ortho_hash{$id_entry}->{'ortho'});
		# find the length of the array from map_hash to be used in loop
		my $max_index = @inter_array;

		# run for each ortholog for reference gene
		foreach my $gene_ortho (@gene_ortho_array) {
			# for each element of the arrays inter_array and gene_map_array
			for (my $i=0; $i<$max_index; $i++) {
				my $gene_map_ortho=$gene_map_array[$i];
				# check to make sure that an ortholog for the interaction gene is defined
				if(exists($ortho_hash{$gene_map_ortho})) {
					# create array of orthologs for interaction genes
					my @gene_map_ortho_array = split("\t", $ortho_hash{$gene_map_ortho}->{'ortho'});
					# loop through each interaction gene
					foreach my $gene_map_inter (@gene_map_ortho_array) {
						# if using paralogs, loop through paralogs for both genes
						if ($para) {
							my @gene_1_array;
							my @gene_2_array;
							#print "$gene_ortho\n";
							if (exists($para_hash{$gene_ortho})) {
								@gene_1_array = split("\t", $para_hash{$gene_ortho}->{'gene2'});
								#print "filling array 1\n";
							}
							# this will add back the original gene to the beginning of the array
							unshift(@gene_1_array, $gene_ortho);
							#print "$gene_ortho\n array @gene_1_array\n";
							#print "$gene_ortho\t$gene_map_inter\t\t\t$gene_1_array[0]\t$gene_1_array[1]\n";
							
							if (exists($para_hash{$gene_map_inter})) {
								@gene_2_array = split("\t", $para_hash{$gene_map_inter}->{'gene2'});
								#print "filling array 2\n";
							}
							unshift(@gene_2_array, $gene_map_inter);
							#print "@gene_2_array\n";
							
							foreach my $gene_1 (@gene_1_array) {
								foreach my $gene_2 (@gene_2_array) {
									if ($gene_1 gt $gene_2) {
										my $temp = $gene_1;
										$gene_1 = $gene_2;
										$gene_2 = $temp;
									}
									printf output_file "$gene_1\t$inter_array[$i]\t$gene_2\n";
									#printf output_file "$gene_1\n";
								}
							}
							
						} else {
							# check to see if the two entries are in alphatical order
							#   if they are not, swap them so they are
							#   Note that if the interaction type depends on direction, this breaks that
							if ($gene_ortho gt $gene_map_inter) { # if inter_type depends on direction, can extend the if statement to account for specific interactions to be skipped
								my $temp = $gene_ortho;
								$gene_ortho = $gene_map_inter;
								$gene_map_inter = $temp;
							}
						
							# print results to output file in sif format
							printf output_file "$gene_ortho\t$inter_array[$i]\t$gene_map_inter\n";
						}
					}

				}
			}
		}
	}

}

print "min_ident = $min_ident\n";

	
close(output_file);
# sort the file and get rid of duplicates
#   Note: this requires replacing the file, hence the mv command
system "sort $out_file | uniq > $out_file.tmp; mv $out_file.tmp $out_file";
# end loop of minimum identities ($min_ident)
} 
	
