#########################################################################
# This script scores all sequences in a fasta file according to a
# motif (PWM), without using nucleosome information.
#
# Parameters:
# 0: the fasta file
# 1: the PWM file (in simple format: 4 lines, probs separated by tabs)
# 2: the output file (one line for each sequence: seqID\tscore)
#
# The background model is hard-coded.
#
#########################################################################

use strict;

if (scalar(@ARGV) < 3) { 
	print "Parameters:\n0: the fasta file\n",
	"1: the PWM file (in simple format: 4 lines, probs separated by tabs)\n",
	"2: the output file (one line for each sequence: seqID\\tscore)\n";
	exit;
}

my $fastafileALL = $ARGV[0]; 
my $PWMfile = $ARGV[1];

# Read the sequences and the motif
my ($hashrefall,$keysrefall) = getfasta_ordered($fastafileALL);
my ($ref, $W) = readMotifFromSimpleFile($PWMfile);
my @matrix = @$ref;

my %bkgr = ("A", 0.249815628838303, "C", 0.253044776990479,
            "G", 0.247065978347663, "T", 0.250073615823555);

# Compute [TF] as the Kd for the most specific site
my $logTF = 0;
for (my $i=0; $i<$W; $i++)
{
   my ($min, $pos) = min2(log($bkgr{"A"})-mylog($matrix[0][$i]), log($bkgr{"C"})-mylog($matrix[1][$i]),
                          log($bkgr{"G"})-mylog($matrix[2][$i]), log($bkgr{"T"})-mylog($matrix[3][$i]));
   $logTF += $min;

}
my $TF = exp($logTF);
#print "\n$TF\n";exit;

my $cnt=0;
my @sorted_probes = ();
my @sorted_probes_names = ();

foreach my $probe (@$keysrefall)
{
   my $seq = uc($$hashrefall{$probe});
   if ($seq =~ /[^ACGT]/)
   {
      die("Error: Sequence $probe is not a DNA sequence.\n";
   }
   
   my ($Wmer,$Wmer_r,$logscoreWmer,$logscoreWmer_r);
   my @logscoresWmer = ();
   my @logscoresWmer_r = ();

   for (my $i=0; $i<length($seq)-$W+1; $i++)
   {
      $Wmer = substr($seq,$i,$W);
      $Wmer_r = get_reverse($Wmer);
      $logscoreWmer = logscore_according_PWM($Wmer) - logscore_according_bkgr($Wmer);
      $logscoreWmer_r = logscore_according_PWM($Wmer_r) - logscore_according_bkgr($Wmer_r);
      push @logscoresWmer, $logscoreWmer;
      push @logscoresWmer_r, $logscoreWmer_r;      
   }
   
   # Compute the probability that no site is bound
   my $logprob = 0;
   foreach my $logscore ( (@logscoresWmer,@logscoresWmer_r) )
   {
      $logprob += mylog( 1/(1+( $TF * exp($logscore) )) );
   }
   my $prob = 1 - exp($logprob);
   
   # Add the probe to the sorted list
   my $i;
   for ($i=0; $i<scalar(@sorted_probes); $i++)
   {
   	if ($sorted_probes[$i] < $prob) { last; }
   }
   splice(@sorted_probes, $i, 0, ($prob));
   splice(@sorted_probes_names, $i, 0, ($probe));

   $cnt++;
}

open(OUT,">$ARGV[2]");
my ($F,$B,$roF,$roB) = (0,0,0,0);
for (my $i=0; $i<scalar(@sorted_probes_names); $i++)
{
	print OUT $sorted_probes_names[$i],"\t",$sorted_probes[$i],"\n";
}




#######################################################################
# Read a motif from a simple text file (4 lines: acgt, W columns)
#######################################################################
sub readMotifFromSimpleFile()
{
   my ($filename) = @_;
   my @matrix = ();
   my $W;
   open(IN, $filename) || die ("$filename does not exist\n");
   for (my $i=0; $i<4; $i++) {
     my $line = <IN>; chomp($line);
     my @array = split("\t", $line); $W = scalar(@array);
     push @matrix, \@array; 
   }
   close(IN);
   return (\@matrix, $W);
}

#######################################################################
# Get the reverse complement of a sequence
#######################################################################
sub get_reverse
{
  my ($string) = @_;
  my $new_string = reverse($string);
  $new_string =~ tr/aAcCgGtT/tTgGcCaA/;
  return $new_string;
}


#######################################################################
# Minimum of the elements in the array
#######################################################################
sub min2{
   my $min = $_[0]; my $pos = 0;
   for (my $i=0; $i<scalar(@_); $i++) { 
      if ($min > $_[$i]) { $min = $_[$i]; $pos=$i; }
   }
   return ($min,$pos);
}

#######################################################################
# 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);
}

#######################################################################
# Computing the log, including log(0)
#######################################################################
sub mylog
{
	if ($_[0]==0) { return 0; }
	else { return log($_[0]); }
}

#######################################################################
# Compute the logscore according to a PWM
#######################################################################
sub logscore_according_PWM
{
   my ($Wmer) = @_;$Wmer =~ tr/ACGTacgt/01230123/;
   my $score = 0;
   for (my $j=0; $j<length($Wmer); $j++) 
   {
	   $score += mylog($matrix[substr($Wmer,$j,1)][$j]);
   }
   return $score;
}

#######################################################################
# Compute the logscore according to the background model
#######################################################################
sub logscore_according_bkgr
{
   my ($Wmer) = @_;$Wmer =~ tr/acgt/ACGT/;
   my $score = 0;
   for (my $j=0; $j<length($Wmer); $j++) 
   {
	   $score += mylog($bkgr{substr($Wmer,$j,1)});
   }
   return $score;
}
