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

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

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


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
open IUTABLE, $iUtable_file or die "Can't open iUtable file $iUtable_file\n";
chomp($header = <IUTABLE>); 
@distances = split/\t/,$header;
shift @distances;   #'gene'
shift @distances;   #'transcript
while(<IUTABLE>){
 chomp;
 my @line = split/\t/;
 my $gid = shift @line;
 my $transcript = shift @line;
 my $gene = "$gid\t$transcript";

 my $EUMA=0;
 for my $i (0..$#distances){
    $d=$distances[$i];
    $pUMA = $line[$i] * ($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)
}
close IUTABLE;