#########################################################################
# This script computes the enrichment of a motif in a set of foreground
# sequences as compared to a set of background sequences
#
# Parameters:
# 0: fasta file with the positive sequences
# 1: fasta file with the negative sequences
# 2: file containing the SORTED scores for pos&neg sequences
#
#########################################################################

use strict;
if (scalar(@ARGV) < 3) { 
	print "Parameters:\n",
	" 0: fasta file with the positive sequences\n",
	" 1: fasta file with the negative sequences\n",
	" 2: file containing the SORTED scores for pos&neg sequences\n";
	exit;
}

my ($hashref,$keysref) = getfasta_ordered($ARGV[0]);
my ($hashref_neg,$keysref_neg) = getfasta_ordered($ARGV[1]);
my $scoreFile = $ARGV[2];

my $enr = compute_enrich($scoreFile, $hashref, $hashref_neg);

print sprintf("%.6f\n",$enr+0.5);


#######################################################################
# Given a file with all the probes sorted according to the score
# and a set of probes (the foreground set), compute the area enrichment
#######################################################################
sub compute_enrich
{
   my ($scoresfile, $hashref, $hashref_neg) = @_;
   open(IN,$scoresfile) || die ("$scoresfile error");
   my @sorted = ();
   while (my $line=<IN>)
   {
       my @temp = split("\t",$line); pop @temp;
       push @sorted, join("\t",@temp);
   }
   
   my ($F,$B,$roF,$roB) = (0,0,0,0);
   my $ind = 0;
   my @pos = ();
   my @neg = ();
   for (my $i=0; $i<scalar(@sorted); $i++)
   {
       if (exists $$hashref{ $sorted[$i] })
       { $F++; $roF += ($ind+1); $ind++; push @pos, length($$hashref{ $sorted[$i] }); }
       elsif (exists $$hashref_neg{ $sorted[$i] })
       { $B++; $roB += ($ind+1); $ind++; push @neg, length($$hashref_neg{ $sorted[$i] }); }
   }
   #print "Pos: ", sum(@pos)/scalar(@pos),"\n";
   #print "Neg: ", sum(@neg)/scalar(@neg),"\n";   
   return 1/($B+$F) * ($roB/$B - $roF/$F);
   close(IN);
}

#######################################################################
# Sum the elements of an array
#######################################################################
sub sum{
  my $sum = 0;
  foreach (@_) { $sum += $_; }
  return $sum;
}

#######################################################################
# Given a fasta input file, return a reference to a hash with the
# name of each sequence as a key for each actual sequence. Thus, if you
# call '($hashref,$keysref) = getfasta($file);' you can see the sequence 
# 'read1' (given that is exists) by typing 'print $$hashref{"read1"};'.
# $keysref is an ordered list of the actual keys.
####################################################################### 
sub getfasta_ordered
{
  my($filename) = @_;
  chomp($filename);
  my $trying = "Unable to open $filename.\n";
  # Create a new file
  open(INPUTFILE, $filename) or die($trying);
  my @filedata = <INPUTFILE>;
  close INPUTFILE;
  my %seqs;
  my @keys = ();
  my $seqname;
  my $seq = '';
  foreach my $line (@filedata) {
    chomp($line);
    if($line =~/^\s*$/) { next; }     # Ignore blank lines
    elsif($line =~/^\s*#/) { next; }  # Ignore comments
    elsif($line =~ /^\>\s*(.*)\s*$/) {
      my $temp = $1;
      #$temp =~ s/\s//sg; 
      if(length($seq) > 0) {
        $seq =~ s/\s//sg;
        $seqs{$seqname} = $seq;
        push(@keys, $seqname);
        $seq = '';
      }
      $seqname = $temp;
      if (defined $seqs{$seqname}) {
         ##print "Error: Duplicate sequence $seqname in file $filename\n"; 
		}
      next;
    } else { $seq .= $line; }
  }
  $seq =~ s/\s//sg; 
  $seqs{$seqname} = $seq;
  push(@keys, $seqname);
  return (\%seqs, \@keys);
}