#!/usr/bin/perl
# Compute gEUMA for each gene from the gUtable file and insertlendis file and print out only genes whose gEUMA is larger than 0.

if(@ARGV<3) { print "usage: $0 gUtable_file insertlendis_file READLENGTH\n"; exit; }

my ($gUtable_file,$insertlendis_file,$READLENGTH) = @ARGV;

open GUTABLE, $gUtable_file or die "Can't open gUtable file $gUtable_file\n";
chomp($header = <GUTABLE>); 
@distances = split/\t/,$header;
shift @distances;   #'gene'
while(<GUTABLE>){
 chomp;
 my @line = split/\t/;
 my $gene = shift @line;
 for my $i (0..$#distances){
   $gUtable{$gene}{"$distances[$i]"}=$line[$i];
 }
}
close GUTABLE;


my $totalfreq=0;
open LENDIS, $insertlendis_file or die "Can't open insertlendis file $insertlendis_file\n";
while(<LENDIS>){
 chomp;
 my ($insert_size, $freq) = split/\t/;
 $insertlendis{"$insert_size"}=$freq;
 $totalfreq+=$freq;
}
close LENDIS;

print "gene\tgEUMA\n";  #header
for my $gene (sort keys %gUtable){
  my $EUMA=0;
  for my $d (keys %{$gUtable{$gene}}){
    $pUMA = $gUtable{$gene}{$d} * ($insertlendis{$d}/$totalfreq) ;
    $EUMA += $pUMA;
  }
  if($EUMA>=0.005) { printf "$gene\t%.2f\n",$EUMA; }  # use >=0.005  instead of >0, because I'm plotting two decimals after point. (in bp)
}