#!/usr/bin/perl

use warnings;
use strict;
use Data::Dumper;


# John Herbert January 2008 (j.m.herbert@bham.ac.uk). Best to be used on a Unix platform. 
# Perl wrapper to run expr_est from Susko and Roger, Bioinformatics. 2004 Sep 22;20(14):2279-87.

# Please edit these variables to point to the expr_est executable on your system. ==================
my $EXPR_EST = './expr_est';

#============================================================================================


my $usage = "# ================================================================================= #
Usage: $0 data_file output_format (one or multi)
The data file should be a tab delimited file in the format:
GNL3    11      14
PTHR2    2      4 
etc
Where GNL3 is the gene name, 11 is the number of hits in library A or library pool A and 14
is the number of hits to library B or library pool B.
Differentially expressed genes are found 
using Susko and Roger's expr_est statistical software.
Output format can be specified either \"one\" gene per line or \"multi\" genes per line
with the same result.
# ================================================================================= #\n";

# Catch the command line args. 
my $data_file = $ARGV[0] or die $usage;

my $out_argument = $ARGV[1] or die $usage;

#Check you have the correct output args
die "WARNING: The second argument is not one or multi \n$usage\n" if(($out_argument ne 'one') and ($out_argument ne 'multi'));

open(my $data_fh, '<', $data_file) or die $usage;

my $storage={};
# Collect the data and put into data structure. 
while(my $line = <$data_fh>){
  chomp($line);
  # Comments or hearder lines can begin with # and will be ignored. 
  next if($line =~ m/^#/); 
  my ($gene,$disease,$healthy)=split(/\s+/,$line);
  if("$gene\t$disease\t$healthy" !~ m/\S+\t\d+\t\d+/){
    warn "line in File not is not in the correct format, will ignore\t$line\n";
    next
  }
  push(@{$storage->{"$disease $healthy"}},$gene);
  
}

# input string for input file
my $file_string;
# Prepare the input for the expr_est program
foreach my $count ( keys %{$storage}){
  $file_string .= "$count\t".(scalar @{$storage->{$count}})."\n";
}

open(my $fh, '>', "temp_file.$$") or die;
# print the data to the temp file
print $fh $file_string;

close $fh or die;


# Count how many lines were in the file
my $num_lines = scalar keys %{$storage};


# Execute the est
my $susko_string = `$EXPR_EST  $num_lines 2 < temp_file.$$`;


# Choose which output to print out
if($out_argument eq 'one' ){
  print_output_one_gene_perline($susko_string);
}else{
  print_output_multi_gene_perline($susko_string);
}

# Get rid of the intermediate temporary file. 
#system("\\rm -f temp_file.$$");
unlink("temp_file.$$");

# ===========================================================================================
# ===========================================================================================

# Subroutine print_out
# Presents the results to the user with one gene per line.==== 
sub print_output_one_gene_perline{
  my $result = shift;
  
  my @data = split(/\n/,$result);
  foreach my $line (@data){
    $line =~ s/ /\t/g;
    my @inner = split(/\s+/,$line);
    my $search = "$inner[1] $inner[2]";
    foreach my $gene (@{$storage->{$search}}){
      print "$gene\t$line\n";
    }
  }
} 
# Subroutine print_out_multi_gene_perline
# =================================================
sub print_output_multi_gene_perline{
  my $result = shift;
  
  my @data = split(/\n/,$result);
  foreach my $line (@data){
    $line =~ s/ /\t/g;
    my @inner = split(/\s+/,$line);
    my $search = "$inner[1] $inner[2]";
    print "$line\t@{$storage->{$search}}\n";
  }
} 
    
